[cig-commits] r14929 - in short/3D/PyLith/trunk: libsrc/feassemble pylith/problems unittests/libtests/feassemble unittests/libtests/feassemble/data unittests/pytests/topology
brad at geodynamics.org
brad at geodynamics.org
Fri May 8 12:39:09 PDT 2009
Author: brad
Date: 2009-05-08 12:39:08 -0700 (Fri, 08 May 2009)
New Revision: 14929
Modified:
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc
short/3D/PyLith/trunk/pylith/problems/Explicit.py
short/3D/PyLith/trunk/pylith/problems/Implicit.py
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc
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
short/3D/PyLith/trunk/unittests/pytests/topology/TestMeshField.py
short/3D/PyLith/trunk/unittests/pytests/topology/TestSolutionFields.py
Log:
Worked on updating explicit time stepping to use increment in displacement. Added unit test to expose bug in binding for Field::operator+=.
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc 2009-05-08 18:56:06 UTC (rev 14928)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityExplicit.cc 2009-05-08 19:39:08 UTC (rev 14929)
@@ -235,7 +235,7 @@
const double valIJ = valI * basis[iQuad*numBasis+jBasis];
for (int iDim=0; iDim < spaceDim; ++iDim)
_cellVector[iBasis*spaceDim+iDim] +=
- valIJ * (+ 2.0 * dispTCell[jBasis*spaceDim+iDim]
+ valIJ * (dispTCell[jBasis*spaceDim+iDim]
- dispTmdtCell[jBasis*spaceDim+iDim]);
} // for
} // for
Modified: short/3D/PyLith/trunk/pylith/problems/Explicit.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Explicit.py 2009-05-08 18:56:06 UTC (rev 14928)
+++ short/3D/PyLith/trunk/pylith/problems/Explicit.py 2009-05-08 19:39:08 UTC (rev 14929)
@@ -70,25 +70,26 @@
Formulation.initialize(self, dimension, normalizer)
self._info.log("Creating other fields and matrices.")
- self.fields.add("disp(t+dt)", "displacement")
+ self.fields.add("dispIncr(t->t+dt)", "displacement_increment")
self.fields.add("disp(t-dt)", "displacement")
self.fields.add("residual", "residual")
- self.fields.createHistory(["disp(t+dt)", "disp(t)", "disp(t-dt)"])
self.fields.copyLayout("disp(t)")
- self.fields.solveSolnName("disp(t+dt)")
+ self.fields.solveSolnName("dispIncr(t->t+dt)")
self._debug.log(resourceUsageString())
- # Create Petsc vectors for fields involved in solve. Since we
- # shift fields through the time history, all fields need a PETSc
- # vector.
- dispTpdt = self.fields.get("disp(t+dt)")
- dispTpdt.createVector()
+ # Set fields to zero
+ dispIncr = self.fields.get("dispIncr(t->t+dt)")
+ dispIncr.zero()
dispT = self.fields.get("disp(t)")
- dispT.createVector()
+ dispT.zero()
dispTmdt = self.fields.get("disp(t-dt)")
- dispTmdt.createVector()
+ dispTmdt.zero()
residual = self.fields.get("residual")
+ residual.zero()
+ # Create Petsc vectors for fields involved in solve
residual.createVector()
+ dispIncr.createVector()
+ self._debug.log(resourceUsageString())
self._info.log("Creating Jacobian matrix.")
from pylith.topology.Jacobian import Jacobian
@@ -100,11 +101,11 @@
self.solver.initialize(self.fields, self.jacobian, self)
self._debug.log(resourceUsageString())
- # Solve for total displacement field
+ # Solve for increment in displacement field.
for constraint in self.constraints:
- constraint.useSolnIncr(False)
+ constraint.useSolnIncr(True)
for integrator in self.integratorsMesh + self.integratorsSubMesh:
- integrator.useSolnIncr(False)
+ integrator.useSolnIncr(True)
self._logger.eventEnd(logEvent)
return
@@ -119,7 +120,7 @@
dispTpdt = self.fields.get("disp(t+dt)")
for constraint in self.constraints:
- constraint.setField(t+dt, dispTpdt)
+ constraint.setFieldIncr(t, t+dt, dispTpdt)
needNewJacobian = False
for integrator in self.integratorsMesh + self.integratorsSubMesh:
@@ -144,8 +145,8 @@
self._info.log("Solving equations.")
residual = self.fields.get("residual")
- dispTpdt = self.fields.solveSoln()
- self.solver.solve(dispTpdt, self.jacobian, residual)
+ dispIncr = self.fields.get("dispIncr(t->t+dt")
+ self.solver.solve(dispIncr, self.jacobian, residual)
self._logger.eventEnd(logEvent)
return
@@ -158,15 +159,13 @@
logEvent = "%spoststep" % self._loggingPrefix
self._logger.eventBegin(logEvent)
- self.fields.shiftHistory()
+ dispIncr = self.fields.get("dispIncr(t->t+dt)")
+ disp = self.fields.get("disp(t)")
+ dispTmdt = self.fields.get("disp(t-dt)")
- # :KLUDGE: only works for KSP solver
- dispTpdt = self.fields.get("disp(t+dt)")
- if not self.solver.guessZero:
- dispT = self.fields.get("disp(t)")
- dispTpdt.copy(dispT)
- else:
- dispTpdt.zero()
+ dispTmdt.copy(dispT)
+ disp += dispIncr
+ dispIncr.zero()
Formulation.poststep(self, t, dt)
Modified: short/3D/PyLith/trunk/pylith/problems/Implicit.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Implicit.py 2009-05-08 18:56:06 UTC (rev 14928)
+++ short/3D/PyLith/trunk/pylith/problems/Implicit.py 2009-05-08 19:39:08 UTC (rev 14929)
@@ -148,19 +148,9 @@
logEvent = "%sprestep" % self._loggingPrefix
self._logger.eventBegin(logEvent)
- # Set dispT to the BC t time t+dt. Unconstrained DOF are
- # unaffected and will be equal to their values at time t.
- self._info.log("Setting constraints.")
- dispIncr = self.fields.get("dispIncr(t->t+dt)")
- if 0 == self._stepCount:
- for constraint in self.constraints:
- constraint.setField(t+dt, dispIncr)
- else:
- for constraint in self.constraints:
- constraint.setFieldIncr(t, t+dt, dispIncr)
-
# If finishing first time step, then switch from solving for total
# displacements to solving for incremental displacements
+ needNewJacobian = False
if 1 == self._stepCount:
self._info.log("Switching from total field solution to incremental " \
"field solution.")
@@ -168,10 +158,18 @@
constraint.useSolnIncr(True)
for integrator in self.integratorsMesh + self.integratorsSubMesh:
integrator.useSolnIncr(True)
- self._reformJacobian(t, dt)
+ needNewJacobian = True
+ self._info.log("Setting constraints.")
+ dispIncr = self.fields.get("dispIncr(t->t+dt)")
+ if 0 == self._stepCount:
+ for constraint in self.constraints:
+ constraint.setField(t+dt, dispIncr)
+ else:
+ for constraint in self.constraints:
+ constraint.setFieldIncr(t, t+dt, dispIncr)
+
### NONLINEAR: Might want to move logic into IntegrateJacobian() and set a flag instead
- needNewJacobian = False
for integrator in self.integratorsMesh + self.integratorsSubMesh:
integrator.timeStep(dt)
if integrator.needNewJacobian():
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc 2009-05-08 18:56:06 UTC (rev 14928)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc 2009-05-08 19:39:08 UTC (rev 14929)
@@ -366,11 +366,11 @@
// Setup fields
CPPUNIT_ASSERT(0 != fields);
fields->add("residual", "residual");
- fields->add("disp(t+dt)", "displacement");
+ fields->add("dispIncr(t->t+dt)", "displacement_increment");
fields->add("disp(t)", "displacement");
fields->add("disp(t-dt)", "displacement");
- fields->solutionName("disp(t+dt)");
- const char* history[] = { "disp(t+dt)", "disp(t)", "disp(t-dt)" };
+ fields->solutionName("dispIncr(t->t+dt)");
+ const char* history[] = { "dispIncr(t->t+dt)", "disp(t)", "disp(t-dt)" };
const int historySize = 3;
fields->createHistory(history, historySize);
@@ -381,18 +381,18 @@
fields->copyLayout("residual");
const int fieldSize = _data->spaceDim * _data->numVertices;
- topology::Field<topology::Mesh>& dispTpdt = fields->get("disp(t+dt)");
+ topology::Field<topology::Mesh>& dispIncr = fields->get("dispIncr(t->t+dt)");
topology::Field<topology::Mesh>& dispT = fields->get("disp(t)");
topology::Field<topology::Mesh>& dispTmdt = fields->get("disp(t-dt)");
- const ALE::Obj<RealSection>& dispTpdtSection = dispTpdt.section();
+ const ALE::Obj<RealSection>& dispIncrSection = dispIncr.section();
const ALE::Obj<RealSection>& dispTSection = dispT.section();
const ALE::Obj<RealSection>& dispTmdtSection = dispTmdt.section();
- CPPUNIT_ASSERT(!dispTpdtSection.isNull());
+ CPPUNIT_ASSERT(!dispIncrSection.isNull());
CPPUNIT_ASSERT(!dispTSection.isNull());
CPPUNIT_ASSERT(!dispTmdtSection.isNull());
const int offset = _data->numCells;
for (int iVertex=0; iVertex < _data->numVertices; ++iVertex) {
- dispTpdtSection->updatePoint(iVertex+offset,
+ dispIncrSection->updatePoint(iVertex+offset,
&_data->fieldTIncr[iVertex*_data->spaceDim]);
dispTSection->updatePoint(iVertex+offset,
&_data->fieldT[iVertex*_data->spaceDim]);
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicit.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicit.py 2009-05-08 18:56:06 UTC (rev 14928)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicit.py 2009-05-08 19:39:08 UTC (rev 14929)
@@ -50,7 +50,7 @@
K = self._calculateStiffnessMat()
M = self._calculateMassMat()
- dispResult = 2.0*self.fieldT - self.fieldTmdt
+ dispResult = 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 2009-05-08 18:56:06 UTC (rev 14928)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicitData1DLinear.cc 2009-05-08 19:39:08 UTC (rev 14929)
@@ -89,8 +89,8 @@
};
const double pylith::feassemble::ElasticityExplicitData1DLinear::_valsResidual[] = {
- 1.60407812e+10,
- -1.59592188e+10,
+ 1.60042188e+10,
+ -1.59957812e+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 2009-05-08 18:56:06 UTC (rev 14928)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicitData1DQuadratic.cc 2009-05-08 19:39:08 UTC (rev 14929)
@@ -103,9 +103,9 @@
};
const double pylith::feassemble::ElasticityExplicitData1DQuadratic::_valsResidual[] = {
- 5.60178125e+10,
- 1.36007500e+11,
- -1.91949375e+11,
+ 5.60018750e+10,
+ 1.36000938e+11,
+ -1.91994375e+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 2009-05-08 18:56:06 UTC (rev 14928)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicitData2DLinear.cc 2009-05-08 19:39:08 UTC (rev 14929)
@@ -94,9 +94,9 @@
};
const double pylith::feassemble::ElasticityExplicitData2DLinear::_valsResidual[] = {
- -1.66790376e+10, 3.09147184e+10,
- -2.30485101e+09, -2.72112551e+10,
- 1.89834302e+10, -3.70266130e+09,
+ -1.66791140e+10, 3.09147184e+10,
+ -2.30469823e+09, -2.72112551e+10,
+ 1.89835830e+10, -3.70266130e+09,
};
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 2009-05-08 18:56:06 UTC (rev 14928)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicitData2DQuadratic.cc 2009-05-08 19:39:08 UTC (rev 14929)
@@ -134,12 +134,12 @@
};
const double pylith::feassemble::ElasticityExplicitData2DQuadratic::_valsResidual[] = {
- -5.90190847e+09, 4.71767000e+10,
- -1.39395410e+10, -6.36899121e+09,
- 4.56827962e+08, 1.65078960e+10,
- 2.42224578e+10, -8.30513314e+09,
- -2.84017435e+10, -6.85808242e+10,
- 2.35658760e+10, 1.95703525e+10,
+ -5.90114822e+09, 4.71777415e+10,
+ -1.39396421e+10, -6.36913623e+09,
+ 4.56555501e+08, 1.65074214e+10,
+ 2.42225720e+10, -8.30477279e+09,
+ -2.84015238e+10, -6.85803760e+10,
+ 2.35671240e+10, 1.95727783e+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 2009-05-08 18:56:06 UTC (rev 14928)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicitData3DLinear.cc 2009-05-08 19:39:08 UTC (rev 14929)
@@ -100,10 +100,10 @@
};
const double pylith::feassemble::ElasticityExplicitData3DLinear::_valsResidual[] = {
- -6.07456678e+09, 3.62530671e+10, 3.19631803e+09,
- -4.01124487e+09, 6.66859044e+10, 2.19201318e+10,
- 6.67153085e+09, -1.05592410e+11, -3.14992920e+10,
- 3.41075018e+09, 2.65461570e+09, 6.38351469e+09,
+ -6.07565959e+09, 3.62534033e+10, 3.19640209e+09,
+ -4.01015205e+09, 6.66855682e+10, 2.19200477e+10,
+ 6.67262366e+09, -1.05592747e+11, -3.14993761e+10,
+ 3.41184299e+09, 2.65427945e+09, 6.38343063e+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 2009-05-08 18:56:06 UTC (rev 14928)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicitData3DQuadratic.cc 2009-05-08 19:39:08 UTC (rev 14929)
@@ -181,16 +181,16 @@
};
const double pylith::feassemble::ElasticityExplicitData3DQuadratic::_valsResidual[] = {
- 2.17455103e+10, -9.08324612e+09, 2.06925719e+10,
- 6.22334143e+10, -2.47874409e+10, 7.42815158e+09,
- -4.60062771e+10, -5.34075482e+10, -3.15147885e+10,
- -1.00608531e+10, 4.19820053e+10, -3.87394106e+10,
- -8.43421071e+09, 6.33100433e+10, -3.41569760e+09,
- 6.62252212e+10, 1.09238467e+11, -3.72823741e+10,
- -9.59455989e+10, -5.29509634e+10, -2.51207822e+10,
- -1.25344031e+10, -6.06788721e+10, 1.03667465e+10,
- -3.24507112e+10, 5.14071214e+10, 7.32766003e+09,
- 5.52057820e+10, -6.50277434e+10, 9.02597104e+10,
+ 2.17446982e+10, -9.08226289e+09, 2.06945299e+10,
+ 6.22338118e+10, -2.47879968e+10, 7.42326275e+09,
+ -4.60076154e+10, -5.34101030e+10, -3.15162055e+10,
+ -1.00555635e+10, 4.19850265e+10, -3.87348834e+10,
+ -8.43197122e+09, 6.33093784e+10, -3.41846873e+09,
+ 6.62258554e+10, 1.09239192e+11, -3.72807133e+10,
+ -9.59444688e+10, -5.29496681e+10, -2.51201133e+10,
+ -1.25318752e+10, -6.06765547e+10, 1.03701057e+10,
+ -3.24446843e+10, 5.14096428e+10, 7.32828558e+09,
+ 5.52108169e+10, -6.50263642e+10, 9.02623199e+10,
};
const double pylith::feassemble::ElasticityExplicitData3DQuadratic::_valsJacobian[] = {
Modified: short/3D/PyLith/trunk/unittests/pytests/topology/TestMeshField.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/topology/TestMeshField.py 2009-05-08 18:56:06 UTC (rev 14928)
+++ short/3D/PyLith/trunk/unittests/pytests/topology/TestMeshField.py 2009-05-08 19:39:08 UTC (rev 14929)
@@ -147,13 +147,12 @@
# No test of result
return
-
def test_operatorAdd(self):
"""
Test newSection(field).
"""
fieldB = MeshField(self.mesh)
- fieldB += self.field
+ self.field += fieldB
# No test of result
return
Modified: short/3D/PyLith/trunk/unittests/pytests/topology/TestSolutionFields.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/topology/TestSolutionFields.py 2009-05-08 18:56:06 UTC (rev 14928)
+++ short/3D/PyLith/trunk/unittests/pytests/topology/TestSolutionFields.py 2009-05-08 19:39:08 UTC (rev 14929)
@@ -109,4 +109,27 @@
return
+ def test_fieldAdd(self):
+ """
+ Test shiftHistory().
+ """
+ fields = self.fields
+ fields.add("field A", "A");
+ fields.add("field B", "B");
+
+ helper_fieldAdd(fields)
+ fieldA = fields.get("field A")
+ fieldB = fields.get("field B")
+ fieldA.copy(fieldB)
+ return
+
+
+def helper_fieldAdd(fields):
+ fieldA = fields.get("field A")
+ fieldB = fields.get("field B")
+
+ fieldA += fieldB
+ return
+
+
# End of file
More information about the CIG-COMMITS
mailing list