Skip to content

Commit

Permalink
Update viscoelastic documentation and formatting requirements
Browse files Browse the repository at this point in the history
  • Loading branch information
naliboff committed Mar 20, 2018
1 parent de29cd7 commit 0881424
Show file tree
Hide file tree
Showing 8 changed files with 101 additions and 34 deletions.
18 changes: 10 additions & 8 deletions include/aspect/material_model/viscoelastic.h
Expand Up @@ -72,8 +72,8 @@ namespace aspect
* except density, which varies linearly with temperature according to the
* thermal expansivity.
*
* When more than one field is present at a point, they are averaged
* arithmetically. An exception is viscosity, which may be averaged
* When more than one compositional field is present at a point, they are
* averaged arithmetically. An exception is viscosity, which may be averaged
* arithmetically, harmonically, geometrically, or by selecting the
* viscosity of the composition with the greatest volume fraction.
*
Expand All @@ -85,7 +85,6 @@ namespace aspect
class Viscoelastic : public MaterialModel::Interface<dim>, public ::aspect::SimulatorAccess<dim>
{
public:

/**
* Function to compute the material properties in @p out given the
* inputs in @p in. If MaterialModelInputs.strain_rate has the length
Expand Down Expand Up @@ -139,7 +138,6 @@ namespace aspect
*/

private:

/**
* The first 3 (2D) or 6 (3D) compositional fields are assumed
* to be components of the viscoelastic stress tensor and
Expand Down Expand Up @@ -208,14 +206,18 @@ namespace aspect

/**
* Bool indicating whether to use a fixed material time scale in the
* viscoelastic rheology during all time step, read from parameter file.
* viscoelastic rheology for all time steps (if true) or to use the
* actual (variable) advection time step of the model (if false). Read
* from parameter file.
*/
bool use_fixed_elastic_time_step;

/**
* Boolindicating whether to use a stress averaging scheme to account
* for differences between the numerical and fixed elastic time step,
* read from parameter file.
* Bool indicating whether to use a stress averaging scheme to account
* for differences between the numerical and fixed elastic time step
* (if true). When set to false, the viscoelastic stresses are not
* modified to account for differences between the viscoelastic time
* step and the numerical time step. Read from parameter file.
*/
bool use_stress_averaging;

Expand Down
75 changes: 70 additions & 5 deletions source/material_model/viscoelastic.cc
Expand Up @@ -125,6 +125,43 @@ namespace aspect
MaterialModel::MaterialModelOutputs<dim> &out) const
{


// Check whether the compositional fields representing the viscoelastic
// stress tensor are both named correctly and listed in the right order.
if ( dim == 2)
{
AssertThrow(this->introspection().compositional_index_for_name("sxx") == 0,
ExcMessage("Material model Viscoelastic only works if the first "
"compositional field is called sxx."));
AssertThrow(this->introspection().compositional_index_for_name("syy") == 1,
ExcMessage("Material model Viscoelastic only works if the second "
"compositional field is called syy."));
AssertThrow(this->introspection().compositional_index_for_name("sxy") == 2,
ExcMessage("Material model Viscoelastic only works if the third "
"compositional field is called sxy."));
}
if ( dim == 3)
{
AssertThrow(this->introspection().compositional_index_for_name("sxx") == 0,
ExcMessage("Material model Viscoelastic only works if the first "
"compositional field is called sxx."));
AssertThrow(this->introspection().compositional_index_for_name("syy") == 1,
ExcMessage("Material model Viscoelastic only works if the second "
"compositional field is called syy."));
AssertThrow(this->introspection().compositional_index_for_name("szz") == 2,
ExcMessage("Material model Viscoelastic only works if the third "
"compositional field is called szz."));
AssertThrow(this->introspection().compositional_index_for_name("sxy") == 3,
ExcMessage("Material model Viscoelastic only works if the fourth "
"compositional field is called sxy."));
AssertThrow(this->introspection().compositional_index_for_name("sxz") == 4,
ExcMessage("Material model Viscoelastic only works if the fifth "
"compositional field is called sxz."));
AssertThrow(this->introspection().compositional_index_for_name("syz") == 5,
ExcMessage("Material model Viscoelastic only works if the sixth "
"compositional field is called syz."));
}

// Define elastic time step
const double dte = ( ( this->get_timestep_number() > 0 && use_fixed_elastic_time_step == false )
?
Expand Down Expand Up @@ -436,11 +473,11 @@ namespace aspect
"rock type or component of the viscoelastic stress tensor. The stress "
"tensor in 2D and 3D, respectively, contains 3 or 6 components. The "
"compositional fields representing these components must be the first "
"listed compositional fields in the parameter file."
"listed compositional fields in the parameter file. "
"\n\n "
"Expanding the model to include non-linear viscous flow (e.g., "
"diffusion/dislocation creep) and plasticity would produce a "
"of constitutive relationship commonly referred to as partial "
"constitutive relationship commonly referred to as partial "
"elastoviscoplastic (e.g., pEVP) in the geodynamics community. "
"While extensively discussed and applied within the geodynamics "
"literature, notable references include: "
Expand All @@ -462,7 +499,7 @@ namespace aspect
"rate of deformation ($\\hat{D}$) as the sum of elastic "
"(($\\hat{D_{e}}$) and viscous (($\\hat{D_{v}}$)) components: "
"$\\hat{D} = \\hat{D_{e}} + \\hat{D_{v}}$ "
"These terms further decompose into"
"These terms further decompose into "
"$\\hat{D_{v}} = \\frac{\\tau}{2\\eta}$ and "
"$\\hat{D_{e}} = \\frac{\\overset{\\triangledown}{\\tau}}{2\\mu}$, where "
"$\\tau$ is the viscous deviatoric stress, $\\eta$ is the shear viscosity, "
Expand Down Expand Up @@ -496,12 +533,40 @@ namespace aspect
"($ \\alpha = \\frac{\\eta}{\\mu} $): "
"$\\eta_{eff} = \\eta \\frac{\\triangle t^{e}}{\\triangle t^{e} + \\alpha}$ "
"The magnitude of the shear modulus thus controls how much the effective "
"viscosity is reduced relative to the initial viscosity."
"viscosity is reduced relative to the initial viscosity. "
"\n\n"
"Elastic effects are introduced into the governing stokes equations through "
"an elastic force term (eqn. 30) using stresses from the previous time step: "
"$F^{e,t} = -\\frac{\\eta_{eff}}{\\mu \\triangle t^{e}} \\tau^{t}$. "
"This force term is added onto the right-hand side force vector in the "
"system of equations. ")
"system of equations. "
"\n\n"
"Several model parameters (densities, elastic shear moduli, thermal expansivities, "
"thermal conductivies, specific heats) can be defined per-compositional field. "
"For each material parameter the user supplies a comma delimited list of length "
"N+1, where N is the number of compositional fields. The additional field corresponds "
"to the value for background material. They should be ordered ''background, "
"composition1, composition2...''. However, the first 3 (2D) or 6 (3D) composition "
"fields correspond to components of the elastic stress tensor and their material "
"values will not contribute to the volume fractions. If a single value is given, then "
"all the compositional fields are given that value. Other lengths of lists are not "
"allowed. For a given compositional field the material parameters are treated as "
"constant, except density, which varies linearly with temperature according to the "
"thermal expansivity. "
"\n\n"
"When more than one compositional field is present at a point, they are averaged "
"arithmetically. An exception is viscosity, which may be averaged arithmetically, "
"harmonically, geometrically, or by selecting the viscosity of the composition field "
"with the greatest volume fraction. "
"\n\n"
"As mentioned above, the viscoelastic stress tensor is tracked through 3 (2D) or "
"6 (3D) individual components on compositional fields or tracers. When using tracers, "
"corresponding compositional fields are still required for the material to access the "
"tracer values. In either case, the stress tensor components must be named and listed "
"in a very specific format, which is designed to minimize mislabeling stress tensor "
"components as distinct 'compositional rock types' (or vice versa). For 2D models, the "
"first three compositional fields must be labeled sxx, syy and sxy. In 3D, the first "
"six compositional fields must be labeled sxx, syy, szz, sxy, sxz, syz. In both cases, "
"x, y and z correspond to the coordinate axes nomenclature used by the Geometry model.")
}
}
2 changes: 1 addition & 1 deletion tests/viscoelastic_fixed_elastic_time_step.prm
Expand Up @@ -62,7 +62,7 @@ end
# Number and name of compositional fields
subsection Compositional fields
set Number of fields = 3
set Names of fields = s11, s22, s12
set Names of fields = sxx, syy, sxy
end

