[cig-commits] [commit] baagaard/dynrup-new-lagrange: Name change (TractPerturbation -> TractionPerturbation). (0736169)

cig_noreply at geodynamics.org cig_noreply at geodynamics.org
Wed Nov 5 15:44:38 PST 2014 

Repository : https://github.com/geodynamics/pylith

On branch  : baagaard/dynrup-new-lagrange
Link       : https://github.com/geodynamics/pylith/compare/f33c75b19fd60eedb2a3405db76a1fee333bb1d7...5b6d812b1612809fea3bd331c4e5af98c25a536a

>---------------------------------------------------------------

commit 0736169fc2892037b3242762ba7b83f014557a1c
Author: Brad Aagaard <baagaard at usgs.gov>
Date:   Tue Jul 1 12:52:50 2014 -0700

    Name change (TractPerturbation -> TractionPerturbation).
    
    Also small change to FaultCohesiveDyn (add Jacobian from
    FaultCohesiveLagrange); needs reimplementation.


>---------------------------------------------------------------

0736169fc2892037b3242762ba7b83f014557a1c
 libsrc/pylith/faults/FaultCohesiveDyn.cc           | 161 +++++++++++++++++++++
 libsrc/pylith/faults/FaultCohesiveDyn.hh           |  12 ++
 modulesrc/faults/FaultCohesiveDyn.i                |  28 ++--
 modulesrc/faults/faults.i                          |   4 +-
 pylith/Makefile.am                                 |   2 +-
 ...ractPerturbation.py => TractionPerturbation.py} |  16 +-
 unittests/pytests/faults/Makefile.am               |   2 +-
 ...Perturbation.py => TestTractionPerturbation.py} |  16 +-
 unittests/pytests/faults/testfaults.py             |   4 +-
 9 files changed, 210 insertions(+), 35 deletions(-)

diff --git a/libsrc/pylith/faults/FaultCohesiveDyn.cc b/libsrc/pylith/faults/FaultCohesiveDyn.cc
index 0784d50..00b3242 100644
--- a/libsrc/pylith/faults/FaultCohesiveDyn.cc
+++ b/libsrc/pylith/faults/FaultCohesiveDyn.cc
@@ -463,6 +463,167 @@ pylith::faults::FaultCohesiveDyn::integrateResidual(const topology::Field& resid
   PYLITH_METHOD_END;
 } // integrateResidual
 
