diff --git a/doc/modules/changes/20240717e_myhill b/doc/modules/changes/20240717e_myhill new file mode 100644 index 00000000000..61c331814b6 --- /dev/null +++ b/doc/modules/changes/20240717e_myhill @@ -0,0 +1,5 @@ +Changed: The CompositeViscoPlastic rheology in ASPECT +now uses a log-stress Newton scheme, which requires +fewer iterations to converge. +
+(Bob Myhill, 2024/07/17) diff --git a/include/aspect/material_model/rheology/composite_visco_plastic.h b/include/aspect/material_model/rheology/composite_visco_plastic.h index 650bf979132..02b6e107b92 100644 --- a/include/aspect/material_model/rheology/composite_visco_plastic.h +++ b/include/aspect/material_model/rheology/composite_visco_plastic.h @@ -102,18 +102,14 @@ namespace aspect const std::vector &n_phase_transitions_per_composition = std::vector()) const; /** - * Compute the strain rate and first stress derivative - * as a function of stress based on the composite viscous creep law. + * Compute the total strain rate and the first derivative of log strain rate + * with respect to log viscoplastic stress from individual log strain rate components. + * Also updates the partial_strain_rates vector */ std::pair - compute_strain_rate_and_derivative (const double creep_stress, - const double pressure, - const double temperature, - const double grain_size, - const DiffusionCreepParameters diffusion_creep_parameters, - const DislocationCreepParameters dislocation_creep_parameters, - const PeierlsCreepParameters peierls_creep_parameters, - const DruckerPragerParameters drucker_prager_parameters) const; + calculate_log_strain_rate_and_derivative(const std::array, 4> &logarithmic_strain_rates_and_stress_derivatives, + const double viscoplastic_stress, + std::vector &partial_strain_rates) const; private: @@ -125,6 +121,11 @@ namespace aspect bool use_peierls_creep; bool use_drucker_prager; + /** + * Vector storing which flow mechanisms are active + */ + std::vector active_flow_mechanisms; + /** * Pointers to objects for computing deformation mechanism * strain rates and effective viscosities. @@ -134,15 +135,59 @@ namespace aspect std::unique_ptr> peierls_creep; std::unique_ptr> drucker_prager; - DruckerPragerParameters drucker_prager_parameters; + /** + * The expected number of chemical compositional fields. + */ + unsigned int number_of_chemical_compositions; - unsigned int number_of_compositions; - double minimum_viscosity; + /** + * The maximum viscosity, imposed via an isoviscous damper + * in series with the composite viscoplastic element + */ double maximum_viscosity; - double min_strain_rate; - double strain_rate_residual_threshold; + + /** + * The viscosity of an isoviscous damper placed in parallel + * with the flow law components. + * The viscosity is equal to the product of the + * minimum and maximum viscosities + * divided by the difference between the maximum and + * minimum viscosities. + */ + double damper_viscosity; + + /** + * A factor used to scale the viscoplastic strain rate + * up to the total strain rate. + * The scaling factor is equal to the maximum viscosity + * divided by the difference between the maximum and + * minimum viscosities. + */ + double strain_rate_scaling_factor; + + /** + * The minimum strain rate allowed by the rheology. + */ + double minimum_strain_rate; + + /** + * The log strain rate threshold which must be passed + * before successful termination of the Newton iteration + * to determine the creep stress. + */ + double log_strain_rate_residual_threshold; + + /** + * The maximum number of iterations allowed to determine + * the creep stress. + */ unsigned int stress_max_iteration_number; + + /** + * Useful number to aid in adding together exponentials + */ + const double logmin = std::log(std::numeric_limits::min()); }; } } diff --git a/include/aspect/material_model/rheology/peierls_creep.h b/include/aspect/material_model/rheology/peierls_creep.h index a63f863ebaf..762a4bf9fcb 100644 --- a/include/aspect/material_model/rheology/peierls_creep.h +++ b/include/aspect/material_model/rheology/peierls_creep.h @@ -179,6 +179,17 @@ namespace aspect const double temperature, const PeierlsCreepParameters creep_parameters) const; + /** + * Compute the natural logarithm of the strain rate norm and its first + * derivative with respect to the natural logarithm of the stress norm + * based on the approximate Peierls creep law. + */ + std::pair + compute_approximate_log_strain_rate_and_derivative (const double log_stress, + const double pressure, + const double temperature, + const PeierlsCreepParameters creep_parameters) const; + /** * Compute the strain rate and first stress derivative * as a function of stress based on the selected Peierls creep law. diff --git a/source/material_model/rheology/composite_visco_plastic.cc b/source/material_model/rheology/composite_visco_plastic.cc index a30a875fc02..d8537ca5aa8 100644 --- a/source/material_model/rheology/composite_visco_plastic.cc +++ b/source/material_model/rheology/composite_visco_plastic.cc @@ -84,7 +84,7 @@ namespace aspect // Isostrain averaging double total_volume_fraction = 1.; - for (unsigned int composition=0; composition < number_of_compositions; ++composition) + for (unsigned int composition=0; composition < number_of_chemical_compositions; ++composition) { // Only include the contribution to the viscosity // from a given composition if the volume fraction exceeds @@ -134,17 +134,20 @@ namespace aspect // 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.)), - min_strain_rate); + minimum_strain_rate); + const double log_edot_ii = std::log(edot_ii); Rheology::DiffusionCreepParameters diffusion_creep_parameters; Rheology::DislocationCreepParameters dislocation_creep_parameters; Rheology::PeierlsCreepParameters peierls_creep_parameters; Rheology::DruckerPragerParameters drucker_prager_parameters; - // 1) Estimate the stress running through the viscoplastic elements and the maximum viscosity element. - // These are all arranged in series. Taking the minimum viscosity assuming that all the strain - // runs through a single element provides an excellent first approximation to the true viscosity. - // The stress can then be calculated as 2 * eta * edot_ii + // 1) Estimate the stress running through the creep elements and the + // maximum viscosity element, whose viscosity is defined in parse_parameters. + // These are all arranged in series. Taking the minimum viscosity from + // all of these elements provides an excellent first approximation + // to the true viscosity. + // The stress can then be calculated as 2 * viscoplastic_viscosity_guess * edot_ii double viscoplastic_viscosity_guess = maximum_viscosity; if (use_diffusion_creep) @@ -172,40 +175,64 @@ namespace aspect drucker_prager_parameters.angle_internal_friction, pressure, edot_ii, drucker_prager_parameters.max_yield_stress), viscoplastic_viscosity_guess); } - double viscoplastic_stress = 2.*viscoplastic_viscosity_guess*edot_ii; + double viscoplastic_stress = 2. * viscoplastic_viscosity_guess * edot_ii; + double log_viscoplastic_stress = std::log(viscoplastic_stress); - // In this rheology model, the total strain rate is partitioned between - // different flow laws. We do not know how the strain is partitioned - // between these flow laws, nor do we know the viscoplastic stress, which is - // required to calculate the partitioning. + // 2) Calculate the log strain rate and first derivative with respect to + // log stress for each element using the guessed creep stress. + std::array, 4> log_edot_and_deriv; - // The following while loop conducts a Newton iteration to obtain the - // viscoplastic stress, which we need in order to calculate the viscosity. - double strain_rate_residual = 2*strain_rate_residual_threshold; - double strain_rate_deriv = 0; + if (use_diffusion_creep) + log_edot_and_deriv[0] = diffusion_creep->compute_log_strain_rate_and_derivative(log_viscoplastic_stress, pressure, temperature, grain_size, diffusion_creep_parameters); + + if (use_dislocation_creep) + log_edot_and_deriv[1] = dislocation_creep->compute_log_strain_rate_and_derivative (log_viscoplastic_stress, pressure, temperature, dislocation_creep_parameters); + + if (use_peierls_creep) + log_edot_and_deriv[2] = peierls_creep->compute_approximate_log_strain_rate_and_derivative(log_viscoplastic_stress, pressure, temperature, peierls_creep_parameters); + + if (use_drucker_prager) + log_edot_and_deriv[3] = drucker_prager->compute_log_strain_rate_and_derivative (log_viscoplastic_stress, pressure, drucker_prager_parameters); + + // 3) Calculate the total log strain rate from the first estimates + // for the component log strain rates. + // This will generally be more than the input total strain rate because + // the creep stress in Step 1 was been calculated assuming that + // only one mechanism was active, whereas the strain rate + // calculated in Step 2 allowed all the mechanisms to + // accommodate strain at that creep stress. + std::pair log_edot_ii_and_deriv_iterate = calculate_log_strain_rate_and_derivative(log_edot_and_deriv, + viscoplastic_stress, + partial_strain_rates); + double log_strain_rate_residual = log_edot_ii_and_deriv_iterate.first - log_edot_ii; + + // 4) In this rheology model, the total strain rate is partitioned between + // different flow components. We do not know how the strain is partitioned + // between these components. + + // The following while loop contains a Newton iteration to obtain the + // viscoplastic stress that is consistent with the total strain rate. unsigned int stress_iteration = 0; - while (std::abs(strain_rate_residual) > strain_rate_residual_threshold + while (std::abs(log_strain_rate_residual) > log_strain_rate_residual_threshold && stress_iteration < stress_max_iteration_number) { + // Apply the Newton update for the log creep stress using the + // strain-rate residual and strain-rate stress derivative + double delta_log_viscoplastic_stress = log_strain_rate_residual/log_edot_ii_and_deriv_iterate.second; + log_viscoplastic_stress += delta_log_viscoplastic_stress; + viscoplastic_stress = std::exp(log_viscoplastic_stress); - const std::pair viscoplastic_edot_and_deriv = compute_strain_rate_and_derivative (viscoplastic_stress, - pressure, - temperature, - grain_size, - diffusion_creep_parameters, - dislocation_creep_parameters, - peierls_creep_parameters, - drucker_prager_parameters); - - const double strain_rate = viscoplastic_stress/(2.*maximum_viscosity) + (maximum_viscosity/(maximum_viscosity - minimum_viscosity))*viscoplastic_edot_and_deriv.first; - strain_rate_deriv = 1./(2.*maximum_viscosity) + (maximum_viscosity/(maximum_viscosity - minimum_viscosity))*viscoplastic_edot_and_deriv.second; + // Update the strain rates of all mechanisms with the new stress + for (auto &i : active_flow_mechanisms) + log_edot_and_deriv[i].first += log_edot_and_deriv[i].second * delta_log_viscoplastic_stress; - strain_rate_residual = strain_rate - edot_ii; + // Compute the new log strain rate residual and log stress derivative + log_edot_ii_and_deriv_iterate = calculate_log_strain_rate_and_derivative(log_edot_and_deriv, + viscoplastic_stress, + partial_strain_rates); + log_strain_rate_residual = log_edot_ii - log_edot_ii_and_deriv_iterate.first; - // If the strain rate derivative is zero, we catch it below. - if (strain_rate_deriv>std::numeric_limits::min()) - viscoplastic_stress -= strain_rate_residual/strain_rate_deriv; - stress_iteration += 1; + ++stress_iteration; // If anything that would be used in the next iteration is not finite, the // Newton iteration would trigger an exception and we want to abort the @@ -213,58 +240,28 @@ namespace aspect // Currently, we still throw an exception, but if this exception is thrown, // another more robust iterative scheme should be implemented // (similar to that seen in the diffusion-dislocation material model). - const bool abort_newton_iteration = !numbers::is_finite(viscoplastic_stress) - || !numbers::is_finite(strain_rate_residual) - || !numbers::is_finite(strain_rate_deriv) - || strain_rate_deriv < std::numeric_limits::min() + const bool abort_newton_iteration = !numbers::is_finite(log_viscoplastic_stress) + || !numbers::is_finite(log_strain_rate_residual) || stress_iteration == stress_max_iteration_number; AssertThrow(!abort_newton_iteration, ExcMessage("No convergence has been reached in the loop that determines " - "the composite viscoplastic stress. Aborting! " - "Residual is " + Utilities::to_string(strain_rate_residual) + + "the composite viscous creep stress. Aborting! " + "Residual is " + Utilities::to_string(log_strain_rate_residual) + " after " + Utilities::to_string(stress_iteration) + " iterations. " "You can increase the number of iterations by adapting the " "parameter 'Maximum creep strain rate iterations'.")); } - // The viscoplastic stress is not the total stress, so we still need to do a little work to obtain the effective viscosity. - // First, we compute the stress running through the strain rate limiter, and then add that to the viscoplastic stress - // NOTE: The viscosity of the strain rate limiter is equal to (minimum_viscosity*maximum_viscosity)/(maximum_viscosity - minimum_viscosity) - const double lim_stress = 2.*minimum_viscosity*(edot_ii - viscoplastic_stress/(2.*maximum_viscosity)); - const double total_stress = viscoplastic_stress + lim_stress; - - // Compute the strain rate experienced by the different mechanisms - // These should sum to the total strain rate - - // The components of partial_strain_rates must be provided in the order - // dictated by make_strain_rate_additional_outputs_names - if (use_diffusion_creep) - { - const std::pair diff_edot_and_deriv = diffusion_creep->compute_strain_rate_and_derivative(viscoplastic_stress, pressure, temperature, grain_size, diffusion_creep_parameters); - partial_strain_rates[0] = diff_edot_and_deriv.first; - } - - if (use_dislocation_creep) - { - const std::pair disl_edot_and_deriv = dislocation_creep->compute_strain_rate_and_derivative(viscoplastic_stress, pressure, temperature, dislocation_creep_parameters); - partial_strain_rates[1] = disl_edot_and_deriv.first; - } - - if (use_peierls_creep) - { - const std::pair prls_edot_and_deriv = peierls_creep->compute_strain_rate_and_derivative(viscoplastic_stress, pressure, temperature, peierls_creep_parameters); - partial_strain_rates[2] = prls_edot_and_deriv.first; - } - - if (use_drucker_prager) - { - const std::pair drpr_edot_and_deriv = drucker_prager->compute_strain_rate_and_derivative(viscoplastic_stress, pressure, drucker_prager_parameters); - partial_strain_rates[3] = drpr_edot_and_deriv.first; - } - - partial_strain_rates[4] = total_stress/(2.*maximum_viscosity); + // 5) We have now calculated the viscoplastic stress consistent with the total + // strain rate, but the viscoplastic stress is only one component of the total stress, + // 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; + const double total_stress = viscoplastic_stress + damper_stress; - // Now we return the viscosity using the total stress + // 6) Return the effective creep viscosity using the total stress return total_stress/(2.*edot_ii); } @@ -280,30 +277,49 @@ namespace aspect template std::pair - CompositeViscoPlastic::compute_strain_rate_and_derivative (const double viscoplastic_stress, - const double pressure, - const double temperature, - const double grain_size, - const DiffusionCreepParameters diffusion_creep_parameters, - const DislocationCreepParameters dislocation_creep_parameters, - const PeierlsCreepParameters peierls_creep_parameters, - const DruckerPragerParameters drucker_prager_parameters) const + CompositeViscoPlastic::calculate_log_strain_rate_and_derivative(const std::array, 4> &logarithmic_strain_rates_and_stress_derivatives, + const double viscoplastic_stress, + std::vector &partial_strain_rates) const { - std::pair viscoplastic_edot_and_deriv = std::make_pair(0., 0.); + // The total strain rate + double viscoplastic_strain_rate_sum = 0.0; - if (use_diffusion_creep) - viscoplastic_edot_and_deriv = viscoplastic_edot_and_deriv + diffusion_creep->compute_strain_rate_and_derivative(viscoplastic_stress, pressure, temperature, grain_size, diffusion_creep_parameters); + // The sum of the stress derivatives multiplied by the mechanism strain rates + double weighted_stress_derivative_sum = 0.0; - if (use_dislocation_creep) - viscoplastic_edot_and_deriv = viscoplastic_edot_and_deriv + dislocation_creep->compute_strain_rate_and_derivative(viscoplastic_stress, pressure, temperature, dislocation_creep_parameters); + // The first derivative of log(strain rate) with respect to log(stress) + // is computed as sum_i(stress_exponent_i * edot_i) / sum_i(edot_i) + // i.e., the stress exponents weighted by strain rate fraction + // summed over the individual flow mechanisms (i). - if (use_peierls_creep) - viscoplastic_edot_and_deriv = viscoplastic_edot_and_deriv + peierls_creep->compute_strain_rate_and_derivative(viscoplastic_stress, pressure, temperature, peierls_creep_parameters); + // Loop over active flow laws and add their contributions + // to the strain rate and stress derivative + for (auto &i : active_flow_mechanisms) + { + double mechanism_log_strain_rate = logarithmic_strain_rates_and_stress_derivatives[i].first; - if (use_drucker_prager) - viscoplastic_edot_and_deriv = viscoplastic_edot_and_deriv + drucker_prager->compute_strain_rate_and_derivative(viscoplastic_stress, pressure, drucker_prager_parameters); + // Check if the mechanism strain rate is within bounds to prevent underflow + if (mechanism_log_strain_rate >= logmin) + { + const double mechanism_strain_rate = std::exp(mechanism_log_strain_rate); + partial_strain_rates[i] = mechanism_strain_rate; + const double log_stress_derivative = logarithmic_strain_rates_and_stress_derivatives[i].second; + viscoplastic_strain_rate_sum += mechanism_strain_rate; + weighted_stress_derivative_sum += log_stress_derivative * mechanism_strain_rate; + } + } + + const double log_viscoplastic_strain_rate_derivative = weighted_stress_derivative_sum / viscoplastic_strain_rate_sum; - return viscoplastic_edot_and_deriv; + // Some opaque mathmatics 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; + // 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; + + return std::make_pair(std::log(strain_rate), log_strain_rate_derivative); } @@ -351,8 +367,8 @@ namespace aspect "Stabilizes strain dependent viscosity. Units: \\si{\\per\\second}."); // Viscosity iteration parameters - prm.declare_entry ("Strain rate residual tolerance", "1e-22", Patterns::Double(0.), - "Tolerance for correct diffusion/dislocation strain rate ratio."); + prm.declare_entry ("Strain rate residual tolerance", "1e-10", Patterns::Double(0.), + "Tolerance for correct log strain rate residual."); prm.declare_entry ("Maximum creep strain rate iterations", "40", Patterns::Integer(0), "Maximum number of iterations to find the correct " "viscoplastic strain rate."); @@ -374,31 +390,48 @@ namespace aspect CompositeViscoPlastic::parse_parameters (ParameterHandler &prm, const std::unique_ptr> &expected_n_phases_per_composition) { - // Retrieve the list of composition names - const std::vector list_of_composition_names = this->introspection().get_composition_names(); - // A background field is required by the subordinate material models - number_of_compositions = list_of_composition_names.size() + 1; + number_of_chemical_compositions = this->introspection().n_chemical_composition_fields() + 1; - min_strain_rate = prm.get_double("Minimum strain rate"); + minimum_strain_rate = prm.get_double("Minimum strain rate"); // Iteration parameters - strain_rate_residual_threshold = prm.get_double ("Strain rate residual tolerance"); + log_strain_rate_residual_threshold = prm.get_double ("Strain rate residual tolerance"); stress_max_iteration_number = prm.get_integer ("Maximum creep strain rate iterations"); - // Read min and max viscosity parameters - minimum_viscosity = prm.get_double ("Minimum viscosity"); + // Read maximum viscosity parameter maximum_viscosity = prm.get_double ("Maximum viscosity"); - // Rheological parameters + // Process minimum viscosity parameter + // In this rheology model, there are two viscous dampers designed + // to stabilise the solution, a "limiting damper" arranged in parallel + // with the flow law components which stops the viscosity going to zero, + // and a "maximum viscosity damper", which is placed in series with the + // combined flow law components and the limiting damper. + // - If the real creep viscosity is infinite, the total strain rate will + // run through only the "maximum viscosity" damper. + // - If the real creep viscosity is equal to zero, the total strain rate + // will run through the "limiting damper" and the "maximum viscosity" damper, + // with a total viscosity equal to the harmonic sum of these dampers. + // Therefore, the "limiting damper" has a viscosity equal to + // eta_max * eta_min / (eta_max - eta_min). + // 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"); + strain_rate_scaling_factor = maximum_viscosity / (maximum_viscosity - minimum_viscosity); + damper_viscosity = maximum_viscosity * minimum_viscosity / (maximum_viscosity - minimum_viscosity); + // Rheological parameters // Diffusion creep parameters - use_diffusion_creep = prm.get_bool ("Include diffusion creep in composite rheology"); + use_diffusion_creep = prm.get_bool("Include diffusion creep in composite rheology"); if (use_diffusion_creep) { diffusion_creep = std::make_unique>(); diffusion_creep->initialize_simulator (this->get_simulator()); diffusion_creep->parse_parameters(prm, expected_n_phases_per_composition); + + active_flow_mechanisms.push_back(0); } // Dislocation creep parameters @@ -408,6 +441,8 @@ namespace aspect dislocation_creep = std::make_unique>(); dislocation_creep->initialize_simulator (this->get_simulator()); dislocation_creep->parse_parameters(prm, expected_n_phases_per_composition); + + active_flow_mechanisms.push_back(1); } // Peierls creep parameters @@ -417,9 +452,11 @@ namespace aspect peierls_creep = std::make_unique>(); peierls_creep->initialize_simulator (this->get_simulator()); peierls_creep->parse_parameters(prm, expected_n_phases_per_composition); + active_flow_mechanisms.push_back(2); AssertThrow((prm.get ("Peierls creep flow law") == "viscosity approximation"), ExcMessage("The Peierls creep flow law parameter needs to be set to viscosity approximation.")); + } // Drucker Prager parameters @@ -429,9 +466,10 @@ namespace aspect drucker_prager = std::make_unique>(); drucker_prager->initialize_simulator (this->get_simulator()); drucker_prager->parse_parameters(prm, expected_n_phases_per_composition); + active_flow_mechanisms.push_back(3); } - AssertThrow(use_diffusion_creep == true || use_dislocation_creep == true || use_peierls_creep == true || use_drucker_prager == true, + AssertThrow(active_flow_mechanisms.size() > 0, ExcMessage("You need to include at least one deformation mechanism.")); } diff --git a/source/material_model/rheology/peierls_creep.cc b/source/material_model/rheology/peierls_creep.cc index 8b46fd86b72..0aaaf358266 100644 --- a/source/material_model/rheology/peierls_creep.cc +++ b/source/material_model/rheology/peierls_creep.cc @@ -334,6 +334,39 @@ namespace aspect + template + std::pair + PeierlsCreep::compute_approximate_log_strain_rate_and_derivative (const double log_stress, + const double pressure, + const double temperature, + const PeierlsCreepParameters creep_parameters) const + { + /** + * b = (E+P*V)/(R*T) + * c = std::pow(gamma, p) + * d = std::pow(1. - c, q) + * s = b*p*q*c*d/(1. - c) + * log_arrhenius = -b*d + * + * log_strain_rate = log_A + (s + n) * std::log(stress) - s * std::log(gamma*peierls_stress) + log_arrhenius + * total_stress_exponent = (s + n) + */ + const PeierlsCreepParameters p = creep_parameters; + + const double b = (p.activation_energy + pressure*p.activation_volume)/(constants::gas_constant * temperature); + const double c = std::pow(p.fitting_parameter, p.glide_parameter_p); + const double d = std::pow(1. - c, p.glide_parameter_q); + const double s = b*p.glide_parameter_p*p.glide_parameter_q*c*d/(1. - c); + const double log_arrhenius = -b*d; + + const double total_stress_exponent = s + p.stress_exponent; + const double log_strain_rate = std::log(p.prefactor) + total_stress_exponent * log_stress - s * std::log(p.fitting_parameter*p.peierls_stress) + log_arrhenius; + + return std::make_pair(log_strain_rate, total_stress_exponent); + } + + + template std::pair PeierlsCreep::compute_exact_strain_rate_and_derivative (const double stress, diff --git a/tests/composite_viscous_outputs.cc b/tests/composite_viscous_outputs.cc index 7a5729a1f7e..866c7787bdf 100644 --- a/tests/composite_viscous_outputs.cc +++ b/tests/composite_viscous_outputs.cc @@ -47,6 +47,10 @@ void f(const aspect::SimulatorAccess &simulator_access, composite_creep = std::make_unique>(); composite_creep->initialize_simulator (simulator_access.get_simulator()); composite_creep->declare_parameters(prm); + 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("Peierls creep flow law", "viscosity approximation"); prm.set("Maximum yield stress", "5e8"); composite_creep->parse_parameters(prm, n_phases); diff --git a/tests/composite_viscous_outputs/screen-output b/tests/composite_viscous_outputs/screen-output index 4e320789606..77b2499d232 100644 --- a/tests/composite_viscous_outputs/screen-output +++ b/tests/composite_viscous_outputs/screen-output @@ -7,13 +7,13 @@ temperature (K) eta (Pas) creep stress (Pa) edot_ii (/s) edot_ii fractio 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.18168e-11 -1500 4.32083e+17 6.64167e+06 1e-11 0.0538305 0.94617 2.09885e-13 1.4634e-101 4.32083e-11 +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.67379e-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.02964e-11 +2000 1.02964e+17 59280.1 1e-11 0.99654 0.00345962 2.26225e-21 4.97674e-204 1.02965e-11 OK Number of active cells: 100 (on 1 levels) Number of degrees of freedom: 5,534 (3,969+242+1,323) diff --git a/tests/composite_viscous_outputs_no_peierls.cc b/tests/composite_viscous_outputs_no_peierls.cc new file mode 100644 index 00000000000..4eb991845ac --- /dev/null +++ b/tests/composite_viscous_outputs_no_peierls.cc @@ -0,0 +1,199 @@ +/* + 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 + +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 = {1.}; + 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("Include diffusion creep in composite rheology", "true"); + prm.set("Include dislocation creep in composite rheology", "true"); + prm.set("Include Peierls creep in composite rheology", "false"); + prm.set("Include Drucker Prager plasticity in composite rheology", "true"); + prm.set("Peierls creep flow law", "viscosity approximation"); + prm.set("Maximum yield stress", "5e8"); + composite_creep->parse_parameters(prm, n_phases); + + 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, n_phases); + + 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, n_phases); + + 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, n_phases); + 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; + SymmetricTensor<2,dim> strain_rate; + strain_rate[0][0] = -1e-11; + strain_rate[0][1] = 0.; + strain_rate[1][1] = 1e-11; + strain_rate[2][0] = 0.; + 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; + + // 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 drpr_stress; + std::vector partial_strain_rates(5, 0.); + + for (unsigned int i=0; i <= 10; i++) + { + temperature = 1000. + i*100.; + + // Compute the viscosity + viscosity = composite_creep->compute_composition_viscosity(pressure, temperature, grain_size, composition, strain_rate, 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_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); + 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((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 << " 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)); +} + +ASPECT_REGISTER_SIGNALS_CONNECTOR(signal_connector<2>, + signal_connector<3>) diff --git a/tests/composite_viscous_outputs_no_peierls.prm b/tests/composite_viscous_outputs_no_peierls.prm new file mode 100644 index 00000000000..dcd9d24dc6d --- /dev/null +++ b/tests/composite_viscous_outputs_no_peierls.prm @@ -0,0 +1,80 @@ +# This test checks whether the composite viscous rheology +# plugin produces the correct composite viscosity and +# decomposed strain rates. + +set Additional shared libraries = tests/libcomposite_viscous_outputs_no_peierls.so +set Dimension = 3 +set End time = 0 +set Use years in output instead of seconds = true +set Nonlinear solver scheme = no Advection, no Stokes + +# 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. + 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_no_peierls/screen-output b/tests/composite_viscous_outputs_no_peierls/screen-output new file mode 100644 index 00000000000..e1c4155b475 --- /dev/null +++ b/tests/composite_viscous_outputs_no_peierls/screen-output @@ -0,0 +1,28 @@ + +Loading shared library <./libcomposite_viscous_outputs_no_peierls.debug.so> + +* Connecting signals +temperature (K) eta (Pas) creep stress (Pa) edot_ii (/s) edot_ii fractions (diff, disl, prls, drpr, max) +1000 3.46071e+19 6.90142e+08 1e-11 1.30023e-06 0.00365246 0 0.996346 3.46071e-09 +1100 3.13843e+19 6.25686e+08 1e-11 7.59662e-05 0.992522 0 0.00740179 3.13843e-09 +1200 7.7057e+18 1.52114e+08 1e-11 0.000594408 0.999406 0 1.44528e-33 7.7057e-10 +1300 2.39282e+18 4.58564e+07 1e-11 0.00338081 0.996619 0 1.32277e-59 2.39282e-10 +1400 9.18168e+17 1.63634e+07 1e-11 0.0149609 0.985039 0 5.56077e-82 9.18169e-11 +1500 4.32083e+17 6.64167e+06 1e-11 0.0538305 0.94617 0 1.4634e-101 4.32084e-11 +1600 2.47246e+17 2.94493e+06 1e-11 0.161077 0.838923 0 3.2006e-119 2.47246e-11 +1700 1.67379e+17 1.34758e+06 1e-11 0.397345 0.602655 0 3.38184e-136 1.67378e-11 +1800 1.28447e+17 568948 1e-11 0.749975 0.250025 0 6.38432e-155 1.28447e-11 +1900 1.09563e+17 191256 1e-11 0.962698 0.0373021 0 1.35594e-178 1.09563e-11 +2000 1.02964e+17 59280.1 1e-11 0.99654 0.00345962 0 4.97674e-204 1.02965e-11 +OK +Number of active cells: 100 (on 1 levels) +Number of degrees of freedom: 5,534 (3,969+242+1,323) + +*** Timestep 0: t=0 years, dt=0 years + + Postprocessing: + +Termination requested by criterion: end time + + +