Merge branch 'geodynamics:main' into point_source

configure.ac

@@ -109,7 +109,7 @@ AC_CHECK_HEADER([mpi.h], [], [AC_MSG_ERROR([header 'mpi.h' not found])])
 
 dnl PETSC
 AC_LANG(C)
-CIT_PATH_PETSC([3.19.4])
+CIT_PATH_PETSC([3.19.5])
 CIT_HEADER_PETSC
 CIT_CHECK_LIB_PETSC
 
@@ -254,6 +254,7 @@ AC_CONFIG_FILES([Makefile
 		tests/fullscale/linearelasticity/nofaults-3d/Makefile
 		tests/fullscale/linearelasticity/faults-2d/Makefile
 		tests/fullscale/linearelasticity/faults-3d/Makefile
+		tests/fullscale/linearelasticity/faults-3d-buried/Makefile
 		tests/fullscale/linearelasticity/greensfns-2d/Makefile
 		tests/fullscale/incompressibleelasticity/Makefile
 		tests/fullscale/incompressibleelasticity/nofaults-2d/Makefile

libsrc/pylith/faults/FaultCohesive.cc

@@ -375,6 +375,10 @@ pylith::faults::FaultCohesive::createConstraints(const pylith::topology::Field&
             } // if
         } // for
     } // for
+    if (pointIS) {
+        err = ISRestoreIndices(pointIS, &points);PYLITH_CHECK_ERROR(err);
+        err = ISDestroy(&pointIS);PYLITH_CHECK_ERROR(err);
+    } // if
 
     std::vector<pylith::feassemble::Constraint*> constraintArray;
     pylith::feassemble::ConstraintSimple *constraint = new pylith::feassemble::ConstraintSimple(this);assert(constraint);

libsrc/pylith/meshio/DataWriter.cc

@@ -248,4 +248,33 @@ pylith::meshio::DataWriter::getCoordsGlobalVec(PetscVec* coordsGlobalVec,
 } // getCoordsGlobalVec
 
 
+// ----------------------------------------------------------------------
+void
+pylith::meshio::DataWriter::_writeVec(PetscVec vector,
+                                      PetscViewer viewer) {
+    // From plexhdf5.c DMPlexGlobalVectorView_HDF5_Internal
+
+    /* To save vec in where we want, we create a new Vec (temp) with           */
+    /* VecCreate(), wrap the vec data in temp, and call VecView(temp, viewer). */
+    PetscVec temp;
+    const PetscScalar *array;
+    const char* vecName;
+    PetscLayout map;
+    PetscErrorCode err = PETSC_SUCCESS;
+
+    err = VecCreate(PetscObjectComm((PetscObject)vector), &temp);PYLITH_CHECK_ERROR(err);
+    err = PetscObjectGetName((PetscObject)vector, &vecName);PYLITH_CHECK_ERROR(err);
+    err = PetscObjectSetName((PetscObject)temp, vecName);PYLITH_CHECK_ERROR(err);
+    err = VecGetLayout(vector, &map);PYLITH_CHECK_ERROR(err);
+    err = VecSetLayout(temp, map);PYLITH_CHECK_ERROR(err);
+    err = VecSetUp(temp);PYLITH_CHECK_ERROR(err);
+    err = VecGetArrayRead(vector, &array);PYLITH_CHECK_ERROR(err);
+    err = VecPlaceArray(temp, array);PYLITH_CHECK_ERROR(err);
+    err = VecView(temp, viewer);PYLITH_CHECK_ERROR(err);
+    err = VecResetArray(temp);PYLITH_CHECK_ERROR(err);
+    err = VecRestoreArrayRead(vector, &array);PYLITH_CHECK_ERROR(err);
+    err = VecDestroy(&temp);PYLITH_CHECK_ERROR(err);
+} // _writeVec
+
+
 // End of file

libsrc/pylith/meshio/DataWriter.hh

@@ -146,6 +146,18 @@ protected:
     void getCoordsGlobalVec(PetscVec* coordinatesVec,
                             const pylith::topology::Mesh& mesh);
 