+
+// ----------------------------------------------------------------------
+// Compute Jacobian matrix (A) associated with operator.
+void
+pylith::faults::FaultCohesiveDyn::integrateJacobian(topology::Jacobian* jacobian,
+						    const PylithScalar t,
+						    topology::SolutionFields* const fields)
+{ // integrateJacobian
+  PYLITH_METHOD_BEGIN;
+
+  assert(jacobian);
+  assert(fields);
+  assert(_fields);
+  assert(_logger);
+
+  const int setupEvent = _logger->eventId("FaIJ setup");
+  const int geometryEvent = _logger->eventId("FaIJ geometry");
+  const int computeEvent = _logger->eventId("FaIJ compute");
+  const int restrictEvent = _logger->eventId("FaIJ restrict");
+  const int updateEvent = _logger->eventId("FaIJ update");
+
+  _logger->eventBegin(setupEvent);
+
+  // Add constraint information to Jacobian matrix; Entries are
+  // associated with vertices ik, jk, ki, and kj.
+
+  // Get cell geometry information that doesn't depend on cell
+  const int spaceDim = _quadrature->spaceDim();
+
+  // Get fields.
+  topology::Field& area = _fields->get("area");
+  topology::VecVisitorMesh areaVisitor(area);
+  const PetscScalar* areaArray = areaVisitor.localArray();
+  
+  PetscDM solnDM = fields->solution().dmMesh();assert(solnDM);
+  PetscSection solnSection = fields->solution().petscSection();assert(solnSection);
+  PetscSection solnGlobalSection = NULL;
+  PetscErrorCode err = DMGetDefaultGlobalSection(solnDM, &solnGlobalSection);PYLITH_CHECK_ERROR(err);assert(solnGlobalSection);
+
+  // Get fault information
+  PetscDM dmMesh = fields->mesh().dmMesh();assert(dmMesh);
+
+  // Allocate vectors for vertex values
+  scalar_array jacobianVertex(spaceDim*spaceDim);
+  int_array indicesL(spaceDim);
+  int_array indicesN(spaceDim);
+  int_array indicesP(spaceDim);
+  int_array indicesRel(spaceDim);
+  for (int i=0; i < spaceDim; ++i) {
+    indicesRel[i] = i;
+  } // for
+
+  // Get sparse matrix
+  const PetscMat jacobianMatrix = jacobian->matrix();assert(jacobianMatrix);
+
+  _logger->eventEnd(setupEvent);
+#if !defined(DETAILED_EVENT_LOGGING)
+  _logger->eventBegin(computeEvent);
+#endif
+
+  const int numVertices = _cohesiveVertices.size();
+  for (int iVertex=0; iVertex < numVertices; ++iVertex) {
+    const int e_lagrange = _cohesiveVertices[iVertex].lagrange;
+    const int v_fault = _cohesiveVertices[iVertex].fault;
+    const int v_negative = _cohesiveVertices[iVertex].negative;
+    const int v_positive = _cohesiveVertices[iVertex].positive;
+
+    if (e_lagrange < 0) { // Skip clamped edges.
+      continue;
+    } // if
+
+    // Compute contribution only if Lagrange constraint is local.
+    PetscInt gloff = 0;
+    err = PetscSectionGetOffset(solnGlobalSection, e_lagrange, &gloff);PYLITH_CHECK_ERROR(err);
+    if (gloff < 0)
+      continue;
+
+    PetscInt gnoff = 0;
+    err = PetscSectionGetOffset(solnGlobalSection, v_negative, &gnoff);PYLITH_CHECK_ERROR(err);
+    gnoff = gnoff < 0 ? -(gnoff+1) : gnoff;
+
+    PetscInt gpoff = 0;
+    err = PetscSectionGetOffset(solnGlobalSection, v_positive, &gpoff);PYLITH_CHECK_ERROR(err);
+    gpoff = gpoff < 0 ? -(gpoff+1) : gpoff;
+
+#if defined(DETAILED_EVENT_LOGGING)
+    _logger->eventBegin(restrictEvent);
+#endif
+
+    // Get area associated with fault vertex.
+    const PetscInt aoff = areaVisitor.sectionOffset(v_fault);
+    assert(1 == areaVisitor.sectionDof(v_fault));
+
+    // Set global order indices
+    indicesL = indicesRel + gloff;
+    indicesN = indicesRel + gnoff;
+    indicesP = indicesRel + gpoff;
+    PetscInt cdof;
+    err = PetscSectionGetConstraintDof(solnSection, v_negative, &cdof);PYLITH_CHECK_ERROR(err);assert(0 == cdof);
+    err = PetscSectionGetConstraintDof(solnSection, v_positive, &cdof);PYLITH_CHECK_ERROR(err);assert(0 == cdof);
+
+#if defined(DETAILED_EVENT_LOGGING)
+    _logger->eventEnd(restrictEvent);
+    _logger->eventBegin(updateEvent);
+#endif
+
+    // Set diagonal entries of Jacobian at positive vertex to area
+    // associated with vertex.
+    for (int iDim=0; iDim < spaceDim; ++iDim)
+      jacobianVertex[iDim*spaceDim+iDim] = areaArray[aoff];
+
+    // Values at positive vertex, entry L,P in Jacobian
+    err = MatSetValues(jacobianMatrix, 
+		       indicesL.size(), &indicesL[0], 
+		       indicesP.size(), &indicesP[0], 
+		       &jacobianVertex[0], ADD_VALUES);PYLITH_CHECK_ERROR(err);
+
+    // Values at positive vertex, entry P,L in Jacobian
+    err = MatSetValues(jacobianMatrix, 
+		       indicesP.size(), &indicesP[0], 
+		       indicesL.size(), &indicesL[0], 
+		       &jacobianVertex[0], ADD_VALUES);PYLITH_CHECK_ERROR(err);
+    
+    // Values at negative vertex, entry L,N in Jacobian
+    jacobianVertex *= -1.0;
+    err = MatSetValues(jacobianMatrix,
+		       indicesL.size(), &indicesL[0],
+		       indicesN.size(), &indicesN[0],
+		       &jacobianVertex[0], ADD_VALUES);PYLITH_CHECK_ERROR(err);
+
+    // Values at negative vertex, entry N,L in Jacobian
+    err = MatSetValues(jacobianMatrix,
+		       indicesN.size(), &indicesN[0],
+		       indicesL.size(), &indicesL[0],
+		       &jacobianVertex[0], ADD_VALUES);PYLITH_CHECK_ERROR(err);
+
+    // Values at Lagrange vertex, entry L,L in Jacobian
+    // We must have entries on the diagonal.
+    jacobianVertex = 0.0;
+    err = MatSetValues(jacobianMatrix,
+		       indicesL.size(), &indicesL[0],
+		       indicesL.size(), &indicesL[0],
+		       &jacobianVertex[0], ADD_VALUES);PYLITH_CHECK_ERROR(err);
+
+#if defined(DETAILED_EVENT_LOGGING)
+    _logger->eventEnd(updateEvent);
+#endif
+
+  } // for
+  PetscLogFlops(numVertices*spaceDim*2);
+
+#if !defined(DETAILED_EVENT_LOGGING)
+  _logger->eventEnd(computeEvent);
+#endif
+
+  _needNewJacobian = false;
+
+  PYLITH_METHOD_END;
+} // integrateJacobian
+
+
 // ----------------------------------------------------------------------
 // Update state variables as needed.
 void
