[cig-commits] r15224 - short/3D/PyLith/trunk/libsrc/materials
willic3 at geodynamics.org
willic3 at geodynamics.org
Fri Jun 12 18:52:36 PDT 2009
Author: willic3
Date: 2009-06-12 18:52:36 -0700 (Fri, 12 Jun 2009)
New Revision: 15224
Modified:
short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc
Log:
Fixed some minor bugs, plus fixed indexing problems with
_calcElasticConstsTimeDep and fixed updating of viscous strain.
Modified: short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc 2009-06-13 00:47:23 UTC (rev 15223)
+++ short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc 2009-06-13 01:52:36 UTC (rev 15224)
@@ -506,18 +506,18 @@
const double lambda = properties[p_lambda];
const double viscosityCoeff = properties[p_viscosityCoeff];
const double powerLawExp = properties[p_powerLawExponent];
- const double visStrainT[] = {stateVars[s_viscousStrain ],
- stateVars[s_viscousStrain+1],
- stateVars[s_viscousStrain+2],
- stateVars[s_viscousStrain+3],
- stateVars[s_viscousStrain+4],
- stateVars[s_viscousStrain+5]};
- const double stressT[] = {stateVars[s_stress ],
- stateVars[s_stress+1],
- stateVars[s_stress+2],
- stateVars[s_stress+3],
- stateVars[s_stress+4],
- stateVars[s_stress+5]};
+ const double visStrainT[] = {stateVars[s_viscousStrain],
+ stateVars[s_viscousStrain + 1],
+ stateVars[s_viscousStrain + 2],
+ stateVars[s_viscousStrain + 3],
+ stateVars[s_viscousStrain + 4],
+ stateVars[s_viscousStrain + 5]};
+ const double stressT[] = {stateVars[s_stress],
+ stateVars[s_stress + 1],
+ stateVars[s_stress + 2],
+ stateVars[s_stress + 3],
+ stateVars[s_stress + 4],
+ stateVars[s_stress + 5]};
const double mu2 = 2.0 * mu;
const double lamPlusMu = lambda + mu;
@@ -533,9 +533,9 @@
const double timeFac = _dt * (1.0 - alpha);
// Initial stress values
- const double meanStressInitial = ( initialStress[0] +
- initialStress[1] +
- initialStress[2] )/3.0;
+ const double meanStressInitial = (initialStress[0] +
+ initialStress[1] +
+ initialStress[2])/3.0;
const double devStressInitial[] = { initialStress[0] - meanStressInitial,
initialStress[1] - meanStressInitial,
initialStress[2] - meanStressInitial,
@@ -546,9 +546,9 @@
_scalarProduct(devStressInitial, devStressInitial);
// Initial strain values
- const double meanStrainInitial = ( initialStrain[0] +
- initialStrain[1] +
- initialStrain[2] )/3.0;
+ const double meanStrainInitial = (initialStrain[0] +
+ initialStrain[1] +
+ initialStrain[2])/3.0;
// Values for current time step
const double e11 = totalStrain[0];
@@ -570,9 +570,9 @@
_scalarProduct(strainPPTpdt, strainPPTpdt);
// Values for previous time step
- const double meanStressT = ( stressT[0] +
- stressT[1] +
- stressT[2] )/3.0;
+ const double meanStressT = (stressT[0] +
+ stressT[1] +
+ stressT[2])/3.0;
const double devStressT[] = { stressT[0] - meanStressT,
stressT[1] - meanStressT,
stressT[2] - meanStressT,
@@ -654,7 +654,7 @@
const double factor1 = 1.0-alpha;
const double effStressTau = factor1 * effStressT + alpha * effStressTpdt;
const double gammaTau = 0.5 * pow((effStressTau/viscosityCoeff),
- (powerLawExp - 1.0))/viscosityCoeff;
+ (powerLawExp - 1.0))/viscosityCoeff;
const double a = ae + alpha * dt * gammaTau;
const double y = a * a * effStressTpdt * effStressTpdt - b +
c * gammaTau - d * d * gammaTau * gammaTau;
@@ -681,7 +681,7 @@
const double factor1 = 1.0-alpha;
const double effStressTau = factor1 * effStressT + alpha * effStressTpdt;
const double gammaTau = 0.5 * pow((effStressTau/viscosityCoeff),
- (powerLawExp - 1.0))/viscosityCoeff;
+ (powerLawExp - 1.0))/viscosityCoeff;
const double a = ae + alpha * dt * gammaTau;
const double dGammaTau = 0.5 * alpha * (powerLawExp - 1.0) *
pow((effStressTau/viscosityCoeff), (powerLawExp - 2.0))/
@@ -717,7 +717,7 @@
const double factor1 = 1.0-alpha;
const double effStressTau = factor1 * effStressT + alpha * effStressTpdt;
const double gammaTau = 0.5 * pow((effStressTau/viscosityCoeff),
- (powerLawExp - 1.0))/viscosityCoeff;
+ (powerLawExp - 1.0))/viscosityCoeff;
const double dGammaTau = 0.5 * alpha * (powerLawExp - 1.0) *
pow((effStressTau/viscosityCoeff), (powerLawExp - 2.0))/
(viscosityCoeff * viscosityCoeff);
@@ -730,7 +730,7 @@
dGammaTau *
(2.0 * a * alpha * dt * effStressTpdt * effStressTpdt +
c - 2.0 * d * d * gammaTau);
-
+
*func = y;
*dfunc = dy;
@@ -922,18 +922,18 @@
const double viscosityCoeff = properties[p_viscosityCoeff];
const double powerLawExp = properties[p_powerLawExponent];
- const double visStrainT[] = {stateVars[s_viscousStrain ],
- stateVars[s_viscousStrain+1],
- stateVars[s_viscousStrain+2],
- stateVars[s_viscousStrain+3],
- stateVars[s_viscousStrain+4],
- stateVars[s_viscousStrain+5]};
- const double stressT[] = {stateVars[s_stress ],
- stateVars[s_stress+1],
- stateVars[s_stress+2],
- stateVars[s_stress+3],
- stateVars[s_stress+4],
- stateVars[s_stress+5]};
+ const double visStrainT[] = {stateVars[s_viscousStrain],
+ stateVars[s_viscousStrain + 1],
+ stateVars[s_viscousStrain + 2],
+ stateVars[s_viscousStrain + 3],
+ stateVars[s_viscousStrain + 4],
+ stateVars[s_viscousStrain + 5]};
+ const double stressT[] = {stateVars[s_stress],
+ stateVars[s_stress + 1],
+ stateVars[s_stress + 2],
+ stateVars[s_stress + 3],
+ stateVars[s_stress + 4],
+ stateVars[s_stress + 5]};
const double mu2 = 2.0 * mu;
const double lamPlusMu = lambda + mu;
@@ -950,9 +950,9 @@
/// Initial state.
// Initial stress values.
- const double meanStressInitial = ( initialStress[0] +
- initialStress[1] +
- initialStress[2] ) / 3.0;
+ const double meanStressInitial = (initialStress[0] +
+ initialStress[1] +
+ initialStress[2])/3.0;
const double devStressInitial[] = { initialStress[0] - meanStressInitial,
initialStress[1] - meanStressInitial,
initialStress[2] - meanStressInitial,
@@ -963,9 +963,9 @@
_scalarProduct(devStressInitial, devStressInitial);
// Initial strain values.
- const double meanStrainInitial = ( initialStrain[0] +
- initialStrain[1] +
- initialStrain[2] ) / 3.0;
+ const double meanStrainInitial = (initialStrain[0] +
+ initialStrain[1] +
+ initialStrain[2])/3.0;
/// Values for current time step
const double e11 = totalStrain[0];
@@ -988,7 +988,7 @@
_scalarProduct(strainPPTpdt, strainPPTpdt);
// Values for previous time step
- const double meanStressT = ( stressT[0] + stressT[1] + stressT[2] ) / 3.0;
+ const double meanStressT = (stressT[0] + stressT[1] + stressT[2])/3.0;
const double devStressT[] = { stressT[0] - meanStressT,
stressT[1] - meanStressT,
stressT[2] - meanStressT,
@@ -1131,20 +1131,20 @@
- dStress11dStrain22
- dStress11dStrain33)/3.0; // C1111
elasticConsts[ 1] = bulkModulus + (2.0 * dStress11dStrain22
- - dStress11dStrain33
- - dStress11dStrain11)/3.0; // C1122
+ - dStress11dStrain11
+ - dStress11dStrain33)/3.0; // C1122
elasticConsts[ 2] = bulkModulus + (2.0 * dStress11dStrain33
- - dStress11dStrain11
- - dStress11dStrain22)/3.0; // C1133
+ - dStress11dStrain22
+ - dStress11dStrain11)/3.0; // C1133
elasticConsts[ 3] = dStress11dStrain12; // C1112
elasticConsts[ 4] = dStress11dStrain23; // C1123
elasticConsts[ 5] = dStress11dStrain13; // C1113
elasticConsts[ 6] = bulkModulus + (2.0 * dStress22dStrain22
- - dStress22dStrain33
- - dStress22dStrain11)/3.0; // C2222
+ - dStress22dStrain11
+ - dStress22dStrain33)/3.0; // C2222
elasticConsts[ 7] = bulkModulus + (2.0 * dStress22dStrain33
- - dStress22dStrain11
- - dStress22dStrain22)/3.0; // C2233
+ - dStress22dStrain22
+ - dStress22dStrain11)/3.0; // C2233
elasticConsts[ 8] = dStress22dStrain12; // C2212
elasticConsts[ 9] = dStress22dStrain23; // C2223
elasticConsts[10] = dStress22dStrain13; // C2213
@@ -1354,6 +1354,7 @@
const double factor2 = timeFac * gammaTau;
double devStressTpdt = 0.0;
double devStressTau = 0.0;
+ double deltaVisStrain = 0.0;
for (int iComp=0; iComp < _tensorSize; ++iComp) {
devStressTpdt = factor1 *
@@ -1362,11 +1363,12 @@
stateVars[s_stress+iComp] = devStressTpdt + diag[iComp] *
(meanStressTpdt + meanStressInitial);
devStressTau = (1.0 - alpha) * devStressT[iComp] + alpha * devStressTpdt;
- stateVars[s_viscousStrain+iComp] += _dt * gammaTau * devStressTau;
+ deltaVisStrain = _dt * gammaTau * devStressTau;
+ stateVars[s_viscousStrain+iComp] = visStrainT[iComp] + deltaVisStrain;
} // for
_needNewJacobian = true;
- PetscLogFlops(14 + _tensorSize * 14);
+ PetscLogFlops(14 + _tensorSize * 15);
} // _updateStateVarsViscoelastic
More information about the CIG-COMMITS
mailing list