Skip to content

Latest commit

 

History

History

stokes

README.md

stokes/

Copyright 2018--2024 Ed Bueler

Introduction

The Python programs in this directory use advanced, open-source libraries and tools:

  • Firedrake, a Python finite element library
    • PETSc, a solver library typically installed by Firedrake
  • Gmsh, a mesher
  • Paraview, for visualization of mesh functions

The .py programs here are documented by the current README and by doc.pdf in this directory.

Installation

Usage

Activate your Firedrake virtual environment (venv) first:

    $ source ~/firedrake/bin/activate

Generate the domain geometry using domain.py and then mesh it using Gmsh:

    (firedrake) $ ./domain.py -o glacier.geo      # create domain outline
    (firedrake) $ gmsh -2 glacier.geo             # writes glacier.msh

You can, at this point, visualize the generated mesh with Gmsh:

    (firedrake) $ gmsh glacier.msh

Now run the Stokes solver flow.py on the mesh, which will write velocity and pressure fields into glacier.pvd.

    (firedrake) $ ./flow.py -mesh glacier.msh -o glacier.pvd

(Note that ./flow.py is equivalent to python3 flow.py because executable permissions have been set on flow.py.) One may visualize the result using Paraview:

    (firedrake) $ paraview glacier.pvd

An alternate visualization plots surface values into an image file surf.png:

    (firedrake) $ ./flow.py -mesh glacier.msh -osurface surf.png

Slab-on-slope

Set the height of the bedrock step to zero when creating the domain geometry:

    (firedrake) $ ./domain.py -bs 0.0 -o slab.geo
    (firedrake) $ gmsh -2 slab.geo
    (firedrake) $ ./flow.py -slab -mesh slab.msh

Here the numerical error is displayed because the exact solution is known. See doc.pdf for the slab-on-slope solution.

Mesh refinement

The default mesh generated by domain.py, i.e. glacier.msh after meshing, has a typical mesh size of 80 m with grid resolution a factor of 4 finer near the interior corners created by the bedrock step, giving 20 m resolution at these corners.

There are four refinement methods to get finer resolution:

  1. One may use the script domain.py to create a refined (or locally-refined) mesh, before it is read by flow.py. The program domain.py optionally takes a target mesh size (-hmesh H) and/or a uniform refinement factor (-refine X). Another option controls the additional refinement at the interior corner (-refine_corner Y, where Y is a factor). The default case corresponds to -hmesh 80 -refine 1 -refine_corner 4. The following example creates a refined mesh with target mesh size varying from 20 m to about 3 m near the interior corner. The resulting grid has about 15 times as many elements as the default mesh:

     (firedrake) $ ./domain.py -refine 4 -refine_corner 8 -o fine1.geo
 (firedrake) $ gmsh -2 fine1.geo
 (firedrake) $ ./flow.py -mesh fine1.msh

  2. One may use Gmsh to refine an existing mesh in .msh, specifically by splitting each triangular cell (element) into four similar triangles:

     (firedrake) $ ./domain.py -refine 2 -refine_corner 4 -o start.geo
 (firedrake) $ gmsh -2 start.geo
 (firedrake) $ gmsh -refine start.msh -o fine2.msh
 (firedrake) $ ./flow.py -mesh fine2.msh

    As with method 1, this refinement is done before the mesh is read by flow.py.

  3. One may read an initial mesh and ask the flow.py script, via Firedrake, to refine:

     (firedrake) $ ./flow.py -mesh start.msh -refine 1

    This is simpler usage compared to method 2, and with the same result.

  4. One may use grid-sequencing. This solves the problem on the initially-read mesh, interpolates onto the next finer one, and then solves there

     (firedrake) $ ./flow.py -mesh start.msh -sequence 1

    This is slightly faster than methods 2 and 3, but the result should be the same.

Solver performance information

Firedrake calls PETSc to solve the Glen-law Stokes equations. Because this is a nonlinear problem, a PETSc SNES solver object does the Newton iterations. This solver has option prefix -s_. Basic information on solver performance includes reporting on the Newton iterations in each Stokes solve:

    (firedrake) $ ./flow.py -mesh glacier.msh -s_snes_monitor -s_snes_converged_reason

A detailed description of the PETSc solver comes from -s_snes_view. Inside the SNES solver are KSP (Krylov space) and PC (preconditioner) components.

For low resolution cases, namely grids with few nodes, one can write-out the system matrix as a Matlab text file foo.m as follows:

    (firedrake) $ ./flow.py -mesh glacier.msh -s_ksp_view_mat :foo.m:ascii_matlab

Note foo.m is now a command in Matlab/Octave, and one may use spy to see the structure of the matrix. Alternatively, the matrix can be read into Python.

For more on Firedrake and PETSc solvers see their respective online documentation. Or see my book PETSc for Partial Differential Equations.

Getting help

There are several ways to get help:

  • ./flow.py -flowhelp shows options specific to flow.py itself.
  • ./flow.py -help intro shows PETSc version information.
  • ./flow.py -mesh glacier.msh -help lists many PETSc solver options. (Option -mesh is needed so that the program proceeds to where a PETSc solver has been created.)

Testing

Activate your Firedrake venv, as above, and then do:

    $ pytest .