diff --git a/libsrc/pylith/faults/FaultCohesiveDyn.hh b/libsrc/pylith/faults/FaultCohesiveDyn.hh
index b0e9e7d..eefb156 100644
--- a/libsrc/pylith/faults/FaultCohesiveDyn.hh
+++ b/libsrc/pylith/faults/FaultCohesiveDyn.hh
@@ -109,6 +109,18 @@ public :
 			 const PylithScalar t,
 			 topology::SolutionFields* const fields);
 
+  /** Integrate contributions to Jacobian matrix (A) associated with
+   * operator that require assembly across processors.
+   *
+   * @param jacobian Sparse matrix
+   * @param t Current time
+   * @param fields Solution fields
+   */
+  virtual
+  void integrateJacobian(topology::Jacobian* jacobian,
+			 const PylithScalar t,
+			 topology::SolutionFields* const fields);
+
   /** Update state variables as needed.
    *
    * @param t Current time
diff --git a/modulesrc/faults/FaultCohesiveDyn.i b/modulesrc/faults/FaultCohesiveDyn.i
index 0d639c7..8333fc2 100644
--- a/modulesrc/faults/FaultCohesiveDyn.i
+++ b/modulesrc/faults/FaultCohesiveDyn.i
@@ -45,7 +45,7 @@ namespace pylith {
        *
        * @param tract Spatial and temporal variation of tractions.
        */
