[cig-commits] r16198 - in short/3D/PyLith/trunk: examples/bar_shearwave/hex8 examples/bar_shearwave/quad4 examples/bar_shearwave/tet4 examples/bar_shearwave/tri3 libsrc/faults unittests/libtests/bc unittests/libtests/faults unittests/libtests/faults/data unittests/libtests/feassemble
brad at geodynamics.org
brad at geodynamics.org
Fri Jan 29 16:40:48 PST 2010
Author: brad
Date: 2010-01-29 16:40:47 -0800 (Fri, 29 Jan 2010)
New Revision: 16198
Modified:
short/3D/PyLith/trunk/examples/bar_shearwave/hex8/pylithapp.cfg
short/3D/PyLith/trunk/examples/bar_shearwave/quad4/pylithapp.cfg
short/3D/PyLith/trunk/examples/bar_shearwave/tet4/pylithapp.cfg
short/3D/PyLith/trunk/examples/bar_shearwave/tri3/pylithapp.cfg
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc
short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataHex8.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataLine2.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataQuad4.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTet4.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTri3.cc
short/3D/PyLith/trunk/unittests/libtests/faults/data/adjustsoln.py
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc
Log:
Fixed some lumped solver bugs (incorrect signs in adjustSolnLumped().
Modified: short/3D/PyLith/trunk/examples/bar_shearwave/hex8/pylithapp.cfg
===================================================================
--- short/3D/PyLith/trunk/examples/bar_shearwave/hex8/pylithapp.cfg 2010-01-29 22:38:01 UTC (rev 16197)
+++ short/3D/PyLith/trunk/examples/bar_shearwave/hex8/pylithapp.cfg 2010-01-30 00:40:47 UTC (rev 16198)
@@ -40,7 +40,8 @@
dimension = 3
# Change to an explicit time stepping formulation
-formulation = pylith.problems.Explicit
+#formulation = pylith.problems.Explicit
+formulation = pylith.problems.ExplicitLumped
# Nondimensionalize problem using wave propagation parameters.
normalizer = spatialdata.units.NondimElasticDynamic
Modified: short/3D/PyLith/trunk/examples/bar_shearwave/quad4/pylithapp.cfg
===================================================================
--- short/3D/PyLith/trunk/examples/bar_shearwave/quad4/pylithapp.cfg 2010-01-29 22:38:01 UTC (rev 16197)
+++ short/3D/PyLith/trunk/examples/bar_shearwave/quad4/pylithapp.cfg 2010-01-30 00:40:47 UTC (rev 16198)
@@ -41,8 +41,8 @@
dimension = 2
# Change to an explicit time stepping formulation
-formulation = pylith.problems.Explicit
-#formulation = pylith.problems.ExplicitLumped
+#formulation = pylith.problems.Explicit
+formulation = pylith.problems.ExplicitLumped
# Nondimensionalize problem using wave propagation parameters.
normalizer = spatialdata.units.NondimElasticDynamic
Modified: short/3D/PyLith/trunk/examples/bar_shearwave/tet4/pylithapp.cfg
===================================================================
--- short/3D/PyLith/trunk/examples/bar_shearwave/tet4/pylithapp.cfg 2010-01-29 22:38:01 UTC (rev 16197)
+++ short/3D/PyLith/trunk/examples/bar_shearwave/tet4/pylithapp.cfg 2010-01-30 00:40:47 UTC (rev 16198)
@@ -41,7 +41,8 @@
dimension = 3
# Change to an explicit time stepping formulation
-formulation = pylith.problems.Explicit
+#formulation = pylith.problems.Explicit
+formulation = pylith.problems.ExplicitLumped
# Nondimensionalize problem using wave propagation parameters.
normalizer = spatialdata.units.NondimElasticDynamic
Modified: short/3D/PyLith/trunk/examples/bar_shearwave/tri3/pylithapp.cfg
===================================================================
--- short/3D/PyLith/trunk/examples/bar_shearwave/tri3/pylithapp.cfg 2010-01-29 22:38:01 UTC (rev 16197)
+++ short/3D/PyLith/trunk/examples/bar_shearwave/tri3/pylithapp.cfg 2010-01-30 00:40:47 UTC (rev 16198)
@@ -41,7 +41,8 @@
dimension = 2
# Change to an explicit time stepping formulation
-formulation = pylith.problems.Explicit
+#formulation = pylith.problems.Explicit
+formulation = pylith.problems.ExplicitLumped
# Nondimensionalize problem using wave propagation parameters.
normalizer = spatialdata.units.NondimElasticDynamic
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2010-01-29 22:38:01 UTC (rev 16197)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc 2010-01-30 00:40:47 UTC (rev 16198)
@@ -708,12 +708,12 @@
// vertex k require 2 adjustments to the solution:
//
// * DOF k: Compute increment in Lagrange multipliers
- // dl_k = 1/(A_i+A_j) (C_ki A_j r_i - C_kj A_i r_j)
- // - A_i A_j / (A_i+A_j) d_k
+ // dl_k = S^{-1} (-C_ki (A_i^{-1} r_i - C_kj A_j^{-1} r_j + u_i - u_j) - d_k)
+ // S = C_ki (A_i^{-1} + A_j^{-1}) C_ki^T
// * DOF i and j: Adjust displacement increment (solution) to account
// for Lagrange multiplier constraints
- // du_i = -A_i^-1 C_ki^T dlk
- // du_j = +A_j^-1 C_kj^T dlk
+ // du_i = +A_i^-1 C_ki^T dlk
+ // du_j = -A_j^-1 C_kj^T dlk
// Get cell information and setup storage for cell data
const int spaceDim = _quadrature->spaceDim();
@@ -895,28 +895,36 @@
switch(spaceDim)
{ // switch
case 1 : {
+ assert(Ai[0] > 0.0);
+ assert(Aj[0] > 0.0);
+
const double Sinv = Ai[0] * Aj[0] / (Ai[0] + Aj[0]);
// Aru = A_i^{-1} r_i - A_j^{-1} r_j + u_i - u_j
const double Aru = ri[0]/Ai[0] - rj[0]/Aj[0] + ui[0] - uj[0];
- // dl_k = D^{-1}( C_{ki} Aru - d_k)
- const double Aruslip = Aru - slipVertex[0];
- const double dlp = Sinv * Aruslip;
+ // dl_k = D^{-1}( d_k - C_{ki} Aru)
+ const double Aruslip = slipVertex[0] - Aru;
+ const double dlp = wt * Sinv * Aruslip;
// Update displacements at node I
- solutionCell[indexI*spaceDim+0] = wt * -1.0/Ai[0] * dlp;
+ solutionCell[indexI*spaceDim+0] = +1.0/Ai[0] * dlp;
// Update displacements at node J
- solutionCell[indexJ*spaceDim+0] = wt * +1.0/Aj[0] * dlp;
+ solutionCell[indexJ*spaceDim+0] = -1.0/Aj[0] * dlp;
// Update Lagrange multipliers
- solutionCell[indexK*spaceDim+0] = wt * dlp;
+ solutionCell[indexK*spaceDim+0] = dlp;
- PetscLogFlops(17);
+ PetscLogFlops(16);
break;
} // case 1
case 2 : {
+ assert(Ai[0] > 0.0);
+ assert(Ai[1] > 0.0);
+ assert(Aj[0] > 0.0);
+ assert(Aj[1] > 0.0);
+
const double Cpx = orientationVertex[0];
const double Cpy = orientationVertex[1];
const double Cqx = orientationVertex[2];
@@ -933,11 +941,11 @@
const double Arux = ri[0]/Ai[0] - rj[0]/Aj[0] + ui[0] - uj[0];
const double Aruy = ri[1]/Ai[1] - rj[1]/Aj[1] + ui[1] - uj[1];
- // dl_k = D^{-1}( C_{ki} Aru - d_k)
+ // dl_k = S^{-1}( d_k - C_{ki} Aru)
const double Arup = Cpx*Arux + Cpy*Aruy;
const double Aruq = Cqx*Arux + Cqy*Aruy;
- const double Arupslip = Arup - slipVertex[0];
- const double Aruqslip = Aruq - slipVertex[1];
+ const double Arupslip = slipVertex[0] - Arup;
+ const double Aruqslip = slipVertex[1] - Aruq;
const double dlp = Sinv * Arupslip;
const double dlq = Sinv * Aruqslip;
@@ -945,12 +953,12 @@
const double dly = wt * (Cpy*dlp + Cqy*dlq);
// Update displacements at node I
- solutionCell[indexI*spaceDim+0] = -dlx / Ai[0];
- solutionCell[indexI*spaceDim+1] = -dly / Ai[1];
+ solutionCell[indexI*spaceDim+0] = dlx / Ai[0];
+ solutionCell[indexI*spaceDim+1] = dly / Ai[1];
// Update displacements at node J
- solutionCell[indexJ*spaceDim+0] = dlx / Aj[0];
- solutionCell[indexJ*spaceDim+1] = dly / Aj[0];
+ solutionCell[indexJ*spaceDim+0] = -dlx / Aj[0];
+ solutionCell[indexJ*spaceDim+1] = -dly / Aj[0];
// Update Lagrange multipliers
solutionCell[indexK*spaceDim+0] = wt * dlp;
@@ -960,6 +968,13 @@
break;
} // case 2
case 3 : {
+ assert(Ai[0] > 0.0);
+ assert(Ai[1] > 0.0);
+ assert(Ai[2] > 0.0);
+ assert(Aj[0] > 0.0);
+ assert(Aj[1] > 0.0);
+ assert(Aj[2] > 0.0);
+
const double Cpx = orientationVertex[0];
const double Cpy = orientationVertex[1];
const double Cpz = orientationVertex[2];
@@ -982,13 +997,13 @@
const double Aruy = ri[1]/Ai[1] - rj[1]/Aj[1] + ui[1] - uj[1];
const double Aruz = ri[2]/Ai[2] - rj[2]/Aj[2] + ui[2] - uj[2];
- // dl_k = D^{-1}( C_{ki} Aru - d_k)
+ // dl_k = D^{-1}( d_k - C_{ki} Aru)
const double Arup = Cpx*Arux + Cpy*Aruy + Cpz*Aruz;
const double Aruq = Cqx*Arux + Cqy*Aruy + Cqz*Aruz;
const double Arur = Crx*Arux + Cry*Aruy + Crz*Aruz;
- const double Arupslip = Arup - slipVertex[0];
- const double Aruqslip = Aruq - slipVertex[1];
- const double Arurslip = Arur - slipVertex[2];
+ const double Arupslip = slipVertex[0] - Arup;
+ const double Aruqslip = slipVertex[1] - Aruq;
+ const double Arurslip = slipVertex[2] - Arur;
const double dlp = Sinv * Arupslip;
const double dlq = Sinv * Aruqslip;
const double dlr = Sinv * Arurslip;
@@ -998,14 +1013,14 @@
const double dlz = wt * (Cpz*dlp + Cqz*dlq + Crz*dlr);
// Update displacements at node I
- solutionCell[indexI*spaceDim+0] = -dlx / Ai[0];
- solutionCell[indexI*spaceDim+1] = -dly / Ai[1];
- solutionCell[indexI*spaceDim+2] = -dlz / Ai[2];
+ solutionCell[indexI*spaceDim+0] = dlx / Ai[0];
+ solutionCell[indexI*spaceDim+1] = dly / Ai[1];
+ solutionCell[indexI*spaceDim+2] = dlz / Ai[2];
// Update displacements at node J
- solutionCell[indexJ*spaceDim+0] = dlx / Aj[0];
- solutionCell[indexJ*spaceDim+1] = dly / Aj[1];
- solutionCell[indexJ*spaceDim+2] = dlz / Aj[2];
+ solutionCell[indexJ*spaceDim+0] = -dlx / Aj[0];
+ solutionCell[indexJ*spaceDim+1] = -dly / Aj[1];
+ solutionCell[indexJ*spaceDim+2] = -dlz / Aj[2];
// Update Lagrange multipliers
solutionCell[indexK*spaceDim+0] = wt * dlp;
Modified: short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc 2010-01-29 22:38:01 UTC (rev 16197)
+++ short/3D/PyLith/trunk/unittests/libtests/bc/TestAbsorbingDampers.cc 2010-01-30 00:40:47 UTC (rev 16198)
@@ -306,7 +306,7 @@
CPPUNIT_ASSERT(!solutionSection.isNull());
const double t = 1.0;
- bc.integrateJacobian(jacobian, t, &fields);
+ bc.integrateJacobian(&jacobian, t, &fields);
CPPUNIT_ASSERT_EQUAL(false, bc.needNewJacobian());
jacobian.complete();
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc 2010-01-29 22:38:01 UTC (rev 16197)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/TestFaultCohesiveKin.cc 2010-01-30 00:40:47 UTC (rev 16198)
@@ -592,7 +592,7 @@
jacobian.allocate();
const double t = 2.134;
- fault.integrateJacobianAssembled(jacobian, t, &fields);
+ fault.integrateJacobianAssembled(&jacobian, t, &fields);
CPPUNIT_ASSERT_EQUAL(false, fault.needNewJacobian());
jacobian.complete();
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataHex8.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataHex8.cc 2010-01-29 22:38:01 UTC (rev 16197)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataHex8.cc 2010-01-30 00:40:47 UTC (rev 16198)
@@ -1455,10 +1455,10 @@
-3.29153423656, 1.17793324732, -0.323110365593, // 15
-2.47367998662, 1.54678442526, 0.338687018672, // 16
-1.27514497151, 2.36269566532, 1.43857366631, // 17
- 8.08563429421, 12.7441237241, 15.3018598761, // 18
- 8.58310012902, 13.3346655484, 15.7873013548, // 19
- 9.16046647706, 13.9142320683, 16.2952559773, // 20
- 7.8208782359, 12.476710034, 14.8327754459, // 21
+ -1.28563429421, -3.94412372412, -4.50185987615, // 18
+ -1.38310012902, -4.13466554839, -4.58730135484, // 19
+ -1.56046647706, -4.31423206826, -4.69525597725, // 20
+ -1.8208782359, -4.47671003401, -4.83277544586, // 21
};
pylith::faults::CohesiveKinDataHex8::CohesiveKinDataHex8(void)
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataLine2.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataLine2.cc 2010-01-29 22:38:01 UTC (rev 16197)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataLine2.cc 2010-01-30 00:40:47 UTC (rev 16198)
@@ -138,7 +138,7 @@
-1.11580527225, // 3
1.3,
3.52282149957, // 5
- 6.59477159896, // 6
+ -3.59477159896, // 6
};
pylith::faults::CohesiveKinDataLine2::CohesiveKinDataLine2(void)
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataQuad4.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataQuad4.cc 2010-01-29 22:38:01 UTC (rev 16197)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataQuad4.cc 2010-01-30 00:40:47 UTC (rev 16198)
@@ -385,8 +385,8 @@
3.6, 4.6,
9.22726159185, 8.9338057688, // 8
8.38975844529, 8.09426880301, // 9
- 10.997469807, 14.1963447061, // 10
- 9.06911072572, 12.530541046, // 11
+ -3.39746980696, -4.59634470614, // 10
+ -3.06911072572, -4.53054104604, // 11
};
pylith::faults::CohesiveKinDataQuad4::CohesiveKinDataQuad4(void)
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTet4.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTet4.cc 2010-01-29 22:38:01 UTC (rev 16197)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTet4.cc 2010-01-30 00:40:47 UTC (rev 16198)
@@ -542,9 +542,9 @@
9.45705113257, 8.45860462062, 11.1042135006, // 7
9.03102479301, 8.26763801792, 10.7719581756, // 8
12.227136455, 10.448359347, 13.6519239684, // 9
- 9.87376739299, 13.5067416009, 15.0712818121, // 10
- 10.1417484323, 13.8495247161, 15.3158446274, // 11
- 9.548359347, 12.7519239684, 14.327136455, // 12
+ -2.47376739299, -4.1067416009, -3.67128181212, // 10
+ -2.34174843225, -4.04952471613, -3.51584462742, // 11
+ -3.348359347, -4.55192396841, -4.12713645498, // 12
};
pylith::faults::CohesiveKinDataTet4::CohesiveKinDataTet4(void)
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTri3.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTri3.cc 2010-01-29 22:38:01 UTC (rev 16197)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/CohesiveKinDataTri3.cc 2010-01-30 00:40:47 UTC (rev 16198)
@@ -306,8 +306,8 @@
3.4, 4.4,
-2.99670600701, -0.524238161213, // 6
-2.17392825851, 0.119527564532, // 7
- 11.1363572418, 14.3450590105, // 8
- 11.5868031403, 14.7856780395, // 9
+ -3.93635724182, -5.14505901051, // 8
+ -3.9868031403, -5.18567803947, // 9
};
pylith::faults::CohesiveKinDataTri3::CohesiveKinDataTri3(void)
Modified: short/3D/PyLith/trunk/unittests/libtests/faults/data/adjustsoln.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/faults/data/adjustsoln.py 2010-01-29 22:38:01 UTC (rev 16197)
+++ short/3D/PyLith/trunk/unittests/libtests/faults/data/adjustsoln.py 2010-01-30 00:40:47 UTC (rev 16198)
@@ -10,11 +10,11 @@
Si = (Ai * Aj) / (Ai + Aj)
Aru = ri/Ai - rj/Aj + ui - uj
- Aruslip = C*Aru - dk
+ Aruslip = dk - C*Aru
dlp = Si * Aruslip
- ddui = -C / Ai * dlp
- dduj = +C / Aj * dlp
+ ddui = +C / Ai * dlp
+ dduj = -C / Aj * dlp
#print "Aru",Aru
#print "Aruslip",Aruslip
@@ -116,17 +116,17 @@
Arup = Cpx*Arux + Cpy*Aruy
Aruq = Cqx*Arux + Cqy*Aruy
- Arupslip = Arup - dkp
- Aruqslip = Aruq - dkq
+ Arupslip = dkp - Arup
+ Aruqslip = dkq - Aruq
dlp = Sppi * Arupslip
dlq = Sqqi * Aruqslip
- dduix = -1.0/Aix * (Cpx*dlp + Cqx*dlq)
- dduiy = -1.0/Aiy * (Cpy*dlp + Cqy*dlq)
+ dduix = +1.0/Aix * (Cpx*dlp + Cqx*dlq)
+ dduiy = +1.0/Aiy * (Cpy*dlp + Cqy*dlq)
- ddujx = +1.0/Ajx * (Cpx*dlp + Cqx*dlq)
- ddujy = +1.0/Ajy * (Cpy*dlp + Cqy*dlq)
+ ddujx = -1.0/Ajx * (Cpx*dlp + Cqx*dlq)
+ ddujy = -1.0/Ajy * (Cpy*dlp + Cqy*dlq)
print "Sppi",Sppi
print "Sqqi",Sqqi
@@ -314,21 +314,21 @@
Arup = Cpx*Arux + Cpy*Aruy + Cpz*Aruz
Aruq = Cqx*Arux + Cqy*Aruy + Cqz*Aruz
Arur = Crx*Arux + Cry*Aruy + Crz*Aruz
- Arupslip = Arup - dkp
- Aruqslip = Aruq - dkq
- Arurslip = Arur - dkr
+ Arupslip = dkp - Arup
+ Aruqslip = dkq - Aruq
+ Arurslip = dkr - Arur
dlp = Sppi * Arupslip
dlq = Sqqi * Aruqslip
dlr = Srri * Arurslip
- dduix = -1.0/Aix * (Cpx*dlp + Cqx*dlq + Crx*dlr)
- dduiy = -1.0/Aiy * (Cpy*dlp + Cqy*dlq + Cry*dlr)
- dduiz = -1.0/Aiy * (Cpz*dlp + Cqz*dlq + Crz*dlr)
+ dduix = +1.0/Aix * (Cpx*dlp + Cqx*dlq + Crx*dlr)
+ dduiy = +1.0/Aiy * (Cpy*dlp + Cqy*dlq + Cry*dlr)
+ dduiz = +1.0/Aiy * (Cpz*dlp + Cqz*dlq + Crz*dlr)
- ddujx = +1.0/Ajx * (Cpx*dlp + Cqx*dlq + Crx*dlr)
- ddujy = +1.0/Ajy * (Cpy*dlp + Cqy*dlq + Cry*dlr)
- ddujz = +1.0/Ajz * (Cpz*dlp + Cqz*dlq + Crz*dlr)
+ ddujx = -1.0/Ajx * (Cpx*dlp + Cqx*dlq + Crx*dlr)
+ ddujy = -1.0/Ajy * (Cpy*dlp + Cqy*dlq + Cry*dlr)
+ ddujz = -1.0/Ajz * (Cpz*dlp + Cqz*dlq + Crz*dlr)
print "Sppi",Sppi
print "Sqqi",Sqqi
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc 2010-01-29 22:38:01 UTC (rev 16197)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc 2010-01-30 00:40:47 UTC (rev 16198)
@@ -274,7 +274,7 @@
jacobian.allocate();
const double t = 1.0;
- integrator.integrateJacobian(jacobian, t, &fields);
+ integrator.integrateJacobian(&jacobian, t, &fields);
CPPUNIT_ASSERT_EQUAL(false, integrator.needNewJacobian());
jacobian.complete();
More information about the CIG-COMMITS
mailing list