From f621469f052e5119b342e30142e137fcc123de6b Mon Sep 17 00:00:00 2001 From: danieldouglas92 Date: Wed, 9 Jul 2025 17:34:30 -0600 Subject: [PATCH 1/2] add plastic prefactors --- .../compositional_plasticity_prefactors.h | 118 +++++++++ .../material_model/rheology/visco_plastic.h | 6 + .../compositional_plasticity_prefactors.cc | 224 ++++++++++++++++++ .../material_model/rheology/visco_plastic.cc | 14 ++ tests/test_plastic_prefactor.prm | 161 +++++++++++++ 5 files changed, 523 insertions(+) create mode 100644 include/aspect/material_model/rheology/compositional_plasticity_prefactors.h create mode 100644 source/material_model/rheology/compositional_plasticity_prefactors.cc create mode 100644 tests/test_plastic_prefactor.prm diff --git a/include/aspect/material_model/rheology/compositional_plasticity_prefactors.h b/include/aspect/material_model/rheology/compositional_plasticity_prefactors.h new file mode 100644 index 00000000000..84eef292ae9 --- /dev/null +++ b/include/aspect/material_model/rheology/compositional_plasticity_prefactors.h @@ -0,0 +1,118 @@ +/* + Copyright (C) 2024 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 + . +*/ + +#ifndef _aspect_material_model_rheology_compositional_plasticity_prefactors_h +#define _aspect_material_model_rheology_compositional_plasticity_prefactors_h + +#include +#include +#include +#include +#include + +namespace aspect +{ + namespace MaterialModel + { + namespace Rheology + { + /** + * A class that handles multiplication of viscosity for a given compositional + * field. The multiplication factors for each composition (viscosity + * prefactors) are also declared, parsed, and in some cases calculated in this class. + */ + template + class CompositionalPlasticityPrefactors : public ::aspect::SimulatorAccess + { + public: + /** + * Constructor. + */ + CompositionalPlasticityPrefactors(); + + /** + * Declare the parameters this function takes through input files. + */ + static + void + declare_parameters (ParameterHandler &prm); + + /** + * Read the parameters from the parameter file. + */ + void + parse_parameters (ParameterHandler &prm); + + /** + * Compute the viscosity. + */ + std::pair + compute_weakening_factor (const MaterialModel::MaterialModelInputs &in, + const unsigned int volume_fraction_index, + const unsigned int i, + const double static_cohesion, + const double static_friction_angle, + const Point &position) const; + + private: + /** + * The viscosity prefactors or terms used to calculate the viscosity + * prefactors, which are read in from the input file by the + * parse_parameters() function. Users can choose between different schemes. + * none: no viscosity change + * hk04_olivine_hydration: calculate the viscosity change due to hydrogen + * incorporation into olivine using Hirth & Kohlstaedt 2004 10.1029/138GM06. + * This method requires a composition called 'bound_fluid' which tracks the wt% + * water in the solid, which is used to compute an atomic ratio of H/Si ppm + * assuming 90 mol% forsterite and 10 mol% fayalite, and finally calculates + * a water fugacity. + * The prefactor for a given compositional field is multiplied with a + * base_viscosity value provided by the material model, which is then returned + * to the material model. + */ + enum WeakeningMechanism + { + none, + porosity, + function + } prefactor_mechanism; + + std::vector max_porosities_for_cohesion_prefactors; + std::vector max_porosities_for_friction_angle_prefactors; + + std::vector minimum_cohesions; + std::vector minimum_friction_angles; + /** + * Parsed functions that specify the friction angle which must be + * given in the input file using the function method. + */ + std::unique_ptr> prefactor_function; + + /** + * The coordinate representation to evaluate the function for the friction angle. + * Possible choices are depth, cartesian and spherical. + */ + Utilities::Coordinates::CoordinateSystem coordinate_system_prefactor_function; + + }; + } + } +} +#endif diff --git a/include/aspect/material_model/rheology/visco_plastic.h b/include/aspect/material_model/rheology/visco_plastic.h index dc9e36d05df..a46462aaedc 100644 --- a/include/aspect/material_model/rheology/visco_plastic.h +++ b/include/aspect/material_model/rheology/visco_plastic.h @@ -25,6 +25,7 @@ #include #include #include +#include #include #include #include @@ -352,6 +353,11 @@ namespace aspect */ Rheology::CompositionalViscosityPrefactors compositional_viscosity_prefactors; + /** + * Object for computing the viscosity multiplied by a given prefactor term. + */ + Rheology::CompositionalPlasticityPrefactors compositional_plasticity_prefactors; + /* * Object for computing plastic stresses, viscosities, and additional outputs */ diff --git a/source/material_model/rheology/compositional_plasticity_prefactors.cc b/source/material_model/rheology/compositional_plasticity_prefactors.cc new file mode 100644 index 00000000000..80786fc85d7 --- /dev/null +++ b/source/material_model/rheology/compositional_plasticity_prefactors.cc @@ -0,0 +1,224 @@ +/* + Copyright (C) 2024 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 +#include + + +#include +#include +#include +#include + +namespace aspect +{ + namespace MaterialModel + { + namespace Rheology + { + template + CompositionalPlasticityPrefactors::CompositionalPlasticityPrefactors () + = default; + + template + std::pair + CompositionalPlasticityPrefactors::compute_weakening_factor (const MaterialModel::MaterialModelInputs &in, + const unsigned int volume_fraction_index, + const unsigned int i, + const double input_cohesion, + const double input_friction_angle, + const Point &position) const + { + std::pair plasticity_parameters (input_cohesion, input_friction_angle); + switch (prefactor_mechanism) + { + case none: + { + return plasticity_parameters; + } + case porosity: + { + // Use the porosity to compute the cohesion and friction angle + // The porosity is computed in the MaterialModel::MaterialModelInputs + // and is available as in.porosity[i] + const unsigned int porosity_idx = this->introspection().compositional_index_for_name("porosity"); + plasticity_parameters.first *= std::max((max_porosities_for_cohesion_prefactors[volume_fraction_index] - in.composition[i][porosity_idx]) / max_porosities_for_cohesion_prefactors[volume_fraction_index], 0.0); + plasticity_parameters.second *= std::max((max_porosities_for_friction_angle_prefactors[volume_fraction_index] - in.composition[i][porosity_idx]) / max_porosities_for_friction_angle_prefactors[volume_fraction_index], 0.0); + + plasticity_parameters.first = std::max(plasticity_parameters.first, minimum_cohesions[volume_fraction_index]); + plasticity_parameters.second = std::max(plasticity_parameters.second, minimum_friction_angles[volume_fraction_index] * constants::degree_to_radians); + return plasticity_parameters; + } + case function: + { + // Use a given function input per composition to get the friction angle + Utilities::NaturalCoordinate point = + this->get_geometry_model().cartesian_to_other_coordinates(position, coordinate_system_prefactor_function); + + // we get time passed as seconds (always) but may want + // to reinterpret it in years + if (this->convert_output_to_years()) + prefactor_function->set_time (this->get_time() / year_in_seconds); + else + prefactor_function->set_time (this->get_time()); + + // determine the friction angle based on position and composition + // The volume_fraction_index is based on the number of chemical compositional fields. + // However, this plugin reads a function for every compositional field, regardless of + // its type. Therefore we have to get the correct index. + // If no fields or no chemical fields are present, but only background material, the index is zero. + // If chemical fields are present, volume_fractions will be of size 1+n_chemical_composition_fields. + // The size of chemical_composition_field_indices will be one less. + unsigned int index = 0; + if (this->introspection().composition_type_exists(CompositionalFieldDescription::chemical_composition)) + index = this->introspection().chemical_composition_field_indices()[volume_fraction_index-1]; + double prefactor_from_function = + prefactor_function->value(Utilities::convert_array_to_point(point.get_coordinates()),index); + + plasticity_parameters.first *= prefactor_from_function; + plasticity_parameters.second *= prefactor_from_function * constants::degree_to_radians; + + return plasticity_parameters; + } + } + } + + + template + void + CompositionalPlasticityPrefactors::declare_parameters (ParameterHandler &prm) + { + prm.declare_entry ("Maximum porosity for cohesion prefactor", "1.0", + Patterns::List(Patterns::Double(0., 1.0)), + "List of the maximum porosity value for calculating the prefactor for the cohesion. " + "entered as a volume fraction. If chosen to be 0.1 (10 percent porosity), the " + "cohesion will be linearly reduced from its default value when the model porosity is 0" + "to a maximum reduction when the model porosity is 0.1. Units: none."); + + prm.declare_entry ("Maximum porosity for friction angle prefactor", "1.0", + Patterns::List(Patterns::Double(0., 1.0)), + "List of the maximum porosity value for calculating the prefactor for the friction " + "angle, entered as a volume fraction. If chosen to be 0.1 (10 percent porosity), the " + "friction angle will be linearly reduced from its default value when the model porosity is 0" + "to a maximum reduction when the model porosity is 0.1. Units: none."); + + prm.declare_entry ("Plasticity prefactor scheme", "none", + Patterns::Selection("none|function|porosity"), + "Select what type of plasticity multiplicative prefactor scheme to apply. " + "Allowed entries are 'none', 'function', and 'porosity'. 'function' allows " + "the cohesion and friction angle prefactors to be specified using a function. " + "'porosity' determines a prefactor based on the modeled porosity, this scheme " + "requires that the user defines a compositional field called 'porosity'. 'none' " + "does not modify the cohesion or the friction angle. Units: none."); + + prm.declare_entry ("Minimum cohesions", "1e6", + Patterns::List(Patterns::Double(0.)), + "List of the minimum cohesions for limiting the reduction to the cohesion. Units: MPa."); + + prm.declare_entry ("Minimum friction angles", "1", + Patterns::List(Patterns::Double(0.)), + "List of the minimum friction angles for limiting the reduction to the friction angle. Units: none."); + + } + + + + template + void + CompositionalPlasticityPrefactors::parse_parameters (ParameterHandler &prm) + { + + // if friction is specified as a function + const unsigned int n_fields = this->n_compositional_fields() + 1; + if (prm.get ("Plasticity prefactor scheme") == "none") + prefactor_mechanism = none; + if (prm.get ("Plasticity prefactor scheme") == "porosity") + { + prefactor_mechanism = porosity; + + // Retrieve the list of chemical names + std::vector chemical_field_names = this->introspection().chemical_composition_field_names(); + std::vector compositional_field_names = this->introspection().get_composition_names(); + + // Establish that a background field is required here + compositional_field_names.insert(compositional_field_names.begin(), "background"); + chemical_field_names.insert(chemical_field_names.begin(),"background"); + + Utilities::MapParsing::Options options(chemical_field_names, "Maximum porosity for cohesion prefactor"); + + options.list_of_allowed_keys = compositional_field_names; + max_porosities_for_cohesion_prefactors = Utilities::MapParsing::parse_map_to_double_array (prm.get("Maximum porosity for cohesion prefactor"), + options); + max_porosities_for_friction_angle_prefactors = Utilities::MapParsing::parse_map_to_double_array (prm.get("Maximum porosity for friction angle prefactor"), + options); + + minimum_cohesions = Utilities::MapParsing::parse_map_to_double_array (prm.get("Minimum cohesions"), + options); + minimum_friction_angles = Utilities::MapParsing::parse_map_to_double_array (prm.get("Minimum friction angles"), + options); + } + if (prm.get ("Plasticity prefactor scheme") == "function") + { + prefactor_mechanism = function; + prm.enter_subsection("Prefactor function"); + { + coordinate_system_prefactor_function = Utilities::Coordinates::string_to_coordinate_system(prm.get("Coordinate system")); + try + { + prefactor_function + = std::make_unique>(n_fields); + prefactor_function->parse_parameters (prm); + } + catch (...) + { + std::cerr << "FunctionParser failed to parse\n" + << "\t friction function\n" + << "with expression \n" + << "\t' " << prm.get("Function expression") << "'"; + throw; + } + } + prm.leave_subsection(); + } + } + } + } +} + +// explicit instantiations +namespace aspect +{ + namespace MaterialModel + { + namespace Rheology + { +#define INSTANTIATE(dim) \ + template class CompositionalPlasticityPrefactors; + + ASPECT_INSTANTIATE(INSTANTIATE) + +#undef INSTANTIATE + } + } +} diff --git a/source/material_model/rheology/visco_plastic.cc b/source/material_model/rheology/visco_plastic.cc index 94f91c241ea..106ebec2a16 100644 --- a/source/material_model/rheology/visco_plastic.cc +++ b/source/material_model/rheology/visco_plastic.cc @@ -351,6 +351,14 @@ namespace aspect output_parameters.drucker_prager_parameters[j].angle_internal_friction, in.position[i]); + // Step 4c: Calculate weakening factors for the cohesion and friction angles based on mechanisms + // other than the strain. + std::pair plasticity_prefactors = compositional_plasticity_prefactors.compute_weakening_factor(in, j, i, + output_parameters.drucker_prager_parameters[j].cohesion, + output_parameters.drucker_prager_parameters[j].angle_internal_friction, + in.position[i]); + output_parameters.drucker_prager_parameters[j].cohesion = plasticity_prefactors.first; + output_parameters.drucker_prager_parameters[j].angle_internal_friction = plasticity_prefactors.second; // Step 5: plastic yielding // Determine if the pressure used in Drucker Prager plasticity will be capped at 0 (default). @@ -737,6 +745,9 @@ namespace aspect // Variable viscosity prefactor parameters Rheology::CompositionalViscosityPrefactors::declare_parameters(prm); + // Variable plasticity prefactor parameters + Rheology::CompositionalPlasticityPrefactors::declare_parameters(prm); + // Drucker Prager plasticity parameters Rheology::DruckerPrager::declare_parameters(prm); @@ -770,6 +781,9 @@ namespace aspect strain_rheology.initialize_simulator (this->get_simulator()); strain_rheology.parse_parameters(prm); + compositional_plasticity_prefactors.initialize_simulator (this->get_simulator()); + compositional_plasticity_prefactors.parse_parameters(prm); + friction_models.initialize_simulator (this->get_simulator()); friction_models.parse_parameters(prm); diff --git a/tests/test_plastic_prefactor.prm b/tests/test_plastic_prefactor.prm new file mode 100644 index 00000000000..1debccb89fe --- /dev/null +++ b/tests/test_plastic_prefactor.prm @@ -0,0 +1,161 @@ +# A test case that checks the water fugacity viscous prefactor +# multiplication scheme in combination with the visco plastic material model. +# Following the flow law for wet olivine composite (diffusion and dislocation) +# creep from Hirth & Kohlstaedt 2004 (10.1029/138GM06), the exact +# viscosity is calculated in python through: +# import numpy as np +# R = 8.3144621 +# P = 1e9 Pa +# T = 1173 K +# edot_ii = 1e-17 +# n_disl = 3.5 +# E_disl = 520e3 J/mol +# V_disl = 22e-6 m^3/mol +# r_disl = 1.2 +# A_disl = 1600 / 1e6 / 1e6**n / 1e6**r 1/Pa^(4.7)/s +# +# n_diff = 1 +# E_diff = 375e3 J/mol +# V_diff = 10e-6 m^3/mol +# r_diff = 0.7 +# d = 1e-3 m +# m = 3 +# A_diff = 2.5e7 / 1e6**n_diff / 1e6**m / 1e6**r_diff + +# A_H2O = 2.6e-5 1/Pa +# activation_energy_H2O = 40e3 J/mol/K +# activation_volume_H2O = 10e-6 m^3/mol +# +# M_H2O = 0.01801528 kg/mol +# M_olivine = 0.1470027 kg/mol +# +# weight_fraction_h2o = 0.01 +# weight_fraction_olivine = 1 - weight_fraction_h2o +# COH = (wt_fraction_h2o/M_H2O) / (weight_fraction_olivine/M_olivine) * 1e6 +# +# fH2O = COH / A_H2O * np.exp((activation_energy_H2O + P*activation_volume_H2O)/(R * T)) +# +# viscosity_disl = 0.5 * A_disl**(-1/n_disl) * edot_ii^((1-n_disl)/n_disl) * np.exp((E_disl + P*V_disl) / (n_disl*R*T)) * fH2O**(-r_disl/n_disl) +# viscosity_diff = 0.5 * A_diff**(-1/n_diff) * d**m * np.exp((E_diff + P*V_diff) / (n_diff*R*T)) * fH2O**(-r_diff/n_diff) +# +# viscosity_comp = (viscosity_disl * viscosity_diff) / (viscosity_disl + viscosity_diff) = 2.62814004323392e20 Pa s +# The viscosity should be equal to 2.6281e20 Pa s, and the material model should return this viscosity. + +# Global parameters +set Dimension = 2 +set Start time = 0 +set End time = 10 +set Maximum time step = 5 +set Use years in output instead of seconds = true +set Nonlinear solver scheme = single Advection, iterated Stokes +set Max nonlinear iterations = 1 +set Surface pressure = 1.e9 + +subsection Compositional fields + set Names of fields = porosity + set Number of fields = 1 +end + +subsection Initial composition model + set Model name = function + + subsection Function + set Function expression = if(y>=50e3, 0.05, 0) + end +end + +# Model geometry (100x100 km, 10 km spacing) +subsection Geometry model + set Model name = box + + subsection Box + set X repetitions = 10 + set Y repetitions = 10 + set X extent = 100e3 + set Y extent = 100e3 + end +end + +# Mesh refinement specifications +subsection Mesh refinement + set Initial adaptive refinement = 0 + set Initial global refinement = 0 + set Time steps between mesh refinement = 0 +end + +# No boundary-driven deformation, but the strain-rate used by the material +# model is set in the material model section through the parameters +# "Reference strain rate" and "Minimum strain rate". +subsection Boundary velocity model + set Tangential velocity boundary indicators = bottom, top, left, right +end + +# Initial temperature conditions +subsection Initial temperature model + set Model name = function + + subsection Function + set Function expression = 1173 + end +end + +subsection Material model + set Model name = visco plastic + + subsection Visco Plastic + # Composite creep, and choose the Hirth & Kohlstaedt 2004 olivine hydration + # viscosity prefactor scheme. + set Viscous flow law = composite + set Reference strain rate = 1e-17 + set Minimum strain rate = 1e-17 + + # Dislocation creep parameters + # The water fugacity exponent is 1.2, but we need to divide by + # the stress exponent, so we input r/n=0.34285714285714286. + set Prefactors for dislocation creep = 1.0095317511683098e-25 + set Stress exponents for dislocation creep = 3.5 + set Activation energies for dislocation creep = 520e3 + set Activation volumes for dislocation creep = 22e-6 + set Water fugacity exponents for dislocation creep = 0.34285714285714286 + + # Diffusion creep parameters + # The same approach is required as for the dislocation parameters, + # but the stress exponent for diffusion creep is 1. + set Prefactors for diffusion creep = 1.5773933612004842e-21 + set Stress exponents for diffusion creep = 1.0 + set Activation energies for diffusion creep = 375e3 + set Activation volumes for diffusion creep = 10e-6 + set Grain size = 1e-3 + set Water fugacity exponents for diffusion creep = 0.7 + + set Cohesions = 100e6 + set Angles of internal friction = 30 + set Plasticity prefactor scheme = porosity + set Maximum porosity for friction angle prefactor = 0.2 + set Maximum porosity for cohesion prefactor = 0.1 + + set Minimum cohesions = 10e6 + set Minimum friction angles = 25 + end +end + +# Gravity model +subsection Gravity model + set Model name = vertical + + subsection Vertical + set Magnitude = 0.0 + end +end + +# Post processing +subsection Postprocess + set List of postprocessors = material statistics, visualization + + subsection Visualization + set Interpolate output = false + set Time between graphical output = 0 + set List of output variables = material properties, named additional outputs + set Output format = vtu + end +end From c05cc7d0b3416e76942741b002eded3f2cad3578 Mon Sep 17 00:00:00 2001 From: danieldouglas92 Date: Thu, 10 Jul 2025 14:09:50 -0600 Subject: [PATCH 2/2] add some asserts --- .../compositional_plasticity_prefactors.h | 26 +++++++++---------- .../compositional_plasticity_prefactors.cc | 20 ++++++++------ tests/test_plastic_prefactor.prm | 2 +- 3 files changed, 25 insertions(+), 23 deletions(-) diff --git a/include/aspect/material_model/rheology/compositional_plasticity_prefactors.h b/include/aspect/material_model/rheology/compositional_plasticity_prefactors.h index 84eef292ae9..4bff5d0d8da 100644 --- a/include/aspect/material_model/rheology/compositional_plasticity_prefactors.h +++ b/include/aspect/material_model/rheology/compositional_plasticity_prefactors.h @@ -34,8 +34,8 @@ namespace aspect namespace Rheology { /** - * A class that handles multiplication of viscosity for a given compositional - * field. The multiplication factors for each composition (viscosity + * A class that handles multiplication of the plastic parameters for a given compositional + * field. The multiplication factors for each composition (plasticity * prefactors) are also declared, parsed, and in some cases calculated in this class. */ template @@ -61,7 +61,7 @@ namespace aspect parse_parameters (ParameterHandler &prm); /** - * Compute the viscosity. + * Compute the plasticity weakening factor. */ std::pair compute_weakening_factor (const MaterialModel::MaterialModelInputs &in, @@ -73,18 +73,16 @@ namespace aspect private: /** - * The viscosity prefactors or terms used to calculate the viscosity + * The plasticity prefactors or terms used to calculate the plasticity * prefactors, which are read in from the input file by the * parse_parameters() function. Users can choose between different schemes. - * none: no viscosity change - * hk04_olivine_hydration: calculate the viscosity change due to hydrogen - * incorporation into olivine using Hirth & Kohlstaedt 2004 10.1029/138GM06. - * This method requires a composition called 'bound_fluid' which tracks the wt% - * water in the solid, which is used to compute an atomic ratio of H/Si ppm - * assuming 90 mol% forsterite and 10 mol% fayalite, and finally calculates - * a water fugacity. + * none: no change in the plasticity parameters + * porosity: calculate the change in the plasticity parameters due to the amount + * of porosity at each quadrature point. + * function: calculate the change in the plasticity parameters by specifying the + * prefactors via a function. * The prefactor for a given compositional field is multiplied with a - * base_viscosity value provided by the material model, which is then returned + * cohesions and the friction angles provided by the material model, which is then returned * to the material model. */ enum WeakeningMechanism @@ -100,13 +98,13 @@ namespace aspect std::vector minimum_cohesions; std::vector minimum_friction_angles; /** - * Parsed functions that specify the friction angle which must be + * Parsed functions that specify the plasticity prefactors which must be * given in the input file using the function method. */ std::unique_ptr> prefactor_function; /** - * The coordinate representation to evaluate the function for the friction angle. + * The coordinate representation to evaluate the function for the plasticity prefactors. * Possible choices are depth, cartesian and spherical. */ Utilities::Coordinates::CoordinateSystem coordinate_system_prefactor_function; diff --git a/source/material_model/rheology/compositional_plasticity_prefactors.cc b/source/material_model/rheology/compositional_plasticity_prefactors.cc index 80786fc85d7..05de926b1e4 100644 --- a/source/material_model/rheology/compositional_plasticity_prefactors.cc +++ b/source/material_model/rheology/compositional_plasticity_prefactors.cc @@ -61,7 +61,6 @@ namespace aspect { // Use the porosity to compute the cohesion and friction angle // The porosity is computed in the MaterialModel::MaterialModelInputs - // and is available as in.porosity[i] const unsigned int porosity_idx = this->introspection().compositional_index_for_name("porosity"); plasticity_parameters.first *= std::max((max_porosities_for_cohesion_prefactors[volume_fraction_index] - in.composition[i][porosity_idx]) / max_porosities_for_cohesion_prefactors[volume_fraction_index], 0.0); plasticity_parameters.second *= std::max((max_porosities_for_friction_angle_prefactors[volume_fraction_index] - in.composition[i][porosity_idx]) / max_porosities_for_friction_angle_prefactors[volume_fraction_index], 0.0); @@ -72,7 +71,7 @@ namespace aspect } case function: { - // Use a given function input per composition to get the friction angle + // Use a given function input per composition to get the prefactors Utilities::NaturalCoordinate point = this->get_geometry_model().cartesian_to_other_coordinates(position, coordinate_system_prefactor_function); @@ -83,7 +82,7 @@ namespace aspect else prefactor_function->set_time (this->get_time()); - // determine the friction angle based on position and composition + // determine the prefactors based on position and composition // The volume_fraction_index is based on the number of chemical compositional fields. // However, this plugin reads a function for every compositional field, regardless of // its type. Therefore we have to get the correct index. @@ -96,6 +95,7 @@ namespace aspect double prefactor_from_function = prefactor_function->value(Utilities::convert_array_to_point(point.get_coordinates()),index); + // Multiply the cohesion and the friction angle by the prefactor. plasticity_parameters.first *= prefactor_from_function; plasticity_parameters.second *= prefactor_from_function * constants::degree_to_radians; @@ -110,14 +110,14 @@ namespace aspect CompositionalPlasticityPrefactors::declare_parameters (ParameterHandler &prm) { prm.declare_entry ("Maximum porosity for cohesion prefactor", "1.0", - Patterns::List(Patterns::Double(0., 1.0)), + Patterns::List(Patterns::Double(std::numeric_limits::epsilon(), 1.0)), "List of the maximum porosity value for calculating the prefactor for the cohesion. " "entered as a volume fraction. If chosen to be 0.1 (10 percent porosity), the " "cohesion will be linearly reduced from its default value when the model porosity is 0" "to a maximum reduction when the model porosity is 0.1. Units: none."); prm.declare_entry ("Maximum porosity for friction angle prefactor", "1.0", - Patterns::List(Patterns::Double(0., 1.0)), + Patterns::List(Patterns::Double(std::numeric_limits::epsilon(), 1.0)), "List of the maximum porosity value for calculating the prefactor for the friction " "angle, entered as a volume fraction. If chosen to be 0.1 (10 percent porosity), the " "friction angle will be linearly reduced from its default value when the model porosity is 0" @@ -149,12 +149,14 @@ namespace aspect CompositionalPlasticityPrefactors::parse_parameters (ParameterHandler &prm) { - // if friction is specified as a function const unsigned int n_fields = this->n_compositional_fields() + 1; if (prm.get ("Plasticity prefactor scheme") == "none") prefactor_mechanism = none; - if (prm.get ("Plasticity prefactor scheme") == "porosity") + else if (prm.get ("Plasticity prefactor scheme") == "porosity") { + AssertThrow(this->introspection().compositional_name_exists("porosity"), + ExcMessage("The 'porosity' Plasticity prefactor scheme work only if " + "is a compositional field called porosity.")); prefactor_mechanism = porosity; // Retrieve the list of chemical names @@ -178,7 +180,7 @@ namespace aspect minimum_friction_angles = Utilities::MapParsing::parse_map_to_double_array (prm.get("Minimum friction angles"), options); } - if (prm.get ("Plasticity prefactor scheme") == "function") + else if (prm.get ("Plasticity prefactor scheme") == "function") { prefactor_mechanism = function; prm.enter_subsection("Prefactor function"); @@ -201,6 +203,8 @@ namespace aspect } prm.leave_subsection(); } + else + AssertThrow(false, ExcMessage("Not a valid plasticity prefactor scheme")); } } } diff --git a/tests/test_plastic_prefactor.prm b/tests/test_plastic_prefactor.prm index 1debccb89fe..c60630be9b6 100644 --- a/tests/test_plastic_prefactor.prm +++ b/tests/test_plastic_prefactor.prm @@ -131,7 +131,7 @@ subsection Material model set Cohesions = 100e6 set Angles of internal friction = 30 set Plasticity prefactor scheme = porosity - set Maximum porosity for friction angle prefactor = 0.2 + set Maximum porosity for friction angle prefactor = 0. set Maximum porosity for cohesion prefactor = 0.1 set Minimum cohesions = 10e6