# Spatial domain of different compositional fields
Expand Down
12 changes: 6 additions & 6 deletions tests/viscoelastic_fixed_elastic_time_step/screen-output
Expand Up @@ -4,9 +4,9 @@ Number of degrees of freedom: 9,287 (3,362+441+441+1,681+1,681+1,681)

*** Timestep 0: t=0 years
Solving temperature system... 0 iterations.
Skipping s11 composition solve because RHS is zero.
Skipping s22 composition solve because RHS is zero.
Skipping s12 composition solve because RHS is zero.
Skipping sxx composition solve because RHS is zero.
Skipping syy composition solve because RHS is zero.
Skipping sxy composition solve because RHS is zero.
Rebuilding Stokes preconditioner...
Solving Stokes system... 27+0 iterations.
Relative nonlinear residual (Stokes system) after nonlinear iteration 1: 1
Expand All @@ -24,9 +24,9 @@ Number of degrees of freedom: 9,287 (3,362+441+441+1,681+1,681+1,681)

*** Timestep 1: t=1000 years
Solving temperature system... 0 iterations.
Solving s11 system ... 11 iterations.
Solving s22 system ... 11 iterations.
Solving s12 system ... 12 iterations.
Solving sxx system ... 11 iterations.
Solving syy system ... 11 iterations.
Solving sxy system ... 12 iterations.
Rebuilding Stokes preconditioner...
Solving Stokes system... 17+0 iterations.
Relative nonlinear residual (Stokes system) after nonlinear iteration 1: 0.00345201
Expand Down
2 changes: 1 addition & 1 deletion tests/viscoelastic_numerical_elastic_time_step.prm
Expand Up @@ -61,7 +61,7 @@ end
# Number and name of compositional fields
subsection Compositional fields
set Number of fields = 3
set Names of fields = s11, s22, s12
set Names of fields = sxx, syy, sxy
end

