From 066becb5b3c58610acd93048e2cb8e396604678a Mon Sep 17 00:00:00 2001 From: Bob Myhill Date: Fri, 26 Jul 2024 11:38:19 +0100 Subject: [PATCH] added elasticity to CompositeViscoPlastic --- doc/modules/changes/20240726_myhill | 4 + .../rheology/composite_visco_plastic.h | 117 ++++- .../material_model/rheology/elasticity.h | 6 + .../rheology/composite_visco_plastic.cc | 412 +++++++++++++++++- source/material_model/rheology/elasticity.cc | 9 + tests/composite_viscous_outputs.cc | 3 +- tests/composite_viscous_outputs_isostress.cc | 9 +- .../screen-output | 24 +- tests/composite_viscous_outputs_no_peierls.cc | 5 +- .../composite_viscous_outputs_w_elasticity.cc | 267 ++++++++++++ ...composite_viscous_outputs_w_elasticity.prm | 85 ++++ .../screen-output | 42 ++ 12 files changed, 941 insertions(+), 42 deletions(-) create mode 100644 doc/modules/changes/20240726_myhill create mode 100644 tests/composite_viscous_outputs_w_elasticity.cc create mode 100644 tests/composite_viscous_outputs_w_elasticity.prm create mode 100644 tests/composite_viscous_outputs_w_elasticity/screen-output 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 2c227308676..a13d094f22a 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; @@ -342,8 +379,12 @@ namespace aspect // because this material model also includes a viscosity damper // 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 - partial_strain_rates[4]); + // damper stress. + 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; @@ -599,7 +656,7 @@ namespace aspect // because this material model also includes a viscosity damper // arranged in parallel with the viscoplastic elements. // The total stress is equal to the sum of the viscoplastic stress and - // minimum stress. + // damper stress. const double damper_stress = 2. * damper_viscosity * (edot_ii - partial_strain_rates[4]); const double total_stress = viscoplastic_stress + damper_stress; @@ -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"); AssertThrow(minimum_viscosity > 0, ExcMessage("Minimum viscosity needs to be larger than zero.")); @@ -832,7 +1183,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 72e1d8d8b8a..0172ab81fa7 100644 --- a/source/material_model/rheology/elasticity.cc +++ b/source/material_model/rheology/elasticity.cc @@ -470,6 +470,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 + + +