Guide to the Meshing for Subsurface Flows in PETSc#

This 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.

Reading the ASCII Output#

For 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,

$ make -f ./gmakefile test globsearch="dm_impls_plex_tutorials-ex10_0"

which outputs

DM Object: Mesh 1 MPI process
  type: plex
Mesh in 3 dimensions:
  Number of 0-cells per rank: 8
  Number of 1-cells per rank: 14 (4)
  Number of 2-cells per rank: 9 (5)
  Number of 3-cells per rank: 2
Labels:
  celltype: 6 strata with value/size (3 (4), 9 (2), 0 (8), 2 (4), 1 (10), 5 (5))
  depth: 4 strata with value/size (0 (8), 1 (14), 2 (9), 3 (2))
  marker: 1 strata with value/size (1 (24))
  Face Sets: 1 strata with value/size (1 (12))

We 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.

Refinement of Simplex Meshes#

In PETSc, we can refine meshes uniformly and adaptively. Adaptive refinement can be controlled using tagging of cells, as well as specifying a target metric for the refined mesh. We will focus on the first type of adaptivity in this tutorial.

Adaptive Refinement#

Adaptive 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,

static PetscErrorCode CreateDomainLabel(DM dm)
{
  DMLabel  label;
  PetscInt cStart, cEnd, c;

  PetscFunctionBeginUser;
  PetscCall(DMGetCoordinatesLocalSetUp(dm));
  PetscCall(DMCreateLabel(dm, "Cell Sets"));
  PetscCall(DMGetLabel(dm, "Cell Sets", &label));
  PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
  for (c = cStart; c < cEnd; ++c) {
    PetscReal centroid[3], volume, x, y;

    PetscCall(DMPlexComputeCellGeometryFVM(dm, c, &volume, centroid, NULL));
    x = centroid[0];
    y = centroid[1];
    /* Headwaters are (0.0,0.25)--(0.1,0.75) */
    if ((x >= 0.0 && x < 0.1) && (y >= 0.25 && y <= 0.75)) {
      PetscCall(DMLabelSetValue(label, c, 1));
      continue;
    }
    /* River channel is (0.1,0.45)--(1.0,0.55) */
    if ((x >= 0.1 && x <= 1.0) && (y >= 0.45 && y <= 0.55)) {
      PetscCall(DMLabelSetValue(label, c, 2));
      continue;
    }
  }
  PetscFunctionReturn(0);
}

or you can label the mesh using a GUI, such as GMsh, and PETSc will read the label values from the input file.

We 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,

  while (adapt) {
    DM       dmAdapt;
    DMLabel  adaptLabel;
    PetscInt nAdaptLoc[2], nAdapt[2];

    adapt        = PETSC_FALSE;
    nAdaptLoc[0] = nAdaptLoc[1] = 0;
    nAdapt[0] = nAdapt[1] = 0;
    /* Adaptation is not preserving the domain label */
    PetscCall(DMHasLabel(dmCur, "Cell Sets", &hasLabel));
    if (!hasLabel) PetscCall(CreateDomainLabel(dmCur));
    PetscCall(DMGetLabel(dmCur, "Cell Sets", &label));
    PetscCall(DMLabelGetValueIS(label, &vIS));
    PetscCall(ISDuplicate(vIS, &valueIS));
    PetscCall(ISDestroy(&vIS));
    /* Sorting directly the label's value IS would corrupt the label so we duplicate the IS first */
    PetscCall(ISSort(valueIS));
    PetscCall(ISGetLocalSize(valueIS, &Nv));
    PetscCall(ISGetIndices(valueIS, &values));
    /* Construct adaptation label */
    PetscCall(DMLabelCreate(PETSC_COMM_SELF, "adapt", &adaptLabel));
    PetscCall(DMPlexGetHeightStratum(dmCur, 0, &cStart, &cEnd));
    for (c = cStart; c < cEnd; ++c) {
      PetscReal volume, centroid[3];
      PetscInt  value, vidx;

      PetscCall(DMPlexComputeCellGeometryFVM(dmCur, c, &volume, centroid, NULL));
      PetscCall(DMLabelGetValue(label, c, &value));
      if (value < 0) continue;
      PetscCall(PetscFindInt(value, Nv, values, &vidx));
      PetscCheck(vidx >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Value %" PetscInt_FMT " for cell %" PetscInt_FMT " does not exist in label", value, c);
      if (volume > volConst[vidx]) {
        PetscCall(DMLabelSetValue(adaptLabel, c, DM_ADAPT_REFINE));
        ++nAdaptLoc[0];
      }
      if (volume < volConst[vidx] * ratio) {
        PetscCall(DMLabelSetValue(adaptLabel, c, DM_ADAPT_COARSEN));
        ++nAdaptLoc[1];
      }
    }
    PetscCall(ISRestoreIndices(valueIS, &values));
    PetscCall(ISDestroy(&valueIS));
    PetscCallMPI(MPI_Allreduce(&nAdaptLoc, &nAdapt, 2, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)dmCur)));
    if (nAdapt[0]) {
      PetscCall(PetscInfo(dmCur, "Adapted mesh, marking %" PetscInt_FMT " cells for refinement, and %" PetscInt_FMT " cells for coarsening\n", nAdapt[0], nAdapt[1]));
      PetscCall(DMAdaptLabel(dmCur, adaptLabel, &dmAdapt));
      PetscCall(DMDestroy(&dmCur));
      PetscCall(DMViewFromOptions(dmAdapt, NULL, "-adapt_dm_view"));
      dmCur = dmAdapt;
      adapt = PETSC_TRUE;
    }
    PetscCall(DMLabelDestroy(&adaptLabel));
}

Test 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

    suffix: 3
    requires: triangle
    args: -init_dm_plex_dim 2 -init_dm_plex_box_faces 5,5 -adapt -volume_constraint_1 0.01 -volume_constraint_2 0.000625 -dm_extrude 10

We can look at a parallel run using extra options for the test system

$ make -f ./gmakefile test globsearch="dm_impls_plex_tutorials-ex10_3" EXTRA_OPTIONS="-dm_view hdf5:$PETSC_DIR/mesh.h5 -dm_partition_view -petscpartitioner_type parmetis" NP=5
$ $PETSC_DIR/lib/petsc/bin/petsc_gen_xmdf.py mesh.h5

../../_images/extrusionAdaptiveParallel.png — Fig. 25 **Parallel extruded mesh with adaptively refined surface**#

By turning on PetscInfo, we can see what decisions the refiner is making

$ make -f ./gmakefile test globsearch="dm_impls_plex_tutorials-ex10_3" EXTRA_OPTIONS="-info :dm"
#       > [0] AdaptMesh(): Adapted mesh, marking 12 cells for refinement, and 0 cells for coarsening
#       > [0] AdaptMesh(): Adapted mesh, marking 29 cells for refinement, and 0 cells for coarsening
#       > [0] AdaptMesh(): Adapted mesh, marking 84 cells for refinement, and 0 cells for coarsening
#       > [0] AdaptMesh(): Adapted mesh, marking 10 cells for refinement, and 0 cells for coarsening

Guide to the Meshing for Subsurface Flows in PETSc#

Reading the ASCII Output#

Refinement of Simplex Meshes#

Regular Refinement#

Adaptive Refinement#