# Spatial domain of different compositional fields
Expand Down
12 changes: 6 additions & 6 deletions tests/viscoelastic_numerical_elastic_time_step/screen-output
Expand Up @@ -4,9 +4,9 @@ Number of degrees of freedom: 9,287 (3,362+441+441+1,681+1,681+1,681)

*** Timestep 0: t=0 years
Solving temperature system... 0 iterations.
Skipping s11 composition solve because RHS is zero.
Skipping s22 composition solve because RHS is zero.
Skipping s12 composition solve because RHS is zero.
Skipping sxx composition solve because RHS is zero.
Skipping syy composition solve because RHS is zero.
Skipping sxy composition solve because RHS is zero.
Rebuilding Stokes preconditioner...
Solving Stokes system... 27+0 iterations.
Relative nonlinear residual (Stokes system) after nonlinear iteration 1: 1
Expand All @@ -24,9 +24,9 @@ Number of degrees of freedom: 9,287 (3,362+441+441+1,681+1,681+1,681)

*** Timestep 1: t=1000 years
Solving temperature system... 0 iterations.
Solving s11 system ... 11 iterations.
Solving s22 system ... 11 iterations.
Solving s12 system ... 12 iterations.
Solving sxx system ... 11 iterations.
Solving syy system ... 11 iterations.
Solving sxy system ... 12 iterations.
Rebuilding Stokes preconditioner...
Solving Stokes system... 17+0 iterations.
Relative nonlinear residual (Stokes system) after nonlinear iteration 1: 0.00275809
Expand Down
2 changes: 1 addition & 1 deletion tests/viscoelastic_stress_averaging.prm
Expand Up @@ -64,7 +64,7 @@ end
# Number and name of compositional fields
subsection Compositional fields
set Number of fields = 3
set Names of fields = s11, s22, s12
set Names of fields = sxx, syy, sxy
end

# Spatial domain of different compositional fields
Expand Down
12 changes: 6 additions & 6 deletions tests/viscoelastic_stress_averaging/screen-output
Expand Up @@ -4,9 +4,9 @@ Number of degrees of freedom: 9,287 (3,362+441+441+1,681+1,681+1,681)

*** Timestep 0: t=0 years
Solving temperature system... 0 iterations.
Skipping s11 composition solve because RHS is zero.
Skipping s22 composition solve because RHS is zero.
Skipping s12 composition solve because RHS is zero.
Skipping sxx composition solve because RHS is zero.
Skipping syy composition solve because RHS is zero.
Skipping sxy composition solve because RHS is zero.
Rebuilding Stokes preconditioner...
Solving Stokes system... 27+0 iterations.
Relative nonlinear residual (Stokes system) after nonlinear iteration 1: 1
Expand All @@ -24,9 +24,9 @@ Number of degrees of freedom: 9,287 (3,362+441+441+1,681+1,681+1,681)

*** Timestep 1: t=1000 years
Solving temperature system... 0 iterations.
Solving s11 system ... 11 iterations.
Solving s22 system ... 11 iterations.
Solving s12 system ... 12 iterations.
Solving sxx system ... 11 iterations.
Solving syy system ... 11 iterations.
Solving sxy system ... 12 iterations.
Rebuilding Stokes preconditioner...
Solving Stokes system... 16+0 iterations.
Relative nonlinear residual (Stokes system) after nonlinear iteration 1: 0.00172532
Expand Down

0 comments on commit 0881424

Please sign in to comment.