[cig-commits] r15873 - in short/3D/PyLith/branches/pylith-friction: . libsrc libsrc/bc libsrc/faults libsrc/feassemble libsrc/problems pylith pylith/problems
brad at geodynamics.org
brad at geodynamics.org
Fri Oct 23 17:39:42 PDT 2009
Author: brad
Date: 2009-10-23 17:39:41 -0700 (Fri, 23 Oct 2009)
New Revision: 15873
Added:
short/3D/PyLith/branches/pylith-friction/libsrc/problems/SolverLumped.cc
short/3D/PyLith/branches/pylith-friction/libsrc/problems/SolverLumped.hh
short/3D/PyLith/branches/pylith-friction/pylith/problems/ExplicitLumped.py
short/3D/PyLith/branches/pylith-friction/pylith/problems/SolverLumped.py
Modified:
short/3D/PyLith/branches/pylith-friction/TODO
short/3D/PyLith/branches/pylith-friction/libsrc/Makefile.am
short/3D/PyLith/branches/pylith-friction/libsrc/bc/Neumann.cc
short/3D/PyLith/branches/pylith-friction/libsrc/bc/Neumann.hh
short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDynL.cc
short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/Integrator.hh
short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/Integrator.icc
short/3D/PyLith/branches/pylith-friction/libsrc/problems/Formulation.cc
short/3D/PyLith/branches/pylith-friction/libsrc/problems/Formulation.hh
short/3D/PyLith/branches/pylith-friction/libsrc/problems/Makefile.am
short/3D/PyLith/branches/pylith-friction/libsrc/problems/problemsfwd.hh
short/3D/PyLith/branches/pylith-friction/pylith/Makefile.am
Log:
Started work on infrastructure for solving system with lumped Jacobian.
Modified: short/3D/PyLith/branches/pylith-friction/TODO
===================================================================
--- short/3D/PyLith/branches/pylith-friction/TODO 2009-10-24 00:18:37 UTC (rev 15872)
+++ short/3D/PyLith/branches/pylith-friction/TODO 2009-10-24 00:39:41 UTC (rev 15873)
@@ -6,39 +6,7 @@
FRICTION
----------------------------------------------------------------------
-1. Update header files for documentation via doxygen. [Brad]
-2. FaultCohesiveDyn
-
- a. Fix implementation so that it works for implicit case [Brad and Surendra]
- Use Lagrange multipliers?
- Enforce interpenetration in strong sense?
- Enforce stick-slip in weak sense?
-
- Create axial and shear test cases in playpen/friction. [Surendra]
-
- b. Implement output of desired fields [Brad]
- orientation
- initial tractions
- slip
- traction
- slip
- fault constitutive model parameters
- fault constitutive model state variables
-
- c. Unit tests [Brad and Surendra]
- Brad will create base class and template.
- Surendra will complete for different cell types (quad, tri, tet, hex)
-
- Use slip-weakening friction (don't use trivial constant friction)
- testInitialize()
- testIntegrateResidual()
- testIntegrateJacobian()??
- testOrientation()
- testInitialTractions()
-
- c. cleanup/optimization
-
3. pylith::friction::FrictionModel
(base class for fault constitutive models) [Brad]
similar to Material + ElasticMaterial
@@ -54,6 +22,43 @@
Rate- and state-friction with slip law
----------------------------------------------------------------------
+LUMPED SOLVER
+----------------------------------------------------------------------
+
+Unit tests
+ FaultCohesiveKin
+ integrateJacobian (lumped)
+ adjustSolnLumped
+
+ FaultCohesiveDynL
+ integrateJacobian (lumped)
+ adjustSolnLumped
+
+ ElasticityExplicit
+ integrateJacobian (lumped)
+
+ AbsorbingDampers
+ integrateJacobian (lumped)
+
+
+FaultCohesiveKin (C++)
+ integrateJacobian (lumped)
+ adjustSolnLumped()
+
+ElasticityExplicit
+ integrateJacobian (lumped)
+
+AbsorbingDampers
+ integrateJacobian (lumped)
+
+FaultCohesiveDynL (C++)
+ integrateJacobian (lumped)
+ adjustSolnLumped()
+ compute Lagrange multipliers (zero increment of slip)
+ calls constrainSolnSpace()
+ adjust displacements
+
+----------------------------------------------------------------------
POST RELEASE 1.4
----------------------------------------------------------------------
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/Makefile.am 2009-10-24 00:18:37 UTC (rev 15872)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/Makefile.am 2009-10-24 00:39:41 UTC (rev 15873)
@@ -102,6 +102,7 @@
problems/Solver.cc \
problems/SolverLinear.cc \
problems/SolverNonlinear.cc \
+ problems/SolverLumped.cc \
topology/FieldBase.cc \
topology/Jacobian.cc \
topology/Mesh.cc \
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/bc/Neumann.cc
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/bc/Neumann.cc 2009-10-24 00:18:37 UTC (rev 15872)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/bc/Neumann.cc 2009-10-24 00:39:41 UTC (rev 15873)
@@ -165,16 +165,6 @@
} // integrateResidual
// ----------------------------------------------------------------------
-// Integrate contributions to Jacobian matrix (A) associated with
-void
-pylith::bc::Neumann::integrateJacobian(topology::Jacobian* jacobian,
- const double t,
- topology::SolutionFields* const fields)
-{ // integrateJacobian
- _needNewJacobian = false;
-} // integrateJacobian
-
-// ----------------------------------------------------------------------
// Verify configuration is acceptable.
void
pylith::bc::Neumann::verifyConfiguration(const topology::Mesh& mesh) const
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/bc/Neumann.hh
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/bc/Neumann.hh 2009-10-24 00:18:37 UTC (rev 15872)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/bc/Neumann.hh 2009-10-24 00:39:41 UTC (rev 15873)
@@ -62,17 +62,6 @@
const double t,
topology::SolutionFields* const fields);
- /** Integrate contributions to Jacobian matrix (A) associated with
- * operator.
- *
- * @param jacobian Sparse matrix for Jacobian of system.
- * @param t Current time
- * @param fields Solution fields
- */
- void integrateJacobian(topology::Jacobian* jacobian,
- const double t,
- topology::SolutionFields* const fields);
-
/** Verify configuration is acceptable.
*
* @param mesh Finite-element mesh
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDynL.cc
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDynL.cc 2009-10-24 00:18:37 UTC (rev 15872)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDynL.cc 2009-10-24 00:39:41 UTC (rev 15873)
@@ -140,6 +140,7 @@
// Create field for diagonal entries of Jacobian at conventional
// vertices i and j associated with Lagrange vertex k
+ _fields->add("Jacobian diagonal", "jacobian_diagonal");
topology::Field<topology::SubMesh>& jacobianDiag =
_fields->get("Jacobian diagonal");
jacobianDiag.newSection(slip, 2*cs->spaceDim());
@@ -675,7 +676,7 @@
double_array jacobianVertex(2*spaceDim);
const ALE::Obj<RealSection>& jacobianSection =
- fields->get("Jacobian diagonal").section();
+ _fields->get("Jacobian diagonal").section();
assert(!jacobianSection.isNull());
_updateJacobianDiagonal(*fields, jacobian);
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/Integrator.hh
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/Integrator.hh 2009-10-24 00:18:37 UTC (rev 15872)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/Integrator.hh 2009-10-24 00:39:41 UTC (rev 15873)
@@ -141,6 +141,18 @@
const double t,
topology::SolutionFields* const fields);
+ /** 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
+ */
+ virtual
+ void integrateResidualAssembled(const topology::Field<topology::Mesh>& residual,
+ const double t,
+ topology::SolutionFields* const fields);
+
/** Integrate contributions to Jacobian matrix (A) associated with
* operator.
*
@@ -153,28 +165,41 @@
const double t,
topology::SolutionFields* const fields);
- /** Integrate contributions to residual term (r) for operator that
- * do not require assembly over cells, vertices, or processors.
+ /** Integrate contributions to Jacobian matrix (A) associated with
+ * operator that do not require assembly over cells, vertices, or
+ * processors
*
- * @param residual Field containing values for residual
+ * @param jacobian Sparse matrix for Jacobian of system.
* @param t Current time
* @param fields Solution fields
*/
- virtual
- void integrateResidualAssembled(const topology::Field<topology::Mesh>& residual,
+ virtual
+ void integrateJacobianAssembled(topology::Jacobian* jacobian,
const double t,
topology::SolutionFields* const fields);
/** Integrate contributions to Jacobian matrix (A) associated with
+ * operator.
+ *
+ * @param jacobian Diagonal matrix (as field) for Jacobian of system.
+ * @param t Current time
+ * @param fields Solution fields
+ */
+ virtual
+ void integrateJacobian(const topology::Field<topology::Mesh>& jacobian,
+ const double t,
+ topology::SolutionFields* const fields);
+
+ /** Integrate contributions to Jacobian matrix (A) associated with
* operator that do not require assembly over cells, vertices, or
* processors
*
- * @param jacobian Sparse matrix for Jacobian of system.
+ * @param jacobian Diagonal matrix (as field) for Jacobian of system.
* @param t Current time
* @param fields Solution fields
*/
virtual
- void integrateJacobianAssembled(topology::Jacobian* jacobian,
+ void integrateJacobianAssembled(const topology::Field<topology::Mesh>& jacobian,
const double t,
topology::SolutionFields* const fields);
@@ -199,6 +224,18 @@
const double t,
const topology::Jacobian& jacobian);
+ /** Adjust solution from solver with lumped Jacobian to match Lagrange
+ * multiplier constraints.
+ *
+ * @param solution Solution field.
+ * @param jacobian Jacobian of the system.
+ * @param residual Residual field.
+ */
+ virtual
+ void adjustSolnLumped(topology::Field<topology::Mesh>* solution,
+ const topology::Field<topology::Mesh>& jacobian,
+ const topology::Field<topology::Mesh>& residuale);
+
/** Verify configuration is acceptable.
*
* @param mesh Finite-element mesh
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/Integrator.icc
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/Integrator.icc 2009-10-24 00:18:37 UTC (rev 15872)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/feassemble/Integrator.icc 2009-10-24 00:39:41 UTC (rev 15873)
@@ -62,6 +62,18 @@
topology::SolutionFields* const fields) {
} // integrateResidual
+// Integrate contributions to residual term (r) for operator that
+// do not require assembly over cells, vertices, or processors.
+template<typename quadrature_type>
+inline
+void
+pylith::feassemble::Integrator<quadrature_type>::integrateResidualAssembled(
+ const topology::Field<topology::Mesh>& residual,
+ const double t,
+ topology::SolutionFields* const fields) {
+ _needNewJacobian = false;
+} // integrateResidualAssembled
+
// Integrate contributions to Jacobian matrix (A) associated with
// operator.
template<typename quadrature_type>
@@ -74,17 +86,29 @@
_needNewJacobian = false;
} // integrateJacobian
-// Integrate contributions to residual term (r) for operator that
-// do not require assembly over cells, vertices, or processors.
+// Integrate contributions to Jacobian matrix (A) associated with
+// operator that do not require assembly over cells, vertices, or
+// processors
template<typename quadrature_type>
inline
void
-pylith::feassemble::Integrator<quadrature_type>::integrateResidualAssembled(
- const topology::Field<topology::Mesh>& residual,
- const double t,
- topology::SolutionFields* const fields) {
+pylith::feassemble::Integrator<quadrature_type>::integrateJacobianAssembled(
+ topology::Jacobian* jacobian,
+ const double t,
+ topology::SolutionFields* const fields) {
+} // integrateJacobianAssembled
+
+// Integrate contributions to Jacobian matrix (A) associated with
+// operator.
+template<typename quadrature_type>
+inline
+void
+pylith::feassemble::Integrator<quadrature_type>::integrateJacobian(
+ const topology::Field<topology::Mesh>& jacobian,
+ const double t,
+ topology::SolutionFields* const fields) {
_needNewJacobian = false;
-} // integrateResidualAssembled
+} // integrateJacobian
// Integrate contributions to Jacobian matrix (A) associated with
// operator that do not require assembly over cells, vertices, or
@@ -93,9 +117,9 @@
inline
void
pylith::feassemble::Integrator<quadrature_type>::integrateJacobianAssembled(
- topology::Jacobian* jacobian,
- const double t,
- topology::SolutionFields* const fields) {
+ const topology::Field<topology::Mesh>& jacobian,
+ const double t,
+ topology::SolutionFields* const fields) {
} // integrateJacobianAssembled
// Update state variables as needed.
@@ -117,7 +141,19 @@
const topology::Jacobian& jacobian) {
} // constrainSolnSpace
+// Adjust solution from solver with lumped Jacobian to match Lagrange
+// multiplier constraints.
+template<typename quadrature_type>
+inline
+void
+pylith::feassemble::Integrator<quadrature_type>::adjustSolnLumped(
+ topology::Field<topology::Mesh>* solution,
+ const topology::Field<topology::Mesh>& jacobian,
+ const topology::Field<topology::Mesh>& residuale) {
+} // adjustSolnLumped
+
+
#endif
// End of file
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/problems/Formulation.cc
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/problems/Formulation.cc 2009-10-24 00:18:37 UTC (rev 15872)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/problems/Formulation.cc 2009-10-24 00:39:41 UTC (rev 15873)
@@ -30,6 +30,7 @@
_t(0.0),
_dt(0.0),
_jacobian(0),
+ _jacobianLumped(0),
_fields(0)
{ // constructor
} // constructor
@@ -47,6 +48,7 @@
pylith::problems::Formulation::deallocate(void)
{ // deallocate
_jacobian = 0; // :TODO: Use shared pointer.
+ _jacobianLumped = 0; // :TODO: Use shared pointer.
_fields = 0; // :TODO: Use shared pointer.
} // deallocate
@@ -104,6 +106,25 @@
} // updateSettings
// ----------------------------------------------------------------------
+// Update handles and parameters for reforming the Jacobian and
+// residual.
+void
+pylith::problems::Formulation::updateSettings(topology::Field<topology::Mesh>* jacobian,
+ topology::SolutionFields* fields,
+ const double t,
+ const double dt)
+{ // updateSettings
+ assert(0 != jacobian);
+ assert(0 != fields);
+ assert(dt > 0.0);
+
+ _jacobianLumped = jacobian;
+ _fields = fields;
+ _t = t;
+ _dt = dt;
+} // updateSettings
+
+// ----------------------------------------------------------------------
// Reform system residual.
void
pylith::problems::Formulation::reformResidual(const PetscVec* tmpResidualVec,
@@ -119,6 +140,8 @@
// ===================================================================
// :KLUDGE: WAITING FOR MATT TO IMPLEMENT DIFFERENT LINE SEARCH
+ // Need a customized replacement routine for line search that will be
+ // set via SNES_SetLineSearch().
int numIntegrators = _meshIntegrators.size();
assert(numIntegrators > 0); // must have at least 1 bulk integrator
for (int i=0; i < numIntegrators; ++i) {
@@ -213,5 +236,59 @@
_jacobian->assemble("final_assembly");
} // reformJacobian
+// ----------------------------------------------------------------------
+// Reform system Jacobian.
+void
+pylith::problems::Formulation::reformJacobianLumped(void)
+{ // reformJacobian
+ assert(0 != _jacobianLumped);
+ assert(0 != _fields);
+ // Set jacobian to zero.
+ _jacobianLumped->zero();
+
+ // Add in contributions that require assembly.
+ int numIntegrators = _meshIntegrators.size();
+ for (int i=0; i < numIntegrators; ++i)
+ _meshIntegrators[i]->integrateJacobian(*_jacobianLumped, _t, _fields);
+ numIntegrators = _submeshIntegrators.size();
+ for (int i=0; i < numIntegrators; ++i)
+ _submeshIntegrators[i]->integrateJacobian(*_jacobianLumped, _t, _fields);
+
+ // Assemble jacbian.
+ _jacobianLumped->complete();
+
+ // Add in contributions that do not require assembly.
+ numIntegrators = _meshIntegrators.size();
+ for (int i=0; i < numIntegrators; ++i)
+ _meshIntegrators[i]->integrateJacobianAssembled(*_jacobianLumped,
+ _t, _fields);
+ numIntegrators = _submeshIntegrators.size();
+ for (int i=0; i < numIntegrators; ++i)
+ _submeshIntegrators[i]->integrateJacobianAssembled(*_jacobianLumped,
+ _t, _fields);
+
+} // reformJacobian
+
+// ----------------------------------------------------------------------
+// Adjust solution from solver with lumped Jacobian to match Lagrange
+// multiplier constraints.
+void
+pylith::problems::Formulation::adjustSolnLumped(
+ topology::Field<topology::Mesh>* solution,
+ const topology::Field<topology::Mesh>& jacobian,
+ const topology::Field<topology::Mesh>& residual)
+{ // adjustSolnLumped
+ assert(0 != solution);
+
+ int numIntegrators = _meshIntegrators.size();
+ for (int i=0; i < numIntegrators; ++i)
+ _meshIntegrators[i]->adjustSolnLumped(solution, jacobian, residual);
+
+ numIntegrators = _submeshIntegrators.size();
+ for (int i=0; i < numIntegrators; ++i)
+ _submeshIntegrators[i]->adjustSolnLumped(solution, jacobian, residual);
+} // adjustSolnLumped
+
+
// End of file
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/problems/Formulation.hh
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/problems/Formulation.hh 2009-10-24 00:18:37 UTC (rev 15872)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/problems/Formulation.hh 2009-10-24 00:39:41 UTC (rev 15873)
@@ -55,6 +55,12 @@
/// Deallocate PETSc and local data structures.
void deallocate(void);
+ /** Get solution fields.
+ *
+ * @returns solution fields.
+ */
+ const topology::SolutionFields& fields(void) const;
+
/** Set handles to integrators over the mesh.
*
* @param integrators Integrators over the mesh.
@@ -84,6 +90,19 @@
const double t,
const double dt);
+ /** Update handles and parameters for reforming the Jacobian and
+ * residual.
+ *
+ * @param jacobian Handle to diagonal matrix (as Field) for system Jacobian.
+ * @param fields Handle to solution fields.
+ * @param t Current time (nondimensional).
+ * @param dt Time step (nondimension).
+ */
+ void updateSettings(topology::Field<topology::Mesh>* jacobian,
+ topology::SolutionFields* fields,
+ const double t,
+ const double dt);
+
/** Reform system residual.
*
* @param tmpResidualVec Temporary PETSc vector for residual.
@@ -98,14 +117,28 @@
*/
void reformJacobian(const PetscVec* tmpSolveSolnVec =0);
- const topology::SolutionFields& fields() const;
+ /* Reform system Jacobian.
+ */
+ void reformJacobianLumped(void);
+ /** Adjust solution from solver with lumped Jacobian to match Lagrange
+ * multiplier constraints.
+ *
+ * @param solution Solution field.
+ * @param jacobian Jacobian of the system.
+ * @param residual Residual field.
+ */
+ void adjustSolnLumped(topology::Field<topology::Mesh>* solution,
+ const topology::Field<topology::Mesh>& jacobian,
+ const topology::Field<topology::Mesh>& residual);
+
// PRIVATE MEMBERS //////////////////////////////////////////////////////
private :
double _t; ///< Current time (nondimensional).
double _dt; ///< Current time step (nondimensional).
topology::Jacobian* _jacobian; ///< Handle to Jacobian of system.
+ topology::Field<topology::Mesh>* _jacobianLumped; ///< Handle to lumped Jacobian of system.
topology::SolutionFields* _fields; ///< Handle to solution fields for system.
/// Integrators over subdomains of the mesh.
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/problems/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/problems/Makefile.am 2009-10-24 00:18:37 UTC (rev 15872)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/problems/Makefile.am 2009-10-24 00:39:41 UTC (rev 15873)
@@ -18,6 +18,7 @@
Solver.hh \
SolverLinear.hh \
SolverNonlinear.hh \
+ SolverLumped.hh \
problemsfwd.hh
noinst_HEADERS =
Added: short/3D/PyLith/branches/pylith-friction/libsrc/problems/SolverLumped.cc
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/problems/SolverLumped.cc (rev 0)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/problems/SolverLumped.cc 2009-10-24 00:39:41 UTC (rev 15873)
@@ -0,0 +1,119 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#include <portinfo>
+
+#include "SolverLumped.hh" // implementation of class methods
+
+#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
+#include "pylith/topology/Jacobian.hh" // USES Jacobian
+#include "pylith/problems/Formulation.hh" // USES Formulation
+
+// ----------------------------------------------------------------------
+typedef pylith::topology::Mesh::SieveMesh SieveMesh;
+typedef pylith::topology::Mesh::RealSection RealSection;
+
+// ----------------------------------------------------------------------
+// Constructor
+pylith::problems::SolverLumped::SolverLumped(void)
+{ // constructor
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor
+pylith::problems::SolverLumped::~SolverLumped(void)
+{ // destructor
+ deallocate();
+} // destructor
+
+// ----------------------------------------------------------------------
+// Deallocate data structures.
+void
+pylith::problems::SolverLumped::deallocate(void)
+{ // deallocate
+ Solver::deallocate();
+} // deallocate
+
+// ----------------------------------------------------------------------
+// Initialize solver.
+void
+pylith::problems::SolverLumped::initialize(
+ const topology::SolutionFields& fields,
+ const topology::Field<topology::Mesh>& jacobian,
+ Formulation* formulation)
+{ // initialize
+ assert(0 != formulation);
+
+ _formulation = formulation;
+} // initialize
+
+// ----------------------------------------------------------------------
+// Solve the system.
+void
+pylith::problems::SolverLumped::solve(
+ topology::Field<topology::Mesh>* solution,
+ const topology::Field<topology::Mesh>& jacobian,
+ const topology::Field<topology::Mesh>& residual)
+{ // solve
+ assert(0 != solution);
+ assert(0 != _formulation);
+
+ // solution = residual / jacobian
+
+ const spatialdata::geocoords::CoordSys* cs = solution->mesh().coordsys();
+ assert(0 != cs);
+ const int spaceDim = cs->spaceDim();
+
+ // Get mesh vertices.
+ const ALE::Obj<SieveMesh>& sieveMesh = solution->mesh().sieveMesh();
+ assert(!sieveMesh.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& vertices =
+ sieveMesh->depthStratum(0);
+ assert(!vertices.isNull());
+ const SieveMesh::label_sequence::iterator verticesBegin =
+ vertices->begin();
+ const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
+
+ // Get sections.
+ double_array solutionVertex(spaceDim);
+ const ALE::Obj<RealSection>& solutionSection = solution->section();
+ assert(!solutionSection.isNull());
+
+ double_array jacobianVertex(spaceDim);
+ const ALE::Obj<RealSection>& jacobianSection = jacobian.section();
+ assert(!jacobianSection.isNull());
+
+ double_array residualVertex(spaceDim);
+ const ALE::Obj<RealSection>& residualSection = residual.section();
+ assert(!residualSection.isNull());
+
+ for (SieveMesh::label_sequence::iterator v_iter=verticesBegin;
+ v_iter != verticesEnd;
+ ++v_iter) {
+ jacobianSection->restrictPoint(*v_iter, &jacobianVertex[0],
+ jacobianVertex.size());
+ residualSection->restrictPoint(*v_iter, &residualVertex[0],
+ residualVertex.size());
+ solutionVertex = residualVertex / jacobianVertex;
+
+ assert(solutionSection->getFiberDimension(*v_iter) == spaceDim);
+ solutionSection->updatePoint(*v_iter, &solutionVertex[0]);
+ } // for
+
+ // Adjust solution to match constraints
+ _formulation->adjustSolnLumped(solution, jacobian, residual);
+
+ PetscLogFlops(vertices->size());
+} // solve
+
+
+// End of file
Added: short/3D/PyLith/branches/pylith-friction/libsrc/problems/SolverLumped.hh
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/problems/SolverLumped.hh (rev 0)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/problems/SolverLumped.hh 2009-10-24 00:39:41 UTC (rev 15873)
@@ -0,0 +1,79 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+/**
+ * @file libsrc/problems/SolverLumped.hh
+ *
+ * @brief Object for using PETSc scalable linear equation solvers (KSP).
+ */
+
+#if !defined(pylith_problems_solverlumped_hh)
+#define pylith_problems_solverlumped_hh
+
+// Include directives ---------------------------------------------------
+#include "Solver.hh" // ISA Solver
+
+#include "pylith/utils/petscfwd.h" // HASA PetscKSP
+
+// SolverLumped ---------------------------------------------------------
+/** @brief Object for using simple solver to solver system with lumped Jacobian.
+ */
+
+class pylith::problems::SolverLumped : Solver
+{ // SolverLumped
+ friend class TestSolverLumped; // unit testing
+
+// PUBLIC MEMBERS ///////////////////////////////////////////////////////
+public :
+
+ /// Constructor
+ SolverLumped(void);
+
+ /// Destructor
+ ~SolverLumped(void);
+
+ /// Deallocate PETSc and local data structures.
+ void deallocate(void);
+
+ /** Initialize solver.
+ *
+ * @param fields Solution fields.
+ * @param jacobian Jacobian of system.
+ * @param formulation Formulation of system of equations.
+ */
+ void
+ initialize(const topology::SolutionFields& fields,
+ const topology::Field<topology::Mesh>& jacobian,
+ Formulation* const formulation);
+
+ /** Solve the system.
+ *
+ * @param solution Solution field.
+ * @param jacobian Jacobian of the system.
+ * @param residual Residual field.
+ */
+ void solve(topology::Field<topology::Mesh>* solution,
+ const topology::Field<topology::Mesh>& jacobian,
+ const topology::Field<topology::Mesh>& residual);
+
+// NOT IMPLEMENTED //////////////////////////////////////////////////////
+private :
+
+ SolverLumped(const SolverLumped&); ///< Not implemented
+ const SolverLumped& operator=(const SolverLumped&); ///< Not implemented
+
+}; // SolverLumped
+
+#endif // pylith_problems_solverlumped_hh
+
+
+// End of file
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/problems/problemsfwd.hh
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/problems/problemsfwd.hh 2009-10-24 00:18:37 UTC (rev 15872)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/problems/problemsfwd.hh 2009-10-24 00:39:41 UTC (rev 15873)
@@ -29,6 +29,7 @@
class Solver;
class SolverLinear;
class SolverNonlinear;
+ class SolverLumped;
} // problems
} // pylith
Modified: short/3D/PyLith/branches/pylith-friction/pylith/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-friction/pylith/Makefile.am 2009-10-24 00:18:37 UTC (rev 15872)
+++ short/3D/PyLith/branches/pylith-friction/pylith/Makefile.am 2009-10-24 00:39:41 UTC (rev 15873)
@@ -104,12 +104,14 @@
perf/Jacobian.py \
problems/__init__.py \
problems/Explicit.py \
+ problems/ExplicitLumped.py \
problems/Formulation.py \
problems/Implicit.py \
problems/Problem.py \
problems/Solver.py \
problems/SolverLinear.py \
problems/SolverNonlinear.py \
+ problems/SolverLumped.py \
problems/TimeDependent.py \
problems/TimeStep.py \
problems/TimeStepAdapt.py \
Added: short/3D/PyLith/branches/pylith-friction/pylith/problems/ExplicitLumped.py
===================================================================
--- short/3D/PyLith/branches/pylith-friction/pylith/problems/ExplicitLumped.py (rev 0)
+++ short/3D/PyLith/branches/pylith-friction/pylith/problems/ExplicitLumped.py 2009-10-24 00:39:41 UTC (rev 15873)
@@ -0,0 +1,246 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file pylith/problems/ExplicitLumped.py
+##
+## @brief Python ExplicitLumped object for solving equations using an
+## explicit formulation with a lumped Jacobian matrix that is stored
+## as a Field.
+##
+## Factory: pde_formulation
+
+from Formulation import Formulation
+from pylith.utils.profiling import resourceUsageString
+
+# ExplicitLumped class
+class ExplicitLumped(Formulation):
+ """
+ Python ExplicitLumped object for solving equations using an explicit
+ formulation.
+
+ The formulation has the general form, [A(t)] {u(t+dt)} = {b(t)},
+ where we want to solve for {u(t+dt)}, A(t) is usually constant
+ (i.e., independent of time), and {b(t)} usually depends on {u(t)}
+ and {u(t-dt)}.
+
+ Jacobian: A(t)
+ solution: u(t+dt)
+ residual: b(t) - A(t) \hat u(t+dt)
+ constant: b(t)
+
+ Factory: pde_formulation.
+ """
+
+ # INVENTORY //////////////////////////////////////////////////////////
+
+ class Inventory(Formulation.Inventory):
+ """
+ Python object for managing ExplicitLumped facilities and properties.
+
+ Provide appropriate solver for lumped Jacobian as the default.
+ """
+
+ ## @class Inventory
+ ## Python object for managing ExplicitLumped facilities and properties.
+ ##
+ ## \b Properties
+ ## @li None
+ ##
+ ## \b Facilities
+ ## @li \b solver Algebraic solver.
+
+ import pyre.inventory
+
+ from SolverLumped import SolverLumped
+ solver = pyre.inventory.facility("solver", family="solver",
+ factory=SolverLumped)
+ solver.meta['tip'] = "Algebraic solver."
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="explicit"):
+ """
+ Constructor.
+ """
+ Formulation.__init__(self, name)
+ self._loggingPrefix = "TSEx "
+ return
+
+
+ def elasticityIntegrator(self):
+ """
+ Get integrator for elastic material.
+ """
+ from pylith.feassemble.ElasticityExplicit import ElasticityExplicit
+ return ElasticityExplicit()
+
+
+ def initialize(self, dimension, normalizer):
+ """
+ Initialize problem for explicit time integration.
+ """
+ logEvent = "%sinit" % self._loggingPrefix
+ self._eventLogger.eventBegin(logEvent)
+
+ Formulation.initialize(self, dimension, normalizer)
+
+ from pylith.utils.petsc import MemoryLogger
+ logger = MemoryLogger.singleton()
+ logger.setDebug(0)
+ logger.stagePush("Problem")
+
+ # Allocate other fields, reusing layout from dispIncr
+ self._info.log("Creating other fields.")
+ self.fields.add("disp(t-dt)", "displacement")
+ self.fields.copyLayout("dispIncr(t->t+dt)")
+ self._debug.log(resourceUsageString())
+
+ # Setup fields and set to zero
+ dispTmdt = self.fields.get("disp(t-dt)")
+ dispTmdt.zero()
+ dispT = self.fields.get("disp(t)")
+ dispT.zero()
+ residual = self.fields.get("residual")
+ residual.zero()
+ residual.createVector()
+ self._debug.log(resourceUsageString())
+ logger.stagePop()
+
+ self._info.log("Creating lumped Jacobian matrix.")
+ from pylith.topology.topology.MeshField import MeshField
+ jacobian = MeshField(self.mesh)
+ jacobian.label("jacobian")
+ solution.vectorFieldType(solution.VECTOR)
+ jacobian.newSection(jacobian.VERTICES_FIELD, dimension)
+ self.jacobian = jacobian
+ self._debug.log(resourceUsageString())
+
+ logger.stagePush("Problem")
+ self._info.log("Initializing solver.")
+ self.solver.initialize(self.fields, self.jacobian, self)
+ self._debug.log(resourceUsageString())
+
+ # Solve for increment in displacement field.
+ for constraint in self.constraints:
+ constraint.useSolnIncr(True)
+ for integrator in self.integratorsMesh + self.integratorsSubMesh:
+ integrator.useSolnIncr(True)
+
+ logger.stagePop()
+ logger.setDebug(0)
+ self._eventLogger.eventEnd(logEvent)
+ return
+
+
+ def prestep(self, t, dt):
+ """
+ Hook for doing stuff before advancing time step.
+ """
+ logEvent = "%sprestep" % self._loggingPrefix
+ self._eventLogger.eventBegin(logEvent)
+
+ dispIncr = self.fields.get("dispIncr(t->t+dt)")
+ for constraint in self.constraints:
+ constraint.setFieldIncr(t, t+dt, dispIncr)
+
+ needNewJacobian = False
+ for integrator in self.integratorsMesh + self.integratorsSubMesh:
+ integrator.timeStep(dt)
+ if integrator.needNewJacobian():
+ needNewJacobian = True
+ if needNewJacobian:
+ self._reformJacobian(t, dt)
+
+ self._eventLogger.eventEnd(logEvent)
+ return
+
+
+ def step(self, t, dt):
+ """
+ Advance to next time step.
+ """
+ logEvent = "%sstep" % self._loggingPrefix
+ self._eventLogger.eventBegin(logEvent)
+
+ self._reformResidual(t, dt)
+
+ self._info.log("Solving equations.")
+ residual = self.fields.get("residual")
+ dispIncr = self.fields.get("dispIncr(t->t+dt)")
+ self.solver.solve(dispIncr, self.jacobian, residual)
+
+ self._eventLogger.eventEnd(logEvent)
+ return
+
+
+ def poststep(self, t, dt):
+ """
+ Hook for doing stuff after advancing time step.
+ """
+ logEvent = "%spoststep" % self._loggingPrefix
+ self._eventLogger.eventBegin(logEvent)
+
+ dispIncr = self.fields.get("dispIncr(t->t+dt)")
+ dispT = self.fields.get("disp(t)")
+ dispTmdt = self.fields.get("disp(t-dt)")
+
+ dispTmdt.copy(dispT)
+ dispT += dispIncr
+ dispIncr.zero()
+
+ Formulation.poststep(self, t, dt)
+
+ self._eventLogger.eventEnd(logEvent)
+ return
+
+
+ # PRIVATE METHODS ////////////////////////////////////////////////////
+
+ def _configure(self):
+ """
+ Set members based using inventory.
+ """
+ Formulation._configure(self)
+ self.solver = self.inventory.solver
+ return
+
+
+ def _reformJacobian(self, t, dt):
+ """
+ Reform Jacobian matrix for operator.
+ """
+ self._debug.log(resourceUsageString())
+ self._info.log("Integrating Jacobian operator.")
+ self._eventLogger.stagePush("Reform Jacobian")
+
+ self.updateSettings(self.jacobian, fields, t, dt)
+ self.reformJacobianLumped()
+
+ self._eventLogger.stagePop()
+
+ if self.viewJacobian:
+ self.jacobian.view("Lumped Jacobian")
+
+ self._debug.log(resourceUsageString())
+ return
+
+
+# FACTORIES ////////////////////////////////////////////////////////////
+
+def pde_formulation():
+ """
+ Factory associated with ExplicitLumped.
+ """
+ return ExplicitLumped()
+
+
+# End of file
Added: short/3D/PyLith/branches/pylith-friction/pylith/problems/SolverLumped.py
===================================================================
--- short/3D/PyLith/branches/pylith-friction/pylith/problems/SolverLumped.py (rev 0)
+++ short/3D/PyLith/branches/pylith-friction/pylith/problems/SolverLumped.py 2009-10-24 00:39:41 UTC (rev 15873)
@@ -0,0 +1,91 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file pylith/solver/SolverLumped.py
+
+## @brief Python PyLith simple solver for system with a lumped (i.e.,
+## diagonal) Jacobian matrix.
+
+from Solver import Solver
+from problems import SolverLumped as ModuleSolverLumped
+
+# SolverLumped class
+class SolverLumped(Solver, ModuleSolverLumped):
+ """
+ Python PyLith linear algebraic solver.
+ """
+
+ # INVENTORY //////////////////////////////////////////////////////////
+
+ class Inventory(Solver.Inventory):
+ """
+ Python object for managing SolverLumped facilities and properties.
+ """
+
+ ## @class Inventory
+ ## Python object for managing SolverLumped facilities and properties.
+ ##
+ ## \b Properties
+ ## @li None
+ ##
+ ## \b Facilities
+ ## @li None
+
+ import pyre.inventory
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="solverlinear"):
+ """
+ Constructor.
+ """
+ Solver.__init__(self, name)
+ ModuleSolverLumped.__init__(self)
+ return
+
+
+ def initialize(self, fields, jacobian, formulation):
+ """
+ Initialize linear solver.
+ """
+ ModuleSolverLumped.initialize(self, fields, jacobian, formulation)
+ return
+
+
+ # PRIVATE METHODS /////////////////////////////////////////////////////
+
+ def _configure(self):
+ """
+ Set members based using inventory.
+ """
+ Solver._configure(self)
+ return
+
+
+ def _cleanup(self):
+ """
+ Deallocate PETSc and local data structures.
+ """
+ self.deallocate()
+ return
+
+
+# FACTORIES ////////////////////////////////////////////////////////////
+
+def solver():
+ """
+ Factory associated with Solver.
+ """
+ return SolverLumped()
+
+
+# End of file
More information about the CIG-COMMITS
mailing list