Update OutputSubfield to allow refinement of mesh with higher order discretizations

[baagaard-usgs]

Goal: Allow output of fields at a spatial resolution that represents the solution and basis order compatible with visualization tools.

Workflow

• Initialization
  1. Refine submesh to get desired resolution
• Output
  1. Project field to "resolution" basis order on submesh
  2. Project field to output basis order on refined submesh

This will be relatively easy to implement because it is confined to adding just a little bit of functionality to OutputSubfield and OutputObserver (Pyre settings).

[baagaard-usgs]

Updated algorithm

1. Create refined version of PETSc DM during initialization. Refer to DMPlexCreateHighOrderSurrogate_Internal
2. Project output field to "original resolution" PETSc DM.
3. Use interpolation to interpolate field from "original resolution" PETSc DM to "refined" PETSc DM.
   • Need to create interpolation for PETSc DM as well and then interpolate field.
   • Get PETSc Vec from DM (for field), so we really need to just manage the DMs and matrices.

[baagaard-usgs]

@knepley I am getting a segmentation fault with DMFilter() to remove the hanging cells. It is almost certainly an uninitialized value as it won't fail when running in the debugger; it only segfaults if not using the debugger.

Stacktrace

[0]PETSC ERROR: #1 PetscFindInt() at /software/baagaard/petsc-dev/src/sys/utils/sorti.c:506
[0]PETSC ERROR: #2 DMPlexFilterLabels_Internal() at /software/baagaard/petsc-dev/src/dm/impls/plex/plexsubmesh.c:3361
[0]PETSC ERROR: #3 DMPlexCreateSubmeshGeneric_Interpolated() at /software/baagaard/petsc-dev/src/dm/impls/plex/plexsubmesh.c:3805
[0]PETSC ERROR: #4 DMPlexFilter() at /software/baagaard/petsc-dev/src/dm/impls/plex/plexsubmesh.c:4213
[0]PETSC ERROR: #5 static _p_DM* pylith::topology::MeshOps::removeHangingCells(_p_DM* const&)() at /home/baagaard/src/cig/pylith/libsrc/pylith/topolo
gy/MeshOps.cc:350

Valgrind report

==1479870== Invalid read of size 4
==1479870==    at 0xA73CA67: PetscFindInt (sorti.c:519)
==1479870==    by 0xC7127C6: DMPlexFilterPoint_Internal (plexsubmesh.c:3301)
==1479870==    by 0xC713D25: DMPlexFilterLabels_Internal (plexsubmesh.c:3362)
==1479870==    by 0xC720893: DMPlexCreateSubmeshGeneric_Interpolated (plexsubmesh.c:3805)
==1479870==    by 0xC72AE34: DMPlexFilter (plexsubmesh.c:4213)
...
==1479870== Conditional jump or move depends on uninitialised value(s)
==1479870==    at 0xA73CA4D: PetscFindInt (sorti.c:513)
==1479870==    by 0xC7127C6: DMPlexFilterPoint_Internal (plexsubmesh.c:3301)
==1479870==    by 0xC713D25: DMPlexFilterLabels_Internal (plexsubmesh.c:3362)
==1479870==    by 0xC720893: DMPlexCreateSubmeshGeneric_Interpolated (plexsubmesh.c:3805)
==1479870==    by 0xC72AE34: DMPlexFilter (plexsubmesh.c:4213)
...
==1479870== Use of uninitialised value of size 8
==1479870==    at 0xA73CA67: PetscFindInt (sorti.c:519)
==1479870==    by 0xC7127C6: DMPlexFilterPoint_Internal (plexsubmesh.c:3301)
==1479870==    by 0xC713D25: DMPlexFilterLabels_Internal (plexsubmesh.c:3362)
==1479870==    by 0xC720893: DMPlexCreateSubmeshGeneric_Interpolated (plexsubmesh.c:3805)
==1479870==    by 0xC72AE34: DMPlexFilter (plexsubmesh.c:4213)

At plexsubmesh.c line 3437, dim=1, so numSubPoints has size 2. At plexsubmesh.c line 3361, we get a depth of 2 and pass numSubPoints[2] to DMPlexFilterPoint which is an invalid read. It looks like it is trying to do something with points corresponding to depths/heights which are not in the filtered mesh. Should the numSubPoints have a size matching the depth/height of the unfiltered mesh?

[baagaard-usgs]

@knepley Here is a solution:

diff --git a/src/dm/impls/plex/plexsubmesh.c b/src/dm/impls/plex/plexsubmesh.c
index 0d5cb6cb1d9..dc5a4f2a5ab 100644
--- a/src/dm/impls/plex/plexsubmesh.c
+++ b/src/dm/impls/plex/plexsubmesh.c
@@ -3356,9 +3356,11 @@ static PetscErrorCode DMPlexFilterLabels_Internal(DM dm, const PetscInt numSubPo
       PetscCall(ISGetIndices(pointIS, &points));
       for (p = 0; p < Np; ++p) {
         const PetscInt point = points[p];
-        PetscInt       subp;
+        PetscInt       subp, subdepth;

         PetscCall(DMPlexGetPointDepth(dm, point, &d));
+        PetscCall(DMPlexGetDepth(subdm, &subdepth));
+        if (d > subdepth)  continue;
         subp = DMPlexFilterPoint_Internal(point, firstSubPoint[d], numSubPoints[d], subpoints[d]);
         if (subp >= 0) PetscCall(DMLabelSetValue(newlabel, subp, values[v]));
       }

[baagaard-usgs]

Remaining issues

• [ ] Switch to output basis order at finest interpolation
• [ ] Add refinement to OutputSoln
• [ ] Update documentation
• [ ] Update examples/strikeslip-2d