[cig-commits] r8274 - in short/3D/PyLith/trunk: . libsrc/feassemble
tests/1d/line2 unittests/libtests/feassemble/data
brad at geodynamics.org
brad at geodynamics.org
Fri Nov 9 21:05:10 PST 2007
Author: brad
Date: 2007-11-09 21:05:10 -0800 (Fri, 09 Nov 2007)
New Revision: 8274
Modified:
short/3D/PyLith/trunk/TODO
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc
short/3D/PyLith/trunk/tests/1d/line2/dislocation_dyn.cfg
short/3D/PyLith/trunk/tests/1d/line2/dislocation_dyn_sliprate.spatialdb
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicit.py
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicitData1DLinear.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicitData1DQuadratic.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicitData2DLinear.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicitData2DQuadratic.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicitData3DLinear.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicitData3DQuadratic.cc
Log:
Fixed bugs in integrating residual for elastic explicit time stepping (don't use disp(t+dt)).
Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO 2007-11-10 00:50:16 UTC (rev 8273)
+++ short/3D/PyLith/trunk/TODO 2007-11-10 05:05:10 UTC (rev 8274)
@@ -4,7 +4,7 @@
Release 1.0.2
- 1. Fix memory use on proc 0 when partitioning.
+ 1. Fix memory imbalance associated with distribution.
2. Fix problems with buildbot on darwin.
@@ -15,6 +15,11 @@
Release 1.1
+ 0. Test dynamic time stepping.
+ Switch to more efficient restrict(). Need tags in ElasticityExplicit.
+ Can we refactor some of the elasticity integration stuff that is
+ common to both implicit and explicit time stepping?
+
1. Finish Neumann BC.
2. Finish velocity Dirichlet BC.
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc 2007-11-10 00:50:16 UTC (rev 8273)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc 2007-11-10 05:05:10 UTC (rev 8274)
@@ -190,12 +190,38 @@
totalStrain[iQuad] = 0.0;
} // for
- int c = 0;
+#ifdef FASTER
+ if (_dispTpdtTags.find(_material->id()) == _dispTpdtTags.end()) {
+ _dispTpdtTags[_material->id()] =
+ mesh->calculateCustomAtlas(dispTpdt, cells);
+ } // if
+ const int dispTpdtTag = _dispTpdtTags[_material->id()];
+
+ if (_dispTTags.find(_material->id()) == _dispTTags.end()) {
+ _dispTTags[_material->id()] =
+ mesh->copyCustomAtlas(dispTpdt, _dispTpdtTags[_material->id()]);
+ } // if
+ const int dispTTag = _dispTTags[_material->id()];
+
+ if (_dispTmdtTags.find(_material->id()) == _dispTmdtTags.end()) {
+ _dispTmdtTags[_material->id()] =
+ mesh->copyCustomAtlas(dispTpdt, _dispTpdtTags[_material->id()]);
+ } // if
+ const int dispTTag = _dispTTags[_material->id()];
+
+ if (_residualTags.find(_material->id()) == _residualTags.end()) {
+ _residualTags[_material->id()] =
+ mesh->copyCustomAtlas(dispTpdt, _dispTpdtTags[_material->id()]);
+ } // if
+ const int residualTag = _residualTags[_material->id()];
+#endif
+
+ int c_index = 0;
for (Mesh::label_sequence::iterator c_iter=cells->begin();
c_iter != cellsEnd;
- ++c_iter) {
+ ++c_iter, ++c_index) {
// Compute geometry information for current cell
- _quadrature->retrieveGeometry(mesh, coordinates, *c_iter, c);
+ _quadrature->retrieveGeometry(mesh, coordinates, *c_iter, c_index);
// Get state variables for cell.
_material->getStateVarsCell(*c_iter, numQuadPts);
@@ -203,12 +229,18 @@
// Reset element vector to zero
_resetCellVector();
- // Restrict input fields to cell
- // :TODO: Replace with optimized restricts().
- // NEED to add optimization tags to ElasticityExplicit object.
+#ifdef FASTER
+ mesh->restrict(dispTpdt, dispTpdtTag, c_index, &dispTpdtCell[0],
+ cellVecSize);
+ mesh->restrict(dispTpdt, dispTTag, c_index, &dispTCell[0],
+ cellVecSize);
+ mesh->restrict(dispTpdt, dispTmdtTag, c_index, &dispTmdtCell[0],
+ cellVecSize);
+#else
mesh->restrict(dispTpdt, *c_iter, &dispTpdtCell[0], cellVecSize);
mesh->restrict(dispT, *c_iter, &dispTCell[0], cellVecSize);
mesh->restrict(dispTmdt, *c_iter, &dispTmdtCell[0], cellVecSize);
+#endif
// Get cell geometry information that depends on cell
const double_array& basis = _quadrature->basis();
@@ -227,8 +259,7 @@
const double valIJ = valI * basis[iQuad*numBasis+jBasis];
for (int iDim=0; iDim < spaceDim; ++iDim)
_cellVector[iBasis*spaceDim+iDim] +=
- valIJ * (- dispTpdtCell[jBasis*spaceDim+iDim]
- + 2.0 * dispTCell[jBasis*spaceDim+iDim]
+ valIJ * (+ 2.0 * dispTCell[jBasis*spaceDim+iDim]
- dispTmdtCell[jBasis*spaceDim+iDim]);
} // for
} // for
@@ -242,7 +273,11 @@
CALL_MEMBER_FN(*this, elasticityResidualFn)(stress);
// Assemble cell contribution into field
+#ifdef FASTER
+ mesh->updateAdd(residual, residualTag, c_index, _cellVector);
+#else
mesh->updateAdd(residual, *c_iter, _cellVector);
+#endif
} // for
} // integrateResidual
@@ -291,13 +326,13 @@
// Allocate vector for cell values (if necessary)
_initCellMatrix();
- int c = 0;
+ int c_index = 0;
for (Mesh::label_sequence::iterator c_iter=cells->begin();
c_iter != cellsEnd;
- ++c_iter) {
+ ++c_iter, ++c_index) {
// Compute geometry information for current cell
- _quadrature->retrieveGeometry(mesh, coordinates, *c_iter, c);
+ _quadrature->retrieveGeometry(mesh, coordinates, *c_iter, c_index);
// Get state variables for cell.
_material->getStateVarsCell(*c_iter, numQuadPts);
@@ -403,15 +438,23 @@
totalStrain[iQuad] = 0.0;
} // for
- int c = 0;
+#ifdef FASTER
+ const int dispTTag = _dispTTags[_material->id()];
+#endif
+
+ int c_index = 0;
for (Mesh::label_sequence::iterator c_iter=cells->begin();
c_iter != cellsEnd;
- ++c_iter) {
+ ++c_iter, ++c_index) {
// Compute geometry information for current cell
- _quadrature->retrieveGeometry(mesh, coordinates, *c_iter, c);
+ _quadrature->retrieveGeometry(mesh, coordinates, *c_iter, c_index);
// Restrict input fields to cell
+#ifdef FASTER
+ mesh->restrict(disp, dispTTag, c_index, &dispCell[0], cellVecSize);
+#else
mesh->restrict(disp, *c_iter, &dispCell[0], cellVecSize);
+#endif
// Get cell geometry information that depends on cell
const double_array& basisDeriv = _quadrature->basisDeriv();
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc 2007-11-10 00:50:16 UTC (rev 8273)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc 2007-11-10 05:05:10 UTC (rev 8274)
@@ -176,7 +176,7 @@
const int spaceDim = _quadrature->spaceDim();
#ifdef FASTER
- int c = 0;
+ int c_index = 0;
if (_dTags.find(_material->id()) == _dTags.end()) {
_dTags[_material->id()] = mesh->calculateCustomAtlas(dispTBctpdt, cells);
@@ -207,10 +207,10 @@
// Loop over cells
for (Mesh::label_sequence::iterator c_iter=cells->begin();
c_iter != cellsEnd;
- ++c_iter) {
+ ++c_iter, ++c_index) {
// Compute geometry information for current cell
PetscLogEventBegin(cellGeomEvent,0,0,0,0);
- _quadrature->retrieveGeometry(mesh, coordinates, *c_iter, c);
+ _quadrature->retrieveGeometry(mesh, coordinates, *c_iter, c_index);
PetscLogEventEnd(cellGeomEvent,0,0,0,0);
// Get state variables for cell.
@@ -224,7 +224,8 @@
// Restrict input fields to cell
PetscLogEventBegin(restrictEvent,0,0,0,0);
#ifdef FASTER
- mesh->restrict(dispTBctpdt, dTag, c, &dispTBctpdtCell[0], cellVecSize);
+ mesh->restrict(dispTBctpdt, dTag, c_index, &dispTBctpdtCell[0],
+ cellVecSize);
#else
mesh->restrict(dispTBctpdt, *c_iter, &dispTBctpdtCell[0], cellVecSize);
#endif
@@ -282,7 +283,7 @@
// Assemble cell contribution into field
PetscLogEventBegin(updateEvent,0,0,0,0);
#ifdef FASTER
- mesh->updateAdd(residual, rTag, c++, _cellVector);
+ mesh->updateAdd(residual, rTag, c_index, _cellVector);
#else
mesh->updateAdd(residual, *c_iter, _cellVector);
#endif
@@ -385,12 +386,12 @@
#endif
// Loop over cells
- int c = 0;
+ int c_index = 0;
for (Mesh::label_sequence::iterator c_iter=cells->begin();
c_iter != cellsEnd;
- ++c_iter, ++c) {
+ ++c_iter, ++c_index) {
// Compute geometry information for current cell
- _quadrature->retrieveGeometry(mesh, coordinates, *c_iter, c);
+ _quadrature->retrieveGeometry(mesh, coordinates, *c_iter, c_index);
// Get state variables for cell.
_material->getStateVarsCell(*c_iter, numQuadPts);
@@ -400,7 +401,8 @@
// Restrict input fields to cell
#ifdef FASTER
- mesh->restrict(dispTBctpdt, dTag, c, &dispTBctpdtCell[0], cellVecSize);
+ mesh->restrict(dispTBctpdt, dTag, c_index, &dispTBctpdtCell[0],
+ cellVecSize);
#else
mesh->restrict(dispTBctpdt, *c_iter, &dispTBctpdtCell[0], cellVecSize);
#endif
@@ -497,16 +499,16 @@
const int dTag = _dTags[_material->id()];
#endif
- int c = 0;
+ int c_index = 0;
for (Mesh::label_sequence::iterator c_iter=cells->begin();
c_iter != cellsEnd;
- ++c_iter, ++c) {
+ ++c_iter, ++c_index) {
// Compute geometry information for current cell
- _quadrature->retrieveGeometry(mesh, coordinates, *c_iter, c);
+ _quadrature->retrieveGeometry(mesh, coordinates, *c_iter, c_index);
// Restrict input fields to cell
#ifdef FASTER
- mesh->restrict(disp, dTag, c, &dispCell[0], cellVecSize);
+ mesh->restrict(disp, dTag, c_index, &dispCell[0], cellVecSize);
#else
mesh->restrict(disp, *c_iter, &dispCell[0], cellVecSize);
#endif
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc 2007-11-10 00:50:16 UTC (rev 8273)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc 2007-11-10 05:05:10 UTC (rev 8274)
@@ -326,9 +326,9 @@
void
pylith::feassemble::Quadrature::retrieveGeometry(
const ALE::Obj<Mesh>& mesh,
- const ALE::Obj<real_section_type>& coordinates,
+ const ALE::Obj<real_section_type>& coordinates,
const Mesh::point_type& cell,
- const int c)
+ const int c)
{ // retrieveGeometry
const real_section_type::value_type* values =
mesh->restrict(_quadPtsPre, _qTag, c);
Modified: short/3D/PyLith/trunk/tests/1d/line2/dislocation_dyn.cfg
===================================================================
--- short/3D/PyLith/trunk/tests/1d/line2/dislocation_dyn.cfg 2007-11-10 00:50:16 UTC (rev 8273)
+++ short/3D/PyLith/trunk/tests/1d/line2/dislocation_dyn.cfg 2007-11-10 05:05:10 UTC (rev 8274)
@@ -76,8 +76,8 @@
[pylithapp.timedependent.interfaces.fault.eq_src.slip_function]
slip.iohandler.filename = dislocation_slip.spatialdb
-slip_rate.iohandler.filename = dislocation_sliprate.spatialdb
-slip_time.iohandler.filename = dislocation_sliptime.spatialdb
+slip_rate.iohandler.filename = dislocation_dyn_sliprate.spatialdb
+slip_time.iohandler.filename = dislocation_dyn_sliptime.spatialdb
# ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/tests/1d/line2/dislocation_dyn_sliprate.spatialdb
===================================================================
--- short/3D/PyLith/trunk/tests/1d/line2/dislocation_dyn_sliprate.spatialdb 2007-11-10 00:50:16 UTC (rev 8273)
+++ short/3D/PyLith/trunk/tests/1d/line2/dislocation_dyn_sliprate.spatialdb 2007-11-10 05:05:10 UTC (rev 8274)
@@ -11,4 +11,4 @@
space-dim = 1
}
}
-0.0 2.0e+3
+0.0 1.0e+3
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicit.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicit.py 2007-11-10 00:50:16 UTC (rev 8273)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicit.py 2007-11-10 05:05:10 UTC (rev 8274)
@@ -50,7 +50,7 @@
K = self._calculateStiffnessMat()
M = self._calculateMassMat()
- dispResult = -self.fieldTpdt + 2.0*self.fieldT - self.fieldTmdt
+ dispResult = 2.0*self.fieldT - self.fieldTmdt
self.valsResidual = 1.0/self.dt**2 * numpy.dot(M, dispResult) - \
numpy.dot(K, self.fieldT)
return
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicitData1DLinear.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicitData1DLinear.cc 2007-11-10 00:50:16 UTC (rev 8273)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicitData1DLinear.cc 2007-11-10 05:05:10 UTC (rev 8274)
@@ -85,8 +85,8 @@
};
const double pylith::feassemble::ElasticityExplicitData1DLinear::_valsResidual[] = {
- 1.60000000e+10,
- -1.60000000e+10,
+ 1.60407812e+10,
+ -1.59592188e+10,
};
const double pylith::feassemble::ElasticityExplicitData1DLinear::_valsJacobian[] = {
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicitData1DQuadratic.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicitData1DQuadratic.cc 2007-11-10 00:50:16 UTC (rev 8273)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicitData1DQuadratic.cc 2007-11-10 05:05:10 UTC (rev 8274)
@@ -99,9 +99,9 @@
};
const double pylith::feassemble::ElasticityExplicitData1DQuadratic::_valsResidual[] = {
- -1.11999375e+11,
- 2.56002500e+11,
- -1.43999375e+11,
+ -1.11997188e+11,
+ 2.56020625e+11,
+ -1.43992500e+11,
};
const double pylith::feassemble::ElasticityExplicitData1DQuadratic::_valsJacobian[] = {
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicitData2DLinear.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicitData2DLinear.cc 2007-11-10 00:50:16 UTC (rev 8273)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicitData2DLinear.cc 2007-11-10 05:05:10 UTC (rev 8274)
@@ -90,9 +90,9 @@
};
const double pylith::feassemble::ElasticityExplicitData2DLinear::_valsResidual[] = {
- -6.72542399e+10, 1.24660734e+11,
- -9.57242172e+09, -1.09884721e+11,
- 7.68275783e+10, -1.47847210e+10,
+ -6.72493510e+10, 1.24660275e+11,
+ -9.56753283e+09, -1.09885179e+11,
+ 7.68324672e+10, -1.47851793e+10,
};
const double pylith::feassemble::ElasticityExplicitData2DLinear::_valsJacobian[] = {
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicitData2DQuadratic.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicitData2DQuadratic.cc 2007-11-10 00:50:16 UTC (rev 8273)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicitData2DQuadratic.cc 2007-11-10 05:05:10 UTC (rev 8274)
@@ -130,12 +130,12 @@
};
const double pylith::feassemble::ElasticityExplicitData2DQuadratic::_valsResidual[] = {
- 2.04455943e+09, 9.48741748e+09,
- -1.86265045e+10, 3.65998688e+09,
- -2.49512114e+09, -3.26311648e+10,
- 1.39608918e+10, 6.18995066e+10,
- -1.90442781e+10, -8.28253071e+10,
- 2.41727473e+10, 4.04175780e+10,
+ 2.04288333e+09, 9.48429385e+09,
+ -1.86251338e+10, 3.66245172e+09,
+ -2.49481575e+09, -3.26305061e+10,
+ 1.39575797e+10, 6.18963488e+10,
+ -1.90456087e+10, -8.28246824e+10,
+ 2.41683699e+10, 4.04126141e+10,
};
const double pylith::feassemble::ElasticityExplicitData2DQuadratic::_valsJacobian[] = {
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicitData3DLinear.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicitData3DLinear.cc 2007-11-10 00:50:16 UTC (rev 8273)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicitData3DLinear.cc 2007-11-10 05:05:10 UTC (rev 8274)
@@ -96,10 +96,10 @@
};
const double pylith::feassemble::ElasticityExplicitData3DLinear::_valsResidual[] = {
- -6.07523928e+09, 3.62541599e+10, 3.19640209e+09,
- -4.00838674e+09, 6.66858204e+10, 2.19195434e+10,
- 6.67438897e+09, -1.05592494e+11, -3.14998804e+10,
- 3.41360830e+09, 2.65453164e+09, 6.38292625e+09,
+ -6.07557553e+09, 3.62544121e+10, 3.19547740e+09,
+ -4.00872299e+09, 6.66860725e+10, 2.19186187e+10,
+ 6.67405272e+09, -1.05592242e+11, -3.15008051e+10,
+ 3.41327205e+09, 2.65478382e+09, 6.38200156e+09,
};
const double pylith::feassemble::ElasticityExplicitData3DLinear::_valsJacobian[] = {
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicitData3DQuadratic.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicitData3DQuadratic.cc 2007-11-10 00:50:16 UTC (rev 8273)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicitData3DQuadratic.cc 2007-11-10 05:05:10 UTC (rev 8274)
@@ -177,16 +177,16 @@
};
const double pylith::feassemble::ElasticityExplicitData3DQuadratic::_valsResidual[] = {
- 1.21810568e+11, -3.29478688e+10, 2.43738946e+11,
- -2.90354241e+10, -7.34518237e+10, 2.34529293e+10,
- -2.49211185e+11, -4.63822584e+11, 2.54633336e+11,
- 3.54901735e+10, -5.62151785e+10, 3.13375453e+11,
- -1.24954972e+11, -1.83164844e+11, -1.74376488e+11,
- 2.33420251e+11, 5.42368499e+11, -3.17019958e+11,
- -9.29980081e+11, -4.08320625e+11, -2.12610662e+11,
- 1.17953752e+11, 4.72743819e+11, -2.44313374e+11,
- 2.27857411e+11, 6.71949765e+11, -5.10996159e+11,
- 5.96615974e+11, -4.69104960e+11, 6.24134995e+11,
+ 1.21815001e+11, -3.29529855e+10, 2.43735803e+11,
+ -2.90438333e+10, -7.34446206e+10, 2.34534559e+10,
+ -2.49209347e+11, -4.63827131e+11, 2.54640449e+11,
+ 3.54857040e+10, -5.62185336e+10, 3.13364516e+11,
+ -1.24960815e+11, -1.83167194e+11, -1.74381553e+11,
+ 2.33419532e+11, 5.42360274e+11, -3.17021730e+11,
+ -9.29983954e+11, -4.08328253e+11, -2.12621459e+11,
+ 1.17946612e+11, 4.72741754e+11, -2.44313311e+11,
+ 2.27852241e+11, 6.71942421e+11, -5.11001827e+11,
+ 5.96605680e+11, -4.69106429e+11, 6.24126033e+11,
};
const double pylith::feassemble::ElasticityExplicitData3DQuadratic::_valsJacobian[] = {
More information about the cig-commits
mailing list