[cig-commits] r8005 - in short/3D/PyLith/trunk: .
examples/twofaults/2d libsrc/feassemble libsrc/materials
knepley at geodynamics.org
knepley at geodynamics.org
Sat Sep 22 14:14:06 PDT 2007
Author: knepley
Date: 2007-09-22 14:14:05 -0700 (Sat, 22 Sep 2007)
New Revision: 8005
Modified:
short/3D/PyLith/trunk/TODO
short/3D/PyLith/trunk/examples/twofaults/2d/pylithapp.cfg
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature0D.cc
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.cc
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din2D.cc
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din3D.cc
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2D.cc
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2Din3D.cc
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature3D.cc
short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc
short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc
Log:
Work to improve speed
Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO 2007-09-22 00:11:56 UTC (rev 8004)
+++ short/3D/PyLith/trunk/TODO 2007-09-22 21:14:05 UTC (rev 8005)
@@ -192,16 +192,22 @@
There is another version that takes an array of maxVolumes, one fo each cell in oldMesh.
-4. Run the benchmarks
+4. Check in nonlinear solver
-5. Check in nonlinear solver
+5. Throw out cohesive midnodes from censored depth (not easy since they are mixed in with new vertices, maybe change numbering)
-6. Look at more than one fault tests
+6. Coarsening for PyLith meshes
-7. Throw out cohesive midnodes from censored depth, fix it for multiple faults
+7. Time fixes
-8. Coarsening for PyLith meshes
+ a. Remove extra time in uninterpolated restrict(). May need special method for P1.
+ Done b. Precompute element geometry
+
+ c. getStateVars() is way to expensive. It looks like _getParameters() allocates an int_array. Try just using _numParamValues.
+
+ d. Same as a) for update
+
======================================================================
QUESTIONS FOR LEIF
======================================================================
Modified: short/3D/PyLith/trunk/examples/twofaults/2d/pylithapp.cfg
===================================================================
--- short/3D/PyLith/trunk/examples/twofaults/2d/pylithapp.cfg 2007-09-22 00:11:56 UTC (rev 8004)
+++ short/3D/PyLith/trunk/examples/twofaults/2d/pylithapp.cfg 2007-09-22 21:14:05 UTC (rev 8005)
@@ -35,7 +35,7 @@
# The settings below control the mesh generation (importing mesh info).
# Turn on debugging output for mesh generation.
[pylithapp.mesh_generator]
-debug = 1
+#debug = 1
# This component specification means we are using PyLith ASCII format,
# and we then specify the filename and number of space dimensions for
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc 2007-09-22 00:11:56 UTC (rev 8004)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc 2007-09-22 21:14:05 UTC (rev 8005)
@@ -90,6 +90,7 @@
_material->useElasticBehavior(!_useSolnIncr);
} // useSolnIncr
+#define FASTER
// ----------------------------------------------------------------------
// Integrate constributions to residual term (r) for operator.
void
@@ -109,6 +110,29 @@
assert(0 != fields);
assert(!mesh.isNull());
+ static PetscEvent setupEvent = 0, cellGeomEvent = 0, stateVarsEvent = 0, restrictEvent = 0, computeEvent = 0, updateEvent = 0;
+
+ if (!setupEvent) {
+ PetscLogEventRegister(&setupEvent, "IRSetup", 0);
+ }
+ if (!cellGeomEvent) {
+ PetscLogEventRegister(&cellGeomEvent, "IRCellGeom", 0);
+ }
+ if (!stateVarsEvent) {
+ PetscLogEventRegister(&stateVarsEvent, "IRStateVars", 0);
+ }
+ if (!restrictEvent) {
+ PetscLogEventRegister(&restrictEvent, "IRRestrict", 0);
+ }
+ if (!computeEvent) {
+ PetscLogEventRegister(&computeEvent, "IRCompute", 0);
+ }
+ if (!updateEvent) {
+ PetscLogEventRegister(&updateEvent, "IRUpdate", 0);
+ }
+ const Obj<sieve_type>& sieve = mesh->getSieve();
+
+ PetscLogEventBegin(setupEvent,0,0,0,0);
// Set variables dependent on dimension of cell
const int cellDim = _quadrature->cellDim();
int tensorSize = 0;
@@ -153,6 +177,9 @@
const int numBasis = _quadrature->numBasis();
const int spaceDim = _quadrature->spaceDim();
+ // Precompute the geometric and function space information
+ _quadrature->precomputeGeometry(mesh, coordinates, cells);
+
// Allocate vector for cell values.
_initCellVector();
const int cellVecSize = numBasis*spaceDim;
@@ -165,22 +192,48 @@
totalStrain[iQuad].resize(tensorSize);
totalStrain[iQuad] = 0.0;
} // for
+ PetscLogEventEnd(setupEvent,0,0,0,0);
// Loop over cells
for (Mesh::label_sequence::iterator c_iter=cells->begin();
c_iter != cellsEnd;
++c_iter) {
// Compute geometry information for current cell
+ PetscLogEventBegin(cellGeomEvent,0,0,0,0);
+#ifdef FASTER
+ _quadrature->retrieveGeometry(mesh, coordinates, *c_iter);
+#else
_quadrature->computeGeometry(mesh, coordinates, *c_iter);
+#endif
+ PetscLogEventEnd(cellGeomEvent,0,0,0,0);
// Get state variables for cell.
+ PetscLogEventBegin(stateVarsEvent,0,0,0,0);
_material->getStateVarsCell(*c_iter, numQuadPts);
+ PetscLogEventEnd(stateVarsEvent,0,0,0,0);
// Reset element vector to zero
_resetCellVector();
// Restrict input fields to cell
+ PetscLogEventBegin(restrictEvent,0,0,0,0);
+#ifdef FASTER
+ const Obj<Mesh::sieve_type::coneSequence>& cone = sieve->cone(*c_iter);
+ Mesh::sieve_type::coneSequence::iterator end = cone->end();
+ int j = -1;
+
+ for(Mesh::sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != end; ++p_iter) {
+ const real_section_type::value_type *array = dispTBctpdt->restrictPoint(*p_iter);
+ const int& dim = dispTBctpdt->getFiberDimension(*p_iter);
+
+ for(int i = 0; i < dim; ++i) {
+ dispTBctpdtCell[++j] = array[i];
+ }
+ }
+#else
mesh->restrict(dispTBctpdt, *c_iter, &dispTBctpdtCell[0], cellVecSize);
+#endif
+ PetscLogEventEnd(restrictEvent,0,0,0,0);
// Get cell geometry information that depends on cell
const double_array& basis = _quadrature->basis();
@@ -218,10 +271,12 @@
#endif
// Compute B(transpose) * sigma, first computing strains
+ PetscLogEventBegin(computeEvent,0,0,0,0);
calcTotalStrainFn(&totalStrain, basisDeriv, dispTBctpdtCell, numBasis);
const std::vector<double_array>& stress =
_material->calcStress(totalStrain);
CALL_MEMBER_FN(*this, elasticityResidualFn)(stress);
+ PetscLogEventEnd(computeEvent,0,0,0,0);
#if 0
std::cout << "Updating residual for cell " << *c_iter << std::endl;
@@ -230,7 +285,17 @@
}
#endif
// Assemble cell contribution into field
+ PetscLogEventBegin(updateEvent,0,0,0,0);
+#ifdef FASTER
+ j = 0;
+ for(Mesh::sieve_type::coneSequence::iterator p_iter = cone->begin(); p_iter != end; ++p_iter) {
+ residual->updateAddPoint(*p_iter, &_cellVector[j]);
+ j += residual->getFiberDimension(*p_iter);
+ }
+#else
mesh->updateAdd(residual, *c_iter, _cellVector);
+#endif
+ PetscLogEventEnd(updateEvent,0,0,0,0);
} // for
} // integrateResidual
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc 2007-09-22 00:11:56 UTC (rev 8004)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc 2007-09-22 21:14:05 UTC (rev 8005)
@@ -29,8 +29,14 @@
_numBasis(0),
_numQuadPts(0),
_spaceDim(0),
- _geometry(0)
+ _geometry(0),
+ _precomputed(false)
{ // constructor
+ _quadPtsPre = new real_section_type(PETSC_COMM_WORLD);
+ _jacobianPre = new real_section_type(PETSC_COMM_WORLD);
+ _jacobianDetPre = new real_section_type(PETSC_COMM_WORLD);
+ _jacobianInvPre = new real_section_type(PETSC_COMM_WORLD);
+ _basisDerivPre = new real_section_type(PETSC_COMM_WORLD);
} // constructor
// ----------------------------------------------------------------------
@@ -57,10 +63,16 @@
_numBasis(q._numBasis),
_numQuadPts(q._numQuadPts),
_spaceDim(q._spaceDim),
- _geometry(0)
+ _geometry(0),
+ _precomputed(q._precomputed)
{ // copy constructor
if (0 != q._geometry)
_geometry = q._geometry->clone();
+ _quadPtsPre = q._quadPtsPre;
+ _jacobianPre = q._jacobianPre;
+ _jacobianDetPre = q._jacobianDetPre;
+ _jacobianInvPre = q._jacobianInvPre;
+ _basisDerivPre = q._basisDerivPre;
} // copy constructor
// ----------------------------------------------------------------------
@@ -245,5 +257,67 @@
} // if
} // _checkJacobianDet
+void
+pylith::feassemble::Quadrature::resetPrecomputation()
+{
+ _precomputed = false;
+ _quadPtsPre->clear();
+ _jacobianPre->clear();
+ _jacobianDetPre->clear();
+ _jacobianInvPre->clear();
+ _basisDerivPre->clear();
+}
+void
+pylith::feassemble::Quadrature::precomputeGeometry(const ALE::Obj<Mesh>& mesh,
+ const ALE::Obj<real_section_type>& coordinates,
+ const ALE::Obj<Mesh::label_sequence>& cells)
+{
+ if (_precomputed) return;
+ const Mesh::label_sequence::iterator end = cells->end();
+
+ _quadPtsPre->setFiberDimension(cells, _numQuadPts * _spaceDim);
+ _quadPtsPre->allocatePoint();
+ _jacobianPre->setFiberDimension(cells, _numQuadPts * _cellDim * _spaceDim);
+ _jacobianPre->allocatePoint();
+ _jacobianDetPre->setFiberDimension(cells, _numQuadPts);
+ _jacobianDetPre->allocatePoint();
+ _jacobianInvPre->setFiberDimension(cells, _numQuadPts * _cellDim * _spaceDim);
+ _jacobianInvPre->allocatePoint();
+ _basisDerivPre->setFiberDimension(cells, _numQuadPts * _numBasis * _spaceDim);
+ _basisDerivPre->allocatePoint();
+ for(Mesh::label_sequence::iterator c_iter = cells->begin(); c_iter != end; ++c_iter) {
+ this->computeGeometry(mesh, coordinates, *c_iter);
+ // Set coordinates of quadrature points in cell
+ _quadPtsPre->updatePoint(*c_iter, &_quadPts[0]);
+ // Set Jacobian at quadrature points in cell
+ _jacobianPre->updatePoint(*c_iter, &_jacobian[0]);
+ // Set determinant of Jacobian at quadrature points in cell
+ _jacobianDetPre->updatePoint(*c_iter, &_jacobianDet[0]);
+ // Set inverse of Jacobian at quadrature points in cell
+ _jacobianInvPre->updatePoint(*c_iter, &_jacobianInv[0]);
+ // Set derivatives of basis functions with respect to global
+ _basisDerivPre->updatePoint(*c_iter, &_basisDeriv[0]);
+ }
+ _precomputed = true;
+}
+
+void
+pylith::feassemble::Quadrature::retrieveGeometry(const ALE::Obj<Mesh>& mesh,
+ const ALE::Obj<real_section_type>& coordinates,
+ const Mesh::point_type& cell)
+{
+ // Do they have a fast copy?
+ const real_section_type::value_type *values = _quadPtsPre->restrictPoint(cell);
+ for(int i = 0; i < _numQuadPts * _spaceDim; ++i) _quadPts[i] = values[i];
+ values = _jacobianPre->restrictPoint(cell);
+ for(int i = 0; i < _numQuadPts * _cellDim * _spaceDim; ++i) _jacobian[i] = values[i];
+ values = _jacobianDetPre->restrictPoint(cell);
+ for(int i = 0; i < _numQuadPts; ++i) _jacobianDet[i] = values[i];
+ values = _jacobianInvPre->restrictPoint(cell);
+ for(int i = 0; i < _numQuadPts * _cellDim * _spaceDim; ++i) _jacobianInv[i] = values[i];
+ values = _basisDerivPre->restrictPoint(cell);
+ for(int i = 0; i < _numQuadPts * _numBasis * _spaceDim; ++i) _basisDeriv[i] = values[i];
+}
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh 2007-09-22 00:11:56 UTC (rev 8004)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh 2007-09-22 21:14:05 UTC (rev 8005)
@@ -213,6 +213,28 @@
const ALE::Obj<real_section_type>& coordinates,
const Mesh::point_type& cell) = 0;
+ /** Reset the precomputation structures. */
+ void resetPrecomputation(void);
+
+ /** Precompute geometric quantities for each cell.
+ *
+ * @param mesh Finite-element mesh
+ * @param coordinates Section containing vertex coordinates
+ */
+ void precomputeGeometry(const ALE::Obj<Mesh>& mesh,
+ const ALE::Obj<real_section_type>& coordinates,
+ const ALE::Obj<Mesh::label_sequence>& cells);
+
+ /** Retrieve precomputed geometric quantities for a cell.
+ *
+ * @param mesh Finite-element mesh
+ * @param coordinates Section containing vertex coordinates
+ * @param cell Finite-element cell
+ */
+ void retrieveGeometry(const ALE::Obj<Mesh>& mesh,
+ const ALE::Obj<real_section_type>& coordinates,
+ const Mesh::point_type& cell);
+
// PROTECTED METHODS ////////////////////////////////////////////////////
protected :
@@ -344,6 +366,13 @@
CellGeometry* _geometry; ///< Geometry of reference cell
+ bool _precomputed;
+ /* Precomputation sections */
+ Obj<real_section_type> _quadPtsPre;
+ Obj<real_section_type> _jacobianPre;
+ Obj<real_section_type> _jacobianDetPre;
+ Obj<real_section_type> _jacobianInvPre;
+ Obj<real_section_type> _basisDerivPre;
}; // Quadrature
#include "Quadrature.icc" // inline methods
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature0D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature0D.cc 2007-09-22 00:11:56 UTC (rev 8004)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature0D.cc 2007-09-22 21:14:05 UTC (rev 8005)
@@ -20,7 +20,7 @@
// ----------------------------------------------------------------------
// Constructor
-pylith::feassemble::Quadrature0D::Quadrature0D(void)
+pylith::feassemble::Quadrature0D::Quadrature0D(void) : Quadrature()
{ // constructor
} // constructor
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.cc 2007-09-22 00:11:56 UTC (rev 8004)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.cc 2007-09-22 21:14:05 UTC (rev 8005)
@@ -20,7 +20,7 @@
// ----------------------------------------------------------------------
// Constructor
-pylith::feassemble::Quadrature1D::Quadrature1D(void)
+pylith::feassemble::Quadrature1D::Quadrature1D(void) : Quadrature()
{ // constructor
} // constructor
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din2D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din2D.cc 2007-09-22 00:11:56 UTC (rev 8004)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din2D.cc 2007-09-22 21:14:05 UTC (rev 8005)
@@ -18,7 +18,7 @@
// ----------------------------------------------------------------------
// Constructor
-pylith::feassemble::Quadrature1Din2D::Quadrature1Din2D(void)
+pylith::feassemble::Quadrature1Din2D::Quadrature1Din2D(void) : Quadrature()
{ // constructor
} // constructor
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din3D.cc 2007-09-22 00:11:56 UTC (rev 8004)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1Din3D.cc 2007-09-22 21:14:05 UTC (rev 8005)
@@ -18,7 +18,7 @@
// ----------------------------------------------------------------------
// Constructor
-pylith::feassemble::Quadrature1Din3D::Quadrature1Din3D(void)
+pylith::feassemble::Quadrature1Din3D::Quadrature1Din3D(void) : Quadrature()
{ // constructor
} // constructor
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2D.cc 2007-09-22 00:11:56 UTC (rev 8004)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2D.cc 2007-09-22 21:14:05 UTC (rev 8005)
@@ -18,7 +18,7 @@
// ----------------------------------------------------------------------
// Constructor
-pylith::feassemble::Quadrature2D::Quadrature2D(void)
+pylith::feassemble::Quadrature2D::Quadrature2D(void) : Quadrature()
{ // constructor
} // constructor
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2Din3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2Din3D.cc 2007-09-22 00:11:56 UTC (rev 8004)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature2Din3D.cc 2007-09-22 21:14:05 UTC (rev 8005)
@@ -21,7 +21,7 @@
// ----------------------------------------------------------------------
// Constructor
-pylith::feassemble::Quadrature2Din3D::Quadrature2Din3D(void)
+pylith::feassemble::Quadrature2Din3D::Quadrature2Din3D(void) : Quadrature()
{ // constructor
} // constructor
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature3D.cc 2007-09-22 00:11:56 UTC (rev 8004)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature3D.cc 2007-09-22 21:14:05 UTC (rev 8005)
@@ -18,7 +18,7 @@
// ----------------------------------------------------------------------
// Constructor
-pylith::feassemble::Quadrature3D::Quadrature3D(void)
+pylith::feassemble::Quadrature3D::Quadrature3D(void) : pylith::feassemble::Quadrature::Quadrature()
{ // constructor
} // constructor
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc 2007-09-22 00:11:56 UTC (rev 8004)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticMaterial.cc 2007-09-22 21:14:05 UTC (rev 8005)
@@ -164,7 +164,7 @@
const char** paramNames = _parameterNames();
for (int iParam=0; iParam < numParams; ++iParam) {
- const ALE::Obj<real_section_type> parameter =
+ const ALE::Obj<real_section_type>& parameter =
_parameters->getReal(paramNames[iParam]);
assert(!parameter.isNull());
Modified: short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc 2007-09-22 00:11:56 UTC (rev 8004)
+++ short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc 2007-09-22 21:14:05 UTC (rev 8005)
@@ -462,14 +462,14 @@
(*parameters)[_MaxwellIsotropic3D::pidVisStrain][iComp] =
totalStrain[iComp] - diag[iComp]*meanStrainTpdt;
} // for
- std::cout << std::endl;
- std::cout << " updateStateElastic: "<< std::endl;
- std::cout << " StrainT VisStrain Stress: " << std::endl;
- for (int iComp=0; iComp < _MaxwellIsotropic3D::tensorSize; ++iComp)
- std::cout << " " << (*parameters)[_MaxwellIsotropic3D::pidStrainT][iComp]
- << " " << (*parameters)[_MaxwellIsotropic3D::pidVisStrain][iComp]
- << " " << stress[iComp]
- << std::endl;
+// std::cout << std::endl;
+// std::cout << " updateStateElastic: "<< std::endl;
+// std::cout << " StrainT VisStrain Stress: " << std::endl;
+// for (int iComp=0; iComp < _MaxwellIsotropic3D::tensorSize; ++iComp)
+// std::cout << " " << (*parameters)[_MaxwellIsotropic3D::pidStrainT][iComp]
+// << " " << (*parameters)[_MaxwellIsotropic3D::pidVisStrain][iComp]
+// << " " << stress[iComp]
+// << std::endl;
_needNewJacobian = true;
} // _updateStateElastic
@@ -546,14 +546,14 @@
_needNewJacobian = false;
- std::cout << std::endl;
- std::cout << " updateStateViscoelastic: "<< std::endl;
- std::cout << " StrainT VisStrain Stress: " << std::endl;
- for (int iComp=0; iComp < _MaxwellIsotropic3D::tensorSize; ++iComp)
- std::cout << " " << (*parameters)[_MaxwellIsotropic3D::pidStrainT][iComp]
- << " " << (*parameters)[_MaxwellIsotropic3D::pidVisStrain][iComp]
- << " " << stress[iComp]
- << std::endl;
+// std::cout << std::endl;
+// std::cout << " updateStateViscoelastic: "<< std::endl;
+// std::cout << " StrainT VisStrain Stress: " << std::endl;
+// for (int iComp=0; iComp < _MaxwellIsotropic3D::tensorSize; ++iComp)
+// std::cout << " " << (*parameters)[_MaxwellIsotropic3D::pidStrainT][iComp]
+// << " " << (*parameters)[_MaxwellIsotropic3D::pidVisStrain][iComp]
+// << " " << stress[iComp]
+// << std::endl;
} // _updateStateViscoelastic
More information about the cig-commits
mailing list