From b1e55a069c0234abee7673aa7515e0a254a99689 Mon Sep 17 00:00:00 2001 From: Bob Myhill Date: Thu, 25 Jul 2024 16:58:15 +0100 Subject: [PATCH 1/8] added isostress averaging --- .../rheology/composite_visco_plastic.h | 37 ++ .../rheology/composite_visco_plastic.cc | 315 ++++++++++++++++-- 2 files changed, 333 insertions(+), 19 deletions(-) diff --git a/include/aspect/material_model/rheology/composite_visco_plastic.h b/include/aspect/material_model/rheology/composite_visco_plastic.h index 02b6e107b92..6ee3a6c764d 100644 --- a/include/aspect/material_model/rheology/composite_visco_plastic.h +++ b/include/aspect/material_model/rheology/composite_visco_plastic.h @@ -83,6 +83,25 @@ namespace aspect const std::vector &phase_function_values = std::vector(), const std::vector &n_phase_transitions_per_composition = std::vector()) const; + /** + * Compute the isostress viscosity over all compositional fields + * based on the composite viscous creep law. + * If @p expected_n_phases_per_composition points to a vector of + * unsigned integers this is considered the number of phase transitions + * for each compositional field and viscosity will be first computed on + * each phase and then averaged for each compositional field. + */ + double + compute_isostress_viscosity (const double pressure, + const double temperature, + const double grain_size, + const std::vector &volume_fractions, + const SymmetricTensor<2,dim> &strain_rate, + std::vector &partial_strain_rates, + const std::vector &phase_function_values = std::vector(), + const std::vector &n_phase_transitions_per_composition = std::vector()) const; + + /** * Compute the compositional field viscosity * based on the composite viscous creep law. @@ -111,7 +130,25 @@ namespace aspect const double viscoplastic_stress, std::vector &partial_strain_rates) const; + /** + * 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 + calculate_log_strain_rate_and_derivative(const std::vector, 4>> &logarithmic_strain_rates_and_stress_derivatives, + const double viscoplastic_stress, + std::vector &partial_strain_rates) const; + private: + /** + * Enumeration for selecting which type of viscosity averaging to use. + */ + enum ViscosityAveragingScheme + { + isostrain, + isostress + } viscosity_averaging_scheme; /** * Whether to use different deformation mechanisms diff --git a/source/material_model/rheology/composite_visco_plastic.cc b/source/material_model/rheology/composite_visco_plastic.cc index d8537ca5aa8..5ee4f4e2acf 100644 --- a/source/material_model/rheology/composite_visco_plastic.cc +++ b/source/material_model/rheology/composite_visco_plastic.cc @@ -82,38 +82,243 @@ namespace aspect double viscosity = 0.; partial_strain_rates.resize(5, 0.); + // Isostress averaging + if (viscosity_averaging_scheme == isostress) + { + viscosity += compute_isostress_viscosity (pressure, + temperature, + grain_size, + volume_fractions, + strain_rate, + partial_strain_rates, + phase_function_values, + n_phase_transitions_per_composition); + } + // Isostrain averaging + if (viscosity_averaging_scheme == isostrain) + { + double total_volume_fraction = 1.; + 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 + // a certain (small) fraction. + if (volume_fractions[composition] > 2.*std::numeric_limits::epsilon()) + { + std::vector partial_strain_rates_composition(5, 0.); + viscosity += (volume_fractions[composition] + * compute_composition_viscosity (pressure, + temperature, + grain_size, + composition, + strain_rate, + partial_strain_rates_composition, + phase_function_values, + n_phase_transitions_per_composition)); + for (unsigned int j=0; j < 5; ++j) + partial_strain_rates[j] += volume_fractions[composition] * partial_strain_rates_composition[j]; + } + else + { + total_volume_fraction -= volume_fractions[composition]; + } + + viscosity /= total_volume_fraction; + for (unsigned int j=0; j < 5; ++j) + partial_strain_rates[j] /= total_volume_fraction; + } + } + return viscosity; + } + + + + + template + double + CompositeViscoPlastic::compute_isostress_viscosity (const double pressure, + const double temperature, + const double grain_size, + const std::vector &volume_fractions, + const SymmetricTensor<2,dim> &strain_rate, + std::vector &partial_strain_rates, + const std::vector &phase_function_values, + const std::vector &n_phase_transitions_per_composition) const + { + // If strain rate is zero (like during the first time step) set it to some very small number + // 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.)), + minimum_strain_rate); + const double log_edot_ii = std::log(edot_ii); + + std::vector diffusion_creep_parameters; + std::vector dislocation_creep_parameters; + std::vector peierls_creep_parameters; + std::vector drucker_prager_parameters; + + // 1) Estimate the stress running through the creep elements and the + // maximum viscosity element, whose viscosity is defined in parse_parameters. + // The creep elements for each composition are arranged in series. + // Taking the minimum viscosity from + // all of these elements provides an excellent first approximation + // to the true viscosity of that composition. + // For a starting guess, we further assume that one composition takes + // up all the strain. + + // The stress can then be calculated as 2 * viscoplastic_viscosity_guess * edot_ii + double viscoplastic_viscosity_guess = maximum_viscosity; + double total_volume_fraction = 1.; - for (unsigned int composition=0; composition < number_of_chemical_compositions; ++composition) + std::vector active_compositions; + 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 // a certain (small) fraction. - if (volume_fractions[composition] > 2.*std::numeric_limits::epsilon()) + const double volume_fraction = volume_fractions[composition]; + if (volume_fraction > 2. * std::numeric_limits::epsilon()) { - std::vector partial_strain_rates_composition(5, 0.); - viscosity += (volume_fractions[composition] - * compute_composition_viscosity (pressure, - temperature, - grain_size, - composition, - strain_rate, - partial_strain_rates_composition, - phase_function_values, - n_phase_transitions_per_composition)); - for (unsigned int j=0; j < 5; ++j) - partial_strain_rates[j] += volume_fractions[composition] * partial_strain_rates_composition[j]; + active_compositions.push_back(composition); + + if (use_diffusion_creep) + { + diffusion_creep_parameters.push_back(diffusion_creep->compute_creep_parameters(composition, phase_function_values, n_phase_transitions_per_composition)); + viscoplastic_viscosity_guess = std::min(diffusion_creep->compute_viscosity(pressure, temperature, grain_size, composition, phase_function_values, n_phase_transitions_per_composition)/volume_fraction, viscoplastic_viscosity_guess); + } + + if (use_dislocation_creep) + { + dislocation_creep_parameters.push_back(dislocation_creep->compute_creep_parameters(composition, phase_function_values, n_phase_transitions_per_composition)); + viscoplastic_viscosity_guess = std::min(dislocation_creep->compute_viscosity(edot_ii/volume_fraction, pressure, temperature, composition, phase_function_values, n_phase_transitions_per_composition)/volume_fraction, viscoplastic_viscosity_guess); + } + + if (use_peierls_creep) + { + peierls_creep_parameters.push_back(peierls_creep->compute_creep_parameters(composition)); + viscoplastic_viscosity_guess = std::min(peierls_creep->compute_approximate_viscosity(edot_ii/volume_fraction, pressure, temperature, composition)/volume_fraction, viscoplastic_viscosity_guess); + } + + if (use_drucker_prager) + { + Rheology::DruckerPragerParameters drucker_prager_parameters_c = drucker_prager->compute_drucker_prager_parameters(composition, phase_function_values, n_phase_transitions_per_composition); + drucker_prager_parameters.push_back(drucker_prager_parameters_c); + viscoplastic_viscosity_guess = std::min(drucker_prager->compute_viscosity(drucker_prager_parameters_c.cohesion, + drucker_prager_parameters_c.angle_internal_friction, pressure, edot_ii/volume_fraction, drucker_prager_parameters_c.max_yield_stress)/volume_fraction, + viscoplastic_viscosity_guess); + } } else { total_volume_fraction -= volume_fractions[composition]; } + } + + double viscoplastic_stress = 2. * viscoplastic_viscosity_guess * edot_ii / total_volume_fraction; + double log_viscoplastic_stress = std::log(viscoplastic_stress); + + // 2) Calculate the log strain rate and first derivative with respect to + // log stress for each element using the guessed creep stress. + std::vector, 4>> log_edot_and_deriv(active_compositions.size()); - viscosity /= total_volume_fraction; - for (unsigned int j=0; j < 5; ++j) - partial_strain_rates[j] /= total_volume_fraction; + for (unsigned int i = 0; i < active_compositions.size(); ++i) + { + const double log_volume_fraction = std::log(volume_fractions[active_compositions[i]]); + + if (use_diffusion_creep) + { + log_edot_and_deriv[i][0] = diffusion_creep->compute_log_strain_rate_and_derivative(log_viscoplastic_stress, pressure, temperature, grain_size, diffusion_creep_parameters[i]); + log_edot_and_deriv[i][0].first += log_volume_fraction; + } + if (use_dislocation_creep) + { + log_edot_and_deriv[i][1] = dislocation_creep->compute_log_strain_rate_and_derivative(log_viscoplastic_stress, pressure, temperature, dislocation_creep_parameters[i]); + log_edot_and_deriv[i][1].first += log_volume_fraction; + } + if (use_peierls_creep) + { + log_edot_and_deriv[i][2] = peierls_creep->compute_approximate_log_strain_rate_and_derivative(log_viscoplastic_stress, pressure, temperature, peierls_creep_parameters[i]); + log_edot_and_deriv[i][2].first += log_volume_fraction; + } + if (use_drucker_prager) + { + log_edot_and_deriv[i][3] = drucker_prager->compute_log_strain_rate_and_derivative(log_viscoplastic_stress, pressure, drucker_prager_parameters[i]); + log_edot_and_deriv[i][3].first += log_volume_fraction; + } } - return viscosity; + + // 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(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); + + // Update the strain rates of all mechanisms with the new stress + for (auto &log_edot_and_deriv_c : log_edot_and_deriv) + { + for (auto &i : active_flow_mechanisms) + log_edot_and_deriv_c[i].first += log_edot_and_deriv_c[i].second * delta_log_viscoplastic_stress; + } + + // 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; + + ++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 + // iteration instead. + // 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(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 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'.")); + } + + // 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; + + // 6) Return the effective creep viscosity using the total stress + return total_stress / (2. * edot_ii); } @@ -311,7 +516,63 @@ namespace aspect const double log_viscoplastic_strain_rate_derivative = weighted_stress_derivative_sum / viscoplastic_strain_rate_sum; - // Some opaque mathmatics converts the viscoplastic strain rate to the total strain rate. + // 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; + // 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); + } + + + + template + std::pair + CompositeViscoPlastic::calculate_log_strain_rate_and_derivative(const std::vector, 4>> &logarithmic_strain_rates_and_stress_derivatives, + const double viscoplastic_stress, + std::vector &partial_strain_rates) const + { + // The total strain rate + double viscoplastic_strain_rate_sum = 0.0; + + // The sum of the stress derivatives multiplied by the mechanism strain rates + double weighted_stress_derivative_sum = 0.0; + + // 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). + // Loop over active flow laws and add their contributions + // to the strain rate and stress derivative + + // First, make sure that the active partial_strain_rates are reset to zero. + for (auto &i : active_flow_mechanisms) + partial_strain_rates[i] = 0; + + for (auto &logarithmic_strain_rates_and_stress_derivatives_c : logarithmic_strain_rates_and_stress_derivatives) + { + for (auto &i : active_flow_mechanisms) + { + double mechanism_log_strain_rate = logarithmic_strain_rates_and_stress_derivatives_c[i].first; + + // 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_c[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; + + // 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; @@ -328,6 +589,13 @@ namespace aspect void CompositeViscoPlastic::declare_parameters (ParameterHandler &prm) { + prm.declare_entry ("Viscosity averaging scheme", "isostress", + Patterns::Selection("isostress|isostrain|unspecified"), + "Determines the relationship between the conditions experienced by " + "each chemical compositional field in a composite material. " + "Select either isostress (the default option, Reuss averaging) " + "or isostrain (Voigt averaging)."); + prm.declare_entry ("Include diffusion creep in composite rheology", "true", Patterns::Bool (), "Whether to include diffusion creep in the composite rheology formulation."); @@ -390,6 +658,15 @@ namespace aspect CompositeViscoPlastic::parse_parameters (ParameterHandler &prm, const std::unique_ptr> &expected_n_phases_per_composition) { + if (prm.get ("Viscosity averaging scheme") == "isostress") + viscosity_averaging_scheme = isostress; + else if (prm.get ("Viscosity averaging scheme") == "isostrain") + viscosity_averaging_scheme = isostrain; + else + { + AssertThrow(false, ExcMessage("Not a valid viscosity averaging scheme. Choose either isostress or isostrain.")); + } + // A background field is required by the subordinate material models number_of_chemical_compositions = this->introspection().n_chemical_composition_fields() + 1; From 4cc572a417f4e6b00a78e58f38cc1b2120fc19cc Mon Sep 17 00:00:00 2001 From: Bob Myhill Date: Thu, 25 Jul 2024 23:27:24 +0100 Subject: [PATCH 2/8] added test --- tests/composite_viscous_outputs_isostress.cc | 239 ++++++++++++++++++ tests/composite_viscous_outputs_isostress.prm | 80 ++++++ .../screen-output | 28 ++ 3 files changed, 347 insertions(+) create mode 100644 tests/composite_viscous_outputs_isostress.cc create mode 100644 tests/composite_viscous_outputs_isostress.prm create mode 100644 tests/composite_viscous_outputs_isostress/screen-output diff --git a/tests/composite_viscous_outputs_isostress.cc b/tests/composite_viscous_outputs_isostress.cc new file mode 100644 index 00000000000..d7be34de6aa --- /dev/null +++ b/tests/composite_viscous_outputs_isostress.cc @@ -0,0 +1,239 @@ +/* + 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("Peierls creep flow law", "viscosity approximation"); + prm.set("Maximum yield stress", "5e8"); + 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; + 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 prls_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_isostress_viscosity(pressure, temperature, grain_size, volume_fractions, 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); + 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("Compositional fields"); + { + prm.declare_entry("Number of fields","1", Patterns::Integer()); + } + 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_isostress.prm b/tests/composite_viscous_outputs_isostress.prm new file mode 100644 index 00000000000..bc3f9c01c5e --- /dev/null +++ b/tests/composite_viscous_outputs_isostress.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_isostress.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_isostress/screen-output b/tests/composite_viscous_outputs_isostress/screen-output new file mode 100644 index 00000000000..545da023131 --- /dev/null +++ b/tests/composite_viscous_outputs_isostress/screen-output @@ -0,0 +1,28 @@ + +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 +OK +Number of active cells: 100 (on 1 levels) +Number of degrees of freedom: 6,857 (3,969+242+1,323+1,323) + +*** Timestep 0: t=0 years, dt=0 years + + Postprocessing: + +Termination requested by criterion: end time + + + From 82c760785fe03537464f04b2efb60d9e4656bd7f Mon Sep 17 00:00:00 2001 From: Bob Myhill Date: Thu, 25 Jul 2024 23:31:06 +0100 Subject: [PATCH 3/8] added changelog --- doc/modules/changes/20240725_myhill | 7 +++++++ 1 file changed, 7 insertions(+) create mode 100644 doc/modules/changes/20240725_myhill diff --git a/doc/modules/changes/20240725_myhill b/doc/modules/changes/20240725_myhill new file mode 100644 index 00000000000..6f7adba213c --- /dev/null +++ b/doc/modules/changes/20240725_myhill @@ -0,0 +1,7 @@ +New: The CompositeViscoPlastic rheology in ASPECT now +allows the user to select isostress (Reuss) or +isostrain (Voigt) averaging of viscosities for +multicomponent materials via the parameter +"Viscosity averaging scheme". +
+(Bob Myhill, 2024/07/25) From f41a86eb2a787190c630c9b257d5fc1d66690bf2 Mon Sep 17 00:00:00 2001 From: Bob Myhill Date: Thu, 25 Jul 2024 23:49:36 +0100 Subject: [PATCH 4/8] set averaging schemes in tests --- tests/composite_viscous_outputs.cc | 1 + tests/composite_viscous_outputs_isostress.cc | 1 - tests/composite_viscous_outputs_no_peierls.cc | 1 + 3 files changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/composite_viscous_outputs.cc b/tests/composite_viscous_outputs.cc index 866c7787bdf..0cec200924d 100644 --- a/tests/composite_viscous_outputs.cc +++ b/tests/composite_viscous_outputs.cc @@ -47,6 +47,7 @@ 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("Viscosity averaging scheme", "isostrain"); 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"); diff --git a/tests/composite_viscous_outputs_isostress.cc b/tests/composite_viscous_outputs_isostress.cc index d7be34de6aa..a1dd2ad1697 100644 --- a/tests/composite_viscous_outputs_isostress.cc +++ b/tests/composite_viscous_outputs_isostress.cc @@ -50,7 +50,6 @@ 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("Viscosity averaging scheme", "isostress"); prm.set("Include diffusion creep in composite rheology", "true"); prm.set("Include dislocation creep in composite rheology", "true"); diff --git a/tests/composite_viscous_outputs_no_peierls.cc b/tests/composite_viscous_outputs_no_peierls.cc index 4eb991845ac..6d3078eabb0 100644 --- a/tests/composite_viscous_outputs_no_peierls.cc +++ b/tests/composite_viscous_outputs_no_peierls.cc @@ -47,6 +47,7 @@ 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("Viscosity averaging scheme", "isostrain"); 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"); From 791ea5535c0fbe0ff5eadbfbf525178fa0dcc16b Mon Sep 17 00:00:00 2001 From: Bob Myhill Date: Fri, 26 Jul 2024 12:12:04 +0100 Subject: [PATCH 5/8] sanitize function order --- .../rheology/composite_visco_plastic.h | 25 ++-- .../rheology/composite_visco_plastic.cc | 140 +++++++++--------- 2 files changed, 82 insertions(+), 83 deletions(-) diff --git a/include/aspect/material_model/rheology/composite_visco_plastic.h b/include/aspect/material_model/rheology/composite_visco_plastic.h index 6ee3a6c764d..4ffd7a5a97e 100644 --- a/include/aspect/material_model/rheology/composite_visco_plastic.h +++ b/include/aspect/material_model/rheology/composite_visco_plastic.h @@ -101,6 +101,15 @@ namespace aspect const std::vector &phase_function_values = std::vector(), const std::vector &n_phase_transitions_per_composition = std::vector()) const; + /** + * 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 + * over all compositional fields. Also updates the partial_strain_rates vector + */ + std::pair + calculate_isostress_log_strain_rate_and_derivative(const std::vector, 4>> &logarithmic_strain_rates_and_stress_derivatives, + const double viscoplastic_stress, + std::vector &partial_strain_rates) const; /** * Compute the compositional field viscosity @@ -126,19 +135,9 @@ namespace aspect * Also updates the partial_strain_rates vector */ std::pair - 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; - - /** - * 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 - calculate_log_strain_rate_and_derivative(const std::vector, 4>> &logarithmic_strain_rates_and_stress_derivatives, - const double viscoplastic_stress, - std::vector &partial_strain_rates) const; + calculate_composition_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: /** diff --git a/source/material_model/rheology/composite_visco_plastic.cc b/source/material_model/rheology/composite_visco_plastic.cc index 5ee4f4e2acf..6aa77928c7d 100644 --- a/source/material_model/rheology/composite_visco_plastic.cc +++ b/source/material_model/rheology/composite_visco_plastic.cc @@ -256,7 +256,7 @@ namespace aspect // 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, + std::pair log_edot_ii_and_deriv_iterate = calculate_isostress_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; @@ -284,9 +284,9 @@ namespace aspect } // 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_edot_ii_and_deriv_iterate = calculate_isostress_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; ++stress_iteration; @@ -323,6 +323,62 @@ namespace aspect + template + std::pair + CompositeViscoPlastic::calculate_isostress_log_strain_rate_and_derivative(const std::vector, 4>> &logarithmic_strain_rates_and_stress_derivatives, + const double viscoplastic_stress, + std::vector &partial_strain_rates) const + { + // The total strain rate + double viscoplastic_strain_rate_sum = 0.0; + + // The sum of the stress derivatives multiplied by the mechanism strain rates + double weighted_stress_derivative_sum = 0.0; + + // 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). + // Loop over active flow laws and add their contributions + // to the strain rate and stress derivative + + // First, make sure that the active partial_strain_rates are reset to zero. + for (auto &i : active_flow_mechanisms) + partial_strain_rates[i] = 0; + + for (auto &logarithmic_strain_rates_and_stress_derivatives_c : logarithmic_strain_rates_and_stress_derivatives) + { + for (auto &i : active_flow_mechanisms) + { + double mechanism_log_strain_rate = logarithmic_strain_rates_and_stress_derivatives_c[i].first; + + // 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_c[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; + + // 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; + // 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); + } + + + template double CompositeViscoPlastic::compute_composition_viscosity (const double pressure, @@ -406,7 +462,7 @@ namespace aspect // 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, + std::pair log_edot_ii_and_deriv_iterate = calculate_composition_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; @@ -432,9 +488,9 @@ namespace aspect log_edot_and_deriv[i].first += log_edot_and_deriv[i].second * delta_log_viscoplastic_stress; // 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_edot_ii_and_deriv_iterate = calculate_composition_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; ++stress_iteration; @@ -472,19 +528,11 @@ namespace aspect - // Overload the + operator to act on two pairs of doubles. - std::pair operator+(const std::pair &x, const std::pair &y) - { - return std::make_pair(x.first+y.first, x.second+y.second); - } - - - template std::pair - 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 + CompositeViscoPlastic::calculate_composition_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 { // The total strain rate double viscoplastic_strain_rate_sum = 0.0; @@ -529,58 +577,10 @@ namespace aspect - template - std::pair - CompositeViscoPlastic::calculate_log_strain_rate_and_derivative(const std::vector, 4>> &logarithmic_strain_rates_and_stress_derivatives, - const double viscoplastic_stress, - std::vector &partial_strain_rates) const + // Overload the + operator to act on two pairs of doubles. + std::pair operator+(const std::pair &x, const std::pair &y) { - // The total strain rate - double viscoplastic_strain_rate_sum = 0.0; - - // The sum of the stress derivatives multiplied by the mechanism strain rates - double weighted_stress_derivative_sum = 0.0; - - // 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). - // Loop over active flow laws and add their contributions - // to the strain rate and stress derivative - - // First, make sure that the active partial_strain_rates are reset to zero. - for (auto &i : active_flow_mechanisms) - partial_strain_rates[i] = 0; - - for (auto &logarithmic_strain_rates_and_stress_derivatives_c : logarithmic_strain_rates_and_stress_derivatives) - { - for (auto &i : active_flow_mechanisms) - { - double mechanism_log_strain_rate = logarithmic_strain_rates_and_stress_derivatives_c[i].first; - - // 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_c[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; - - // 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; - // 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); + return std::make_pair(x.first+y.first, x.second+y.second); } From 84ab5dbfcda1f169b7640149c7b840ec1fc4205c Mon Sep 17 00:00:00 2001 From: Bob Myhill Date: Fri, 26 Jul 2024 19:06:30 +0100 Subject: [PATCH 6/8] improve starting guess --- .../rheology/composite_visco_plastic.cc | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/source/material_model/rheology/composite_visco_plastic.cc b/source/material_model/rheology/composite_visco_plastic.cc index 6aa77928c7d..810486999b4 100644 --- a/source/material_model/rheology/composite_visco_plastic.cc +++ b/source/material_model/rheology/composite_visco_plastic.cc @@ -169,7 +169,7 @@ namespace aspect // up all the strain. // The stress can then be calculated as 2 * viscoplastic_viscosity_guess * edot_ii - double viscoplastic_viscosity_guess = maximum_viscosity; + double inverse_viscoplastic_viscosity_guess = 1./maximum_viscosity; double total_volume_fraction = 1.; std::vector active_compositions; @@ -178,6 +178,7 @@ namespace aspect // Only include the contribution to the viscosity // from a given composition if the volume fraction exceeds // a certain (small) fraction. + double viscoplastic_viscosity_guess_c = maximum_viscosity; const double volume_fraction = volume_fractions[composition]; if (volume_fraction > 2. * std::numeric_limits::epsilon()) { @@ -186,29 +187,30 @@ namespace aspect if (use_diffusion_creep) { diffusion_creep_parameters.push_back(diffusion_creep->compute_creep_parameters(composition, phase_function_values, n_phase_transitions_per_composition)); - viscoplastic_viscosity_guess = std::min(diffusion_creep->compute_viscosity(pressure, temperature, grain_size, composition, phase_function_values, n_phase_transitions_per_composition)/volume_fraction, viscoplastic_viscosity_guess); + viscoplastic_viscosity_guess_c = std::min(diffusion_creep->compute_viscosity(pressure, temperature, grain_size, composition, phase_function_values, n_phase_transitions_per_composition), viscoplastic_viscosity_guess_c); } if (use_dislocation_creep) { dislocation_creep_parameters.push_back(dislocation_creep->compute_creep_parameters(composition, phase_function_values, n_phase_transitions_per_composition)); - viscoplastic_viscosity_guess = std::min(dislocation_creep->compute_viscosity(edot_ii/volume_fraction, pressure, temperature, composition, phase_function_values, n_phase_transitions_per_composition)/volume_fraction, viscoplastic_viscosity_guess); + viscoplastic_viscosity_guess_c = std::min(dislocation_creep->compute_viscosity(edot_ii, pressure, temperature, composition, phase_function_values, n_phase_transitions_per_composition), viscoplastic_viscosity_guess_c); } if (use_peierls_creep) { peierls_creep_parameters.push_back(peierls_creep->compute_creep_parameters(composition)); - viscoplastic_viscosity_guess = std::min(peierls_creep->compute_approximate_viscosity(edot_ii/volume_fraction, pressure, temperature, composition)/volume_fraction, viscoplastic_viscosity_guess); + viscoplastic_viscosity_guess_c = std::min(peierls_creep->compute_approximate_viscosity(edot_ii, pressure, temperature, composition), viscoplastic_viscosity_guess_c); } if (use_drucker_prager) { Rheology::DruckerPragerParameters drucker_prager_parameters_c = drucker_prager->compute_drucker_prager_parameters(composition, phase_function_values, n_phase_transitions_per_composition); drucker_prager_parameters.push_back(drucker_prager_parameters_c); - viscoplastic_viscosity_guess = std::min(drucker_prager->compute_viscosity(drucker_prager_parameters_c.cohesion, - drucker_prager_parameters_c.angle_internal_friction, pressure, edot_ii/volume_fraction, drucker_prager_parameters_c.max_yield_stress)/volume_fraction, - viscoplastic_viscosity_guess); + viscoplastic_viscosity_guess_c = std::min(drucker_prager->compute_viscosity(drucker_prager_parameters_c.cohesion, + drucker_prager_parameters_c.angle_internal_friction, pressure, edot_ii, drucker_prager_parameters_c.max_yield_stress), + viscoplastic_viscosity_guess_c); } + inverse_viscoplastic_viscosity_guess += volume_fraction / viscoplastic_viscosity_guess_c; } else { @@ -216,6 +218,7 @@ namespace aspect } } + const double viscoplastic_viscosity_guess = 1. / inverse_viscoplastic_viscosity_guess; double viscoplastic_stress = 2. * viscoplastic_viscosity_guess * edot_ii / total_volume_fraction; double log_viscoplastic_stress = std::log(viscoplastic_stress); From 56a8c6aedc9d47e17161cf66a8e7bf3479afd621 Mon Sep 17 00:00:00 2001 From: Bob Myhill Date: Fri, 26 Jul 2024 19:08:52 +0100 Subject: [PATCH 7/8] fix comment --- source/material_model/rheology/composite_visco_plastic.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/source/material_model/rheology/composite_visco_plastic.cc b/source/material_model/rheology/composite_visco_plastic.cc index 810486999b4..4bffd3763c7 100644 --- a/source/material_model/rheology/composite_visco_plastic.cc +++ b/source/material_model/rheology/composite_visco_plastic.cc @@ -165,8 +165,8 @@ namespace aspect // Taking the minimum viscosity from // all of these elements provides an excellent first approximation // to the true viscosity of that composition. - // For a starting guess, we further assume that one composition takes - // up all the strain. + // For a starting guess, we further assume that each composition + // experiences the same strain rate. // The stress can then be calculated as 2 * viscoplastic_viscosity_guess * edot_ii double inverse_viscoplastic_viscosity_guess = 1./maximum_viscosity; From 2f09e153aac8e54424fb3200cdcf53b21b28c451 Mon Sep 17 00:00:00 2001 From: Bob Myhill Date: Tue, 30 Jul 2024 23:45:41 +0100 Subject: [PATCH 8/8] addressed comments --- .../rheology/composite_visco_plastic.h | 77 ++++++-- .../rheology/composite_visco_plastic.cc | 187 ++++++++++++------ tests/composite_viscous_outputs.cc | 13 +- tests/composite_viscous_outputs_isostress.cc | 2 +- tests/composite_viscous_outputs_no_peierls.cc | 11 +- 5 files changed, 207 insertions(+), 83 deletions(-) diff --git a/include/aspect/material_model/rheology/composite_visco_plastic.h b/include/aspect/material_model/rheology/composite_visco_plastic.h index 4ffd7a5a97e..8889479a907 100644 --- a/include/aspect/material_model/rheology/composite_visco_plastic.h +++ b/include/aspect/material_model/rheology/composite_visco_plastic.h @@ -38,6 +38,31 @@ namespace aspect namespace Rheology { + /** + * A namespace for selecting how to average the viscosity + * over multiple chemical compositions. + * + * Current options are: + * 'isostress': Reuss (harmonic) averaging of viscosities. + * 'isostrain': Voigt (arithmetic) averaging of viscosities. + */ + namespace ViscosityAveraging + { + enum Kind + { + isostress, + isostrain + }; + + /** + * Read the viscosity averaging scheme from the parameter file + * using the parameter name given in @p parameter_name and return + * the enum that corresponds to this operation. + */ + ViscosityAveraging::Kind + parse (const std::string ¶meter_name, + const ParameterHandler &prm); + } template class CompositeViscoPlastic : public ::aspect::SimulatorAccess @@ -83,6 +108,7 @@ namespace aspect const std::vector &phase_function_values = std::vector(), const std::vector &n_phase_transitions_per_composition = std::vector()) const; + private: /** * Compute the isostress viscosity over all compositional fields * based on the composite viscous creep law. @@ -104,7 +130,7 @@ namespace aspect /** * 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 - * over all compositional fields. Also updates the partial_strain_rates vector + * over all compositional fields. Also updates the partial_strain_rates vector. */ std::pair calculate_isostress_log_strain_rate_and_derivative(const std::vector, 4>> &logarithmic_strain_rates_and_stress_derivatives, @@ -112,7 +138,26 @@ namespace aspect std::vector &partial_strain_rates) const; /** - * Compute the compositional field viscosity + * Compute the isostrain viscosity over all compositional fields + * based on the composite viscous creep law. + * If @p expected_n_phases_per_composition points to a vector of + * unsigned integers this is considered the number of phase transitions + * for each compositional field and viscosity will be first computed on + * each phase and then averaged for each compositional field. + */ + double + compute_isostrain_viscosity (const double pressure, + const double temperature, + const double grain_size, + const std::vector &volume_fractions, + const SymmetricTensor<2,dim> &strain_rate, + std::vector &partial_strain_rates, + const std::vector &phase_function_values = std::vector(), + const std::vector &n_phase_transitions_per_composition = std::vector()) const; + + + /** + * Compute the viscosity for a single compositional field * based on the composite viscous creep law. * If @p n_phase_transitions_per_composition points to a vector of * unsigned integers this is considered the number of phase transitions @@ -131,26 +176,22 @@ namespace aspect /** * 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 + * with respect to log viscoplastic stress from individual log strain rate components + * for a single compositional field. + * Also updates the partial_strain_rates vector. */ std::pair calculate_composition_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: /** * Enumeration for selecting which type of viscosity averaging to use. */ - enum ViscosityAveragingScheme - { - isostrain, - isostress - } viscosity_averaging_scheme; + ViscosityAveraging::Kind viscosity_averaging_scheme; /** - * Whether to use different deformation mechanisms + * Whether to use different deformation mechanisms. */ bool use_diffusion_creep; bool use_dislocation_creep; @@ -158,10 +199,18 @@ namespace aspect bool use_drucker_prager; /** - * Vector storing which flow mechanisms are active + * Vector storing which flow mechanisms are active. */ std::vector active_flow_mechanisms; + /** + * Total number of decomposed strain rates. This includes the maximum + * viscosity limiter but not the minimum viscosity limiter, + * 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; + /** * Pointers to objects for computing deformation mechanism * strain rates and effective viscosities. @@ -178,7 +227,7 @@ namespace aspect /** * The maximum viscosity, imposed via an isoviscous damper - * in series with the composite viscoplastic element + * in series with the composite viscoplastic element. */ double maximum_viscosity; @@ -221,7 +270,7 @@ namespace aspect unsigned int stress_max_iteration_number; /** - * Useful number to aid in adding together exponentials + * Useful number to aid in adding together exponentials. */ const double logmin = std::log(std::numeric_limits::min()); }; diff --git a/source/material_model/rheology/composite_visco_plastic.cc b/source/material_model/rheology/composite_visco_plastic.cc index 4bffd3763c7..ab020ed431b 100644 --- a/source/material_model/rheology/composite_visco_plastic.cc +++ b/source/material_model/rheology/composite_visco_plastic.cc @@ -62,6 +62,29 @@ namespace aspect //} + + namespace ViscosityAveraging + { + ViscosityAveraging::Kind + parse (const std::string ¶meter_name, + const ParameterHandler &prm) + { + ViscosityAveraging::Kind scheme; + if (prm.get (parameter_name) == "isostress") + scheme = isostress; + else if (prm.get (parameter_name) == "isostrain") + scheme = isostrain; + else + { + AssertThrow(false, ExcMessage("Not a valid viscosity averaging scheme. Choose either isostress or isostrain.")); + } + + return scheme; + } + } + + + template CompositeViscoPlastic::CompositeViscoPlastic () = default; @@ -80,61 +103,49 @@ namespace aspect const std::vector &n_phase_transitions_per_composition) const { double viscosity = 0.; - partial_strain_rates.resize(5, 0.); - // Isostress averaging - if (viscosity_averaging_scheme == isostress) - { - viscosity += compute_isostress_viscosity (pressure, - temperature, - grain_size, - volume_fractions, - strain_rate, - partial_strain_rates, - phase_function_values, - n_phase_transitions_per_composition); - } + // Make sure partial_strain_rates is filled with zeros and is the right length + std::fill(partial_strain_rates.begin(), partial_strain_rates.end(), 0); + partial_strain_rates.resize(n_decomposed_strain_rates, 0.); - // Isostrain averaging - if (viscosity_averaging_scheme == isostrain) + // Compute the viscosity and the partial strain rates + // according to the isostress or isostrain viscosity averaging scheme. + switch (viscosity_averaging_scheme) { - double total_volume_fraction = 1.; - 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 - // a certain (small) fraction. - if (volume_fractions[composition] > 2.*std::numeric_limits::epsilon()) - { - std::vector partial_strain_rates_composition(5, 0.); - viscosity += (volume_fractions[composition] - * compute_composition_viscosity (pressure, - temperature, - grain_size, - composition, - strain_rate, - partial_strain_rates_composition, - phase_function_values, - n_phase_transitions_per_composition)); - for (unsigned int j=0; j < 5; ++j) - partial_strain_rates[j] += volume_fractions[composition] * partial_strain_rates_composition[j]; - } - else - { - total_volume_fraction -= volume_fractions[composition]; - } - - viscosity /= total_volume_fraction; - for (unsigned int j=0; j < 5; ++j) - partial_strain_rates[j] /= total_volume_fraction; - } + case ViscosityAveraging::Kind::isostress: + { + viscosity = compute_isostress_viscosity (pressure, + temperature, + grain_size, + volume_fractions, + strain_rate, + partial_strain_rates, + phase_function_values, + n_phase_transitions_per_composition); + break; + } + case ViscosityAveraging::Kind::isostrain: + { + viscosity = compute_isostrain_viscosity (pressure, + temperature, + grain_size, + volume_fractions, + strain_rate, + partial_strain_rates, + phase_function_values, + n_phase_transitions_per_composition); + break; + } + default: + { + Assert (false, ExcNotImplemented()); + } } return viscosity; } - template double CompositeViscoPlastic::compute_isostress_viscosity (const double pressure, @@ -228,6 +239,19 @@ namespace aspect for (unsigned int i = 0; i < active_compositions.size(); ++i) { + // In the isostress viscosity averaging formalism, + // the total strain rate is the sum of the strain rates for each + // composition i multiplied by the volume fraction of that + // composition. The strain rate for each composition is the sum + // of strain rates for each deformation mechanism j: + // edot_total = sum_i (volume_fraction_i * (sum_j (edot_ij))) + // This can be refactored: + // edot_total = sum_i (sum_j (volume_fraction_i * (edot_ij))) + // Therefore the logarithm of the contribution of each + // component-deformation mechanism to the total strain rate + // is: + // ln(edot_total_contrib_ij) = ln(volume_fraction_i) + ln(edot_ij) + const double log_volume_fraction = std::log(volume_fractions[active_compositions[i]]); if (use_diffusion_creep) @@ -269,7 +293,9 @@ namespace aspect // between these components. // The following while loop contains a Newton iteration to obtain the - // viscoplastic stress that is consistent with the total strain rate. + // viscoplastic stress that is consistent with the total strain rate + // and results in the correct partitioning of strain rate amongst + // the different flow components. unsigned int stress_iteration = 0; while (std::abs(log_strain_rate_residual) > log_strain_rate_residual_threshold && stress_iteration < stress_max_iteration_number) { @@ -382,6 +408,57 @@ namespace aspect + template + double + CompositeViscoPlastic::compute_isostrain_viscosity (const double pressure, + const double temperature, + const double grain_size, + const std::vector &volume_fractions, + const SymmetricTensor<2,dim> &strain_rate, + std::vector &partial_strain_rates, + const std::vector &phase_function_values, + const std::vector &n_phase_transitions_per_composition) const + { + // The isostrain viscosity is calculated by summing the viscosities + //of each composition weighted by their volume fractions. + double viscosity = 0.; + double total_volume_fraction = 1.; + 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 + // a certain (small) fraction. + if (volume_fractions[composition] > 2. * std::numeric_limits::epsilon()) + { + std::vector partial_strain_rates_composition(n_decomposed_strain_rates, 0.); + viscosity += (volume_fractions[composition] + * compute_composition_viscosity (pressure, + temperature, + grain_size, + composition, + strain_rate, + partial_strain_rates_composition, + phase_function_values, + n_phase_transitions_per_composition)); + for (unsigned int j=0; j < n_decomposed_strain_rates; ++j) + partial_strain_rates[j] += volume_fractions[composition] * partial_strain_rates_composition[j]; + } + else + { + total_volume_fraction -= volume_fractions[composition]; + } + } + + viscosity /= total_volume_fraction; + + for (unsigned int j=0; j < n_decomposed_strain_rates; ++j) + partial_strain_rates[j] /= total_volume_fraction; + + return viscosity; + } + + + template double CompositeViscoPlastic::compute_composition_viscosity (const double pressure, @@ -393,6 +470,7 @@ namespace aspect const std::vector &phase_function_values, const std::vector &n_phase_transitions_per_composition) const { + // Calculate the viscosity for a single compositional field based on the local state. // If strain rate is zero (like during the first time step) set it to some very small number // 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- @@ -537,7 +615,8 @@ namespace aspect const double viscoplastic_stress, std::vector &partial_strain_rates) const { - // The total strain rate + // Initialize the double containing the total strain rate + // for the single compositional field of interest. double viscoplastic_strain_rate_sum = 0.0; // The sum of the stress derivatives multiplied by the mechanism strain rates @@ -593,7 +672,7 @@ namespace aspect CompositeViscoPlastic::declare_parameters (ParameterHandler &prm) { prm.declare_entry ("Viscosity averaging scheme", "isostress", - Patterns::Selection("isostress|isostrain|unspecified"), + Patterns::Selection("isostress|isostrain"), "Determines the relationship between the conditions experienced by " "each chemical compositional field in a composite material. " "Select either isostress (the default option, Reuss averaging) " @@ -661,14 +740,8 @@ namespace aspect CompositeViscoPlastic::parse_parameters (ParameterHandler &prm, const std::unique_ptr> &expected_n_phases_per_composition) { - if (prm.get ("Viscosity averaging scheme") == "isostress") - viscosity_averaging_scheme = isostress; - else if (prm.get ("Viscosity averaging scheme") == "isostrain") - viscosity_averaging_scheme = isostrain; - else - { - AssertThrow(false, ExcMessage("Not a valid viscosity averaging scheme. Choose either isostress or isostrain.")); - } + viscosity_averaging_scheme = ViscosityAveraging::parse("Viscosity averaging scheme", + prm); // A background field is required by the subordinate material models number_of_chemical_compositions = this->introspection().n_chemical_composition_fields() + 1; diff --git a/tests/composite_viscous_outputs.cc b/tests/composite_viscous_outputs.cc index 0cec200924d..7b5efbcc0b7 100644 --- a/tests/composite_viscous_outputs.cc +++ b/tests/composite_viscous_outputs.cc @@ -20,6 +20,7 @@ #include #include +#include template void f(const aspect::SimulatorAccess &simulator_access, @@ -54,32 +55,32 @@ void f(const aspect::SimulatorAccess &simulator_access, 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); + 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, n_phases); + 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, n_phases); + 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, n_phases); + 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, n_phases); + 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. @@ -125,7 +126,7 @@ void f(const aspect::SimulatorAccess &simulator_access, temperature = 1000. + i*100.; // Compute the viscosity - viscosity = composite_creep->compute_composition_viscosity(pressure, temperature, grain_size, composition, strain_rate, partial_strain_rates); + viscosity = composite_creep->compute_viscosity(pressure, temperature, grain_size, volume_fractions, 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 diff --git a/tests/composite_viscous_outputs_isostress.cc b/tests/composite_viscous_outputs_isostress.cc index a1dd2ad1697..fcd433f9247 100644 --- a/tests/composite_viscous_outputs_isostress.cc +++ b/tests/composite_viscous_outputs_isostress.cc @@ -128,7 +128,7 @@ void f(const aspect::SimulatorAccess &simulator_access, temperature = 1000. + i*100.; // Compute the viscosity - viscosity = composite_creep->compute_isostress_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, 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_no_peierls.cc b/tests/composite_viscous_outputs_no_peierls.cc index 6d3078eabb0..4b15fb65b54 100644 --- a/tests/composite_viscous_outputs_no_peierls.cc +++ b/tests/composite_viscous_outputs_no_peierls.cc @@ -34,6 +34,7 @@ void f(const aspect::SimulatorAccess &simulator_access, // 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; @@ -54,26 +55,26 @@ void f(const aspect::SimulatorAccess &simulator_access, 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); + 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, n_phases); + 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, n_phases); + dislocation_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, n_phases); + 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. @@ -118,7 +119,7 @@ void f(const aspect::SimulatorAccess &simulator_access, temperature = 1000. + i*100.; // Compute the viscosity - viscosity = composite_creep->compute_composition_viscosity(pressure, temperature, grain_size, composition, strain_rate, partial_strain_rates); + viscosity = composite_creep->compute_viscosity(pressure, temperature, grain_size, volume_fractions, 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