[cig-commits] r14923 - short/3D/PyLith/trunk/libsrc/materials
willic3 at geodynamics.org
willic3 at geodynamics.org
Thu May 7 20:56:50 PDT 2009
Author: willic3
Date: 2009-05-07 20:56:50 -0700 (Thu, 07 May 2009)
New Revision: 14923
Modified:
short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc
Log:
Started updating PowerLaw3D for new SWIG version.
Also made some other fixes.
Code is still incomplete.
Modified: short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc 2009-05-08 01:26:58 UTC (rev 14922)
+++ short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc 2009-05-08 03:56:50 UTC (rev 14923)
@@ -30,6 +30,9 @@
namespace materials {
namespace _PowerLaw3D{
+ /// Dimension of material.
+ const int dimension = 3;
+
/// Number of entries in stress/strain tensors.
const int tensorSize = 6;
@@ -37,72 +40,117 @@
const int numElasticConsts = 21;
/// Number of physical properties.
- const int numProperties = 8;
+ const int numProperties = 5;
/// Physical properties.
- // We are including Maxwell time for now, even though it is not used in
- // computations. It is used to determine the stable time step size.
- const Material::PropMetaData properties[] = {
- { "density", 1, OTHER_FIELD },
- { "mu", 1, OTHER_FIELD },
- { "lambda", 1, OTHER_FIELD },
- { "viscosity_coeff", 1, OTHER_FIELD },
- { "power_law_exponent", 1, OTHER_FIELD },
- { "maxwell_time", 1, OTHER_FIELD },
- { "total_strain", tensorSize, OTHER_FIELD },
- { "viscous_strain_t", tensorSize, OTHER_FIELD },
- { "stress_t", tensorSize, OTHER_FIELD },
+ const MetaData::ParamDescription properties[] = {
+ { "density", 1, pylith::topology::FieldBase::SCALAR },
+ { "mu", 1, pylith::topology::FieldBase::SCALAR },
+ { "lambda", 1, pylith::topology::FieldBase::SCALAR },
+ { "viscosity_coeff", 1, pylith::topology::FieldBase::SCALAR },
+ { "power_law_exponent", 1, pylith::topology::FieldBase::SCALAR }
};
- /// Indices (order) of properties.
- const int pidDensity = 0;
- const int pidMu = pidDensity + 1;
- const int pidLambda = pidMu + 1;
- const int pidViscosityCoeff = pidLambda + 1;
- const int pidPowerLawExp = pidViscosityCoeff + 1;
- const int pidMaxwellTime = pidPowerLawExp + 1;
- const int pidStrainT = pidMaxwellTime + 1;
- const int pidVisStrainT = pidStrainT + tensorSize;
- const int pidStressT = pidVisStrainT + tensorSize;
- /// Values expected in spatial database
- const int numDBValues = 5;
- const char* namesDBValues[] = {"density", "vs", "vp" , "viscosity_coeff",
- "power_law_exponent"};
+ // Values expected in properties spatial database
+ const int numDBProperties = 5;
+ const char* dbProperties[] = {"density", "vs", "vp" ,
+ "viscosity_coeff",
+ "power_law_exponent"};
- /// Indices (order) of database values
- const int didDensity = 0;
- const int didVs = 1;
- const int didVp = 2;
- const int didViscosityCoeff = 3;
- const int didPowerLawExp = 4;
+ /// Number of state variables.
+ const int numStateVars = 2;
- /// Initial state values expected in spatial database
- const int numInitialStateDBValues = tensorSize;
- const char* namesInitialStateDBValues[] = { "stress_xx", "stress_yy",
- "stress_zz", "stress_xy",
- "stress_yz", "stress_xz" };
+ /// State variables.
+ const Metadata::ParamDescription stateVars[] = {
+ { "viscous_strain", tensorSize, pylith::topology::FieldBase::TENSOR },
+ { "stress", tensorSize, pylith::topology::FieldBase::TENSOR }
+ };
+ // Values expected in state variables spatial database.
+ const int numDBStateVars = 12;
+ const char* dbStateVars[] = { "viscous-strain-xx",
+ "viscous-strain-yy",
+ "viscous-strain-zz",
+ "viscous-strain-xy",
+ "viscous-strain-yz",
+ "viscous-strain-xz",
+ "stress-xx",
+ "stress-yy",
+ "stress-zz",
+ "stress-xy",
+ "stress-yz",
+ "stress-xz"
+ };
+
} // _PowerLaw3D
} // materials
} // pylith
+// Indices of physical properties.
+const int pylith::materials::PowerLaw3D::p_density = 0;
+
+const int pylith::materials::PowerLaw3D::p_mu =
+ pylith::materials::PowerLaw3D::p_density + 1;
+
+const int pylith::materials::PowerLaw3D::p_lambda =
+ pylith::materials::PowerLaw3D::p_mu + 1;
+
+const int pylith::materials::PowerLaw3D::p_viscosityCoeff =
+ pylith::materials::PowerLaw3D::p_lambda + 1;
+
+const int pylith::materials::PowerLaw3D::p_powerLawExponent =
+ pylith::materials::PowerLaw3D::p_viscosityCoeff + 1;
+
+// Indices of property database values (order must match dbProperties).
+const int pylith::materials::PowerLaw3D::db_density = 0;
+
+const int pylith::materials::PowerLaw3D::db_vs =
+ pylith::materials::PowerLaw3D::db_density + 1;
+
+const int pylith::materials::PowerLaw3D::db_vp =
+ pylith::materials::PowerLaw3D::db_vs + 1;
+
+const int pylith::materials::PowerLaw3D::db_viscosityCoeff =
+ pylith::materials::PowerLaw3D::db_vp + 1;
+
+const int pylith::materials::PowerLaw3D::db_powerLawExponent =
+ pylith::materials::PowerLaw3D::db_viscosityCoeff + 1;
+
+// Indices of state variables.
+const int pylith::materials::PowerLaw3D::s_viscousStrain = 0;
+
+const int pylith::materials::PowerLaw3D::s_stress =
+ pylith::materials::PowerLaw3D::s_viscousStrain +
+ pylith::materials::PowerLaw3D::tensorSize;
+
+// Indices of state variable database values (order must match dbStateVars).
+const int pylith::materials::PowerLaw3D::db_viscousStrain = 0;
+
+const int pylith::materials::PowerLaw3D::db_stress =
+ pylith::materials::PowerLaw3D::db_viscousStrain +
+ pylith::materials::PowerLaw3D::tensorSize;
+
// ----------------------------------------------------------------------
// Default constructor.
pylith::materials::PowerLaw3D::PowerLaw3D(void) :
- ElasticMaterial(_PowerLaw3D::tensorSize,
+ ElasticMaterial(_PowerLaw3D::dimension,
+ _PowerLaw3D::tensorSize,
_PowerLaw3D::numElasticConsts,
- _PowerLaw3D::namesDBValues,
- _PowerLaw3D::namesInitialStateDBValues,
- _PowerLaw3D::numDBValues,
- _PowerLaw3D::properties,
- _PowerLaw3D::numProperties),
- _calcElasticConstsFn(&pylith::materials::PowerLaw3D::_calcElasticConstsElastic),
- _calcStressFn(&pylith::materials::PowerLaw3D::_calcStressElastic),
- _updatePropertiesFn(&pylith::materials::PowerLaw3D::_updatePropertiesElastic)
+ Metadata(_PowerLaw3D::properties,
+ _PowerLaw3D::numProperties,
+ _PowerLaw3D::dbProperties,
+ _PowerLaw3D::numDBProperties,
+ _PowerLaw3D::stateVars,
+ _PowerLaw3D::numStateVars,
+ _PowerLaw3D::dbStateVars,
+ _PowerLaw3D::numDBStateVars)),
+ _calcElasticConstsFn(0),
+ _calcStressFn(0),
+ _updateStateVarsFn(0)
{ // constructor
- _dimension = 3;
- _visStrainT.resize(_PowerLaw3D::tensorSize);
- _stressT.resize(_PowerLaw3D::tensorSize);
+ useElasticBehavior(true);
+ _viscousStrain.resize(_tensorSize);
+ _stress.resize(_tensorSize);
} // constructor
// ----------------------------------------------------------------------
@@ -112,6 +160,29 @@
} // destructor
// ----------------------------------------------------------------------
+// Set whether elastic or inelastic constitutive relations are used.
+void
+pylith::materials::PowerLaw3D::useElasticBehavior(const bool flag)
+{ // useElasticBehavior
+ if (flag) {
+ _calcStressFn =
+ &pylith::materials::PowerLaw3D::_calcStressElastic;
+ _calcElasticConstsFn =
+ &pylith::materials::PowerLaw3D::_calcElasticConstsElastic;
+ _updateStateVarsFn =
+ &pylith::materials::PowerLaw3D::_updateStateVarsElastic;
+
+ } else {
+ _calcStressFn =
+ &pylith::materials::PowerLaw3D::_calcStressViscoelastic;
+ _calcElasticConstsFn =
+ &pylith::materials::PowerLaw3D::_calcElasticConstsViscoelastic;
+ _updateStateVarsFn =
+ &pylith::materials::PowerLaw3D::_updateStateVarsViscoelastic;
+ } // if/else
+} // useElasticBehavior
+
+// ----------------------------------------------------------------------
// Compute properties from values in spatial database.
void
pylith::materials::PowerLaw3D::_dbToProperties(
@@ -122,14 +193,14 @@
const int numDBValues = dbValues.size();
assert(_PowerLaw3D::numDBValues == numDBValues);
- const double density = dbValues[_PowerLaw3D::didDensity];
- const double vs = dbValues[_PowerLaw3D::didVs];
- const double vp = dbValues[_PowerLaw3D::didVp];
- const double viscosityCoeff = dbValues[_PowerLaw3D::didViscosityCoeff];
- const double powerLawExp = dbValues[_PowerLaw3D::didPowerLawExp];
+ const double density = dbValues[db_density];
+ const double vs = dbValues[db_vs];
+ const double vp = dbValues[db_vp];
+ const double viscosityCoeff = dbValues[db_viscosityCoeff];
+ const double powerLawExponent = dbValues[db_powerLawExponent];
if (density <= 0.0 || vs <= 0.0 || vp <= 0.0 || viscosityCoeff <= 0.0
- || powerLawExp < 1.0) {
+ || powerLawExponent < 1.0) {
std::ostringstream msg;
msg << "Spatial database returned illegal value for physical "
<< "properties.\n"
@@ -137,7 +208,7 @@
<< "vp: " << vp << "\n"
<< "vs: " << vs << "\n"
<< "viscosityCoeff: " << viscosityCoeff << "\n"
- << "powerLawExp: " << powerLawExp << "\n";
+ << "powerLawExponent: " << powerLawExponent << "\n";
throw std::runtime_error(msg.str());
} // if
@@ -154,11 +225,11 @@
} // if
assert(mu > 0);
- propValues[_PowerLaw3D::pidDensity] = density;
- propValues[_PowerLaw3D::pidMu] = mu;
- propValues[_PowerLaw3D::pidLambda] = lambda;
- propValues[_PowerLaw3D::pidViscosityCoeff] = viscosityCoeff;
- propValues[_PowerLaw3D::pidPowerLawExp] = powerLawExp;
+ propValues[p_density] = density;
+ propValues[p_mu] = mu;
+ propValues[p_lambda] = lambda;
+ propValues[p_viscosityCoeff] = viscosityCoeff;
+ propValues[p_powerLawExponent] = powerLawExponent;
PetscLogFlops(6);
} // _dbToProperties
@@ -171,38 +242,26 @@
{ // _nondimProperties
assert(0 != _normalizer);
assert(0 != values);
- assert(nvalues == _totalPropsQuadPt);
+ assert(nvalues == _numPropsQuadPt);
const double densityScale = _normalizer->densityScale();
const double pressureScale = _normalizer->pressureScale();
const double timeScale = _normalizer->timeScale();
// **** NOTE: Make sure scaling is correct for viscosity coefficient.
- const double powerLawExp = _PowerLaw3D::pidPowerLawExp;
- const double viscosityCoeffScale = (pressureScale^(1.0/powerLawExp))/timeScale;
- const double powerLawExpScale = 1.0;
- values[_PowerLaw3D::pidDensity] =
- _normalizer->nondimensionalize(values[_PowerLaw3D::pidDensity],
- densityScale);
- values[_PowerLaw3D::pidMu] =
- _normalizer->nondimensionalize(values[_PowerLaw3D::pidMu],
- pressureScale);
- values[_PowerLaw3D::pidLambda] =
- _normalizer->nondimensionalize(values[_PowerLaw3D::pidLambda],
- pressureScale);
- values[_PowerLaw3D::pidViscosityCoeff] =
- _normalizer->nondimensionalize(values[_PowerLaw3D::pidViscosityCoeff],
+ const double powerLawExponent = values[p_powerLawExponent];
+ const double viscosityCoeffScale =
+ (pressureScale^(1.0/powerLawExponent))/timeScale;
+ values[p_density] =
+ _normalizer->nondimensionalize(values[p_density], densityScale);
+ values[p_mu] =
+ _normalizer->nondimensionalize(values[p_mu], pressureScale);
+ values[p_lambda] =
+ _normalizer->nondimensionalize(values[p_lambda], pressureScale);
+ values[p_viscosityCoeff] =
+ _normalizer->nondimensionalize(values[p_viscosityCoeff],
viscosityCoeffScale);
- values[_PowerLaw3D::pidPowerLawExp] =
- _normalizer->nondimensionalize(values[_PowerLaw3D::pidPowerLawExp],
- powerLawExpScale);
- values[_PowerLaw3D::pidMaxwellTime] =
- _normalizer->nondimensionalize(values[_PowerLaw3D::pidMaxwellTime],
- timeScale);
- values[_PowerLaw3D::pidStressT] =
- _normalizer->nondimensionalize(values[_PowerLaw3D::pidStressT],
- pressureScale);
- PetscLogFlops(9 + _PowerLaw3D::tensorSize);
+ PetscLogFlops(7);
} // _nondimProperties
// ----------------------------------------------------------------------
@@ -219,81 +278,93 @@
const double pressureScale = _normalizer->pressureScale();
const double timeScale = _normalizer->timeScale();
// **** NOTE: Make sure scaling is correct for viscosity coefficient.
- const double powerLawExp = _PowerLaw3D::pidPowerLawExp;
- const double viscosityCoeffScale = (pressureScale^(1.0/powerLawExp))/timeScale;
- const double powerLawExpScale = 1.0;
- values[_PowerLaw3D::pidDensity] =
- _normalizer->dimensionalize(values[_PowerLaw3D::pidDensity],
- densityScale);
- values[_PowerLaw3D::pidMu] =
- _normalizer->dimensionalize(values[_PowerLaw3D::pidMu],
- pressureScale);
- values[_PowerLaw3D::pidLambda] =
- _normalizer->dimensionalize(values[_PowerLaw3D::pidLambda],
- pressureScale);
- values[_PowerLaw3D::pidViscosityCoeff] =
- _normalizer->dimensionalize(values[_PowerLaw3D::pidViscosityCoeff],
+ const double powerLawExponent = values[p_powerLawExponent];
+ const double viscosityCoeffScale =
+ (pressureScale^(1.0/powerLawExponent))/timeScale;
+ values[p_density] =
+ _normalizer->dimensionalize(values[p_density], densityScale);
+ values[p_mu] =
+ _normalizer->dimensionalize(values[p_mu], pressureScale);
+ values[p_lambda] =
+ _normalizer->dimensionalize(values[p_lambda], pressureScale);
+ values[p_viscosityCoeff] =
+ _normalizer->dimensionalize(values[p_viscosityCoeff],
viscosityCoeffScale);
- values[_PowerLaw3D::pidPowerLawExp] =
- _normalizer->dimensionalize(values[_PowerLaw3D::pidPowerLawExp],
- powerLawExpScale);
- values[_PowerLaw3D::pidMaxwellTime] =
- _normalizer->dimensionalize(values[_PowerLaw3D::pidMaxwellTime],
- timeScale);
- values[_PowerLaw3D::pidStressT] =
- _normalizer->dimensionalize(values[_PowerLaw3D::pidStressT],
- pressureScale);
- PetscLogFlops(9 + _PowerLaw3D::tensorSize);
+ PetscLogFlops(7);
} // _dimProperties
// ----------------------------------------------------------------------
-// Nondimensionalize initial state.
+// Compute initial state variables from values in spatial database.
void
-pylith::materials::PowerLaw3D::_nondimInitState(double* const values,
- const int nvalues) const
-{ // _nondimInitState
- assert(0 != _normalizer);
- assert(0 != values);
- assert(nvalues == _PowerLaw3D::numInitialStateDBValues);
+pylith::materials::PowerLaw3D::_dbToStateVars(
+ double* const stateValues,
+ const double_array& dbValues) const
+{ // _dbToStateVars
+ assert(0 != stateValues);
+ const int numDBValues = dbValues.size();
+ assert(_PowerLaw3D::numDBStateVars == numDBValues);
- const double pressureScale = _normalizer->pressureScale();
- _normalizer->nondimensionalize(values, nvalues, pressureScale);
+ const int totalSize = 2 * _tensorSize;
+ assert(totalSize == _numVarsQuadPt);
+ assert(totalSize == numDBValues);
+ memcpy(stateValues[s_viscousStrain], &dbValues[db_viscousStrain],
+ _tensorSize*sizeof(double));
+ memcpy(stateValues[s_stress], &dbValues[db_stress],
+ _tensorSize*sizeof(double));
- PetscLogFlops(nvalues);
-} // _nondimInitState
+ PetscLogFlops(0);
+} // _dbToStateVars
// ----------------------------------------------------------------------
-// Dimensionalize initial state.
-void
-pylith::materials::PowerLaw3D::_dimInitState(double* const values,
- const int nvalues) const
-{ // _dimInitState
- assert(0 != _normalizer);
- assert(0 != values);
- assert(nvalues == _PowerLaw3D::numInitialStateDBValues);
-
- const double pressureScale = _normalizer->pressureScale();
- _normalizer->dimensionalize(values, nvalues, pressureScale);
-
- PetscLogFlops(nvalues);
-} // _dimInitState
-
-// ----------------------------------------------------------------------
// Compute density at location from properties.
void
pylith::materials::PowerLaw3D::_calcDensity(double* const density,
- const double* properties,
- const int numProperties)
+ const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars)
{ // _calcDensity
assert(0 != density);
assert(0 != properties);
assert(_totalPropsQuadPt == numProperties);
- density[0] = properties[_PowerLaw3D::pidDensity];
+ density[0] = properties[p_density];
} // _calcDensity
// ----------------------------------------------------------------------
+// Get stable time step for implicit time integration.
+double
+pylith::materials::PowerLaw3D::_stableTimeStepImplicit(
+ const double* properties,
+ const int numProperties,
+ const double* stateVars,
+ const int numStateVars) const
+{ // _stableTimeStepImplicit
+ assert(0 != properties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(0 != stateVars);
+ assert(_numVarsQuadPt == numStateVars);
+
+ memcpy(&stress[0], &stateVars[s_stress],
+ tensorSize * sizeof(double));
+ const double meanStress = (stress[0] + stress[1] + stress[2])/3.0;
+ const double devStress[] = {stress[0] - meanStress,
+ stress[1] - meanStress,
+ stress[2] - meanStress,
+ stress[3],
+ stress[4],
+ stress[5] };
+ const double effStress = sqrt(0.5 * _scalarProduct(devStress, devStress));
+ dtStable = 1.0;
+ if (effStress != 0.0)
+ dtStable = 0.1 * ((viscosityCoeff/effStressTpdt) ^ (powerLawExp - 1.0)) *
+ (viscosityCoeff/mu);
+
+ return dtStable;
+} // _stableTimeStepImplicit
+
+// ----------------------------------------------------------------------
// Compute stress tensor at location from properties as an elastic
// material.
void
@@ -302,155 +373,54 @@
const int stressSize,
const double* properties,
const int numProperties,
+ const double* stateVars,
+ const int numStateVars,
const double* totalStrain,
const int strainSize,
- const double* initialState,
- const int initialStateSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize,
const bool computeStateVars)
{ // _calcStressElastic
assert(0 != stress);
assert(_PowerLaw3D::tensorSize == stressSize);
assert(0 != properties);
- assert(_totalPropsQuadPt == numProperties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(0 != stateVars);
+ assert(_numVarsQuadPt == numStateVars);
assert(0 != totalStrain);
assert(_PowerLaw3D::tensorSize == strainSize);
- assert(_PowerLaw3D::tensorSize == initialStateSize);
+ assert(0 != initialStress);
+ assert(_PowerLaw3D::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_PowerLaw3D::tensorSize == initialStrainSize);
- const double mu = properties[_PowerLaw3D::pidMu];
- const double lambda = properties[_PowerLaw3D::pidLambda];
- const double lambda = properties[_PowerLaw3D::pidLambda];
- const double viscosityCoeff = properties[_PowerLaw3D::pidViscosityCoeff];
- const double powerLawExp = properties[_PowerLaw3D::pidPowerLawExp];
+ const double mu = properties[p_mu];
+ const double lambda = properties[p_lambda];
const double mu2 = 2.0 * mu;
- const double e11 = totalStrain[0];
- const double e22 = totalStrain[1];
- const double e33 = totalStrain[2];
- const double e12 = totalStrain[3];
- const double e23 = totalStrain[4];
- const double e13 = totalStrain[5];
+ const double e11 = totalStrain[0] - initialStrain[0];
+ const double e22 = totalStrain[1] - initialStrain[1];
+ const double e33 = totalStrain[2] - initialStrain[2];
+ const double e12 = totalStrain[3] - initialStrain[3];
+ const double e23 = totalStrain[4] - initialStrain[4];
+ const double e13 = totalStrain[5] - initialStrain[5];
const double traceStrainTpdt = e11 + e22 + e33;
const double s123 = lambda * traceStrainTpdt;
- stress[0] = s123 + mu2*e11 + initialState[0];
- stress[1] = s123 + mu2*e22 + initialState[1];
- stress[2] = s123 + mu2*e33 + initialState[2];
- stress[3] = mu2 * e12 + initialState[3];
- stress[4] = mu2 * e23 + initialState[4];
- stress[5] = mu2 * e13 + initialState[5];
+ stress[0] = s123 + mu2*e11 + initialStress[0];
+ stress[1] = s123 + mu2*e22 + initialStress[1];
+ stress[2] = s123 + mu2*e33 + initialStress[2];
+ stress[3] = mu2 * e12 + initialStress[3];
+ stress[4] = mu2 * e23 + initialStress[4];
+ stress[5] = mu2 * e13 + initialStress[5];
- // The following section is put in here for now just to compute the Maxwell
- // time based on the elastic solution.
- const double meanStressTpdt = (stress[0] + stress[1] + stress[2])/3.0;
- double devStressTpdt[] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
- const double diag[] = { 1.0, 1.0, 1.0, 0.0, 0.0, 0.0};
- for (int iComp=0; iComp < stressSize; ++iComp) {
- devStressTpdt[iComp] = stress[iComp] - diag[iComp] * meanStressTpdt;
- } // for
- const double effStressTpdt = sqrt(0.5 *
- _scalarProduct(devStressTpdt,
- devStressTpdt));
- properties[_PowerLaw3D::pidMaxwellTime] = 1.0e30;
- if (effStressTpdt != 0.0)
- properties[_PowerLaw3D::pidMaxwellTime] = ((viscosityCoeff/effStressTpdt) ^
- (powerLawExp - 1.0)) *
- (viscosityCoeff/mu);
-
- PetscLogFlops(29 + stressSize * 2);
+ PetscLogFlops(25);
} // _calcStressElastic
// ----------------------------------------------------------------------
-// Effective stress function that computes effective stress function only
-// (no derivative).
-double
-pylith::materials::PowerLaw3D::_effStressFunc(double effStressTpdt,
- double *params)
-{ // _effStressFunc
- double ae = params[0];
- double b = params[1];
- double c = params[2];
- double d = params[3];
- double alpha = params[4];
- double dt = params[5];
- double effStressT = params[6];
- double powerLawExp = params[7];
- double viscosityCoeff = params[8];
- double factor1 = 1.0-alpha;
- double effStressTau = factor1 * effStressT + alpha * effStressTpdt;
- double gammaTau = 0.5 * ((effStressTau/viscosityCoeff)^
- (powerLawExp - 1.0))/viscosityCoeff;
- double a = ae + alpha * dt * gammaTau;
- double y = a * a * effStressTpdt * effStressTpdt - b +
- c * gammaTau - d * d * gammaTau * gammaTau;
- PetscLogFlops(21);
- return y;
-} // _effStressFunc
-
-// ----------------------------------------------------------------------
-// Effective stress function that computes effective stress function
-// derivative only (no function value).
-double
-pylith::materials::PowerLaw3D::_effStressDFunc(double effStressTpdt,
- double *params)
-{ // _effStressDFunc
- double ae = params[0];
- double c = params[2];
- double d = params[3];
- double alpha = params[4];
- double dt = params[5];
- double effStressT = params[6];
- double powerLawExp = params[7];
- double viscosityCoeff = params[8];
- double factor1 = 1.0-alpha;
- double effStressTau = factor1 * effStressT + alpha * effStressTpdt;
- double gammaTau = 0.5 * ((effStressTau/viscosityCoeff)^
- (powerLawExp - 1.0))/viscosityCoeff;
- double a = ae + alpha * dt * gammaTau;
- double dGammaTau = 0.5 * alpha * (powerLawExp - 1.0) *
- ((effStressTau/viscosityCoeff) ^ (powerLawExp - 2.0))/
- (viscosityCoeff * viscosityCoeff);
- double dy = 2.0 * a * a * effStressTpdt + dGammaTau *
- (2.0 * a * alpha * dt * effStressTpdt * effStressTpdt +
- c - 2.0 * d * d * gammaTau);
- PetscLogFlops(36);
- return dy;
-} // _effStressDFunc
-
-// ----------------------------------------------------------------------
-// Effective stress function that computes effective stress function
-// and derivative.
-void
-pylith::materials::PowerLaw3D::_effStressFuncDFunc(double effStressTpdt,
- double *params,
- double *y,
- double *dy)
-{ // _effStressFunc
- double ae = params[0];
- double b = params[1];
- double c = params[2];
- double d = params[3];
- double alpha = params[4];
- double dt = params[5];
- double effStressT = params[6];
- double powerLawExp = params[7];
- double viscosityCoeff = params[8];
- double factor1 = 1.0-alpha;
- double effStressTau = factor1 * effStressT + alpha * effStressTpdt;
- double gammaTau = 0.5 * ((effStressTau/viscosityCoeff)^
- (powerLawExp - 1.0))/viscosityCoeff;
- double dGammaTau = 0.5 * alpha * (powerLawExp - 1.0) *
- ((effStressTau/viscosityCoeff) ^ (powerLawExp - 2.0))/
- (viscosityCoeff * viscosityCoeff);
- double a = ae + alpha * dt * gammaTau;
- double *y = a * a * effStressTpdt * effStressTpdt - b +
- c * gammaTau - d * d * gammaTau * gammaTau;
- double *dy = 2.0 * a * a * effStressTpdt + dGammaTau *
- (2.0 * a * alpha * dt * effStressTpdt * effStressTpdt +
- c - 2.0 * d * d * gammaTau);
- PetscLogFlops(46);
-} // _effStressFuncDFunc
-// ----------------------------------------------------------------------
// Compute stress tensor at location from properties as a viscoelastic
// material.
void
@@ -459,19 +429,28 @@
const int stressSize,
const double* properties,
const int numProperties,
+ const double* stateVars,
+ const int numStateVars,
const double* totalStrain,
const int strainSize,
- const double* initialState,
- const int initialStateSize,
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize,
const bool computeStateVars)
{ // _calcStressViscoelastic
assert(0 != stress);
assert(_PowerLaw3D::tensorSize == stressSize);
assert(0 != properties);
- assert(_totalPropsQuadPt == numProperties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(0 != stateVars);
+ assert(_numPropsQuadPt == numStateVars);
assert(0 != totalStrain);
assert(_PowerLaw3D::tensorSize == strainSize);
- assert(_PowerLaw3D::tensorSize == initialStateSize);
+ assert(0 != initialStress);
+ assert(_PowerLaw3D::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_PowerLaw3D::tensorSize == initialStrainSize);
const int tensorSize = _PowerLaw3D::tensorSize;
@@ -479,13 +458,13 @@
// time step.
if (computeStateVars) {
- const double mu = properties[_PowerLaw3D::pidMu];
- const double lambda = properties[_PowerLaw3D::pidLambda];
- const double viscosityCoeff = properties[_PowerLaw3D::pidViscosityCoeff];
- const double powerLawExp = properties[_PowerLaw3D::pidPowerLawExp];
- memcpy(&stressT[0], &properties[_PowerLaw3D::pidStressT],
+ const double mu = properties[p_mu];
+ const double lambda = properties[p_lambda];
+ const double viscosityCoeff = properties[p_viscosityCoeff];
+ const double powerLawExp = properties[p_powerLawExponent];
+ memcpy(&visStrainT[0], &stateVars[s_viscousStrain],
tensorSize * sizeof(double));
- memcpy(&visStrainT[0], &properties[_PowerLaw3D::pidVisStrainT],
+ memcpy(&stressT[0], &stateVars[s_stress],
tensorSize * sizeof(double));
const double mu2 = 2.0 * mu;
@@ -501,45 +480,50 @@
const double alpha = 0.5;
const double timeFac = _dt * (1.0 - alpha);
+ // Initial stress values
+ const double meanStressInitial = (initialStress[0] + initialStress[1] +
+ initialStress[2])/3.0;
+ const double devStressInitial[] = { initialStress[0] - meanStressInitial,
+ initialStress[1] - meanStressInitial,
+ initialStress[2] - meanStressInitial,
+ initialStress[3],
+ initialStress[4],
+ initialStress[5] };
+ const double stressInvar2Initial = 0.5 *
+ _scalarProduct(devStressInitial, devStressInitial);
+
+ // Initial strain values
+ const double meanStrainInitial = (initialStrain[0] + initialStrain[1] +
+ initialStrain[2])/3.0;
+
// Values for current time step
const double e11 = totalStrain[0];
const double e22 = totalStrain[1];
const double e33 = totalStrain[2];
const double traceStrainTpdt = e11 + e22 + e33;
- const double meanStrainTpdt = traceStrainTpdt/3.0;
+ const double meanStrainTpdt = traceStrainTpdt/3.0 - meanStrainInitial;
const double meanStressTpdt = bulkModulus * traceStrainTpdt;
- // Initial stress values
- const double meanStressInitial = (_stressInitial[0] + _stressInitial[1] +
- _stressInitial[2])/3.0;
- const double devStressInitial[] = { _stressInitial[0] - meanStressInitial,
- _stressInitial[1] - meanStressInitial,
- _stressInitial[2] - meanStressInitial,
- _stressInitial[3],
- _stressInitial[4],
- _stressInitial[5] };
- const double stressInvar2Initial = 0.5 *
- _scalarProduct(devStressInitial, devStressInitial);
-
- // Values for current time step
+ // Note that I use the initial strain rather than the deviatoric initial
+ // strain since otherwise the initial mean strain would get used twice.
const double strainPPTpdt[] =
- { _totalStrain[0] - meanStrainTpdt - visStrainT[0],
- _totalStrain[1] - meanStrainTpdt - visStrainT[1],
- _totalStrain[2] - meanStrainTpdt - visStrainT[2],
- _totalStrain[3] - visStrainT[3],
- _totalStrain[4] - visStrainT[4],
- _totalStrain[5] - visStrainT[5] };
+ { totalStrain[0] - meanStrainTpdt - visStrainT[0] - initialStrain[0],
+ totalStrain[1] - meanStrainTpdt - visStrainT[1] - initialStrain[1],
+ totalStrain[2] - meanStrainTpdt - visStrainT[2] - initialStrain[2],
+ totalStrain[3] - visStrainT[3] - initialStrain[3],
+ totalStrain[4] - visStrainT[4] - initialStrain[4],
+ totalStrain[5] - visStrainT[5] - initialStrain[5] };
const double strainPPInvar2Tpdt = 0.5 *
_scalarProduct(strainPPTpdt, strainPPTpdt);
// Values for previous time step
- const double meanStressT = (_stressT[0] + _stressT[1] + _stressT[2])/3.0;
- const double devStressT[] = { _stressT[0] - meanStressT,
- _stressT[1] - meanStressT,
- _stressT[2] - meanStressT,
- _stressT[3],
- _stressT[4],
- _stressT[5] };
+ const double meanStressT = (stressT[0] + stressT[1] + stressT[2])/3.0;
+ const double devStressT[] = { stressT[0] - meanStressT,
+ stressT[1] - meanStressT,
+ stressT[2] - meanStressT,
+ stressT[3],
+ stressT[4],
+ stressT[5] };
const double stressInvar2T = 0.5 * _scalarProduct(devStressT, devStressT);
const double effStressT = sqrt(stressInvar2T);
@@ -552,7 +536,7 @@
timeFac;
const double d = timeFac * effStressT;
- PetscLogFlops(45);
+ PetscLogFlops(55);
// Put parameters into a vector and call root-finding algorithm.
// This could also be a struct.
const double effStressParams[] = {ae,
@@ -573,15 +557,6 @@
pylith::materials::PowerLaw3D::_effStressFunc,
pylith::materials::PowerLaw3D::_effStressFuncDFunc);
- // Compute Maxwell time
- properties[_PowerLaw3D::pidMaxwellTime] = 1.0e30;
- if (effStressTpdt != 0.0) {
- properties[_PowerLaw3D::pidMaxwellTime] =
- ((viscosityCoeff/effStressTpdt) ^ (powerLawExp - 1.0)) *
- (viscosityCoeff/mu);
- PetscLogFlops(5);
- } // if
-
// Compute stresses from effective stress.
const double effStressTau = (1.0 - alpha) * effStressT +
alpha * effStressTpdt;
@@ -601,15 +576,105 @@
PetscLogFlops(14 + 8 * tensorSize);
// If state variables have already been updated, current stress is already
- // contained in stressT.
+ // contained in stress.
} else {
- memcpy(&stress[0], &properties[_PowerLaw3D::pidStressT],
- tensorSize * sizeof(double));
+ memcpy(&stress[0], &stateVars[s_stress], tensorSize * sizeof(double));
} // else
} // _calcStressViscoelastic
// ----------------------------------------------------------------------
+// Effective stress function that computes effective stress function only
+// (no derivative).
+double
+pylith::materials::PowerLaw3D::_effStressFunc(double effStressTpdt,
+ double *params)
+{ // _effStressFunc
+ double ae = params[0];
+ double b = params[1];
+ double c = params[2];
+ double d = params[3];
+ double alpha = params[4];
+ double dt = params[5];
+ double effStressT = params[6];
+ double powerLawExp = params[7];
+ double viscosityCoeff = params[8];
+ double factor1 = 1.0-alpha;
+ double effStressTau = factor1 * effStressT + alpha * effStressTpdt;
+ double gammaTau = 0.5 * ((effStressTau/viscosityCoeff)^
+ (powerLawExp - 1.0))/viscosityCoeff;
+ double a = ae + alpha * dt * gammaTau;
+ double y = a * a * effStressTpdt * effStressTpdt - b +
+ c * gammaTau - d * d * gammaTau * gammaTau;
+ PetscLogFlops(21);
+ return y;
+} // _effStressFunc
+
+// ----------------------------------------------------------------------
+// Effective stress function that computes effective stress function
+// derivative only (no function value).
+double
+pylith::materials::PowerLaw3D::_effStressDFunc(double effStressTpdt,
+ double *params)
+{ // _effStressDFunc
+ double ae = params[0];
+ double c = params[2];
+ double d = params[3];
+ double alpha = params[4];
+ double dt = params[5];
+ double effStressT = params[6];
+ double powerLawExp = params[7];
+ double viscosityCoeff = params[8];
+ double factor1 = 1.0-alpha;
+ double effStressTau = factor1 * effStressT + alpha * effStressTpdt;
+ double gammaTau = 0.5 * ((effStressTau/viscosityCoeff)^
+ (powerLawExp - 1.0))/viscosityCoeff;
+ double a = ae + alpha * dt * gammaTau;
+ double dGammaTau = 0.5 * alpha * (powerLawExp - 1.0) *
+ ((effStressTau/viscosityCoeff) ^ (powerLawExp - 2.0))/
+ (viscosityCoeff * viscosityCoeff);
+ double dy = 2.0 * a * a * effStressTpdt + dGammaTau *
+ (2.0 * a * alpha * dt * effStressTpdt * effStressTpdt +
+ c - 2.0 * d * d * gammaTau);
+ PetscLogFlops(36);
+ return dy;
+} // _effStressDFunc
+
+// ----------------------------------------------------------------------
+// Effective stress function that computes effective stress function
+// and derivative.
+void
+pylith::materials::PowerLaw3D::_effStressFuncDFunc(double effStressTpdt,
+ double *params,
+ double *y,
+ double *dy)
+{ // _effStressFunc
+ double ae = params[0];
+ double b = params[1];
+ double c = params[2];
+ double d = params[3];
+ double alpha = params[4];
+ double dt = params[5];
+ double effStressT = params[6];
+ double powerLawExp = params[7];
+ double viscosityCoeff = params[8];
+ double factor1 = 1.0-alpha;
+ double effStressTau = factor1 * effStressT + alpha * effStressTpdt;
+ double gammaTau = 0.5 * ((effStressTau/viscosityCoeff)^
+ (powerLawExp - 1.0))/viscosityCoeff;
+ double dGammaTau = 0.5 * alpha * (powerLawExp - 1.0) *
+ ((effStressTau/viscosityCoeff) ^ (powerLawExp - 2.0))/
+ (viscosityCoeff * viscosityCoeff);
+ double a = ae + alpha * dt * gammaTau;
+ double *y = a * a * effStressTpdt * effStressTpdt - b +
+ c * gammaTau - d * d * gammaTau * gammaTau;
+ double *dy = 2.0 * a * a * effStressTpdt + dGammaTau *
+ (2.0 * a * alpha * dt * effStressTpdt * effStressTpdt +
+ c - 2.0 * d * d * gammaTau);
+ PetscLogFlops(46);
+} // _effStressFuncDFunc
+
+// ----------------------------------------------------------------------
// Compute derivative of elasticity matrix at location from properties.
void
pylith::materials::PowerLaw3D::_calcElasticConstsElastic(
@@ -617,21 +682,30 @@
const int numElasticConsts,
const double* properties,
const int numProperties,
+ const double* stateVars,
+ const int numStateVars,
const double* totalStrain,
const int strainSize,
- const double* initialState,
- const int initialStateSize)
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize)
{ // _calcElasticConstsElastic
assert(0 != elasticConsts);
assert(_PowerLaw3D::numElasticConsts == numElasticConsts);
assert(0 != properties);
- assert(_totalPropsQuadPt == numProperties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(0 != stateVars);
+ assert(_numVarsQuadPt == numStateVars);
assert(0 != totalStrain);
assert(_PowerLaw3D::tensorSize == strainSize);
- assert(_PowerLaw3D::tensorSize == initialStateSize);
+ assert(0 != initialStress);
+ assert(_PowerLaw3D::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_PowerLaw3D::tensorSize == initialStrainSize);
- const double mu = properties[_PowerLaw3D::pidMu];
- const double lambda = properties[_PowerLaw3D::pidLambda];
+ const double mu = properties[p_mu];
+ const double lambda = properties[p_lambda];
const double mu2 = 2.0 * mu;
const double lambda2mu = lambda + mu2;
@@ -658,7 +732,7 @@
elasticConsts[19] = 0; // C2313
elasticConsts[20] = mu2; // C1313
- PetscLogFlops(4);
+ PetscLogFlops(2);
} // _calcElasticConstsElastic
// ----------------------------------------------------------------------
@@ -671,26 +745,34 @@
const int numElasticConsts,
const double* properties,
const int numProperties,
+ const double* stateVars,
+ const int numStateVars,
const double* totalStrain,
const int strainSize,
- const double* initialState,
- const int initialStateSize)
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize)
{ // _calcElasticConstsViscoelasticInitial
assert(0 != elasticConsts);
assert(_PowerLaw3D::numElasticConsts == numElasticConsts);
assert(0 != properties);
- assert(_totalPropsQuadPt == numProperties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(0 != stateVars);
+ assert(_numPropsQuadPt == numStateVars);
assert(0 != totalStrain);
assert(_PowerLaw3D::tensorSize == strainSize);
- assert(_PowerLaw3D::tensorSize == initialStateSize);
+ assert(0 != initialStress);
+ assert(_PowerLaw3D::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_PowerLaw3D::tensorSize == initialStrainSize);
const int tensorSize = _PowerLaw3D::tensorSize;
- const double mu = properties[_PowerLaw3D::pidMu];
- const double lambda = properties[_PowerLaw3D::pidLambda];
- const double viscosityCoeff = properties[_PowerLaw3D::pidViscosityCoeff];
- const double powerLawExp = properties[_PowerLaw3D::pidPowerLawExp];
- memcpy(&stress[0], &properties[_PowerLaw3D::pidStressT],
- tensorSize * sizeof(double));
+ const double mu = properties[p_mu];
+ const double lambda = properties[p_lambda];
+ const double viscosityCoeff = properties[p_viscosityCoeff];
+ const double powerLawExp = properties[p_powerLawExponent];
+ memcpy(&stress[0], &stateVars[s_stress], tensorSize * sizeof(double));
const double mu2 = 2.0 * mu;
const double ae = 1.0/mu2;
@@ -743,20 +825,35 @@
const int numElasticConsts,
const double* properties,
const int numProperties,
+ const double* stateVars,
+ const int numStateVars,
const double* totalStrain,
const int strainSize,
- const double* initialState,
- const int initialStateSize)
+ const double* initialStress,
+ const int initialStressSize,
+ const double* initialStrain,
+ const int initialStrainSize)
{ // _calcElasticConstsViscoelastic
+ assert(0 != elasticConsts);
+ assert(_PowerLaw3D::numElasticConsts == numElasticConsts);
+ assert(0 != properties);
+ assert(_numPropsQuadPt == numProperties);
+ assert(0 != stateVars);
+ assert(_numPropsQuadPt == numStateVars);
+ assert(0 != totalStrain);
+ assert(_PowerLaw3D::tensorSize == strainSize);
+ assert(0 != initialStress);
+ assert(_PowerLaw3D::tensorSize == initialStressSize);
+ assert(0 != initialStrain);
+ assert(_PowerLaw3D::tensorSize == initialStrainSize);
- const double mu = properties[_PowerLaw3D::pidMu];
- const double lambda = properties[_PowerLaw3D::pidLambda];
- const double viscosityCoeff = properties[_PowerLaw3D::pidViscosityCoeff];
- const double powerLawExp = properties[_PowerLaw3D::pidPowerLawExp];
- memcpy(&stressT[0], &properties[_PowerLaw3D::pidStressT],
+ const double mu = properties[p_mu];
+ const double lambda = properties[p_lambda];
+ const double viscosityCoeff = properties[p_viscosityCoeff];
+ const double powerLawExp = properties[p_powerLawExponent];
+ memcpy(&visStrainT[0], &stateVars[s_viscousStrain],
tensorSize * sizeof(double));
- memcpy(&visStrainT[0], &properties[_PowerLaw3D::pidVisStrainT],
- tensorSize * sizeof(double));
+ memcpy(&stressT[0], &stateVars[s_stress], tensorSize * sizeof(double));
const double mu2 = 2.0 * mu;
const double lamPlusMu = lambda + mu;
@@ -913,22 +1010,6 @@
} // _calcElasticConstsViscoelastic
// ----------------------------------------------------------------------
-// Get stable time step for implicit time integration.
-double
-pylith::materials::PowerLaw3D::_stableTimeStepImplicit(const double* properties,
- const int numProperties) const
-{ // _stableTimeStepImplicit
- assert(0 != properties);
- assert(_totalPropsQuadPt == numProperties);
-
- const double maxwellTime =
- properties[_PowerLaw3D::pidMaxwellTime];
- const double dtStable = 0.1*maxwellTime;
-
- return dtStable;
-} // _stableTimeStepImplicit
-
-// ----------------------------------------------------------------------
// Update state variables.
void
pylith::materials::PowerLaw3D::_updatePropertiesElastic(
More information about the CIG-COMMITS
mailing list