diff --git a/doc/modules/changes/20240726_myhill b/doc/modules/changes/20240726_myhill
new file mode 100644
index 00000000000..6671eb3214e
--- /dev/null
+++ b/doc/modules/changes/20240726_myhill
@@ -0,0 +1,4 @@
+New: The CompositeViscoPlastic rheology in ASPECT now
+includes an implementation of elasticity.
+
+(Bob Myhill, 2024/07/26)
diff --git a/include/aspect/material_model/rheology/composite_visco_plastic.h b/include/aspect/material_model/rheology/composite_visco_plastic.h
index 8889479a907..790f34c1f9c 100644
--- a/include/aspect/material_model/rheology/composite_visco_plastic.h
+++ b/include/aspect/material_model/rheology/composite_visco_plastic.h
@@ -28,6 +28,7 @@
#include
#include
#include
+#include
#include
namespace aspect
@@ -91,6 +92,24 @@ namespace aspect
parse_parameters (ParameterHandler &prm,
const std::unique_ptr> &expected_n_phases_per_composition = nullptr);
+ /**
+ * Compute the inverse of the scalar elastic viscosity
+ * obtained from the elasticity rheology. The required scalar shear modulus is
+ * calculated by harmonically averaging the individual component shear moduli
+ * weighted by the @p volume_fractions of the components.
+ */
+ double
+ compute_inverse_kelvin_viscosity(const std::vector &volume_fractions) const;
+
+ /**
+ * Compute the effective viscoelastic strain rate used to calculate the
+ * viscosity.
+ */
+ SymmetricTensor<2, dim>
+ compute_effective_strain_rate(const SymmetricTensor<2, dim> &strain_rate,
+ const SymmetricTensor<2, dim> &elastic_stress,
+ const double inverse_kelvin_viscosity) const;
+
/**
* Compute the viscosity based on the composite viscous creep law.
* If @p n_phase_transitions_per_composition points to a vector of
@@ -103,7 +122,8 @@ namespace aspect
const double temperature,
const double grain_size,
const std::vector &volume_fractions,
- const SymmetricTensor<2,dim> &strain_rate,
+ const SymmetricTensor<2,dim> &effective_strain_rate,
+ const double inverse_kelvin_viscosity,
std::vector &partial_strain_rates,
const std::vector &phase_function_values = std::vector(),
const std::vector &n_phase_transitions_per_composition = std::vector()) const;
@@ -122,7 +142,8 @@ namespace aspect
const double temperature,
const double grain_size,
const std::vector &volume_fractions,
- const SymmetricTensor<2,dim> &strain_rate,
+ const SymmetricTensor<2,dim> &effective_strain_rate,
+ const double inverse_kelvin_viscosity,
std::vector &partial_strain_rates,
const std::vector &phase_function_values = std::vector(),
const std::vector &n_phase_transitions_per_composition = std::vector()) const;
@@ -135,6 +156,7 @@ namespace aspect
std::pair
calculate_isostress_log_strain_rate_and_derivative(const std::vector, 4>> &logarithmic_strain_rates_and_stress_derivatives,
const double viscoplastic_stress,
+ const double inverse_kelvin_viscosity,
std::vector &partial_strain_rates) const;
/**
@@ -185,6 +207,69 @@ namespace aspect
const double viscoplastic_stress,
std::vector &partial_strain_rates) const;
+ /**
+ * 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;
+ */
+
+ /**
+ * Given the stress of the previous time step in the material model inputs @p in,
+ * the elastic shear moduli @p average_elastic_shear_moduli at each point,
+ * and the (viscous) viscosities given in the material model outputs object @p out,
+ * fill a material model outputs objects with the elastic force terms, viscoelastic
+ * strain rate and viscous dissipation.
+ */
+ /*
+ void
+ fill_elastic_outputs (const MaterialModel::MaterialModelInputs &in,
+ const std::vector &average_elastic_shear_moduli,
+ MaterialModel::MaterialModelOutputs &out) const;
+ */
+
+ /**
+ * Given the stress of the previous time step in the material model inputs @p in,
+ * the elastic shear moduli @p average_elastic_shear_moduli at each point,
+ * and the (viscous) viscosities given in the material model outputs object @p out,
+ * fill a material model outputs (ElasticAdditionalOutputs) object with the
+ * average shear modulus, elastic viscosity, and the deviatoric stress of the current timestep.
+ */
+ /*
+ void
+ fill_elastic_additional_outputs (const MaterialModel::MaterialModelInputs &in,
+ const std::vector &average_elastic_shear_moduli,
+ MaterialModel::MaterialModelOutputs &out) const;
+ */
+
+ /**
+ * Given the stress of the previous time step in the material model inputs @p in,
+ * the elastic shear moduli @p average_elastic_shear_moduli at each point,
+ * and the (viscous) viscosities given in the material model outputs object @p out,
+ * compute an update to the elastic stresses and use it to fill the reaction terms
+ * material model output property.
+ */
+ void
+ fill_reaction_outputs (const MaterialModel::MaterialModelInputs &in,
+ const std::vector &average_elastic_shear_moduli,
+ MaterialModel::MaterialModelOutputs &out) const;
+
+ /**
+ * Given the stress of the previous time step in the material model inputs @p in,
+ * the elastic shear moduli @p average_elastic_shear_moduli at each point,
+ * and the (viscous) viscosities given in the material model outputs object @p out,
+ * compute the update to the elastic stresses of the previous timestep and use it
+ * to fill the reaction rates material model output property.
+ */
+ void
+ fill_reaction_rates (const MaterialModel::MaterialModelInputs &in,
+ const std::vector &average_elastic_shear_moduli,
+ MaterialModel::MaterialModelOutputs &out) const;
+
+ private:
/**
* Enumeration for selecting which type of viscosity averaging to use.
*/
@@ -197,6 +282,7 @@ namespace aspect
bool use_dislocation_creep;
bool use_peierls_creep;
bool use_drucker_prager;
+ bool use_elasticity;
/**
* Vector storing which flow mechanisms are active.
@@ -209,16 +295,31 @@ namespace aspect
* which is arranged in parallel with the viscoplastic elements and
* therefore does not contribute to the total strain rate.
*/
- static constexpr unsigned int n_decomposed_strain_rates = 5;
+ static constexpr unsigned int n_decomposed_strain_rates = 6;
+
+ /**
+ * The index of the hard damper in the decomposed strain rates.
+ * This is always the last element.
+ */
+ static constexpr unsigned int damper_strain_rate_index = 5;
+ static constexpr unsigned int isostrain_damper_strain_rate_index = 4;
+
+ /**
+ * The index of the elastic element in the decomposed strain rates.
+ * This is always the penultimate element.
+ */
+ static constexpr unsigned int elastic_strain_rate_index = 4;
/**
* Pointers to objects for computing deformation mechanism
* strain rates and effective viscosities.
*/
- std::unique_ptr> diffusion_creep;
+ std::unique_ptr>
+ diffusion_creep;
std::unique_ptr> dislocation_creep;
std::unique_ptr> peierls_creep;
std::unique_ptr> drucker_prager;
+ std::unique_ptr> elasticity;
/**
* The expected number of chemical compositional fields.
@@ -229,8 +330,14 @@ namespace aspect
* The maximum viscosity, imposed via an isoviscous damper
* in series with the composite viscoplastic element.
*/
+ double inverse_maximum_viscosity;
double maximum_viscosity;
+ /**
+ * The minimum viscosity, imposed via an isoviscous damper
+ * in parallel with the flow law components
+ */
+ double minimum_viscosity;
/**
* The viscosity of an isoviscous damper placed in parallel
@@ -273,6 +380,8 @@ namespace aspect
* Useful number to aid in adding together exponentials.
*/
const double logmin = std::log(std::numeric_limits::min());
+
+ static constexpr unsigned int n_independent_components = SymmetricTensor<2, dim>::n_independent_components;
};
}
}
diff --git a/include/aspect/material_model/rheology/elasticity.h b/include/aspect/material_model/rheology/elasticity.h
index f3f3698c8f6..79a8b5af1d5 100644
--- a/include/aspect/material_model/rheology/elasticity.h
+++ b/include/aspect/material_model/rheology/elasticity.h
@@ -112,6 +112,12 @@ namespace aspect
const std::vector &
get_elastic_shear_moduli () const;
+ /**
+ * Return the values of the damper viscosity used in the rheology model.
+ */
+ double
+ get_damper_viscosity() const;
+
/**
* Calculates the effective elastic viscosity (this is the equivalent viscosity of
* a material which was unstressed at the end of the previous timestep).
diff --git a/source/material_model/rheology/composite_visco_plastic.cc b/source/material_model/rheology/composite_visco_plastic.cc
index ab020ed431b..c94b302574a 100644
--- a/source/material_model/rheology/composite_visco_plastic.cc
+++ b/source/material_model/rheology/composite_visco_plastic.cc
@@ -20,6 +20,7 @@
#include
+#include
#include
#include
@@ -42,7 +43,8 @@ namespace aspect
// The composite visco plastic rheology calculates the decomposed strain
// rates for each of the following deformation mechanisms:
// diffusion creep, dislocation creep, Peierls creep,
- // Drucker-Prager plasticity and a constant (high) viscosity limiter.
+ // Drucker-Prager plasticity, Kelvin (damped) elasticity
+ // and a constant (high) viscosity limiter.
// The values are provided in this order as a vector of additional
// outputs. If the user declares one or more mechanisms inactive
// (by assigning use_mechanism = False) then the corresponding
@@ -56,6 +58,7 @@ namespace aspect
// names.emplace_back("edot_dislocation");
// names.emplace_back("edot_peierls");
// names.emplace_back("edot_drucker_prager");
+ // names.emplace_back("edot_kelvin");
// names.emplace_back("edot_limiter");
// return names;
// }
@@ -91,13 +94,43 @@ namespace aspect
+ template
+ double
+ CompositeViscoPlastic::compute_inverse_kelvin_viscosity(const std::vector &volume_fractions) const
+ {
+ double inverse_kelvin_viscosity = 0.;
+ if (this->get_parameters().enable_elasticity)
+ {
+ // Take the volume-weighted harmonic average of the individual component
+ // shear moduli, as required for isostress (Reuss) material averaging.
+ const std::vector &elastic_shear_moduli = elasticity->get_elastic_shear_moduli();
+ const double elastic_shear_modulus = MaterialUtilities::average_value(volume_fractions, elastic_shear_moduli, MaterialUtilities::harmonic);
+ inverse_kelvin_viscosity = 1./elasticity->calculate_elastic_viscosity(elastic_shear_modulus);
+ }
+ return inverse_kelvin_viscosity;
+ }
+
+
+
+ template
+ SymmetricTensor<2,dim>
+ CompositeViscoPlastic::compute_effective_strain_rate(const SymmetricTensor<2,dim> &strain_rate,
+ const SymmetricTensor<2,dim> &elastic_stress,
+ const double inverse_kelvin_viscosity) const
+ {
+ return strain_rate + (0.5 * elastic_stress * inverse_kelvin_viscosity);
+ }
+
+
+
template
double
CompositeViscoPlastic::compute_viscosity (const double pressure,
const double temperature,
const double grain_size,
const std::vector &volume_fractions,
- const SymmetricTensor<2,dim> &strain_rate,
+ const SymmetricTensor<2,dim> &effective_strain_rate,
+ const double inverse_kelvin_viscosity,
std::vector &partial_strain_rates,
const std::vector &phase_function_values,
const std::vector &n_phase_transitions_per_composition) const
@@ -118,7 +151,8 @@ namespace aspect
temperature,
grain_size,
volume_fractions,
- strain_rate,
+ effective_strain_rate,
+ inverse_kelvin_viscosity,
partial_strain_rates,
phase_function_values,
n_phase_transitions_per_composition);
@@ -130,7 +164,7 @@ namespace aspect
temperature,
grain_size,
volume_fractions,
- strain_rate,
+ effective_strain_rate,
partial_strain_rates,
phase_function_values,
n_phase_transitions_per_composition);
@@ -152,7 +186,8 @@ namespace aspect
const double temperature,
const double grain_size,
const std::vector &volume_fractions,
- const SymmetricTensor<2,dim> &strain_rate,
+ const SymmetricTensor<2,dim> &effective_strain_rate,
+ const double inverse_kelvin_viscosity,
std::vector &partial_strain_rates,
const std::vector &phase_function_values,
const std::vector &n_phase_transitions_per_composition) const
@@ -161,7 +196,7 @@ namespace aspect
// to prevent a division-by-zero, and a floating point exception.
// Otherwise, calculate the square-root of the norm of the second invariant of the deviatoric-
// strain rate (often simplified as epsilondot_ii)
- const double edot_ii = std::max(std::sqrt(std::max(-second_invariant(deviator(strain_rate)), 0.)),
+ const double edot_ii = std::max(std::sqrt(std::max(-second_invariant(deviator(effective_strain_rate)), 0.)),
minimum_strain_rate);
const double log_edot_ii = std::log(edot_ii);
@@ -285,6 +320,7 @@ namespace aspect
// accommodate strain at that creep stress.
std::pair log_edot_ii_and_deriv_iterate = calculate_isostress_log_strain_rate_and_derivative(log_edot_and_deriv,
viscoplastic_stress,
+ inverse_kelvin_viscosity,
partial_strain_rates);
double log_strain_rate_residual = log_edot_ii_and_deriv_iterate.first - log_edot_ii;
@@ -315,6 +351,7 @@ namespace aspect
// Compute the new log strain rate residual and log stress derivative
log_edot_ii_and_deriv_iterate = calculate_isostress_log_strain_rate_and_derivative(log_edot_and_deriv,
viscoplastic_stress,
+ inverse_kelvin_viscosity,
partial_strain_rates);
log_strain_rate_residual = log_edot_ii - log_edot_ii_and_deriv_iterate.first;
@@ -343,7 +380,11 @@ namespace aspect
// arranged in parallel with the viscoplastic elements.
// The total stress is equal to the sum of the viscoplastic stress and
// minimum stress.
- const double damper_stress = 2. * damper_viscosity * edot_ii;
+ double viscoplastic_strain_rate = 0.;
+ for (auto &i : active_flow_mechanisms)
+ viscoplastic_strain_rate += partial_strain_rates[i];
+
+ const double damper_stress = 2. * damper_viscosity * viscoplastic_strain_rate;
const double total_stress = viscoplastic_stress + damper_stress;
// 6) Return the effective creep viscosity using the total stress
@@ -356,6 +397,7 @@ namespace aspect
std::pair
CompositeViscoPlastic::calculate_isostress_log_strain_rate_and_derivative(const std::vector, 4>> &logarithmic_strain_rates_and_stress_derivatives,
const double viscoplastic_stress,
+ const double inverse_kelvin_viscosity,
std::vector &partial_strain_rates) const
{
// The total strain rate
@@ -396,12 +438,21 @@ namespace aspect
const double log_viscoplastic_strain_rate_derivative = weighted_stress_derivative_sum / viscoplastic_strain_rate_sum;
// Some opaque mathematics converts the viscoplastic strain rate to the total strain rate.
- const double f = viscoplastic_stress / (2. * maximum_viscosity);
- const double strain_rate = (strain_rate_scaling_factor * viscoplastic_strain_rate_sum) + f;
- partial_strain_rates[4] = strain_rate - viscoplastic_strain_rate_sum;
+ const double inverse_hard_and_elastic_viscosity = (inverse_maximum_viscosity + inverse_kelvin_viscosity);
+ const double f = 0.5 * inverse_hard_and_elastic_viscosity * viscoplastic_stress;
+ const double g = 1. + minimum_viscosity * inverse_kelvin_viscosity;
+ const double strain_rate = (strain_rate_scaling_factor * g * viscoplastic_strain_rate_sum) + f;
+ const double strain_rate_hard_and_elastic = strain_rate - viscoplastic_strain_rate_sum;
+
+ // Elastic strain rate
+ partial_strain_rates[elastic_strain_rate_index] = strain_rate_hard_and_elastic * (inverse_kelvin_viscosity / inverse_hard_and_elastic_viscosity);
+
+ // Hard damper strain rate
+ partial_strain_rates[damper_strain_rate_index] = strain_rate_hard_and_elastic - partial_strain_rates[elastic_strain_rate_index];
+
// And the partial derivative of the log *total* strain rate
// with respect to log *viscoplastic* stress follows as
- const double log_strain_rate_derivative = (strain_rate_scaling_factor * viscoplastic_strain_rate_sum * log_viscoplastic_strain_rate_derivative + f) / strain_rate;
+ const double log_strain_rate_derivative = (strain_rate_scaling_factor * g * viscoplastic_strain_rate_sum * log_viscoplastic_strain_rate_derivative + f) / strain_rate;
return std::make_pair(std::log(strain_rate), log_strain_rate_derivative);
}
@@ -430,7 +481,10 @@ namespace aspect
// a certain (small) fraction.
if (volume_fractions[composition] > 2. * std::numeric_limits::epsilon())
{
- std::vector partial_strain_rates_composition(n_decomposed_strain_rates, 0.);
+ // There is no elastic component allowed in isostrain models,
+ // so the size of vector partial_strain_rates_composition is
+ // one smaller than partial_strain_rates.
+ std::vector partial_strain_rates_composition(n_decomposed_strain_rates-1, 0.);
viscosity += (volume_fractions[composition]
* compute_composition_viscosity (pressure,
temperature,
@@ -440,8 +494,11 @@ namespace aspect
partial_strain_rates_composition,
phase_function_values,
n_phase_transitions_per_composition));
- for (unsigned int j=0; j < n_decomposed_strain_rates; ++j)
+ for (auto &j : active_flow_mechanisms)
partial_strain_rates[j] += volume_fractions[composition] * partial_strain_rates_composition[j];
+
+ // Shift the strain rate for the hard viscosity damper into the last position in partial_strain_rates.
+ partial_strain_rates[damper_strain_rate_index] += volume_fractions[composition] * partial_strain_rates_composition[isostrain_damper_strain_rate_index];
}
else
{
@@ -451,7 +508,7 @@ namespace aspect
viscosity /= total_volume_fraction;
- for (unsigned int j=0; j < n_decomposed_strain_rates; ++j)
+ for (auto &j : active_flow_mechanisms)
partial_strain_rates[j] /= total_volume_fraction;
return viscosity;
@@ -649,7 +706,7 @@ namespace aspect
// Some opaque mathematics converts the viscoplastic strain rate to the total strain rate.
const double f = viscoplastic_stress / (2. * maximum_viscosity);
const double strain_rate = (strain_rate_scaling_factor * viscoplastic_strain_rate_sum) + f;
- partial_strain_rates[4] = strain_rate - viscoplastic_strain_rate_sum;
+ partial_strain_rates[isostrain_damper_strain_rate_index] = strain_rate - viscoplastic_strain_rate_sum;
// And the partial derivative of the log *total* strain rate
// with respect to log *viscoplastic* stress follows as
const double log_strain_rate_derivative = (strain_rate_scaling_factor * viscoplastic_strain_rate_sum * log_viscoplastic_strain_rate_derivative + f) / strain_rate;
@@ -659,6 +716,292 @@ namespace aspect
+ // TODO: UNCOMMENT
+ /*
+ template
+ void
+ CompositeViscoPlastic::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();
+ out.additional_outputs.push_back(
+ std::make_unique> (n_points));
+ }
+
+ // We need to modify the shear heating outputs to correctly account for elastic stresses.
+ if (out.template get_additional_output>() == nullptr)
+ {
+ const unsigned int n_points = out.n_evaluation_points();
+ out.additional_outputs.push_back(
+ std::make_unique> (n_points));
+ }
+
+ // Create the ReactionRateOutputs that are necessary for the operator splitting
+ // step (either on the fields or directly on the particles)
+ // that sets both sets of stresses to the total stress of the
+ // previous timestep.
+ if (out.template get_additional_output>() == nullptr &&
+ (this->get_parameters().use_operator_splitting || (this->get_parameters().mapped_particle_properties).count(this->introspection().compositional_index_for_name("ve_stress_xx"))))
+ {
+ const unsigned int n_points = out.n_evaluation_points();
+ out.additional_outputs.push_back(
+ std::make_unique>(n_points, this->n_compositional_fields()));
+ }
+ }
+
+
+
+ template
+ void
+ CompositeViscoPlastic::fill_elastic_outputs (const MaterialModel::MaterialModelInputs &in,
+ const std::vector &inverse_kelvin_viscosities,
+ MaterialModel::MaterialModelOutputs &out) const
+ {
+ // Create a reference to the structure for the elastic outputs.
+ // The structure is created during the Stokes assembly.
+ MaterialModel::ElasticOutputs
+ *elastic_out = out.template get_additional_output>();
+
+ // Create a reference to the structure for the prescribed shear heating outputs.
+ // The structure is created during the advection assembly.
+ HeatingModel::PrescribedShearHeatingOutputs
+ *heating_out = out.template get_additional_output>();
+
+ if (elastic_out == nullptr && heating_out == nullptr)
+ return;
+
+ // TODO should a RHS term be a separate MaterialProperties?
+ if (in.requests_property(MaterialProperties::additional_outputs))
+ {
+ // The viscosity should be averaged if material averaging is applied.
+ std::vector effective_creep_viscosities;
+ if (this->get_parameters().material_averaging != MaterialAveraging::none)
+ {
+ MaterialModelOutputs out_copy(out.n_evaluation_points(),
+ this->introspection().n_compositional_fields);
+ out_copy.viscosities = out.viscosities;
+
+ const MaterialAveraging::AveragingOperation averaging_operation_for_viscosity =
+ get_averaging_operation_for_viscosity(this->get_parameters().material_averaging);
+ MaterialAveraging::average(averaging_operation_for_viscosity,
+ in.current_cell,
+ this->introspection().quadratures.velocities,
+ this->get_mapping(),
+ in.requested_properties,
+ out_copy);
+
+ effective_creep_viscosities = out_copy.viscosities;
+ }
+ else
+ effective_creep_viscosities = out.viscosities;
+
+ const unsigned int stress_start_index = this->introspection().compositional_index_for_name("ve_stress_xx");
+
+ for (unsigned int i=0; i < in.n_evaluation_points(); ++i)
+ {
+ const SymmetricTensor<2, dim> deviatoric_strain_rate = deviator(in.strain_rate[i]);
+
+ // Get stress from timestep $t$ rotated and advected into the current
+ // timestep $t+\Delta t_c$ from the compositional fields.
+ // This function is only evaluated during the assembly of the Stokes equations
+ // (the force term goes into the rhs of the momentum equation).
+ // This happens after the advection equations have been solved, and hence in.composition
+ // contains the rotated and advected stresses $tau^{0adv}$.
+ // 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));
+
+ // Average effective creep viscosity
+ // Use the viscosity corresponding to the stresses selected above.
+ // out.viscosities is computed during the assembly of the Stokes equations
+ // based on the current_linearization_point. This means that it will be updated after every
+ // nonlinear Stokes iteration.
+ // The effective creep viscosity has already been scaled with the timestep ratio dtc/dte.
+ const double effective_creep_viscosity = effective_creep_viscosities[i];
+
+ // The force term is computed as:
+ // $\frac{-\eta_{effcreep} \tau_{0adv}}{\eta_{e}}$, where $\eta_{effcreep}$ is the
+ // current harmonic average of the viscous and elastic viscosity, or the yield stress
+ // divided by two times the second invariant of the deviatoric strain rate.
+ // In case the computational timestep differs from the elastic timestep,
+ // linearly interpolate between the two.
+ // The elastic viscosity has also already been scaled with the timestep ratio.
+ const double viscosity_ratio = effective_creep_viscosity / inverse_kelvin_viscosities[i];
+
+ if (elastic_out != nullptr)
+ {
+ elastic_out->elastic_force[i] = -1. * viscosity_ratio * stress_0_advected;
+
+ // The viscoelastic strain rate is needed only when the Newton method is selected.
+ const typename Parameters::NonlinearSolver::Kind nonlinear_solver = this->get_parameters().nonlinear_solver;
+ if ((nonlinear_solver == Parameters::NonlinearSolver::iterated_Advection_and_Newton_Stokes) ||
+ (nonlinear_solver == Parameters::NonlinearSolver::single_Advection_iterated_Newton_Stokes))
+ elastic_out->viscoelastic_strain_rate[i] = compute_effective_strain_rate(in.strain_rate[i], stress_0_advected, inverse_kelvin_viscosities[i]);
+ }
+
+ // The shear heating rate (used by the heating model) depends not only
+ // on the stress and strain rates from the viscous flow laws, but also from
+ // the three viscous dampers. Usually, the heating from the dampers will be
+ // minor, but for consistency they should be included.
+
+ // 1) The total stress
+ const SymmetricTensor<2, dim> stress = 2. * effective_creep_viscosity * deviatoric_strain_rate + viscosity_ratio * stress_0_advected;
+
+ // 2) The elastic damper
+ const SymmetricTensor<2, dim> kelvin_strain_rate = 0.5 * (stress - stress_0_advected) * inverse_kelvin_viscosities[i];
+ const double damper_viscosity = elasticity->get_damper_viscosity();
+ const double damper_power_density = 2. * damper_viscosity * kelvin_strain_rate * kelvin_strain_rate;
+
+ // 3) Total other work (assuming incompressibility)
+ // This includes the power density from the other two dampers
+ const SymmetricTensor<2, dim> visco_plastic_strain_rate = deviatoric_strain_rate - kelvin_strain_rate;
+ const double viscoplastic_power_density = stress * visco_plastic_strain_rate;
+ // If compressible,
+ // visco_plastic_strain_rate = visco_plastic_strain_rate -
+ // 1. / 3. * trace(visco_plastic_strain_rate) * unit_symmetric_tensor();
+
+ // The shear heating term needs to account for the elastic stress, but only the visco_plastic strain rate.
+ // This is best computed here, and stored for later use by the heating model.
+ if (heating_out != nullptr)
+ heating_out->prescribed_shear_heating_rates[i] = viscoplastic_power_density + damper_power_density;
+ }
+
+ }
+ }
+
+
+
+ template
+ void
+ CompositeViscoPlastic::fill_elastic_additional_outputs (const MaterialModel::MaterialModelInputs &in,
+ const std::vector &inverse_kelvin_viscosities,
+ MaterialModel::MaterialModelOutputs &out) const
+ {
+ // Create a reference to the structure for the elastic additional outputs
+ MaterialModel::ElasticAdditionalOutputs
+ *elastic_additional_out = out.template get_additional_output>();
+
+ if (elastic_additional_out == nullptr || !in.requests_property(MaterialProperties::additional_outputs))
+ return;
+
+ const unsigned int stress_start_index = this->introspection().compositional_index_for_name("ve_stress_xx");
+ const double dtc = this->get_timestep();
+ const double elastic_damper_viscosity = elasticity->get_damper_viscosity();
+
+ for (unsigned int i = 0; i < in.n_evaluation_points(); ++i)
+ {
+ const double effective_viscosity = out.viscosities[i];
+ const SymmetricTensor<2, dim> deviatoric_strain_rate = deviator(in.strain_rate[i]);
+ 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));
+
+ // Apply the stress update to get the total deviatoric stress of timestep t.
+ // This is not the stress passing through the elastic element, because of the elastic damper.
+ elastic_additional_out->deviatoric_stress[i] = 2. * effective_viscosity * deviatoric_strain_rate + effective_viscosity * inverse_kelvin_viscosities[i] * stress_0_advected;
+ elastic_additional_out->elastic_viscosity[i] = 1. / inverse_kelvin_viscosities[i];
+ elastic_additional_out->elastic_shear_moduli[i] = (elastic_additional_out->elastic_viscosity[i] - elastic_damper_viscosity)/dtc;
+ }
+ }
+ */
+
+
+ // Rotate the elastic stresses of the previous timestep $t$ into the current timestep $t+dtc$.
+ template
+ void
+ CompositeViscoPlastic::fill_reaction_outputs (const MaterialModel::MaterialModelInputs &in,
+ const std::vector &,
+ MaterialModel::MaterialModelOutputs &out) const
+ {
+ elasticity->fill_reaction_outputs(in, std::vector(), out);
+ }
+
+
+
+ // The following function computes the reaction rates for the operator
+ // splitting step that at the beginning of the new timestep $t+dtc$ updates the
+ // stored compositions $tau^{0\mathrm{adv}}$ at time $t$ to $tau^{t}$.
+ // This update consists of the stress change resulting from system evolution,
+ // but does not advect or rotate the elastic stress tensor. Advection is done by
+ // solving the advection equation and the elastic stress tensor is rotated through
+ // the source term (reaction_terms) of that same equation.
+ template
+ void
+ CompositeViscoPlastic::fill_reaction_rates (const MaterialModel::MaterialModelInputs &in,
+ const std::vector &inverse_kelvin_viscosities,
+ MaterialModel::MaterialModelOutputs &out) const
+ {
+ ReactionRateOutputs *reaction_rate_out = out.template get_additional_output>();
+
+ if (reaction_rate_out == nullptr)
+ return;
+
+ // At the moment when the reaction rates are required (at the beginning of the timestep),
+ // the solution vector 'solution' holds the stress from the previous timestep,
+ // advected into the new position of the previous timestep, so $\tau^{t}_{0adv}$.
+ // This is the same as the vector 'old_solution' holds. At later moments during the current timestep,
+ // 'solution' will hold the current_linearization_point instead of the solution of the previous timestep.
+ //
+ // In case fields are used to track the stresses, MaterialModelInputs are based on 'solution'
+ // when calling the MaterialModel for the reaction rates. When particles are used, MaterialModelInputs
+ // for this function are filled with the old_solution (including for the strain rate), except for the
+ // compositions that represent the stress tensor components, these are taken directly from the
+ // particles. As the particles are restored to their pre-advection location at the beginning of
+ // each nonlinear iteration, their values and positions correspond to the old solution.
+ // This means that in both cases we can use 'in' to get to the $\tau^{t}_{0adv}$ and velocity/strain rate of the
+ // previous timestep.
+ if (in.current_cell.state() == IteratorState::valid && this->get_timestep_number() > 0 && in.requests_property(MaterialProperties::reaction_rates))
+ {
+ const unsigned int stress_start_index = this->introspection().compositional_index_for_name("ve_stress_xx");
+ const double elastic_damper_viscosity = elasticity->get_damper_viscosity();
+ for (unsigned int i = 0; i < in.n_evaluation_points(); ++i)
+ {
+ // Set all reaction rates to zero
+ for (unsigned int c = 0; c < in.composition[i].size(); ++c)
+ reaction_rate_out->reaction_rates[i][c] = 0.0;
+
+ // Get $\tau^{0adv}$ of the previous timestep t from the compositional fields.
+ // This stress includes the rotation and advection of the previous timestep,
+ // i.e., the reaction term (which prescribes the change in stress due to rotation
+ // over the previous timestep) has already been applied during the previous timestep.
+ const SymmetricTensor<2, dim> stress_0_t (Utilities::Tensors::to_symmetric_tensor(&in.composition[i][stress_start_index],
+ &in.composition[i][stress_start_index]+n_independent_components));
+
+ // $\eta^{t}_{effcreep}$. This viscosity has been calculated with the timestep_ratio dtc/dte.
+ const double effective_creep_viscosity = out.viscosities[i];
+
+ // Compute the total stress at time t.
+ const SymmetricTensor<2, dim>
+ stress_t = 2. * effective_creep_viscosity * deviator(in.strain_rate[i])
+ + effective_creep_viscosity * inverse_kelvin_viscosities[i] * stress_0_t;
+
+ // Fill reaction rates.
+ // During this timestep, the reaction rates will be multiplied
+ // with the current timestep size to turn the rate of change into a change.
+ // However, this update belongs
+ // to the previous timestep. Therefore we divide by the
+ // current timestep and multiply with the previous one.
+ // When multiplied with the current timestep, this will give
+ // (rate * previous_dt / current_dt) * current_dt = rate * previous_dt = previous_change.
+ // previous_change = (1 - eta_d/eta_kel)*(stress_t - stress_0_t).
+ // To compute the rate we should return to the operator splitting scheme,
+ // we therefore divide the change in stress by the current timestep current_dt (= dtc).
+
+ const SymmetricTensor<2, dim> stress_update = (1. - (elastic_damper_viscosity*inverse_kelvin_viscosities[i])) * (stress_t - stress_0_t) / this->get_timestep();
+
+ Utilities::Tensors::unroll_symmetric_tensor_into_array(stress_update,
+ &reaction_rate_out->reaction_rates[i][stress_start_index],
+ &reaction_rate_out->reaction_rates[i][stress_start_index]+n_independent_components);
+ }
+ }
+ }
+
+
+
// Overload the + operator to act on two pairs of doubles.
std::pair operator+(const std::pair &x, const std::pair &y)
{
@@ -694,6 +1037,10 @@ namespace aspect
Patterns::Bool (),
"Whether to include Drucker-Prager plasticity in the composite rheology formulation.");
+ prm.declare_entry ("Include elasticity in composite rheology", "false",
+ Patterns::Bool (),
+ "Whether to include elasticity in the composite rheology formulation.");
+
// Diffusion creep parameters
Rheology::DiffusionCreep::declare_parameters(prm);
@@ -706,6 +1053,9 @@ namespace aspect
// Drucker Prager parameters
Rheology::DruckerPragerPower::declare_parameters(prm);
+ // Elastic parameters
+ Rheology::Elasticity::declare_parameters(prm);
+
// Some of the parameters below are shared with the subordinate
// rheology models (diffusion, dislocation, ...),
// and will already have been declared. This is fine, the deal.II
@@ -754,6 +1104,7 @@ namespace aspect
// Read maximum viscosity parameter
maximum_viscosity = prm.get_double ("Maximum viscosity");
+ inverse_maximum_viscosity = 1. / maximum_viscosity;
// Process minimum viscosity parameter
// In this rheology model, there are two viscous dampers designed
@@ -771,7 +1122,7 @@ namespace aspect
// When scaling the viscoplastic strain up to the total strain,
// eta_max / (eta_max - eta_min) becomes a useful value,
// which we here call the "strain_rate_scaling_factor".
- const double minimum_viscosity = prm.get_double("Minimum viscosity");
+ minimum_viscosity = prm.get_double("Minimum viscosity");
strain_rate_scaling_factor = maximum_viscosity / (maximum_viscosity - minimum_viscosity);
damper_viscosity = maximum_viscosity * minimum_viscosity / (maximum_viscosity - minimum_viscosity);
@@ -823,7 +1174,30 @@ namespace aspect
}
AssertThrow(active_flow_mechanisms.size() > 0,
- ExcMessage("You need to include at least one deformation mechanism."));
+ ExcMessage("You need to include at least one non-elastic deformation mechanism."));
+
+ // Elastic parameters
+ // Elasticity is not treated as a flow mechanism,
+ // so we do not push back the active_flow_mechanisms vector.
+ use_elasticity = prm.get_bool ("Include elasticity in composite rheology");
+ if (use_elasticity)
+ {
+ elasticity = std::make_unique>();
+ elasticity->initialize_simulator (this->get_simulator());
+ elasticity->parse_parameters(prm);
+ AssertThrow(viscosity_averaging_scheme == ViscosityAveraging::isostress,
+ ExcMessage("Elasticity in the CompositeViscoPlastic rheology "
+ "requires that the 'Viscosity averaging scheme' be "
+ "set to isostress."));
+ AssertThrow(prm.get ("Use fixed elastic time step") == "false",
+ ExcMessage("Elasticity in the CompositeViscoPlastic rheology "
+ "requires that 'Use fixed elastic time step' be "
+ "set to false."));
+ AssertThrow(std::abs(prm.get_double ("Stabilization time scale factor") - 1.) < std::numeric_limits::min(),
+ ExcMessage("Elasticity in the CompositeViscoPlastic rheology "
+ "requires that the 'Stabilization time scale factor' be "
+ "set to 1."));
+ }
}
}
diff --git a/source/material_model/rheology/elasticity.cc b/source/material_model/rheology/elasticity.cc
index 57ec63a1f2b..86b2585232e 100644
--- a/source/material_model/rheology/elasticity.cc
+++ b/source/material_model/rheology/elasticity.cc
@@ -466,6 +466,15 @@ namespace aspect
+ template
+ double
+ Elasticity::get_damper_viscosity () const
+ {
+ return elastic_damper_viscosity;
+ }
+
+
+
template
double
Elasticity::
diff --git a/tests/composite_viscous_outputs.cc b/tests/composite_viscous_outputs.cc
index 7b5efbcc0b7..6f8f65301af 100644
--- a/tests/composite_viscous_outputs.cc
+++ b/tests/composite_viscous_outputs.cc
@@ -98,6 +98,7 @@ void f(const aspect::SimulatorAccess &simulator_access,
double temperature;
const double pressure = 1.e9;
const double grain_size = 1.e-3;
+ const double inverse_kelvin_viscosity = 0.;
SymmetricTensor<2,dim> strain_rate;
strain_rate[0][0] = -1e-11;
strain_rate[0][1] = 0.;
@@ -126,7 +127,7 @@ void f(const aspect::SimulatorAccess &simulator_access,
temperature = 1000. + i*100.;
// Compute the viscosity
- viscosity = composite_creep->compute_viscosity(pressure, temperature, grain_size, volume_fractions, strain_rate, partial_strain_rates);
+ viscosity = composite_creep->compute_viscosity(pressure, temperature, grain_size, volume_fractions, strain_rate, inverse_kelvin_viscosity, partial_strain_rates);
total_strain_rate = std::accumulate(partial_strain_rates.begin(), partial_strain_rates.end(), 0.);
// The creep strain rate is calculated by subtracting the strain rate
diff --git a/tests/composite_viscous_outputs_isostress.cc b/tests/composite_viscous_outputs_isostress.cc
index fcd433f9247..a8073a5e2c7 100644
--- a/tests/composite_viscous_outputs_isostress.cc
+++ b/tests/composite_viscous_outputs_isostress.cc
@@ -100,6 +100,7 @@ void f(const aspect::SimulatorAccess &simulator_access,
double temperature;
const double pressure = 1.e9;
const double grain_size = 1.e-3;
+ const double inverse_kelvin_viscosity = 0.;
SymmetricTensor<2,dim> strain_rate;
strain_rate[0][0] = -1e-11;
strain_rate[0][1] = 0.;
@@ -108,7 +109,7 @@ void f(const aspect::SimulatorAccess &simulator_access,
strain_rate[2][1] = 0.;
strain_rate[2][2] = 0.;
- std::cout << "temperature (K) eta (Pas) creep stress (Pa) edot_ii (/s) edot_ii fractions (diff, disl, prls, drpr, max)" << std::endl;
+ std::cout << "temperature (K) eta (Pas) creep stress (Pa) edot_ii (/s) edot_ii fractions (diff, disl, prls, drpr, kel, max)" << std::endl;
// Loop through strain rates, tracking whether there is a discrepancy in
// the decomposed strain rates.
@@ -121,21 +122,21 @@ void f(const aspect::SimulatorAccess &simulator_access,
double disl_stress;
double prls_stress;
double drpr_stress;
- std::vector partial_strain_rates(5, 0.);
+ std::vector partial_strain_rates(6, 0.);
for (unsigned int i=0; i <= 10; i++)
{
temperature = 1000. + i*100.;
// Compute the viscosity
- viscosity = composite_creep->compute_viscosity(pressure, temperature, grain_size, volume_fractions, strain_rate, partial_strain_rates);
+ viscosity = composite_creep->compute_viscosity(pressure, temperature, grain_size, volume_fractions, strain_rate, inverse_kelvin_viscosity, partial_strain_rates);
total_strain_rate = std::accumulate(partial_strain_rates.begin(), partial_strain_rates.end(), 0.);
// The creep strain rate is calculated by subtracting the strain rate
// of the max viscosity dashpot from the total strain rate
// The creep stress is then calculated by subtracting the stress running
// through the strain rate limiter from the total stress
- creep_strain_rate = total_strain_rate - partial_strain_rates[4];
+ creep_strain_rate = total_strain_rate - partial_strain_rates[5];
creep_stress = 2.*(viscosity*total_strain_rate - lim_visc*creep_strain_rate);
// Print the output
diff --git a/tests/composite_viscous_outputs_isostress/screen-output b/tests/composite_viscous_outputs_isostress/screen-output
index 545da023131..c804c667689 100644
--- a/tests/composite_viscous_outputs_isostress/screen-output
+++ b/tests/composite_viscous_outputs_isostress/screen-output
@@ -2,18 +2,18 @@
Loading shared library <./libcomposite_viscous_outputs_isostress.debug.so>
* Connecting signals
-temperature (K) eta (Pas) creep stress (Pa) edot_ii (/s) edot_ii fractions (diff, disl, prls, drpr, max)
-1000 3.45602e+19 6.89205e+08 1e-11 1.29846e-06 0.00363514 0.0654305 0.930933 3.45602e-09
-1100 2.98863e+19 5.95725e+08 1e-11 7.23285e-05 0.835897 0.163395 0.0006365 2.98863e-09
-1200 7.70568e+18 1.52114e+08 1e-11 0.000594407 0.999399 6.44893e-06 1.44515e-33 7.70568e-10
-1300 2.39282e+18 4.58564e+07 1e-11 0.00338081 0.996619 3.19443e-09 1.32277e-59 2.39282e-10
-1400 9.18168e+17 1.63634e+07 1e-11 0.0149609 0.985039 1.26589e-11 5.56077e-82 9.18169e-11
-1500 4.32083e+17 6.64167e+06 1e-11 0.0538305 0.94617 2.09885e-13 1.4634e-101 4.32084e-11
-1600 2.47246e+17 2.94493e+06 1e-11 0.161077 0.838923 8.74429e-15 3.2006e-119 2.47246e-11
-1700 1.67379e+17 1.34758e+06 1e-11 0.397345 0.602655 5.42942e-16 3.38184e-136 1.67378e-11
-1800 1.28447e+17 568948 1e-11 0.749975 0.250025 2.2682e-17 6.38432e-155 1.28447e-11
-1900 1.09563e+17 191256 1e-11 0.962698 0.0373021 2.59166e-19 1.35594e-178 1.09563e-11
-2000 1.02964e+17 59280.1 1e-11 0.99654 0.00345962 2.26225e-21 4.97674e-204 1.02965e-11
+temperature (K) eta (Pas) creep stress (Pa) edot_ii (/s) edot_ii fractions (diff, disl, prls, drpr, kel, max)
+1000 3.45602e+19 6.89205e+08 1e-11 1.29846e-06 0.00363514 0.0654305 0.930933 0 3.45602e-09
+1100 2.98863e+19 5.95725e+08 1e-11 7.23285e-05 0.835897 0.163395 0.0006365 0 2.98863e-09
+1200 7.70568e+18 1.52114e+08 1e-11 0.000594407 0.999399 6.44893e-06 1.44515e-33 0 7.70568e-10
+1300 2.39282e+18 4.58564e+07 1e-11 0.00338081 0.996619 3.19443e-09 1.32277e-59 0 2.39282e-10
+1400 9.18168e+17 1.63634e+07 1e-11 0.0149609 0.985039 1.26589e-11 5.56077e-82 0 9.18169e-11
+1500 4.32083e+17 6.64167e+06 1e-11 0.0538305 0.94617 2.09885e-13 1.4634e-101 0 4.32084e-11
+1600 2.47246e+17 2.94493e+06 1e-11 0.161077 0.838923 8.74429e-15 3.2006e-119 0 2.47246e-11
+1700 1.67379e+17 1.34758e+06 1e-11 0.397345 0.602655 5.42942e-16 3.38184e-136 0 1.67378e-11
+1800 1.28447e+17 568948 1e-11 0.749975 0.250025 2.2682e-17 6.38432e-155 0 1.28447e-11
+1900 1.09563e+17 191256 1e-11 0.962698 0.0373021 2.59166e-19 1.35594e-178 0 1.09563e-11
+2000 1.02964e+17 59280.1 1e-11 0.99654 0.00345962 2.26225e-21 4.97674e-204 0 1.02965e-11
OK
Number of active cells: 100 (on 1 levels)
Number of degrees of freedom: 6,857 (3,969+242+1,323+1,323)
diff --git a/tests/composite_viscous_outputs_no_peierls.cc b/tests/composite_viscous_outputs_no_peierls.cc
index 4b15fb65b54..e9bc6b0a158 100644
--- a/tests/composite_viscous_outputs_no_peierls.cc
+++ b/tests/composite_viscous_outputs_no_peierls.cc
@@ -92,7 +92,8 @@ void f(const aspect::SimulatorAccess &simulator_access,
double temperature;
const double pressure = 1.e9;
const double grain_size = 1.e-3;
- SymmetricTensor<2,dim> strain_rate;
+ const double inverse_kelvin_viscosity = 0.;
+ SymmetricTensor<2, dim> strain_rate;
strain_rate[0][0] = -1e-11;
strain_rate[0][1] = 0.;
strain_rate[1][1] = 1e-11;
@@ -119,7 +120,7 @@ void f(const aspect::SimulatorAccess &simulator_access,
temperature = 1000. + i*100.;
// Compute the viscosity
- viscosity = composite_creep->compute_viscosity(pressure, temperature, grain_size, volume_fractions, strain_rate, partial_strain_rates);
+ viscosity = composite_creep->compute_viscosity(pressure, temperature, grain_size, volume_fractions, strain_rate, inverse_kelvin_viscosity, partial_strain_rates);
total_strain_rate = std::accumulate(partial_strain_rates.begin(), partial_strain_rates.end(), 0.);
// The creep strain rate is calculated by subtracting the strain rate
diff --git a/tests/composite_viscous_outputs_w_elasticity.cc b/tests/composite_viscous_outputs_w_elasticity.cc
new file mode 100644
index 00000000000..bdc85469c37
--- /dev/null
+++ b/tests/composite_viscous_outputs_w_elasticity.cc
@@ -0,0 +1,267 @@
+/*
+ Copyright (C) 2022 - 2023 by the authors of the ASPECT code.
+
+ This file is part of ASPECT.
+
+ ASPECT is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2, or (at your option)
+ any later version.
+
+ ASPECT is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with ASPECT; see the file LICENSE. If not see
+ .
+*/
+
+#include
+#include
+#include
+#include
+
+template
+void f(const aspect::SimulatorAccess &simulator_access,
+ aspect::Assemblers::Manager &)
+{
+ // This function tests whether the composite creep rheology is producing
+ // the correct composite viscosity and partial strain rates corresponding to
+ // the different creep mechanisms incorporated into the rheology.
+ // It is assumed that each individual creep mechanism has already been tested.
+
+ using namespace aspect::MaterialModel;
+
+ // First, we set up a few objects which are used by the rheology model.
+ aspect::ParameterHandler prm;
+
+ const std::vector list_of_composition_names = simulator_access.introspection().get_composition_names();
+ auto n_phases = std::make_unique>(1); // 1 phase per composition
+ const unsigned int composition = 0;
+ const std::vector volume_fractions = {0.6, 0.4};
+ const std::vector phase_function_values = std::vector();
+ const std::vector n_phase_transitions_per_composition = std::vector(1);
+
+ // Next, we initialise instances of the composite rheology and
+ // individual creep mechanisms.
+ std::unique_ptr> composite_creep;
+ composite_creep = std::make_unique>();
+ composite_creep->initialize_simulator (simulator_access.get_simulator());
+ composite_creep->declare_parameters(prm);
+ prm.set("Viscosity averaging scheme", "isostress");
+ prm.set("Include diffusion creep in composite rheology", "true");
+ prm.set("Include dislocation creep in composite rheology", "true");
+ prm.set("Include Peierls creep in composite rheology", "true");
+ prm.set("Include Drucker Prager plasticity in composite rheology", "true");
+ prm.set("Include elasticity in composite rheology", "true");
+ prm.set("Peierls creep flow law", "viscosity approximation");
+ prm.set("Maximum yield stress", "5e8");
+ prm.set("Use fixed elastic time step", "false");
+ composite_creep->parse_parameters(prm);
+
+ std::unique_ptr> diffusion_creep;
+ diffusion_creep = std::make_unique>();
+ diffusion_creep->initialize_simulator (simulator_access.get_simulator());
+ diffusion_creep->declare_parameters(prm);
+ diffusion_creep->parse_parameters(prm);
+
+ std::unique_ptr> dislocation_creep;
+ dislocation_creep = std::make_unique>();
+ dislocation_creep->initialize_simulator (simulator_access.get_simulator());
+ dislocation_creep->declare_parameters(prm);
+ dislocation_creep->parse_parameters(prm);
+
+ std::unique_ptr> peierls_creep;
+ peierls_creep = std::make_unique>();
+ peierls_creep->initialize_simulator (simulator_access.get_simulator());
+ peierls_creep->declare_parameters(prm);
+ peierls_creep->parse_parameters(prm);
+
+ std::unique_ptr> drucker_prager_power;
+ drucker_prager_power = std::make_unique>();
+ drucker_prager_power->initialize_simulator (simulator_access.get_simulator());
+ drucker_prager_power->declare_parameters(prm);
+ prm.set("Maximum yield stress", "5e8");
+ drucker_prager_power->parse_parameters(prm);
+ Rheology::DruckerPragerParameters p = drucker_prager_power->compute_drucker_prager_parameters(composition, phase_function_values, n_phase_transitions_per_composition);
+
+ // The creep components are arranged in series with each other.
+ // This package of components is then arranged in parallel with
+ // a strain rate limiter with a constant viscosity lim_visc.
+ // The whole system is then arranged in series with a viscosity limiter with
+ // viscosity max_visc.
+ // lim_visc is equal to (min_visc*max_visc)/(max_visc - min_visc)
+ double min_visc = prm.get_double("Minimum viscosity");
+ double max_visc = prm.get_double("Maximum viscosity");
+ double lim_visc = (min_visc*max_visc)/(max_visc - min_visc);
+
+ // Assign values to the variables which will be passed to compute_viscosity
+ // The test involves pure shear calculations at 1 GPa and variable temperature
+ double temperature;
+ const double pressure = 1.e9;
+ const double grain_size = 1.e-3;
+ const double inverse_kelvin_viscosity = composite_creep->compute_inverse_kelvin_viscosity(volume_fractions);
+ SymmetricTensor<2,dim> strain_rate;
+ strain_rate[0][0] = -2e-11;
+ strain_rate[0][1] = 0.;
+ strain_rate[1][1] = 2e-11;
+ strain_rate[2][0] = 0.;
+ strain_rate[2][1] = 0.;
+ strain_rate[2][2] = 0.;
+
+ SymmetricTensor<2,dim> elastic_stress;
+ elastic_stress[0][0] = 2e-11 / inverse_kelvin_viscosity;
+ elastic_stress[0][1] = 0.;
+ elastic_stress[1][1] = -2e-11 / inverse_kelvin_viscosity;
+ elastic_stress[2][0] = 0.;
+ elastic_stress[2][1] = 0.;
+ elastic_stress[2][2] = 0.;
+
+ SymmetricTensor<2,dim> effective_strain_rate = composite_creep->compute_effective_strain_rate(strain_rate, elastic_stress, inverse_kelvin_viscosity);
+
+
+ std::cout << "temperature (K) eta (Pas) creep stress (Pa) edot_ii (/s) edot_ii fractions (diff, disl, prls, drpr, kel, max)" << std::endl;
+
+ // Loop through strain rates, tracking whether there is a discrepancy in
+ // the decomposed strain rates.
+ bool error = false;
+ double viscosity;
+ double total_strain_rate;
+ double creep_strain_rate;
+ double creep_stress;
+ double diff_stress;
+ double disl_stress;
+ double prls_stress;
+ double drpr_stress;
+ std::vector partial_strain_rates(6, 0.);
+
+ for (unsigned int i=0; i <= 10; i++)
+ {
+ temperature = 1000. + i*100.;
+
+ // Compute the viscosity
+ viscosity = composite_creep->compute_viscosity(pressure, temperature, grain_size, volume_fractions, effective_strain_rate, inverse_kelvin_viscosity, partial_strain_rates);
+ total_strain_rate = std::accumulate(partial_strain_rates.begin(), partial_strain_rates.end(), 0.);
+
+ // The creep strain rate is calculated by subtracting the strain rate
+ // of the max viscosity dashpot from the total strain rate
+ // The creep stress is then calculated by subtracting the stress running
+ // through the strain rate limiter from the total stress
+ creep_strain_rate = total_strain_rate - partial_strain_rates[4] - partial_strain_rates[5];
+ creep_stress = 2.*(viscosity*total_strain_rate - lim_visc*creep_strain_rate);
+
+ // Print the output
+ std::cout << temperature << ' ' << viscosity << ' ' << creep_stress << ' ' << total_strain_rate;
+ for (unsigned int i=0; i < partial_strain_rates.size(); ++i)
+ {
+ std::cout << ' ' << partial_strain_rates[i]/total_strain_rate;
+ }
+ std::cout << std::endl;
+
+ // The following lines test that each individual creep mechanism
+ // experiences the same creep stress
+
+ // Each creep mechanism should experience the same stress
+ diff_stress = 2.*partial_strain_rates[0]*diffusion_creep->compute_viscosity(pressure, temperature, grain_size, composition);
+ disl_stress = 2.*partial_strain_rates[1]*dislocation_creep->compute_viscosity(partial_strain_rates[1], pressure, temperature, composition);
+ prls_stress = 2.*partial_strain_rates[2]*peierls_creep->compute_viscosity(partial_strain_rates[2], pressure, temperature, composition);
+ if (partial_strain_rates[3] > 0.)
+ {
+ drpr_stress = 2.*partial_strain_rates[3]*drucker_prager_power->compute_viscosity(p.cohesion,
+ p.angle_internal_friction,
+ pressure,
+ partial_strain_rates[3],
+ p.max_yield_stress);
+ }
+ else
+ {
+ drpr_stress = creep_stress;
+ }
+
+ if ((std::fabs((diff_stress - creep_stress)/creep_stress) > 1e-6)
+ || (std::fabs((disl_stress - creep_stress)/creep_stress) > 1e-6)
+ || (std::fabs((prls_stress - creep_stress)/creep_stress) > 1e-6)
+ || (std::fabs((drpr_stress - creep_stress)/creep_stress) > 1e-6))
+ {
+ error = true;
+ std::cout << " creep stress: " << creep_stress;
+ std::cout << " diffusion stress: " << diff_stress;
+ std::cout << " dislocation stress: " << disl_stress;
+ std::cout << " peierls stress: " << prls_stress;
+ std::cout << " drucker prager stress: " << drpr_stress << std::endl;
+ }
+ }
+
+ if (error)
+ {
+ std::cout << " Error: The individual creep stresses differ by more than the required tolerance." << std::endl;
+ std::cout << "Some parts of the test were not successful." << std::endl;
+ }
+ else
+ {
+ std::cout << "OK" << std::endl;
+ }
+
+}
+
+template <>
+void f(const aspect::SimulatorAccess<2> &,
+ aspect::Assemblers::Manager<2> &)
+{
+ AssertThrow(false,dealii::ExcInternalError());
+}
+
+template
+void signal_connector (aspect::SimulatorSignals &signals)
+{
+ using namespace dealii;
+ std::cout << "* Connecting signals" << std::endl;
+ signals.set_assemblers.connect (std::bind(&f,
+ std::placeholders::_1,
+ std::placeholders::_2));
+}
+
+
+using namespace aspect;
+
+
+void declare_parameters(const unsigned int dim,
+ ParameterHandler &prm)
+{
+ prm.enter_subsection("Formulation");
+ {
+ prm.declare_entry("Enable elasticity", "true", Patterns::Bool());
+ }
+ prm.leave_subsection();
+
+ prm.enter_subsection("Compositional fields");
+ {
+ if (dim==2)
+ {
+ prm.declare_entry("Number of fields","4", Patterns::Integer());
+ prm.declare_entry("Names of fields","ve_stress_xx, ve_stress_yy, ve_stress_xy, foreground", Patterns::Anything());
+ }
+ else
+ {
+ prm.declare_entry("Number of fields","7", Patterns::Integer());
+ prm.declare_entry("Names of fields","ve_stress_xx, ve_stress_yy, ve_stress_zz, ve_stress_xy, ve_stress_xz, ve_stress_yz, foreground", Patterns::Anything());
+ }
+ }
+ prm.leave_subsection();
+}
+
+
+
+void parameter_connector ()
+{
+ SimulatorSignals<2>::declare_additional_parameters.connect (&declare_parameters);
+ SimulatorSignals<3>::declare_additional_parameters.connect (&declare_parameters);
+}
+
+
+
+ASPECT_REGISTER_SIGNALS_CONNECTOR(signal_connector<2>,
+ signal_connector<3>)
+ASPECT_REGISTER_SIGNALS_PARAMETER_CONNECTOR(parameter_connector)
diff --git a/tests/composite_viscous_outputs_w_elasticity.prm b/tests/composite_viscous_outputs_w_elasticity.prm
new file mode 100644
index 00000000000..10f467f0318
--- /dev/null
+++ b/tests/composite_viscous_outputs_w_elasticity.prm
@@ -0,0 +1,85 @@
+# This test checks whether the composite viscous rheology
+# plugin produces the correct composite viscosity and
+# decomposed strain rates with elasticity enabled.
+
+set Additional shared libraries = tests/libcomposite_viscous_outputs_w_elasticity.so
+set Dimension = 3
+set End time = 0
+set Use years in output instead of seconds = true
+set Nonlinear solver scheme = single Advection, iterated Stokes
+
+subsection Formulation
+ set Enable elasticity = true
+end
+
+# Model geometry (100x100 km, 10 km spacing)
+subsection Geometry model
+ set Model name = box
+
+ subsection Box
+ set X repetitions = 10
+ set Y repetitions = 10
+ set X extent = 100e3
+ set Y extent = 100e3
+ end
+end
+
+# Mesh refinement specifications
+subsection Mesh refinement
+ set Initial adaptive refinement = 0
+ set Initial global refinement = 0
+ set Time steps between mesh refinement = 0
+end
+
+# Boundary classifications (fixed T boundaries, prescribed velocity)
+
+# Temperature boundary and initial conditions
+subsection Boundary temperature model
+ set Fixed temperature boundary indicators = bottom, top, left, right
+ set List of model names = box
+
+ subsection Box
+ set Bottom temperature = 273
+ set Left temperature = 273
+ set Right temperature = 273
+ set Top temperature = 273
+ end
+end
+
+# Velocity on boundaries characterized by functions
+subsection Boundary velocity model
+ set Prescribed velocity boundary indicators = bottom y: function, top y: function, left x: function, right x: function
+
+ subsection Function
+ set Variable names = x,y,z
+ set Function constants = m=0.0005, year=1
+ set Function expression = if (x<50e3 , -1*m/year, 1*m/year); if (y<50e3 , 1*m/year, -1*m/year);0
+ end
+end
+
+subsection Initial temperature model
+ set Model name = function
+
+ subsection Function
+ set Function expression = 273
+ end
+end
+
+# Material model
+subsection Material model
+ set Model name = visco plastic
+
+ subsection Visco Plastic
+ set Angles of internal friction = 30.
+ set Use fixed elastic time step = false
+ end
+end
+
+# Gravity model
+subsection Gravity model
+ set Model name = vertical
+
+ subsection Vertical
+ set Magnitude = 10.0
+ end
+end
diff --git a/tests/composite_viscous_outputs_w_elasticity/screen-output b/tests/composite_viscous_outputs_w_elasticity/screen-output
new file mode 100644
index 00000000000..1db09e16011
--- /dev/null
+++ b/tests/composite_viscous_outputs_w_elasticity/screen-output
@@ -0,0 +1,42 @@
+
+Loading shared library <./libcomposite_viscous_outputs_w_elasticity.debug.so>
+
+* Connecting signals
+temperature (K) eta (Pas) creep stress (Pa) edot_ii (/s) edot_ii fractions (diff, disl, prls, drpr, kel, max)
+1000 3.45481e+19 6.8899e+08 1e-11 1.29806e-06 0.00363118 0.0652213 0.916549 0.0145971 3.45481e-09
+1100 2.98014e+19 5.94052e+08 1e-11 7.21254e-05 0.827709 0.159074 0.000553006 0.0125916 2.98014e-09
+1200 7.69828e+18 1.51972e+08 1e-11 0.000593854 0.996147 6.39569e-06 1.3794e-33 0.00325265 7.69828e-10
+1300 2.39205e+18 4.58431e+07 1e-11 0.00337984 0.995609 3.18669e-09 1.30376e-59 0.00101068 2.39205e-10
+1400 9.18038e+17 1.63615e+07 1e-11 0.0149592 0.984653 1.26476e-11 5.5297e-82 0.000387886 9.18038e-11
+1500 4.32047e+17 6.64131e+06 1e-11 0.0538275 0.94599 2.09799e-13 1.45943e-101 0.000182547 4.32047e-11
+1600 2.47231e+17 2.94483e+06 1e-11 0.161072 0.838824 8.74218e-15 3.19521e-119 0.000104459 2.47231e-11
+1700 1.6737e+17 1.34754e+06 1e-11 0.397334 0.602596 5.42836e-16 3.37707e-136 7.07165e-05 1.6737e-11
+1800 1.28441e+17 568929 1e-11 0.74995 0.249996 2.2677e-17 6.37367e-155 5.42684e-05 1.28441e-11
+1900 1.09558e+17 191248 1e-11 0.962657 0.0372966 2.59096e-19 1.35307e-178 4.629e-05 1.09558e-11
+2000 1.0296e+17 59277.5 1e-11 0.996497 0.0034591 2.26166e-21 4.96601e-204 4.35021e-05 1.0296e-11
+OK
+Number of active cells: 100 (on 1 levels)
+Number of degrees of freedom: 14,795 (3,969+242+1,323+1,323+1,323+1,323+1,323+1,323+1,323+1,323)
+
+*** Timestep 0: t=0 years, dt=0 years
+ Solving temperature system... 0 iterations.
+ Skipping ve_stress_xx composition solve because RHS is zero.
+ Skipping ve_stress_yy composition solve because RHS is zero.
+ Skipping ve_stress_zz composition solve because RHS is zero.
+ Skipping ve_stress_xy composition solve because RHS is zero.
+ Skipping ve_stress_xz composition solve because RHS is zero.
+ Skipping ve_stress_yz composition solve because RHS is zero.
+ Skipping foreground composition solve because RHS is zero.
+ Solving Stokes system... 43+0 iterations.
+ Relative nonlinear residual (Stokes system) after nonlinear iteration 1: 1
+
+ Solving Stokes system... 0+0 iterations.
+ Relative nonlinear residual (Stokes system) after nonlinear iteration 2: 9.98903e-08
+
+
+ Postprocessing:
+
+Termination requested by criterion: end time
+
+
+