-      void tractionPerturbation(TractPerturbation* tract);
+      void tractionPerturbation(TractionPerturbation* tract);
   
       /** Set the friction (constitutive) model.
        *
@@ -92,8 +92,20 @@ namespace pylith {
        */
       virtual
       void integrateResidual(const pylith::topology::Field& residual,
-				      const PylithScalar t,
-				      pylith::topology::SolutionFields* const fields);
+			     const PylithScalar t,
+			     pylith::topology::SolutionFields* const fields);
+
+      /** Integrate contributions to Jacobian matrix (A) associated with
+       * operator that require assembly across processors.
+       *
+       * @param jacobian Sparse matrix
+       * @param t Current time
+       * @param fields Solution fields
+       */
+      virtual
+      void integrateJacobian(pylith::topology::Jacobian* jacobian,
+			     const PylithScalar t,
+			     pylith::topology::SolutionFields* const fields);
 
       /** Update state variables as needed.
        *
@@ -104,16 +116,6 @@ namespace pylith {
       void updateStateVars(const PylithScalar t,
 			   pylith::topology::SolutionFields* const fields);
       
-      /** Constrain solution space based on friction.
-       *
-       * @param fields Solution fields.
-       * @param t Current time.
-       * @param jacobian Sparse matrix for system Jacobian.
-       */
-      void constrainSolnSpace(pylith::topology::SolutionFields* const fields,
-			      const PylithScalar t,
-			      const pylith::topology::Jacobian& jacobian);
-
       /** Adjust solution from solver with lumped Jacobian to match Lagrange
        *  multiplier constraints.
        *
diff --git a/modulesrc/faults/faults.i b/modulesrc/faults/faults.i
index 098d350..d078f54 100644
--- a/modulesrc/faults/faults.i
+++ b/modulesrc/faults/faults.i
@@ -28,7 +28,7 @@
 #include "pylith/faults/LiuCosSlipFn.hh"
 #include "pylith/faults/TimeHistorySlipFn.hh"
 #include "pylith/faults/EqKinSrc.hh"
-#include "pylith/faults/TractPerturbation.hh"
+#include "pylith/faults/TractionPerturbation.hh"
 #include "pylith/faults/Fault.hh"
 #include "pylith/faults/FaultCohesive.hh"
 #include "pylith/faults/FaultCohesiveLagrange.hh"
@@ -72,7 +72,7 @@ import_array();
 %include "LiuCosSlipFn.i"
 %include "TimeHistorySlipFn.i"
 %include "EqKinSrc.i"
-%include "TractPerturbation.i"
+%include "TractionPerturbation.i"
 %include "Fault.i"
 %include "FaultCohesive.i"
 %include "FaultCohesiveLagrange.i"
diff --git a/pylith/Makefile.am b/pylith/Makefile.am
index bf0ddcc..4ae0759 100644
--- a/pylith/Makefile.am
+++ b/pylith/Makefile.am
@@ -41,7 +41,7 @@ nobase_pkgpyexec_PYTHON = \
 	faults/SlipTimeFn.py \
 	faults/StepSlipFn.py \
 	faults/TimeHistorySlipFn.py \
-	faults/TractPerturbation.py \
+	faults/TractionPerturbation.py \
 	faults/Fault.py \
 	faults/FaultCohesive.py \
 	faults/FaultCohesiveKin.py \
diff --git a/pylith/faults/TractPerturbation.py b/pylith/faults/TractionPerturbation.py
similarity index 84%
rename from pylith/faults/TractPerturbation.py
rename to pylith/faults/TractionPerturbation.py
index dd2b595..c935692 100644
--- a/pylith/faults/TractPerturbation.py
+++ b/pylith/faults/TractionPerturbation.py
@@ -16,24 +16,24 @@
 # ----------------------------------------------------------------------
 #
 
-## @file pylith/faults/TractPerturbation.py
+## @file pylith/faults/TractionPerturbation.py
 ##
 
 ## @brief Python object for managing parameters for a kinematic
 ## earthquake sources.
 ##
-## TractPerturbation is responsible for providing the value of
+## TractionPerturbation is responsible for providing the value of
 ## specified traction at time t over a fault surface.
 ##
 ## Factory: eq_kinematic_src
 
 from pylith.bc.TimeDependent import TimeDependent
-from faults import TractPerturbation as ModuleTractPerturbation
+from faults import TractionPerturbation as ModuleTractionPerturbation
 
 from pylith.utils.NullComponent import NullComponent
 
-# TractPerturbation class
-class TractPerturbation(TimeDependent, ModuleTractPerturbation):
+# TractionPerturbation class
+class TractionPerturbation(TimeDependent, ModuleTractionPerturbation):
   """
   Python object for managing specified tractions on a fault surface.
 
@@ -78,7 +78,7 @@ class TractPerturbation(TimeDependent, ModuleTractPerturbation):
     """
     Create handle to corresponding C++ object.
     """
-    ModuleTractPerturbation.__init__(self)
+    ModuleTractionPerturbation.__init__(self)
     return
   
 
@@ -86,9 +86,9 @@ class TractPerturbation(TimeDependent, ModuleTractPerturbation):
 
 def traction_perturbation():
   """
-  Factory associated with TractPerturbation.
+  Factory associated with TractionPerturbation.
   """
-  return TractPerturbation()
+  return TractionPerturbation()
 
 
 # End of file 
diff --git a/unittests/pytests/faults/Makefile.am b/unittests/pytests/faults/Makefile.am
index 67e3676..94e8000 100644
--- a/unittests/pytests/faults/Makefile.am
+++ b/unittests/pytests/faults/Makefile.am
@@ -32,7 +32,7 @@ noinst_PYTHON = \
 	TestEqKinSrc.py \
 	TestFaultCohesiveKin.py \
 	TestFaultCohesiveDyn.py \
-	TestTractPerturbation.py \
+	TestTractionPerturbation.py \
 	TestFaultCohesiveImpulses.py \
 	TestLiuCosSlipFn.py \
 	TestSingleRupture.py \
diff --git a/unittests/pytests/faults/TestTractPerturbation.py b/unittests/pytests/faults/TestTractionPerturbation.py
similarity index 80%
rename from unittests/pytests/faults/TestTractPerturbation.py
rename to unittests/pytests/faults/TestTractionPerturbation.py
index f6a3032..2b1ac62 100644
--- a/unittests/pytests/faults/TestTractPerturbation.py
+++ b/unittests/pytests/faults/TestTractionPerturbation.py
@@ -16,25 +16,25 @@
 # ======================================================================
 #
 
-## @file unittests/pytests/faults/TestTractPerturbation.py
+## @file unittests/pytests/faults/TestTractionPerturbation.py
 
-## @brief Unit testing of TractPerturbation object.
+## @brief Unit testing of TractionPerturbation object.
 
 import unittest
 
-from pylith.faults.TractPerturbation import TractPerturbation
+from pylith.faults.TractionPerturbation import TractionPerturbation
 
 # ----------------------------------------------------------------------
-class TestTractPerturbation(unittest.TestCase):
+class TestTractionPerturbation(unittest.TestCase):
   """
-  Unit testing of TractPerturbation object.
+  Unit testing of TractionPerturbation object.
   """
 
   def test_constructor(self):
     """
     Test constructor.
     """
-    tract = TractPerturbation()
+    tract = TractionPerturbation()
     return
 
 
@@ -62,7 +62,7 @@ class TestTractPerturbation(unittest.TestCase):
     dbChange.inventory.label = "traction change"
     dbChange._configure()
     
-    tract = TractPerturbation()
+    tract = TractionPerturbation()
     tract.inventory.dbInitial = dbInitial
     tract.inventory.dbChange = dbChange
     tract._configure()
@@ -73,7 +73,7 @@ class TestTractPerturbation(unittest.TestCase):
     """
     Test factory method.
     """
-    from pylith.faults.TractPerturbation import traction_perturbation
+    from pylith.faults.TractionPerturbation import traction_perturbation
     fn = traction_perturbation()
     return
 
diff --git a/unittests/pytests/faults/testfaults.py b/unittests/pytests/faults/testfaults.py
index 4add3b6..b80aa4e 100644
--- a/unittests/pytests/faults/testfaults.py
+++ b/unittests/pytests/faults/testfaults.py
@@ -81,8 +81,8 @@ class TestApp(Script):
     from TestEqKinSrc import TestEqKinSrc
     suite.addTest(unittest.makeSuite(TestEqKinSrc))
 
-    from TestTractPerturbation import TestTractPerturbation
-    suite.addTest(unittest.makeSuite(TestTractPerturbation))
+    from TestTractionPerturbation import TestTractionPerturbation
+    suite.addTest(unittest.makeSuite(TestTractionPerturbation))
 
     from TestFaultCohesiveKin import TestFaultCohesiveKin
     suite.addTest(unittest.makeSuite(TestFaultCohesiveKin))

More information about the CIG-COMMITS mailing list