diff --git a/include/aspect/material_model/rheology/elasticity.h b/include/aspect/material_model/rheology/elasticity.h index 7d48e227a4d..7a3b07e8a4d 100644 --- a/include/aspect/material_model/rheology/elasticity.h +++ b/include/aspect/material_model/rheology/elasticity.h @@ -34,9 +34,9 @@ namespace aspect using namespace dealii; /** - * Additional output fields for the elastic shear modulus to be added to - * the MaterialModel::MaterialModelOutputs structure and filled in the - * MaterialModel::Interface::evaluate() function. + * Additional output fields for the elastic shear modulus and other + * elastic outputs to be added to the MaterialModel::MaterialModelOutputs + * structure and filled in the MaterialModel::Interface::evaluate() function. */ template class ElasticAdditionalOutputs : public NamedAdditionalMaterialOutputs @@ -69,7 +69,9 @@ namespace aspect /** * The deviatoric stress of the current timestep, so including - * the rotation, advection and stress update. + * the rotation, advection and stress update, at the evaluation points + * passed to the instance of MaterialModel::Interface::evaluate() + * that fills the current object. */ std::vector> deviatoric_stress; }; @@ -96,8 +98,9 @@ namespace aspect parse_parameters (ParameterHandler &prm); /** - * Create the additional material model output objects that contain the - * elastic shear moduli, elastic viscosity, timestep ratio and reaction rates. + * Create the two additional material model output objects that contain the + * elastic shear moduli, elastic viscosity, ratio of computational to elastic timestep, + * and deviatoric stress of the current timestep and the reaction rates. */ void create_elastic_additional_outputs (MaterialModel::MaterialModelOutputs &out) const; @@ -119,7 +122,8 @@ namespace aspect * the elastic shear moduli @p average_elastic_shear_moduli a each point, * and the (viscous) viscosities given in the material model outputs object @p out, * fill a material model outputs (ElasticAdditionalOutputs) objects with the - * average shear modulus, elastic viscosity and ratio of computational to elastic timestep. + * average shear modulus, elastic viscosity, ratio of computational to elastic timestep, + * and the deviatoric stress of the current timestep. */ void fill_elastic_additional_outputs (const MaterialModel::MaterialModelInputs &in, @@ -158,7 +162,7 @@ namespace aspect get_elastic_shear_moduli () const; /** - * Calculates the effective elastic viscosity (this is the equivalent viscosity of + * Calculate the effective elastic viscosity (this is the equivalent viscosity of * a material which was unstressed at the end of the previous timestep). */ double @@ -189,8 +193,6 @@ namespace aspect const double viscosity_pre_yield, const double shear_modulus) const; - - /** * Compute the elastic time step. */ @@ -208,7 +210,7 @@ namespace aspect /** * Get the stored stress of the previous timestep. For fields, use a * composition evaluator of the old solution. For particles, get the - * stress directly from the particles that is stored in in.composition. + * stress directly from the particles, which is available from in.composition. */ std::vector> retrieve_stress_previous_timestep (const MaterialModel::MaterialModelInputs &in, @@ -234,7 +236,7 @@ namespace aspect bool use_fixed_elastic_time_step; /** - * Double for fixed elastic time step value, read from parameter file + * Double for fixed elastic time step value, read from parameter file. */ double fixed_elastic_time_step; @@ -249,8 +251,8 @@ namespace aspect /** * We cache the evaluators that are necessary to evaluate the velocity - * gradients and the old compositions. They are required to compute the elastic stresses, but - * are not provided by the material model. + * gradients and the old compositions. They are required to compute the elastic stresses, + * but are not provided by the material model. * By caching the evaluators, we can avoid recreating them every time we need them. */ mutable std::unique_ptr> evaluator; diff --git a/include/aspect/material_model/rheology/visco_plastic.h b/include/aspect/material_model/rheology/visco_plastic.h index 0b221484fb5..4b420d3e62a 100644 --- a/include/aspect/material_model/rheology/visco_plastic.h +++ b/include/aspect/material_model/rheology/visco_plastic.h @@ -71,7 +71,8 @@ namespace aspect std::vector friction_angles; /** - * The plastic yield stress. + * The current plastic yield stress, depending on composition, + * pressure, and strain. */ std::vector yield_stresses; diff --git a/include/aspect/material_model/visco_plastic.h b/include/aspect/material_model/visco_plastic.h index b8921753f59..874e9d6dec9 100644 --- a/include/aspect/material_model/visco_plastic.h +++ b/include/aspect/material_model/visco_plastic.h @@ -77,9 +77,12 @@ namespace aspect * representing these 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 ve_stress_xx, ve_stress_yy - * and ve_stress_xy. In 3D, the first six compositional fields must be labeled - * ve_stress_xx, ve_stress_yy, ve_stress_zz, ve_stress_xy, ve_stress_xz, ve_stress_yz. + * 3+3 consecutive compositional fields must be labeled ve_stress_xx, ve_stress_yy + * ve_stress_xy, ve_stress_xx_old, ve_stress_yy_old, ve_stress_xy_old. + * In 3D, six consecutive compositional fields must be labeled + * ve_stress_xx, ve_stress_yy, ve_stress_zz, ve_stress_xy, ve_stress_xz, ve_stress_yz, + * ve_stress_xx_old, ve_stress_yy_old, ve_stress_zz_old, ve_stress_xy_old, + * ve_stress_xz_old, ve_stress_yz_old. * * Combining this viscoelasticity implementation with non-linear viscous flow * and plasticity produces a constitutive relationship commonly referred to @@ -95,10 +98,11 @@ namespace aspect * * The overview below directly follows Moresi et al. (2003) eqns. 23-32. * However, an important distinction between this material model and - * the studies above is the use of compositional fields, rather than + * the studies above is the option to use of compositional fields, rather than * particles, to track individual components of the viscoelastic stress - * tensor. The material model will be updated when an option to track - * and calculate viscoelastic stresses with particles is implemented. + * tensor. Calculating viscoelastic stresses with particles is also implemented, + * and can be switched on by using particles with the particle property + * 'elastic stress'. * Moresi et al. (2003) begins (eqn. 23) by writing the deviatoric * rate of deformation ($\hat{D}$) as the sum of elastic * ($\hat{D_{e}}$) and viscous ($\hat{D_{v}}$) components: @@ -130,12 +134,11 @@ namespace aspect * W^{t}\tau^{t} + \tau^{t}W^{t}$. * In this material model, the size of the time step above ($\Delta t^{e}$) * can be specified as the numerical time step size or an independent fixed time - * step. If the latter case is selected, the user has an option to apply a - * stress averaging scheme to account for the differences between the numerical - * and fixed elastic time step (eqn. 32). If one selects to use a fixed elastic time - * step throughout the model run, an equal numerical and elastic time step can be - * achieved by using CFL and maximum time step values that restrict the numerical - * time step to the fixed elastic time step. + * step. If the latter case is selected, a stress averaging scheme is applied + * to account for the differences between the numerical + * and fixed elastic time step (eqn. 32). A fixed computational time + * step throughout the model run can be achieved by using a large CFL number and + * smaller maximum time step values that restrict the numerical time step to a specific time. * * The formulation above allows rewriting the total rate of deformation (eqn. 29) as * $\tau^{t + \Delta t^{e}} = \eta_{eff} \left ( @@ -150,8 +153,10 @@ namespace aspect * viscosity is reduced relative to the initial viscosity. * * 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 \Delta t^{e}} \tau^{t}$. + * an elastic force term (eqn. 30 updated to the version in Farrington et al. 2014) + * using stresses from the previous time step rotated and advected into the current + * time step: + * $F^{e,t} = -\frac{\eta_{eff}}{\mu \Delta t^{e}} \tau^{0adv}$. * This force term is added onto the right-hand side force vector in the * system of equations. * diff --git a/include/aspect/material_model/viscoelastic.h b/include/aspect/material_model/viscoelastic.h index beda817fc61..1a37530d0a0 100644 --- a/include/aspect/material_model/viscoelastic.h +++ b/include/aspect/material_model/viscoelastic.h @@ -45,9 +45,12 @@ namespace aspect * representing these 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 ve_stress_xx, ve_stress_yy - * and ve_stress_xy. In 3D, the first six compositional fields must be labeled - * ve_stress_xx, ve_stress_yy, ve_stress_zz, ve_stress_xy, ve_stress_xz, ve_stress_yz. + * 3+3 consecutive compositional fields must be labeled ve_stress_xx, ve_stress_yy, + * ve_stress_xy, ve_stress_xx_old, ve_stress_yy_old, ve_stress_xy_old. + * In 3D, 6+6 consecutive compositional fields must be labeled + * ve_stress_xx, ve_stress_yy, ve_stress_zz, ve_stress_xy, ve_stress_xz, ve_stress_yz, + * ve_stress_xx_old, ve_stress_yy_old, ve_stress_zz_old, ve_stress_xy_old, + * ve_stress_xz_old, ve_stress_yz_old. * * Expanding the model to include non-linear viscous flow (e.g., * diffusion/dislocation creep) and plasticity would produce a constitutive @@ -63,10 +66,11 @@ namespace aspect * * The overview below directly follows Moresi et al. (2003) eqns. 23-32. * However, an important distinction between this material model and - * the studies above is the use of compositional fields, rather than + * the studies above is the option to use of compositional fields, rather than * particles, to track individual components of the viscoelastic stress - * tensor. The material model will be updated when an option to track - * and calculate viscoelastic stresses with particles is implemented. + * tensor. Calculating viscoelastic stresses with particles is also implemented, + * and can be switched on by using particles with the particle property + * 'elastic stress'. * * Moresi et al. (2003) begins (eqn. 23) by writing the deviatoric * rate of deformation ($\hat{D}$) as the sum of elastic @@ -92,11 +96,11 @@ namespace aspect * W^{t}\tau^{t} + \tau^{t}W^{t}$. * In this material model, the size of the time step above ($\Delta t^{e}$) * can be specified as the numerical time step size or an independent fixed time - * step. If the latter case is a selected, the user has an option to apply a - * stress averaging scheme to account for the differences between the numerical - * and fixed elastic time step (eqn. 32). If one selects to use a fixed elastic time - * step throughout the model run, this can still be achieved by using CFL and - * maximum time step values that restrict the numerical time step to a specific time. + * step. If the latter case is selected, a stress averaging scheme is applied + * to account for the differences between the numerical + * and fixed elastic time step (eqn. 32). A fixed computational time + * step throughout the model run can be achieved by using a large CFL number and + * smaller maximum time step values that restrict the numerical time step to a specific time. * * The formulation above allows rewriting the total rate of deformation (eqn. 29) as * $\tau^{t + \Delta t^{e}} = \eta_{eff} \left ( @@ -110,9 +114,11 @@ namespace aspect * The magnitude of the shear modulus thus controls how much the effective * viscosity is reduced relative to the initial viscosity. * - * 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 \Delta t^{e}} \tau^{t}$. + * Elastic effects are introduced into the governing Stokes equations through + * an elastic force term (eqn. 30 updated to the version in Farrington et al. 2014) + * using stresses from the previous time step rotated and advected into the current + * time step: + * $F^{e,t} = -\frac{\eta_{eff}}{\mu \Delta t^{e}} \tau^{0adv}$. * This force term is added onto the right-hand side force vector in the * system of equations. * @@ -128,9 +134,11 @@ namespace aspect * 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 + * ``background, composition1, composition2...''. 3+3 (2D) + * or 6+6 (3D) consecutive composition fields should correspond to components + * of the elastic stress tensor of the previous timestep rotated and advected + * into the current timestep and the tensor of the previous timestep advected + * into the current timestep only. 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, diff --git a/source/material_model/rheology/elasticity.cc b/source/material_model/rheology/elasticity.cc index b060870ecce..96fc80c63ae 100644 --- a/source/material_model/rheology/elasticity.cc +++ b/source/material_model/rheology/elasticity.cc @@ -39,6 +39,8 @@ namespace aspect std::vector make_elastic_additional_outputs_names() { std::vector names; + // We do not output all the additional output to + // visualization output, only the ones listed here. names.emplace_back("elastic_shear_modulus"); names.emplace_back("elastic_viscosity"); return names; @@ -64,8 +66,6 @@ namespace aspect (void)idx; // suppress warning in release mode AssertIndexRange (idx, 2); - // We do not output the timestep_ratio, - // as it is constant over the domain at a given timestep. switch (idx) { case 0: @@ -110,8 +110,9 @@ namespace aspect "the model to run the user must select 'true' or 'false'."); prm.declare_entry ("Fixed elastic time step", "1.e3", Patterns::Double (0.), - "The fixed elastic time step $dte$. Units: years if the " - "'Use years in output instead of seconds' parameter is set; " + "The fixed elastic time step $dte$. It is always used during the first " + "timestep; afterwards on if 'Used fixed elastic time step' is true. " + "Units: years if the 'Use years in output instead of seconds' parameter is set; " "seconds otherwise."); prm.declare_entry ("Stabilization time scale factor", "1.", Patterns::Double (1.), @@ -184,7 +185,7 @@ namespace aspect AssertThrow(elastic_damper_viscosity == 0., ExcMessage("The viscoelastic material model and the visco-plastic material model with elasticity enabled require " "that no elastic damping is applied.")); - // An update of the stored stresses is done in an operator splitting step or by the particle property 'elastic stress'. + // An update of the stored stresses is done in an operator splitting step for fields or by the particle property 'elastic stress'. AssertThrow(this->get_parameters().use_operator_splitting || (this->get_parameters().mapped_particle_properties).count(this->introspection().compositional_index_for_name("ve_stress_xx")), ExcMessage("The viscoelastic material model and the visco-plastic material model with elasticity enabled require " "operator splitting for stresses tracked on compositional fields or the particle property 'elastic stress' " @@ -205,9 +206,9 @@ namespace aspect std::vector stress_field_names = this->introspection().get_names_for_fields_of_type(CompositionalFieldDescription::stress); std::vector stress_field_indices = this->introspection().get_indices_for_fields_of_type(CompositionalFieldDescription::stress); - // As long as the FEPointEvaluation requires a consecutive range of indices + // We require a consecutive range of indices (for example for FEPointEvaluation) // to extract the fields representing the viscoelastic stress tensor components, - // check that they are listed without interruption by other fields. + // so check that they are listed without interruption by other fields. // They do not, however, have to be the first fields listed. AssertThrow(((stress_field_indices[2*n_independent_components-1] - stress_field_indices[0]) == (2*n_independent_components-1)), ExcMessage("Rheology model Elasticity requires that the compositional fields representing stress tensor components are listed in consecutive order.")); @@ -307,6 +308,8 @@ namespace aspect void Elasticity::create_elastic_additional_outputs (MaterialModel::MaterialModelOutputs &out) const { + // Create the ElasticAdditionalOutputs that include the average shear modulus, elastic + // viscosity, timestep ratio and total deviatoric stress of the current timestep. if (out.template get_additional_output>() == nullptr) { const unsigned int n_points = out.n_evaluation_points(); @@ -415,7 +418,6 @@ namespace aspect // Only at the beginning of the next timestep do we add the stress update of the // current timestep to the stress stored in the compositional fields, giving // $\tau{t+\Delta t_c}$ with $t+\Delta t_c$ being the current timestep. - const SymmetricTensor<2,dim> stress_0_advected (Utilities::Tensors::to_symmetric_tensor(&in.composition[i][stress_start_index], &in.composition[i][stress_start_index]+n_independent_components)); @@ -493,7 +495,7 @@ namespace aspect { const double eta = out.viscosities[i]; const SymmetricTensor<2, dim> deviatoric_strain_rate = deviator(in.strain_rate[i]); - const SymmetricTensor<2,dim> stress_0 (Utilities::Tensors::to_symmetric_tensor(&in.composition[i][stress_start_index], + const SymmetricTensor<2,dim> stress_0_advected (Utilities::Tensors::to_symmetric_tensor(&in.composition[i][stress_start_index], &in.composition[i][stress_start_index]+n_independent_components)); const SymmetricTensor<2,dim> stress_old (Utilities::Tensors::to_symmetric_tensor(&in.composition[i][stress_start_index+n_independent_components], &in.composition[i][stress_start_index+n_independent_components]+n_independent_components)); @@ -501,7 +503,7 @@ namespace aspect const double elastic_viscosity = calculate_elastic_viscosity(average_elastic_shear_moduli[i]); // Apply the stress update to get the total deviatoric stress of timestep t. - elastic_additional_out->deviatoric_stress[i] = 2. * eta * deviatoric_strain_rate + eta / elastic_viscosity * stress_0 + (1. - timestep_ratio) * (1. - eta / elastic_viscosity) * stress_old; + elastic_additional_out->deviatoric_stress[i] = 2. * eta * deviatoric_strain_rate + eta / elastic_viscosity * stress_0_advected + (1. - timestep_ratio) * (1. - eta / elastic_viscosity) * stress_old; elastic_additional_out->elastic_shear_moduli[i] = average_elastic_shear_moduli[i]; elastic_additional_out->elastic_viscosity[i] = elastic_viscosity; elastic_additional_out->timestep_ratio[i] = timestep_ratio; @@ -709,6 +711,8 @@ namespace aspect return dte; } + + template double Elasticity::calculate_timestep_ratio() const @@ -805,6 +809,8 @@ namespace aspect return edot_deviator; } + + template std::vector> Elasticity::