Cohesive cells are not forced to have a consistent orientation in v4.0.0, 4.1.0, or v4.1.1

baagaard-usgs commented 3 weeks ago

Describe the bug

Cohesive cells are not forced to have a consistent orientation. In many simple meshes (for example, all of our full-scale tests), the cohesive cells have consistent orientations. In more complex meshes, the orientation of a few cohesive cells is likely to be wrong. This creates errors in the slip, traction change, and stresses.

To Reproduce

PETSc branch: knepley/pylith PyLith fork: https://github.com/baagaard-usgs/pylith.git PyLith branch: fix-cohesive-cells-orientation

Steps to reproduce the behavior. Dump the fault_normal_dir vertex field in output/step01_slip-west_branch_info.h5.

cd examples/crustal-strikeslip-2d
pylith step01_slip.cfg
h5dump -d /vertex_fields/normal_dir output/step01_slip-west_branch_info.h5

BAD (normal dir should not flip)

HDF5 "output/step01_slip-west_branch_info.h5" {
DATASET "/vertex_fields/normal_dir" {
   DATATYPE  H5T_IEEE_F64LE
   DATASPACE  SIMPLE { ( 1, 13, 2 ) / ( H5S_UNLIMITED, 13, 2 ) }
   DATA {
   (0,0,0): -0.790654, 0.612263,
   (0,1,0): -0.767542, 0.640999,
   (0,2,0): 0.688091, -0.725624,
   (0,3,0): 0.636798, -0.771031,
   (0,4,0): 0.625379, -0.780321,
   (0,5,0): 0.634403, -0.773003,
   (0,6,0): -0.653672, 0.756778,
   (0,7,0): -0.663265, 0.748385,
   (0,8,0): 0.667521, -0.744591,
   (0,9,0): -0.668686, 0.743545,
   (0,10,0): -0.667078, 0.744988,
   (0,11,0): -0.667078, 0.744988,
   (0,12,0): -0.812921, 0.582374
   }
   ATTRIBUTE "timestepping" {
      DATATYPE  H5T_STD_I32LE
      DATASPACE  SCALAR
      DATA {
      (0): 1
      }
   }
   ATTRIBUTE "vector_field_type" {
      DATATYPE  H5T_STRING {
         STRSIZE 7;
         STRPAD H5T_STR_NULLTERM;
         CSET H5T_CSET_ASCII;
         CTYPE H5T_C_S1;
      }
      DATASPACE  SCALAR
      DATA {
      (0): "vector"
      }
   }
}
}

Attempted fix

diff --git a/libsrc/pylith/faults/TopologyOps.cc b/libsrc/pylith/faults/TopologyOps.cc
index 29a8cb1aa..742f6c1d1 100644
--- a/libsrc/pylith/faults/TopologyOps.cc
+++ b/libsrc/pylith/faults/TopologyOps.cc
@@ -216,6 +216,7 @@ pylith::faults::TopologyOps::create(pylith::topology::Mesh* mesh,
         }
     }
     err = DMLabelDestroy(&label);PYLITH_CHECK_ERROR(err);
+    err = DMPlexOrient(sdm);PYLITH_CHECK_ERROR(err);

     PetscReal lengthScale = 1.0;
     err = DMPlexGetScale(dm, PETSC_UNIT_LENGTH, &lengthScale);PYLITH_CHECK_ERROR(err);

This generates an error:

 -- Pre-initializing fault 'fault_main'.
[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
[0]PETSC ERROR: Invalid argument
[0]PETSC ERROR: Previously seen cells 3684 and 3685 do not match: Fault mesh is non-orientable

-- snip --

[0]PETSC ERROR: #1 DMPlexCheckFace_Old_Internal() at /Users/baagaard/software/unix/petsc-dev/src/dm/impls/plex/plexorient.c:218
[0]PETSC ERROR: #2 DMPlexOrient() at /Users/baagaard/software/unix/petsc-dev/src/dm/impls/plex/plexorient.c:421
[0]PETSC ERROR: #3 static void pylith::faults::TopologyOps::create(pylith::topology::Mesh *, const pylith::topology::Mesh &, PetscDMLabel, const int, const int)() at /Users/baagaard/src/cig/pylith/libsrc/pylith/faults/TopologyOps.cc:219
Fatal error. Calling MPI_Abort() to abort PyLith application.
Traceback (most recent call last):
  File "/Users/baagaard/software/unix/py3.12-venv/pylith-debug/lib/python3.12/site-packages/pylith/apps/PetscApplication.py", line 55, in onComputeNodes
    self.main(*args, **kwds)
  File "/Users/baagaard/software/unix/py3.12-venv/pylith-debug/lib/python3.12/site-packages/pylith/apps/PyLithApp.py", line 90, in main
    mesh = self.mesher.create(self.problem, interfaces)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/baagaard/software/unix/py3.12-venv/pylith-debug/lib/python3.12/site-packages/pylith/topology/MeshImporter.py", line 95, in create
    self._adjustTopology(mesh, faults, problem)
  File "/Users/baagaard/software/unix/py3.12-venv/pylith-debug/lib/python3.12/site-packages/pylith/topology/MeshGenerator.py", line 59, in _adjustTopology
    interface.adjustTopology(mesh)
  File "/Users/baagaard/software/unix/py3.12-venv/pylith-debug/lib/python3.12/site-packages/pylith/faults/faults.py", line 216, in adjustTopology
    return _faults.FaultCohesive_adjustTopology(self, mesh)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
RuntimeError: Error occurred while adjusting topology to create cohesive cells for fault 'fault_main'.

knepley commented 3 weeks ago

If you use knepley/fix-orientation-input in PETSc and baagaard-usgs/fix-cohesive-cells-orientation in PyLith, then the orientations are consistent

knepley commented 3 weeks ago

@baagaard-usgs I confirm that this also fixes the solve for degree 2 subduction-2d

geodynamics / pylith