Residual for refined mesh is assembled incorrectly

[baagaard-usgs]

Describe the bug

Refinement appears to be working for simulations without faults. However, the residual is wrong for refinement in simulations with a fault. The PyLith Jacobian looks correct but the finite-difference Jacobian does not; this is consistent with the residual being wrong. I don't see anything wrong in the topology or the labels. My wild guess is that we are taking the wrong path in the code when assembling the residual for a hybrid cell for the Lagrange multiplier field.

Update: Originally I thought the bug was only associated with meshes from Gmsh, but this does not appear to be the case as the same ASCII mesh also generates the error.

To Reproduce

PETSc branch: knepley/pylith-4.0
PyLith fork: https://github.com/baagaard-usgs/pylith.git
PyLith branch: add-refine-tests

cd $SRC/tests/fullscale/linearelasticity/faults-2d
pylith pylithapp_minmesh.cfg minmesh_noslip.cfg minmesh_noslip_quad.cfg --mesh_generator.refiner=pylith.topology.RefineUniform --mesh_generator.refiner.levels=1 --petsc.snes_test_jacobian --petsc.snes_test_jacobian_view

Diagram of mesh
refinedmesh.pdf

PyLith Jacobian

row 24: (4, 0.166667)  (5, 0.)  (6, 0.166667)  (7, 0.)  (8, -0.166667)  (9, 0.)  (10, -0.166667)  (11, 0.)  (14, 0.666667)  (15, 0.)  (22, -0.666667)  (23, 0.)  (24, 0.)  (25, 0.)  (26, 0.)  (27, 0.)  (28, 0.)  (29, 0.)
row 25: (4, 0.)  (5, 0.166667)  (6, 0.)  (7, 0.166667)  (8, 0.)  (9, -0.166667)  (10, 0.)  (11, -0.166667)  (14, 0.)  (15, 0.666667)  (22, 0.)  (23, -0.666667)  (24, 0.)  (25, 0.)  (26, 0.)  (27, 0.)  (28, 0.)  (29, 0.)
row 26: (4, 0.333333)  (5, 0.)  (8, -0.333333)  (9, 0.)  (14, 0.166667)  (15, 0.)  (22, -0.166667)  (23, 0.)  (24, 0.)  (25, 0.)  (26, 0.)  (27, 0.)
row 27: (4, 0.)  (5, 0.333333)  (8, 0.)  (9, -0.333333)  (14, 0.)  (15, 0.166667)  (22, 0.)  (23, -0.166667)  (24, 0.)  (25, 0.)  (26, 0.)  (27, 0.)
row 28: (6, 0.333333)  (7, 0.)  (10, -0.333333)  (11, 0.)  (14, 0.166667)  (15, 0.)  (22, -0.166667)  (23, 0.)  (24, 0.)  (25, 0.)  (28, 0.)  (29, 0.)
row 29: (6, 0.)  (7, 0.333333)  (10, 0.)  (11, -0.333333)  (14, 0.)  (15, 0.166667)  (22, 0.)  (23, -0.166667)  (24, 0.)  (25, 0.)  (28, 0.)  (29, 0.)

Finite-difference Jacobian

row 24: (2, -0.333333)  (4, 0.166667)  (6, 0.166667)  (14, 0.666667)  (18, -0.333333)  (20, -0.333333)  (24, 0.)
row 25: (3, -0.333333)  (5, 0.166667)  (7, 0.166667)  (15, 0.666667)  (19, -0.333333)  (21, -0.333333)  (25, 0.)
row 26: (2, -0.333333)  (4, 0.333333)  (14, 0.166667)  (18, -0.166667)  (26, 0.)
row 27: (3, -0.333333)  (5, 0.333333)  (15, 0.166667)  (19, -0.166667)  (27, 0.)
row 28: (2, -0.333333)  (6, 0.333333)  (14, 0.166667)  (20, -0.166667)  (28, 0.)
row 29: (3, -0.333333)  (7, 0.333333)  (15, 0.166667)  (21, -0.166667)  (29, 0.)

[knepley]

It will not let me push to your branch. I am not sure what is wrong. Here is the diff that fixes it

 141add-refine-tests *$:/PETSc3/cig/pylith-git$ git log -1 -p
commit 9da7a4f5e51dfefd3d8139f980d3ef55c91536e0 (HEAD -> add-refine-tests)
Author: Matthew G. Knepley <knepley@gmail.com>
Date:   Sun Feb 4 11:19:41 2024 -0500

    Refine: Must order cohesive supports after refinement


diff --git a/libsrc/pylith/topology/RefineUniform.cc b/libsrc/pylith/topology/RefineUniform.cc
index efd90cda9..be5e45059 100644
--- a/libsrc/pylith/topology/RefineUniform.cc
+++ b/libsrc/pylith/topology/RefineUniform.cc
@@ -79,6 +79,7 @@ pylith::topology::RefineUniform::refine(Mesh* const newMesh,

         err = DMDestroy(&dmCur);PYLITH_CHECK_ERROR(err);
     } // for
+    err = DMPlexReorderCohesiveSupports(dmNew);PYLITH_CHECK_ERROR(err);

     newMesh->setDM(dmNew);

@@ -117,6 +118,7 @@ pylith::topology::RefineUniform::refine(Mesh* const newMesh,

     topology::MeshOps::checkTopology(*newMesh);

     // newMesh->view("REFINED_MESH", "::ascii_info_detail");
+    err = DMViewFromOptions(dmNew, NULL, "-pylith_ref_dm_view");PYLITH_CHECK_ERROR(err);


     PYLITH_METHOD_END;
 } // refine

[baagaard-usgs]

This fix resolves the issue.

• Finish adding refinement tests in tests/fullscale/linearelasticity/faults-2d
• Add static C++ method in MeshOps.cc for computing the centroid of a cell
• Update MMS tests to include refinement in tests/mmstests/linearelasticity/faults-2d