+    /** Write PETSc Vec using viewer.
+     *
+     * This method allows us to circumvent the usual place a vector is written and write it to
+     * a desired location. See DMPlexGlobalVectorView_HDF5_Internal() in plexhdf5.c in PETSc.
+     *
+     * @param[in] vector PETSc Vec to write.
+     * @param[inout] viewer PETSc viewer to use to write vector.
+     */
+    static
+    void _writeVec(PetscVec vector,
+                   PetscViewer viewer);
+
     // NOT IMPLEMENTED //////////////////////////////////////////////////////
 private:
 

libsrc/pylith/meshio/DataWriterHDF5.cc

@@ -40,15 +40,6 @@
 #include <sstream> // USES std::ostringstream
 #include <stdexcept> // USES std::runtime_error
 
-extern "C" {
-    extern PetscErrorCode VecView_Seq(Vec,
-                                      PetscViewer);
-
-    extern PetscErrorCode VecView_MPI(Vec,
-                                      PetscViewer);
-
-}
-
 #if H5_VERS_MAJOR == 1 && H5_VERS_MINOR >= 8
 #define PYLITH_HDF5_USE_API_18
 #endif
@@ -207,13 +198,7 @@ pylith::meshio::DataWriterHDF5::writeVertexField(const PylithScalar t,
         err = PetscViewerHDF5SetTimestep(_viewer, istep);PYLITH_CHECK_ERROR(err);
 
         PetscVec vector = subfield.getVector();assert(vector);
-        PetscBool isseq;
-        err = PetscObjectTypeCompare((PetscObject) vector, VECSEQ, &isseq);PYLITH_CHECK_ERROR(err);
-        if (isseq) {
-            err = VecView_Seq(vector, _viewer);PYLITH_CHECK_ERROR(err);
-        } else {
-            err = VecView_MPI(vector, _viewer);PYLITH_CHECK_ERROR(err);
-        }
+        DataWriter::_writeVec(vector, _viewer);PYLITH_CHECK_ERROR(err);
         err = PetscViewerHDF5PopTimestepping(_viewer);PYLITH_CHECK_ERROR(err);
         err = PetscViewerHDF5PopGroup(_viewer);PYLITH_CHECK_ERROR(err);
 
@@ -276,13 +261,7 @@ pylith::meshio::DataWriterHDF5::writeCellField(const PylithScalar t,
         err = PetscViewerHDF5SetTimestep(_viewer, istep);PYLITH_CHECK_ERROR(err);
 
         PetscVec vector = subfield.getVector();assert(vector);
-        PetscBool isseq;
-        err = PetscObjectTypeCompare((PetscObject) vector, VECSEQ, &isseq);PYLITH_CHECK_ERROR(err);
-        if (isseq) {
-            err = VecView_Seq(vector, _viewer);PYLITH_CHECK_ERROR(err);
-        } else {
-            err = VecView_MPI(vector, _viewer);PYLITH_CHECK_ERROR(err);
-        } // if/else
+        DataWriter::_writeVec(vector, _viewer);
         err = PetscViewerHDF5PopTimestepping(_viewer);PYLITH_CHECK_ERROR(err);
         err = PetscViewerHDF5PopGroup(_viewer);PYLITH_CHECK_ERROR(err);
 

libsrc/pylith/meshio/DataWriterHDF5Ext.cc

@@ -37,15 +37,6 @@
 #include <sstream> // USES std::ostringstream
 #include <stdexcept> // USES std::runtime_error
 
-extern "C" {
-    extern PetscErrorCode VecView_Seq(Vec,
-                                      PetscViewer);
-
-    extern PetscErrorCode VecView_MPI(Vec,
-                                      PetscViewer);
-
-}
-
 // ----------------------------------------------------------------------
 // Constructor
 pylith::meshio::DataWriterHDF5Ext::DataWriterHDF5Ext(void) :
@@ -195,13 +186,7 @@ pylith::meshio::DataWriterHDF5Ext::writeVertexField(const PylithScalar t,
         assert(binaryViewer);
 
         PetscVec vector = subfield.getVector();assert(vector);
-        PetscBool isseq;
-        err = PetscObjectTypeCompare((PetscObject) vector, VECSEQ, &isseq);PYLITH_CHECK_ERROR(err);
-        if (isseq) {
-            err = VecView_Seq(vector, binaryViewer);PYLITH_CHECK_ERROR(err);
-        } else {
-            err = VecView_MPI(vector, binaryViewer);PYLITH_CHECK_ERROR(err);
-        } // if/else
+        DataWriter::_writeVec(vector, binaryViewer);
 
         ExternalDataset& datasetInfo = _datasets[name];
         ++datasetInfo.numTimeSteps;
@@ -339,13 +324,7 @@ pylith::meshio::DataWriterHDF5Ext::writeCellField(const PylithScalar t,
         assert(binaryViewer);
 
         PetscVec vector = subfield.getVector();assert(vector);
-        PetscBool isseq;
-        err = PetscObjectTypeCompare((PetscObject) vector, VECSEQ, &isseq);PYLITH_CHECK_ERROR(err);
-        if (isseq) {
-            err = VecView_Seq(vector, binaryViewer);PYLITH_CHECK_ERROR(err);
-        } else {
-            err = VecView_MPI(vector, binaryViewer);PYLITH_CHECK_ERROR(err);
-        } // if/else
+        DataWriter::_writeVec(vector, binaryViewer);
 
         ExternalDataset& datasetInfo = _datasets[name];
         ++datasetInfo.numTimeSteps;

pylith/faults/KinSrcTimeHistory.py

@@ -61,12 +61,12 @@ def __init__(self, name="kinsrctimehistory"):
         KinSrc.__init__(self, name)
         return
 
-    def preinitialize(self):
+    def preinitialize(self, problem):
         """Do pre-initialization setup.
         """
-        KinSrc.preinitialize(self)
+        KinSrc.preinitialize(self, problem)
 
-        ModuleKinSrc.setTimeHistory(self, self.dbTimeHistory)
+        ModuleKinSrc.setTimeHistoryDB(self, self.dbTimeHistory)
         return
 
     # PRIVATE METHODS ////////////////////////////////////////////////////

tests/fullscale/linearelasticity/Makefile.am

@@ -20,6 +20,7 @@ SUBDIRS = \
 	nofaults-3d \
 	faults-2d \
 	faults-3d \
+	faults-3d-buried \
 	greensfns-2d
 
 

tests/fullscale/linearelasticity/faults-2d/Makefile.am

@@ -49,7 +49,8 @@ dist_noinst_DATA = \
 	threeblocks_ic_tri.cfg \
 	shearnoslip.cfg \
 	shearnoslip_quad.cfg \
-	shearnoslip_tri.cfg
+	shearnoslip_tri.cfg \
+	zeroslipfn.timedb
 
 
 noinst_TMP =

tests/fullscale/linearelasticity/faults-2d/shearnoslip.cfg

@@ -20,12 +20,18 @@ features = [
 # ----------------------------------------------------------------------
 # faults
 # ----------------------------------------------------------------------
+[pylithapp.problem.interfaces.fault_xmid.eq_ruptures]
+rupture = pylith.faults.KinSrcTimeHistory
+
 [pylithapp.problem.interfaces.fault_xmid.eq_ruptures.rupture]
 db_auxiliary_field = spatialdata.spatialdb.UniformDB
 db_auxiliary_field.description = Fault rupture auxiliary field spatial database
 db_auxiliary_field.values = [initiation_time, final_slip_left_lateral, final_slip_opening]
 db_auxiliary_field.data = [0.0*s, 0.0*m, 0.0*m]
 
+time_history.description = Slip time function time history
+time_history.filename = zeroslipfn.timedb
+
 [pylithapp.problem.interfaces.fault_xneg.eq_ruptures.rupture]
 db_auxiliary_field = spatialdata.spatialdb.UniformDB
 db_auxiliary_field.description = Fault rupture auxiliary field spatial database

tests/fullscale/linearelasticity/faults-2d/zeroslipfn.timedb

@@ -0,0 +1,10 @@
+#TIME HISTORY ascii
+TimeHistory {
+  num-points = 5
+  time-units = year
+}
+ 0.0  0.0
+ 1.0  0.0
+ 4.0  0.0
+ 8.0  0.0
+10.0  0.0

tests/fullscale/linearelasticity/faults-3d-buried/Makefile.am

@@ -0,0 +1,46 @@
+# -*- Makefile -*-
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard, U.S. Geological Survey
+# Charles A. Williams, GNS Science
+# Matthew G. Knepley, University at Buffalo
+#
+# This code was developed as part of the Computational Infrastructure
+# for Geodynamics (http://geodynamics.org).
+#
+# Copyright (c) 2010-2022 University of California, Davis
+#
+# See LICENSE.md for license information.
+#
+# ----------------------------------------------------------------------
+
+TESTS = test_pylith.py
+
+dist_check_SCRIPTS = test_pylith.py
+
+dist_noinst_PYTHON = \
+	generate_gmsh.py \
+	meshes.py \
+	TestUniformSlip.py \
+	uniformslip_soln.py
+
+dist_noinst_DATA = \
+	mesh_tet.msh \
+	pylithapp.cfg \
+	uniformslip.cfg
+
+
+noinst_TMP =
+
+
+export_datadir = $(abs_builddir)
+include $(top_srcdir)/tests/data.am
+
+clean-local: clean-local-tmp clean-data
+.PHONY: clean-local-tmp
+clean-local-tmp:
+	$(RM) $(RM_FLAGS) -r output __pycache__
+
+
+# End of file

tests/fullscale/linearelasticity/faults-3d-buried/TestUniformSlip.py

@@ -0,0 +1,93 @@
+#!/usr/bin/env nemesis
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard, U.S. Geological Survey
+# Charles A. Williams, GNS Science
+# Matthew G. Knepley, University at Buffalo
+#
+# This code was developed as part of the Computational Infrastructure
+# for Geodynamics (http://geodynamics.org).
+#
+# Copyright (c) 2010-2022 University of California, Davis
+#
+# See LICENSE.md for license information.
+#
+# ----------------------------------------------------------------------
+
+import unittest
+
+from pylith.testing.FullTestApp import (FullTestCase, Check)
+
+import meshes
+import uniformslip_soln
+
+
+# -------------------------------------------------------------------------------------------------
+class TestCase(FullTestCase):
+
+    def setUp(self):
+        defaults = {
+            "filename": "output/{name}-{mesh_entity}.h5",
+            "exact_soln": uniformslip_soln.AnalyticalSoln(),
+            "mesh": self.mesh,
+        }
+        self.checks = [
+            Check(
+                mesh_entities=["domain"],
+                vertex_fields=[],
+                defaults=defaults,
+            ),
+            Check(
+                mesh_entities=["mat_elastic"],
+                filename="output/{name}-{mesh_entity}_info.h5",
+                cell_fields = ["density", "bulk_modulus", "shear_modulus"],
+                defaults=defaults,
+            ),
+            Check(
+                mesh_entities=["bc_xneg", "bc_xpos", "bc_yneg", "bc_ypos", "bc_zneg"],
+                filename="output/{name}-{mesh_entity}_info.h5",
+                vertex_fields=["initial_amplitude"],
+                defaults=defaults,
+            ),
+            Check(
+                mesh_entities=["fault"],
+                vertex_fields=["slip"],
+                defaults=defaults,
+            ),
+        ]
+
+    def run_pylith(self, testName, args, nprocs=1):
+        FullTestCase.run_pylith(self, testName, args, nprocs=nprocs)
+
+
+# -------------------------------------------------------------------------------------------------
+class TestTetGmsh(TestCase):
+
+    def setUp(self):
+        self.name = "uniformslip"
+        self.mesh = meshes.TetGmsh()
+        super().setUp()
+
+        TestCase.run_pylith(self, self.name, ["uniformslip.cfg"], nprocs=1)
+        return
+
+
+# -------------------------------------------------------------------------------------------------
+def test_cases():
+    return [
+        TestTetGmsh,
+    ]
+
+
+# -------------------------------------------------------------------------------------------------
+if __name__ == '__main__':
+    FullTestCase.parse_args()
+
+    suite = unittest.TestSuite()
+    for test in test_cases():
+        suite.addTest(unittest.makeSuite(test))
+    unittest.TextTestRunner(verbosity=2).run(suite)
+
+
+# End of file