1# Meshing for Subsurface Flows in PETSc 2 3This tutorials guides users in creating meshes for the TDyCore simulation framework for subsurface flows. The user inputs a surface mesh, a refinement prescription, and an extrusion prescription in order to create the simulation mesh. 4 5Reading the ASCII Output 6 7For example, a very simple mesh would start with a square surface mesh divided into two triangles, which is then extruded to form two triangular prisms. This is the first test in the DMPlex tutorial code ex10, 8 9```console 10$ make -f ./gmakefile test search="dm_impls_plex_tutorials-ex10_0" 11``` 12 13which outputs 14 15```{literalinclude} /../src/dm/impls/plex/tutorials/output/ex10_0.out 16``` 17 18We can see that there are two 3-cells, meaning three-dimensional cells, and from the `celltype` label we see that those cells have celltype 9, meaning they are triangular prisms. The original surface mesh had 5 edges, so we would expect 10 edges for the two surfaces and four edges connecting those surfaces. This is exactly what we see, since there are 14 1-cells, but 4 of them noted in parentheses are tensor cells created by extrusion. We can see this another way in the celltype label, where there are ten mesh points of type 1, meaning segments, and four mesh points of type 2, meaning tensor products of a vertex and segment. Similarly, there are 9 2-cells, but 5 of them stretch between the two surfaces, meaning they are tensor products of two segments. 19 20Regular Refinement of Simplex Meshes 21 22We can regularly refine the surface before extrusion using `-dm_refine <k>`, where `k` is the number of refinements, 23 24```console 25$ make -f ./gmakefile test search="dm_impls_plex_tutorials-ex10_1" EXTRA_OPTIONS="-srf_dm_refine 2 -srf_dm_view draw -draw_save $PETSC_DIR/surface.png -draw_save_single_file" 26``` 27 28which produces the following surface 29 30:::{figure} /images/tutorials/meshing/surface.png 31:align: center 32 33**Surface mesh refined twice** 34::: 35 36and the extruded mesh can be visualized using VTK. Here I make the image using Paraview, and give the extrusion 3 layers 37 38```console 39$ make -f ./gmakefile test search="dm_impls_plex_tutorials-ex10_1" EXTRA_OPTIONS="-dm_view hdf5:$PETSC_DIR/mesh.h5 -dm_extrude 3" 40$ $PETSC_DIR/lib/petsc/bin/petsc_gen_xmdf.py mesh.h5 41``` 42 43:::{figure} /images/tutorials/meshing/extrusion.png 44:align: center 45 46**Extruded mesh with refined surface** 47::: 48 49We can similarly look at this in parallel. Test 2 uses three refinements and three extrusion layers on five processes 50 51```console 52$ make -f ./gmakefile test search="dm_impls_plex_tutorials-ex10_2" EXTRA_OPTIONS="-dm_view hdf5:$PETSC_DIR/mesh.h5 -dm_partition_view -petscpartitioner_type parmetis" 53$ $PETSC_DIR/lib/petsc/bin/petsc_gen_xmdf.py mesh.h5 54``` 55 56:::{figure} /images/tutorials/meshing/extrusionParallel.png 57:align: center 58 59**Parallel extruded mesh with refined surface** 60::: 61 62Adaptive Refinement of Simplex Meshes 63 64Adaptive refinement of simplicial meshes is somewhat tricky when we demand that the meshes be conforming, as we do in this case. We would like different grid cells to have different levels of refinement, for example headwaters cells in a watershed be refined twice, while river channel cells be refined four times. In order to differentiate between cells, we first mark the cells on the surface using a `DMLabel`. We can do this programmatically, 65 66```{literalinclude} /../src/dm/impls/plex/tutorials/ex10.c 67:append: '}' 68:end-at: PetscFunctionReturn(PETSC_SUCCESS); 69:start-at: static PetscErrorCode CreateDomainLabel( 70``` 71 72or you can label the mesh using a GUI, such as GMsh, and PETSc will read the label values from the input file. 73 74We next create a label marking each cell in the mesh with an action, such as `DM_ADAPT_REFINE` or `DM_ADAPT_COARSEN`. We do this based on a volume constraint, namely that cells with a certain label value should have a certain volume. You could, of course, choose a more complex strategy, but here we just want a clear criterion. We can give volume constraints for label value `v` using the command line argument `-volume_constraint_<v> <vol>`. The mesh is then refined iteratively, checking the volume constraints each time, 75 76```{literalinclude} /../src/dm/impls/plex/tutorials/ex10.c 77:append: '}' 78:end-at: PetscCall(DMLabelDestroy(&adaptLabel)); 79:start-at: while (adapt) { 80``` 81 82Test 3 from `ex10` constrains the headwater cells (with marker 1) to have volume less than 0.01, and the river channel cells (with marker 2) to be smaller than 0.000625 83 84```{literalinclude} /../src/dm/impls/plex/tutorials/ex10.c 85:lines: 1-3 86:start-at: 'suffix: 3' 87``` 88 89We can look at a parallel run using extra options for the test system 90 91```console 92$ make -f ./gmakefile test search="dm_impls_plex_tutorials-ex10_3" EXTRA_OPTIONS="-dm_view hdf5:$PETSC_DIR/mesh.h5 -dm_partition_view -petscpartitioner_type parmetis" NP=5 93$ $PETSC_DIR/lib/petsc/bin/petsc_gen_xmdf.py mesh.h5 94``` 95 96:::{figure} /images/tutorials/meshing/extrusionAdaptiveParallel.png 97:align: center 98 99**Parallel extruded mesh with adaptively refined surface** 100::: 101 102By turning on `PetscInfo`, we can see what decisions the refiner is making 103 104```console 105$ make -f ./gmakefile test search="dm_impls_plex_tutorials-ex10_3" EXTRA_OPTIONS="-info :dm" 106# > [0] AdaptMesh(): Adapted mesh, marking 12 cells for refinement, and 0 cells for coarsening 107# > [0] AdaptMesh(): Adapted mesh, marking 29 cells for refinement, and 0 cells for coarsening 108# > [0] AdaptMesh(): Adapted mesh, marking 84 cells for refinement, and 0 cells for coarsening 109# > [0] AdaptMesh(): Adapted mesh, marking 10 cells for refinement, and 0 cells for coarsening 110``` 111
