[cig-commits] r12984 - in short/3D/PyLith/trunk: libsrc/faults libsrc/feassemble modulesrc/bc modulesrc/faults modulesrc/feassemble pylith/feassemble pylith/problems unittests/libtests/faults
brad at geodynamics.org
brad at geodynamics.org
Thu Oct 2 17:18:00 PDT 2008
Author: brad
Date: 2008-10-02 17:18:00 -0700 (Thu, 02 Oct 2008)
New Revision: 12984
Modified:
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh
short/3D/PyLith/trunk/libsrc/feassemble/Integrator.icc
short/3D/PyLith/trunk/modulesrc/bc/bc.pyxe.src
short/3D/PyLith/trunk/modulesrc/faults/faults.pyxe.src
short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.pyxe.src
short/3D/PyLith/trunk/pylith/feassemble/Integrator.py
short/3D/PyLith/trunk/pylith/problems/Formulation.py
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.hh
Log:
Added integrateResisualAssembled() and integrateJacobianAssembled() to Integrator to account for cases where terms do not need assembly (e.g., faults using Lagrange multiplier constraints). This should fix the bug for the fault implementation when running in parallel.
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh 2008-10-02 01:07:36 UTC (rev 12983)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh 2008-10-03 00:18:00 UTC (rev 12984)
@@ -109,31 +109,33 @@
const double_array& normalDir,
spatialdata::spatialdb::SpatialDB* matDB);
- /** Integrate contributions to residual term (r) for operator.
+ /** Integrate contributions to residual term (r) for operator that
+ * do not require assembly across cells, vertices, or processors.
*
* @param residual Field containing values for residual
* @param t Current time
* @param fields Solution fields
* @param mesh Finite-element mesh
*/
- void integrateResidual(const ALE::Obj<real_section_type>& residual,
- const double t,
- topology::FieldsManager* const fields,
- const ALE::Obj<Mesh>& mesh,
- const spatialdata::geocoords::CoordSys* cs);
+ void integrateResidualAssembled(const ALE::Obj<real_section_type>& residual,
+ const double t,
+ topology::FieldsManager* const fields,
+ const ALE::Obj<Mesh>& mesh,
+ const spatialdata::geocoords::CoordSys* cs);
/** Integrate contributions to Jacobian matrix (A) associated with
- * operator.
+ * operator that do not require assembly across cells, vertices, or
+ * processors.
*
* @param mat Sparse matrix
* @param t Current time
* @param fields Solution fields
* @param mesh Finite-element mesh
*/
- void integrateJacobian(PetscMat* mat,
- const double t,
- topology::FieldsManager* const fields,
- const ALE::Obj<Mesh>& mesh);
+ void integrateJacobianAssembled(PetscMat* mat,
+ const double t,
+ topology::FieldsManager* const fields,
+ const ALE::Obj<Mesh>& mesh);
/** Verify configuration is acceptable.
*
@@ -194,11 +196,6 @@
void _calcOrientation(const double_array& upDir,
const double_array& normalDir);
- /** Calculate pairing between fault vertices and first cell they
- * appear in to prevent double counting in integrating Jacobian.
- */
- void _calcVertexCellPairs(void);
-
/** Calculate conditioning field.
*
* @param cs Coordinate system for mesh
@@ -264,11 +261,6 @@
/// Field over the fault mesh vertices of vector field of cumulative slip.
ALE::Obj<real_section_type> _cumSlip;
- /// Label of cell used to compute Jacobian for each fault vertex (must
- /// prevent overlap so that only 1 cell will contribute for
- /// each vertex).
- ALE::Obj<int_section_type> _faultVertexCell;
-
std::map<Mesh::point_type, Mesh::point_type> _cohesiveToFault;
/// Scalar field for vertex information over fault mesh.
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh 2008-10-02 01:07:36 UTC (rev 12983)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh 2008-10-03 00:18:00 UTC (rev 12984)
@@ -120,7 +120,7 @@
const double t,
topology::FieldsManager* const fields,
const ALE::Obj<Mesh>& mesh,
- const spatialdata::geocoords::CoordSys* cs) = 0;
+ const spatialdata::geocoords::CoordSys* cs);
/** Integrate contributions to Jacobian matrix (A) associated with
* operator.
@@ -130,12 +130,43 @@
* @param fields Solution fields
* @param mesh Finite-element mesh
*/
- virtual
+ virtual
void integrateJacobian(PetscMat* mat,
const double t,
topology::FieldsManager* const fields,
- const ALE::Obj<Mesh>& mesh) = 0;
+ const ALE::Obj<Mesh>& mesh);
+ /** Integrate contributions to residual term (r) for operator that
+ * do not require assembly over cells, vertices, or processors.
+ *
+ * @param residual Field containing values for residual
+ * @param t Current time
+ * @param fields Solution fields
+ * @param mesh Finite-element mesh
+ * @param cs Mesh coordinate system
+ */
+ virtual
+ void integrateResidualAssembled(const ALE::Obj<real_section_type>& residual,
+ const double t,
+ topology::FieldsManager* const fields,
+ const ALE::Obj<Mesh>& mesh,
+ const spatialdata::geocoords::CoordSys* cs);
+
+ /** Integrate contributions to Jacobian matrix (A) associated with
+ * operator that do not require assembly over cells, vertices, or
+ * processors
+ *
+ * @param mat Sparse matrix
+ * @param t Current time
+ * @param fields Solution fields
+ * @param mesh Finite-element mesh
+ */
+ virtual
+ void integrateJacobianAssembled(PetscMat* mat,
+ const double t,
+ topology::FieldsManager* const fields,
+ const ALE::Obj<Mesh>& mesh);
+
/** Update state variables as needed.
*
* @param t Current time
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Integrator.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Integrator.icc 2008-10-02 01:07:36 UTC (rev 12983)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Integrator.icc 2008-10-03 00:18:00 UTC (rev 12984)
@@ -35,6 +35,52 @@
_useSolnIncr = flag;
} // useSolnIncr
+// Integrate contributions to residual term (r) for operator.
+inline
+void
+pylith::feassemble::Integrator::integrateResidual(
+ const ALE::Obj<real_section_type>& residual,
+ const double t,
+ topology::FieldsManager* const fields,
+ const ALE::Obj<Mesh>& mesh,
+ const spatialdata::geocoords::CoordSys* cs) {
+} // integrateResidual
+
+// Integrate contributions to Jacobian matrix (A) associated with
+// operator.
+inline
+void
+pylith::feassemble::Integrator::integrateJacobian(
+ PetscMat* mat,
+ const double t,
+ topology::FieldsManager* const fields,
+ const ALE::Obj<Mesh>& mesh) {
+} // integrateJacobian
+
+// Integrate contributions to residual term (r) for operator that
+// do not require assembly over cells, vertices, or processors.
+inline
+void
+pylith::feassemble::Integrator::integrateResidualAssembled(
+ const ALE::Obj<real_section_type>& residual,
+ const double t,
+ topology::FieldsManager* const fields,
+ const ALE::Obj<Mesh>& mesh,
+ const spatialdata::geocoords::CoordSys* cs) {
+} // integrateResidualAssembled
+
+// Integrate contributions to Jacobian matrix (A) associated with
+// operator that do not require assembly over cells, vertices, or
+// processors
+inline
+void
+pylith::feassemble::Integrator::integrateJacobianAssembled(
+ PetscMat* mat,
+ const double t,
+ topology::FieldsManager* const fields,
+ const ALE::Obj<Mesh>& mesh) {
+} // integrateJacobianAssembled
+
// Update state variables as needed.
inline
void
Modified: short/3D/PyLith/trunk/modulesrc/bc/bc.pyxe.src
===================================================================
--- short/3D/PyLith/trunk/modulesrc/bc/bc.pyxe.src 2008-10-02 01:07:36 UTC (rev 12983)
+++ short/3D/PyLith/trunk/modulesrc/bc/bc.pyxe.src 2008-10-03 00:18:00 UTC (rev 12984)
@@ -902,6 +902,94 @@
return
+ def integrateResidualAssembled(self, residual, t, fields, mesh, cs):
+ """
+ Integrate contributions to residual term (r) for operator that do
+ not require assembly over cells, vertices, or processors.
+ """
+ # create shim for method 'integrateResidualAssembled'
+ #embed{ void AbsorbingDampers_integrateResidualAssembled(void* objVptr, void* residualVptr, double t, void* fieldsVptr, void* meshVptr, void* csVptr)
+ try {
+ assert(0 != objVptr);
+ assert(0 != residualVptr);
+ assert(0 != fieldsVptr);
+ assert(0 != meshVptr);
+ assert(0 != csVptr);
+ ALE::Obj<pylith::Mesh>* mesh =
+ (ALE::Obj<pylith::Mesh>*) meshVptr;
+ ALE::Obj<pylith::real_section_type>* residual =
+ (ALE::Obj<pylith::real_section_type>*) residualVptr;
+ pylith::topology::FieldsManager* fields =
+ (pylith::topology::FieldsManager*) fieldsVptr;
+ spatialdata::geocoords::CoordSys* cs =
+ (spatialdata::geocoords::CoordSys*) csVptr;
+ ((pylith::bc::AbsorbingDampers*) objVptr)->integrateResidualAssembled(*residual,
+ t, fields, *mesh, cs);
+ } catch (const std::exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.what()));
+ } catch (const ALE::Exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.msg().c_str()));
+ } catch (...) {
+ PyErr_SetString(PyExc_RuntimeError,
+ "Caught unknown C++ exception.");
+ } // try/catch
+ #}embed
+ if mesh.name != "pylith_topology_Mesh":
+ raise TypeError, \
+ "Argument must be extension module type 'Mesh'."
+ AbsorbingDampers_integrateResidualAssembled(self.thisptr,
+ PyCObject_AsVoidPtr(residual),
+ t,
+ ptrFromHandle(fields),
+ ptrFromHandle(mesh),
+ ptrFromHandle(cs))
+ return
+
+
+ def integrateJacobianAssembled(self, mat, t, fields, mesh):
+ """
+ Compute contributions to Jacobian matrix (A) associated with
+ operator that do not require assembly over cells, vertices, or
+ processors.
+ """
+ # create shim for method 'integrateJacobianAssembled'
+ #embed{ void AbsorbingDampers_integrateJacobianAssembled(void* objVptr, void* matVptr, double t, void* fieldsVptr, void* meshVptr)
+ try {
+ assert(0 != objVptr);
+ assert(0 != matVptr);
+ assert(0 != fieldsVptr);
+ assert(0 != meshVptr);
+ ALE::Obj<pylith::Mesh>* mesh =
+ (ALE::Obj<pylith::Mesh>*) meshVptr;
+ PetscMat* mat = (PetscMat*) matVptr;
+ pylith::topology::FieldsManager* fields =
+ (pylith::topology::FieldsManager*) fieldsVptr;
+ ((pylith::bc::AbsorbingDampers*) objVptr)->integrateJacobianAssembled(
+ mat, t, fields, *mesh);
+ } catch (const std::exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.what()));
+ } catch (const ALE::Exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.msg().c_str()));
+ } catch (...) {
+ PyErr_SetString(PyExc_RuntimeError,
+ "Caught unknown C++ exception.");
+ } // try/catch
+ #}embed
+ if mesh.name != "pylith_topology_Mesh":
+ raise TypeError, \
+ "Argument 'mesh' must be extension module type 'Mesh'."
+ AbsorbingDampers_integrateJacobianAssembled(self.thisptr,
+ PyCObject_AsVoidPtr(mat),
+ t,
+ ptrFromHandle(fields),
+ ptrFromHandle(mesh))
+ return
+
+
def updateState(self, t, fields, mesh):
"""
Update state variables as needed.
@@ -1200,6 +1288,94 @@
return
+ def integrateResidualAssembled(self, residual, t, fields, mesh, cs):
+ """
+ Integrate contributions to residual term (r) for operator that do
+ not require assembly over cells, vertices, or processors.
+ """
+ # create shim for method 'integrateResidualAssembled'
+ #embed{ void Neumann_integrateResidualAssembled(void* objVptr, void* residualVptr, double t, void* fieldsVptr, void* meshVptr, void* csVptr)
+ try {
+ assert(0 != objVptr);
+ assert(0 != residualVptr);
+ assert(0 != fieldsVptr);
+ assert(0 != meshVptr);
+ assert(0 != csVptr);
+ ALE::Obj<pylith::Mesh>* mesh =
+ (ALE::Obj<pylith::Mesh>*) meshVptr;
+ ALE::Obj<pylith::real_section_type>* residual =
+ (ALE::Obj<pylith::real_section_type>*) residualVptr;
+ pylith::topology::FieldsManager* fields =
+ (pylith::topology::FieldsManager*) fieldsVptr;
+ spatialdata::geocoords::CoordSys* cs =
+ (spatialdata::geocoords::CoordSys*) csVptr;
+ ((pylith::bc::Neumann*) objVptr)->integrateResidualAssembled(*residual,
+ t, fields, *mesh, cs);
+ } catch (const std::exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.what()));
+ } catch (const ALE::Exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.msg().c_str()));
+ } catch (...) {
+ PyErr_SetString(PyExc_RuntimeError,
+ "Caught unknown C++ exception.");
+ } // try/catch
+ #}embed
+ if mesh.name != "pylith_topology_Mesh":
+ raise TypeError, \
+ "Argument must be extension module type 'Mesh'."
+ Neumann_integrateResidualAssembled(self.thisptr,
+ PyCObject_AsVoidPtr(residual),
+ t,
+ ptrFromHandle(fields),
+ ptrFromHandle(mesh),
+ ptrFromHandle(cs))
+ return
+
+
+ def integrateJacobianAssembled(self, mat, t, fields, mesh):
+ """
+ Compute contributions to Jacobian matrix (A) associated with
+ operator that do not require assembly over cells, vertices, or
+ processors.
+ """
+ # create shim for method 'integrateJacobianAssembled'
+ #embed{ void Neumann_integrateJacobianAssembled(void* objVptr, void* matVptr, double t, void* fieldsVptr, void* meshVptr)
+ try {
+ assert(0 != objVptr);
+ assert(0 != matVptr);
+ assert(0 != fieldsVptr);
+ assert(0 != meshVptr);
+ ALE::Obj<pylith::Mesh>* mesh =
+ (ALE::Obj<pylith::Mesh>*) meshVptr;
+ PetscMat* mat = (PetscMat*) matVptr;
+ pylith::topology::FieldsManager* fields =
+ (pylith::topology::FieldsManager*) fieldsVptr;
+ ((pylith::bc::Neumann*) objVptr)->integrateJacobianAssembled(
+ mat, t, fields, *mesh);
+ } catch (const std::exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.what()));
+ } catch (const ALE::Exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.msg().c_str()));
+ } catch (...) {
+ PyErr_SetString(PyExc_RuntimeError,
+ "Caught unknown C++ exception.");
+ } // try/catch
+ #}embed
+ if mesh.name != "pylith_topology_Mesh":
+ raise TypeError, \
+ "Argument 'mesh' must be extension module type 'Mesh'."
+ Neumann_integrateJacobianAssembled(self.thisptr,
+ PyCObject_AsVoidPtr(mat),
+ t,
+ ptrFromHandle(fields),
+ ptrFromHandle(mesh))
+ return
+
+
def updateState(self, t, fields, mesh):
"""
Update state variables as needed.
Modified: short/3D/PyLith/trunk/modulesrc/faults/faults.pyxe.src
===================================================================
--- short/3D/PyLith/trunk/modulesrc/faults/faults.pyxe.src 2008-10-02 01:07:36 UTC (rev 12983)
+++ short/3D/PyLith/trunk/modulesrc/faults/faults.pyxe.src 2008-10-03 00:18:00 UTC (rev 12984)
@@ -556,6 +556,94 @@
return
+ def integrateResidualAssembled(self, residual, t, fields, mesh, cs):
+ """
+ Integrate contributions to residual term (r) for operator that do
+ not require assembly over cells, vertices, or processors.
+ """
+ # create shim for method 'integrateResidualAssembled'
+ #embed{ void FaultCohesiveKin_integrateResidualAssembled(void* objVptr, void* residualVptr, double t, void* fieldsVptr, void* meshVptr, void* csVptr)
+ try {
+ assert(0 != objVptr);
+ assert(0 != residualVptr);
+ assert(0 != fieldsVptr);
+ assert(0 != meshVptr);
+ assert(0 != csVptr);
+ ALE::Obj<pylith::Mesh>* mesh =
+ (ALE::Obj<pylith::Mesh>*) meshVptr;
+ ALE::Obj<pylith::real_section_type>* residual =
+ (ALE::Obj<pylith::real_section_type>*) residualVptr;
+ pylith::topology::FieldsManager* fields =
+ (pylith::topology::FieldsManager*) fieldsVptr;
+ spatialdata::geocoords::CoordSys* cs =
+ (spatialdata::geocoords::CoordSys*) csVptr;
+ ((pylith::faults::FaultCohesiveKin*) objVptr)->integrateResidualAssembled(*residual,
+ t, fields, *mesh, cs);
+ } catch (const std::exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.what()));
+ } catch (const ALE::Exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.msg().c_str()));
+ } catch (...) {
+ PyErr_SetString(PyExc_RuntimeError,
+ "Caught unknown C++ exception.");
+ } // try/catch
+ #}embed
+ if mesh.name != "pylith_topology_Mesh":
+ raise TypeError, \
+ "Argument must be extension module type 'Mesh'."
+ FaultCohesiveKin_integrateResidualAssembled(self.thisptr,
+ PyCObject_AsVoidPtr(residual),
+ t,
+ ptrFromHandle(fields),
+ ptrFromHandle(mesh),
+ ptrFromHandle(cs))
+ return
+
+
+ def integrateJacobianAssenbked(self, mat, t, fields, mesh):
+ """
+ Compute contributions to Jacobian matrix (A) associated with
+ operator that do not require assembly over cells, vertices, or
+ processors.
+ """
+ # create shim for method 'integrateJacobianAssembled'
+ #embed{ void FaultCohesiveKin_integrateJacobianAssembled(void* objVptr, void* matVptr, double t, void* fieldsVptr, void* meshVptr)
+ try {
+ assert(0 != objVptr);
+ assert(0 != matVptr);
+ assert(0 != fieldsVptr);
+ assert(0 != meshVptr);
+ ALE::Obj<pylith::Mesh>* mesh =
+ (ALE::Obj<pylith::Mesh>*) meshVptr;
+ PetscMat* mat = (PetscMat*) matVptr;
+ pylith::topology::FieldsManager* fields =
+ (pylith::topology::FieldsManager*) fieldsVptr;
+ ((pylith::faults::FaultCohesiveKin*) objVptr)->integrateJacobianAssembled(
+ mat, t, fields, *mesh);
+ } catch (const std::exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.what()));
+ } catch (const ALE::Exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.msg().c_str()));
+ } catch (...) {
+ PyErr_SetString(PyExc_RuntimeError,
+ "Caught unknown C++ exception.");
+ } // try/catch
+ #}embed
+ if mesh.name != "pylith_topology_Mesh":
+ raise TypeError, \
+ "Argument 'mesh' must be extension module type 'Mesh'."
+ FaultCohesiveKin_integrateJacobianAssembled(self.thisptr,
+ PyCObject_AsVoidPtr(mat),
+ t,
+ ptrFromHandle(fields),
+ ptrFromHandle(mesh))
+ return
+
+
def updateState(self, t, fields, mesh):
"""
Update state variables as needed.
Modified: short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.pyxe.src
===================================================================
--- short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.pyxe.src 2008-10-02 01:07:36 UTC (rev 12983)
+++ short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.pyxe.src 2008-10-03 00:18:00 UTC (rev 12984)
@@ -1012,6 +1012,98 @@
return
+ def integrateResidualAssembled(self, residual, t, fields, mesh, cs):
+ """
+ Integrate contributions to residual term (r) for operator that do
+ not require assembly over cells, vertices, or processors.
+ """
+ # create shim for method 'integrateResidualAssembled'
+ #embed{ void Integrator_integrateResidualAssembled(void* objVptr, void* residualVptr, double t, void* fieldsVptr, void* meshVptr, void* csVptr)
+ try {
+ assert(0 != objVptr);
+ assert(0 != residualVptr);
+ assert(0 != fieldsVptr);
+ assert(0 != meshVptr);
+ assert(0 != csVptr);
+ ALE::Obj<pylith::Mesh>* mesh =
+ (ALE::Obj<pylith::Mesh>*) meshVptr;
+ ALE::Obj<pylith::real_section_type>* residual =
+ (ALE::Obj<pylith::real_section_type>*) residualVptr;
+ pylith::topology::FieldsManager* fields =
+ (pylith::topology::FieldsManager*) fieldsVptr;
+ spatialdata::geocoords::CoordSys* cs =
+ (spatialdata::geocoords::CoordSys*) csVptr;
+ ((pylith::feassemble::Integrator*) objVptr)->integrateResidualAssembled(*residual,
+ t, fields, *mesh,
+ cs);
+ } catch (const std::exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.what()));
+ } catch (const ALE::Exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.msg().c_str()));
+ } catch (...) {
+ PyErr_SetString(PyExc_RuntimeError,
+ "Caught unknown C++ exception.");
+ } // try/catch
+ #}embed
+ if mesh.name != "pylith_topology_Mesh":
+ raise TypeError, \
+ "Argument must be extension module type 'Mesh'."
+ Integrator_integrateResidualAssembled(self.thisptr,
+ PyCObject_AsVoidPtr(residual),
+ t,
+ ptrFromHandle(fields),
+ ptrFromHandle(mesh),
+ ptrFromHandle(cs))
+ return
+
+
+ def integrateJacobianAssembled(self, mat, t, fields, mesh):
+ """
+ Compute contributions to Jacobian matrix (A) associated with
+ operator that do not require assembly over cells, vertices, or
+ processors.
+ """
+ # create shim for method 'integrateJacobianAssembled'
+ #embed{ void Integrator_integrateJacobianAssembled(void* objVptr, void* matVptr, double t, void* fieldsVptr, void* meshVptr)
+ try {
+ assert(0 != objVptr);
+ assert(0 != matVptr);
+ assert(0 != fieldsVptr);
+ assert(0 != meshVptr);
+ ALE::Obj<pylith::Mesh>* mesh =
+ (ALE::Obj<pylith::Mesh>*) meshVptr;
+ PetscMat* mat = (PetscMat*) matVptr;
+ pylith::topology::FieldsManager* fields =
+ (pylith::topology::FieldsManager*) fieldsVptr;
+ ((pylith::feassemble::Integrator*) objVptr)->integrateJacobianAssembled(
+ mat, t, fields, *mesh);
+ } catch (const std::exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.what()));
+ } catch (const ALE::Exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.msg().c_str()));
+ } catch (...) {
+ PyErr_SetString(PyExc_RuntimeError,
+ "Caught unknown C++ exception.");
+ } // try/catch
+ #}embed
+ if fields.name != "pylith_topology_FieldsManager":
+ raise TypeError, \
+ "Argument 'fields' must be extension module type 'FieldsManager'."
+ if mesh.name != "pylith_topology_Mesh":
+ raise TypeError, \
+ "Argument 'mesh' must be extension module type 'Mesh'."
+ Integrator_integrateJacobianAssembled(self.thisptr,
+ PyCObject_AsVoidPtr(mat),
+ t,
+ ptrFromHandle(fields),
+ ptrFromHandle(mesh))
+ return
+
+
def updateState(self, t, fields, mesh):
"""
Update state variables as needed.
Modified: short/3D/PyLith/trunk/pylith/feassemble/Integrator.py
===================================================================
--- short/3D/PyLith/trunk/pylith/feassemble/Integrator.py 2008-10-02 01:07:36 UTC (rev 12983)
+++ short/3D/PyLith/trunk/pylith/feassemble/Integrator.py 2008-10-03 00:18:00 UTC (rev 12984)
@@ -28,6 +28,8 @@
"useSolnIncr",
"integrateResidual",
"integrateJacobian",
+ "integrateResidualAssembled",
+ "integrateJacobianAssembled",
"preinitialize",
"verifyConfiguration",
"initialize",
@@ -137,6 +139,20 @@
return
+ def needNewJacobian(self):
+ """
+ Returns true if we need to recompute Jacobian matrix for operator,
+ false otherwise.
+ """
+ logEvent = "%snewJacobian" % self._loggingPrefix
+ self._logger.eventBegin(logEvent)
+
+ assert(None != self.cppHandle)
+ flag = self.cppHandle.needNewJacobian
+ self._logger.eventEnd(logEvent)
+ return flag
+
+
def integrateResidual(self, residual, t, fields):
"""
Integrate contributions to residual term at time t.
@@ -152,34 +168,53 @@
return
- def needNewJacobian(self):
+ def integrateJacobian(self, jacobian, t, fields):
"""
- Returns true if we need to recompute Jacobian matrix for operator,
- false otherwise.
+ Integrate contributions to Jacobian term at time t.
"""
- logEvent = "%snewJacobian" % self._loggingPrefix
+ logEvent = "%sjacobian" % self._loggingPrefix
self._logger.eventBegin(logEvent)
assert(None != self.cppHandle)
- flag = self.cppHandle.needNewJacobian
+ self.cppHandle.integrateJacobian(jacobian, t.value, fields.cppHandle,
+ self.mesh.cppHandle)
self._logger.eventEnd(logEvent)
- return flag
+ return
- def integrateJacobian(self, jacobian, t, fields):
+ def integrateResidualAssembled(self, residual, t, fields):
"""
- Integrate contributions to Jacobian term at time t.
+ Integrate contributions to residual term at time t that do not
+ require assembly over cells, vertices, or processors.
"""
- logEvent = "%sjacobian" % self._loggingPrefix
+ logEvent = "%sresidualAs" % self._loggingPrefix
self._logger.eventBegin(logEvent)
assert(None != self.cppHandle)
- self.cppHandle.integrateJacobian(jacobian, t.value, fields.cppHandle,
- self.mesh.cppHandle)
+ self.cppHandle.integrateResidualAssembled(residual, t.value,
+ fields.cppHandle,
+ self.mesh.cppHandle,
+ self.mesh.coordsys.cppHandle)
self._logger.eventEnd(logEvent)
return
+ def integrateJacobianAssembled(self, jacobian, t, fields):
+ """
+ Integrate contributions to Jacobian term at time t that do not
+ require assembly over cells, vertices, or processors.
+ """
+ logEvent = "%sjacobianAs" % self._loggingPrefix
+ self._logger.eventBegin(logEvent)
+
+ assert(None != self.cppHandle)
+ self.cppHandle.integrateJacobianAssembled(jacobian, t.value,
+ fields.cppHandle,
+ self.mesh.cppHandle)
+ self._logger.eventEnd(logEvent)
+ return
+
+
def poststep(self, t, dt, totalTime, fields):
"""
Hook for doing stuff after advancing time step.
@@ -219,9 +254,11 @@
"init",
"timestep",
"solnIncr",
+ "newJacobian",
"residual",
- "newJacobian",
"jacobian",
+ "residualAs",
+ "jacobianAs",
"poststep",
"finalize"]
for event in events:
Modified: short/3D/PyLith/trunk/pylith/problems/Formulation.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Formulation.py 2008-10-02 01:07:36 UTC (rev 12983)
+++ short/3D/PyLith/trunk/pylith/problems/Formulation.py 2008-10-03 00:18:00 UTC (rev 12984)
@@ -415,10 +415,22 @@
self._debug.log(resourceUsageString())
import pylith.utils.petsc as petsc
petsc.mat_setzero(self.jacobian)
+
+ # Add in contributions that require assembly
for integrator in self.integrators:
integrator.timeStep(dt)
integrator.integrateJacobian(self.jacobian, t+dt, self.fields)
+ self._debug.log(resourceUsageString()) # TEMPORARY
+
+ self._info.log("Assembling Jacobian of operator.")
petsc.mat_assemble(self.jacobian)
+
+ self._info.log("Reforming assembled portion of Jacobian of operator.")
+ # Add in contributions that do not require assembly
+ for integrator in self.integrators:
+ integrator.integrateJacobianAssembled(self.jacobian, t+dt, self.fields)
+ self._debug.log(resourceUsageString()) # TEMPORARY
+
if self.viewJacobian:
filename = self._createJacobianFilename(t+dt)
petsc.mat_view_binary(self.jacobian, filename)
@@ -440,6 +452,7 @@
timeStamp = repr(time).rjust(self.jacobianTimeWidth, '0')
filename = basename + "_t" + timeStamp + ".mat"
return filename
+
def _reformResidual(self, t, dt):
"""
@@ -449,12 +462,20 @@
residual = self.fields.getReal("residual")
import pylith.topology.topology as bindings
bindings.zeroRealSection(residual)
+
+ # Add in contributions that require assembly
for integrator in self.integrators:
integrator.timeStep(dt)
integrator.integrateResidual(residual, t, self.fields)
self._info.log("Completing residual.")
bindings.completeSection(self.mesh.cppHandle, residual)
+
+ self._info.log("Integrating assembled residual term in operator.")
+ # Add in contributions that do not require assembly
+ for integrator in self.integrators:
+ integrator.integrateResidualAssembled(residual, t, self.fields)
+
return
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc 2008-10-02 01:07:36 UTC (rev 12983)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc 2008-10-03 00:18:00 UTC (rev 12984)
@@ -181,19 +181,6 @@
tolerance);
} // for
- // Check pairing of constraint vertices with cells
- iVertex = 0;
- for (Mesh::label_sequence::iterator v_iter=vertices->begin();
- v_iter != verticesEnd;
- ++v_iter, ++iVertex) {
- const int fiberDim = fault._faultVertexCell->getFiberDimension(*v_iter);
- CPPUNIT_ASSERT_EQUAL(1, fiberDim);
- const int_section_type::value_type* vertexCell =
- fault._faultVertexCell->restrictPoint(*v_iter);
- CPPUNIT_ASSERT(0 != vertexCell);
- CPPUNIT_ASSERT_EQUAL(_data->constraintCells[iVertex], vertexCell[0]);
- } // for
-
// Check pseudoStiffness
iVertex = 0;
for (Mesh::label_sequence::iterator v_iter=vertices->begin();
@@ -275,6 +262,202 @@
residual->restrictPoint(*v_iter);
CPPUNIT_ASSERT(0 != vals);
+ for (int i=0; i < fiberDimE; ++i) {
+ const double valE = 0.0; // no contribution
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, vals[i], tolerance);
+ } // for
+ } // for
+ } // Integrate residual with solution (as opposed to solution increment).
+
+ residual->zero();
+ { // Integrate residual with solution increment.
+ fault.useSolnIncr(true);
+ fault.integrateResidual(residual, t, &fields, mesh, &cs);
+
+ //residual->view("RESIDUAL"); // DEBUGGING
+
+ // Check values
+ const double* valsE = _data->valsResidualIncr;
+ iVertex = 0;
+ const int fiberDimE = spaceDim;
+ const double tolerance = 1.0e-06;
+ for (Mesh::label_sequence::iterator v_iter=vBegin;
+ v_iter != vEnd;
+ ++v_iter, ++iVertex) {
+ const int fiberDim = residual->getFiberDimension(*v_iter);
+ CPPUNIT_ASSERT_EQUAL(fiberDimE, fiberDim);
+ const real_section_type::value_type* vals =
+ residual->restrictPoint(*v_iter);
+ CPPUNIT_ASSERT(0 != vals);
+
+ for (int i=0; i < fiberDimE; ++i) {
+ const double valE = 0.0; // no contribution
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, vals[i], tolerance);
+ } // for
+ } // for
+ } // Integrate residual with solution increment.
+} // testIntegrateResidual
+
+// ----------------------------------------------------------------------
+// Test integrateJacobian().
+void
+pylith::faults::TestFaultCohesiveKin::testIntegrateJacobian(void)
+{ // testIntegrateJacobian
+ ALE::Obj<Mesh> mesh;
+ FaultCohesiveKin fault;
+ _initialize(&mesh, &fault);
+
+ // Setup fields
+ topology::FieldsManager fields(mesh);
+ fields.addReal("residual");
+ fields.addReal("solution");
+ fields.solutionField("solution");
+
+ const ALE::Obj<real_section_type>& residual = fields.getReal("residual");
+ CPPUNIT_ASSERT(!residual.isNull());
+ const int spaceDim = _data->spaceDim;
+ residual->setChart(mesh->getSieve()->getChart());
+ residual->setFiberDimension(mesh->depthStratum(0), spaceDim);
+ mesh->allocate(residual);
+ residual->zero();
+ fields.copyLayout("residual");
+
+ const ALE::Obj<real_section_type>& solution = fields.getReal("solution");
+ CPPUNIT_ASSERT(!solution.isNull());
+
+ const ALE::Obj<Mesh::label_sequence>& vertices = mesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+ const Mesh::label_sequence::iterator vBegin = vertices->begin();
+ const Mesh::label_sequence::iterator vEnd = vertices->end();
+ int iVertex = 0;
+ for (Mesh::label_sequence::iterator v_iter=vBegin;
+ v_iter != vEnd;
+ ++v_iter, ++iVertex) {
+ solution->updatePoint(*v_iter, &_data->fieldT[iVertex*spaceDim]);
+ } // for
+
+ PetscMat jacobian;
+ PetscErrorCode err = MeshCreateMatrix(mesh, solution, MATMPIBAIJ, &jacobian);
+ CPPUNIT_ASSERT(0 == err);
+
+ const double t = 2.134;
+ fault.integrateJacobian(&jacobian, t, &fields, mesh);
+ CPPUNIT_ASSERT_EQUAL(false, fault.needNewJacobian());
+
+ err = MatAssemblyBegin(jacobian, MAT_FINAL_ASSEMBLY);
+ CPPUNIT_ASSERT(0 == err);
+ err = MatAssemblyEnd(jacobian, MAT_FINAL_ASSEMBLY);
+ CPPUNIT_ASSERT(0 == err);
+
+ //MatView(jacobian, PETSC_VIEWER_STDOUT_WORLD); // DEBUGGING
+
+ const double* valsE = _data->valsJacobian;
+ const int nrowsE = solution->sizeWithBC();
+ const int ncolsE = nrowsE;
+
+ int nrows = 0;
+ int ncols = 0;
+ MatGetSize(jacobian, &nrows, &ncols);
+ CPPUNIT_ASSERT_EQUAL(nrowsE, nrows);
+ CPPUNIT_ASSERT_EQUAL(ncolsE, ncols);
+
+ PetscMat jDense;
+ PetscMat jSparseAIJ;
+ MatConvert(jacobian, MATSEQAIJ, MAT_INITIAL_MATRIX, &jSparseAIJ);
+ MatConvert(jSparseAIJ, MATSEQDENSE, MAT_INITIAL_MATRIX, &jDense);
+
+ double_array vals(nrows*ncols);
+ int_array rows(nrows);
+ int_array cols(ncols);
+ for (int iRow=0; iRow < nrows; ++iRow)
+ rows[iRow] = iRow;
+ for (int iCol=0; iCol < ncols; ++iCol)
+ cols[iCol] = iCol;
+ MatGetValues(jDense, nrows, &rows[0], ncols, &cols[0], &vals[0]);
+ const double tolerance = 1.0e-06;
+ for (int iRow=0; iRow < nrows; ++iRow)
+ for (int iCol=0; iCol < ncols; ++iCol) {
+ const int index = ncols*iRow+iCol;
+ const double valE = 0.0;
+#if 0 // DEBUGGING
+ if (fabs(valE-vals[index]) > tolerance)
+ std::cout << "ERROR: iRow: " << iRow << ", iCol: " << iCol
+ << "valE: " << valE
+ << ", val: " << vals[index]
+ << std::endl;
+#endif // DEBUGGING
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, vals[index], tolerance);
+ } // for
+ MatDestroy(jDense);
+ MatDestroy(jSparseAIJ);
+ CPPUNIT_ASSERT_EQUAL(false, fault.needNewJacobian());
+} // testIntegrateJacobian
+
+// ----------------------------------------------------------------------
+// Test integrateResidualAssembled().
+void
+pylith::faults::TestFaultCohesiveKin::testIntegrateResidualAssembled(void)
+{ // testIntegrateResidualAssembled
+ ALE::Obj<Mesh> mesh;
+ FaultCohesiveKin fault;
+ _initialize(&mesh, &fault);
+
+ spatialdata::geocoords::CSCart cs;
+ cs.setSpaceDim((mesh)->getDimension());
+ cs.initialize();
+
+ // Setup fields
+ topology::FieldsManager fields(mesh);
+ fields.addReal("residual");
+ fields.addReal("solution");
+ fields.solutionField("solution");
+
+ const ALE::Obj<real_section_type>& residual = fields.getReal("residual");
+ CPPUNIT_ASSERT(!residual.isNull());
+ const int spaceDim = _data->spaceDim;
+ residual->setChart(mesh->getSieve()->getChart());
+ residual->setFiberDimension(mesh->depthStratum(0), spaceDim);
+ mesh->allocate(residual);
+ residual->zero();
+ fields.copyLayout("residual");
+
+ const ALE::Obj<real_section_type>& solution = fields.getReal("solution");
+ CPPUNIT_ASSERT(!solution.isNull());
+
+ const ALE::Obj<Mesh::label_sequence>& vertices = mesh->depthStratum(0);
+ CPPUNIT_ASSERT(!vertices.isNull());
+ const Mesh::label_sequence::iterator vBegin = vertices->begin();
+ const Mesh::label_sequence::iterator vEnd = vertices->end();
+ int iVertex = 0;
+ for (Mesh::label_sequence::iterator v_iter=vBegin;
+ v_iter != vEnd;
+ ++v_iter, ++iVertex) {
+ solution->updatePoint(*v_iter, &_data->fieldT[iVertex*spaceDim]);
+ } // for
+
+ const double t = 2.134;
+ const double dt = 0.01;
+ fault.timeStep(dt);
+ { // Integrate residual with solution (as opposed to solution increment).
+ fault.useSolnIncr(false);
+ fault.integrateResidual(residual, t, &fields, mesh, &cs);
+
+ //residual->view("RESIDUAL"); // DEBUGGING
+
+ // Check values
+ const double* valsE = _data->valsResidual;
+ iVertex = 0;
+ const int fiberDimE = spaceDim;
+ const double tolerance = 1.0e-06;
+ for (Mesh::label_sequence::iterator v_iter=vBegin;
+ v_iter != vEnd;
+ ++v_iter, ++iVertex) {
+ const int fiberDim = residual->getFiberDimension(*v_iter);
+ CPPUNIT_ASSERT_EQUAL(fiberDimE, fiberDim);
+ const real_section_type::value_type* vals =
+ residual->restrictPoint(*v_iter);
+ CPPUNIT_ASSERT(0 != vals);
+
const bool isConstraint = _isConstraintVertex(*v_iter);
const double pseudoStiffness =
(!isConstraint) ? _data->pseudoStiffness : 1.0;
@@ -323,13 +506,13 @@
} // for
} // for
} // Integrate residual with solution increment.
-} // testIntegrateResidual
+} // testIntegrateResidualAssembled
// ----------------------------------------------------------------------
-// Test integrateJacobian().
+// Test integrateJacobianAssembled().
void
-pylith::faults::TestFaultCohesiveKin::testIntegrateJacobian(void)
-{ // testIntegrateJacobian
+pylith::faults::TestFaultCohesiveKin::testIntegrateJacobianAssembled(void)
+{ // testIntegrateJacobianAssembled
ALE::Obj<Mesh> mesh;
FaultCohesiveKin fault;
_initialize(&mesh, &fault);
@@ -423,10 +606,10 @@
MatDestroy(jDense);
MatDestroy(jSparseAIJ);
CPPUNIT_ASSERT_EQUAL(false, fault.needNewJacobian());
-} // testIntegrateJacobian
+} // testIntegrateJacobianAssembled
// ----------------------------------------------------------------------
-// Test integrateJacobian().
+// Test calcTractionsChange().
void
pylith::faults::TestFaultCohesiveKin::testCalcTractionsChange(void)
{ // testCalcTractionsChange
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.hh 2008-10-02 01:07:36 UTC (rev 12983)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.hh 2008-10-03 00:18:00 UTC (rev 12984)
@@ -99,6 +99,12 @@
/// Test integrateJacobian().
void testIntegrateJacobian(void);
+ /// Test integrateResidualAssembled().
+ void testIntegrateResidualAssembled(void);
+
+ /// Test integrateJacobianAssembled().
+ void testIntegrateJacobianAssembled(void);
+
/// Test _calcTractionsChange().
void testCalcTractionsChange(void);
More information about the cig-commits
mailing list