diff --git a/include/aspect/material_model/entropy_model.h b/include/aspect/material_model/entropy_model.h index 6f01646ffc1..47e83f9d73e 100644 --- a/include/aspect/material_model/entropy_model.h +++ b/include/aspect/material_model/entropy_model.h @@ -25,6 +25,7 @@ #include #include +#include #include #include #include @@ -124,31 +125,11 @@ namespace aspect double max_lateral_eta_variation; /** - * The value for thermal conductivity. It can be a constant - * for the whole domain, or P-T dependent. + * The thermal conductivity parametrization to use. This material + * model supports either a constant thermal conductivity or a + * pressure- and temperature-dependent thermal conductivity. */ - double thermal_conductivity_value; - double thermal_conductivity (const double temperature, - const double pressure, - const Point &position) const; - - enum ConductivityFormulation - { - constant, - p_T_dependent - } conductivity_formulation; - - /** - * Parameters for the temperature- and pressure dependence of the - * thermal conductivity. - */ - std::vector conductivity_transition_depths; - std::vector reference_thermal_conductivities; - std::vector conductivity_pressure_dependencies; - std::vector conductivity_reference_temperatures; - std::vector conductivity_exponents; - std::vector saturation_scaling; - double maximum_conductivity; + std::unique_ptr> thermal_conductivity; /** * Information about the location of data files. diff --git a/include/aspect/material_model/steinberger.h b/include/aspect/material_model/steinberger.h index 8b3f1b99ff9..d9e3df2c0e1 100644 --- a/include/aspect/material_model/steinberger.h +++ b/include/aspect/material_model/steinberger.h @@ -23,6 +23,7 @@ #include #include +#include #include #include @@ -116,7 +117,8 @@ namespace aspect * The viscosity of this model is based on the paper * Steinberger & Calderwood 2006: "Models of large-scale viscous flow in the * Earth's mantle with constraints from mineral physics and surface - * observations". The thermal conductivity is constant and the other + * observations". The thermal conductivity is constant or follows a + * pressure-temperature dependent approximation and the other * parameters are provided via lookup tables from the software PERPLEX. * * @ingroup MaterialModels @@ -202,19 +204,6 @@ namespace aspect private: - /** - * Compute the pressure- and temperature-dependent thermal - * conductivity either as a constant value, or based on the - * equation given in Stackhouse et al., 2015: First-principles - * calculations of the lattice thermal conductivity of the - * lower mantle, or based on the equation given in Tosi et al., - * 2013: Mantle dynamics with pressure- and temperature-dependent - * thermal expansivity and conductivity. - */ - double thermal_conductivity (const double temperature, - const double pressure, - const Point &position) const; - /** * Whether the compositional fields representing mass fractions * should be normalized to one when computing their fractions @@ -248,31 +237,11 @@ namespace aspect bool use_lateral_average_temperature; /** - * The value of the thermal conductivity if a constant thermal - * conductivity is used for the whole domain. - */ - double thermal_conductivity_value; - - /** - * Enumeration for selecting which type of conductivity law to use. - */ - enum ConductivityFormulation - { - constant, - p_T_dependent - } conductivity_formulation; - - /** - * Parameters for the temperature- and pressure dependence of the - * thermal conductivity. + * The thermal conductivity parametrization to use. This material + * model supports either a constant thermal conductivity or a + * pressure- and temperature-dependent thermal conductivity. */ - std::vector conductivity_transition_depths; - std::vector reference_thermal_conductivities; - std::vector conductivity_pressure_dependencies; - std::vector conductivity_reference_temperatures; - std::vector conductivity_exponents; - std::vector saturation_scaling; - double maximum_conductivity; + std::unique_ptr> thermal_conductivity; /** * Compositional prefactors with which to multiply the reference viscosity. diff --git a/include/aspect/material_model/thermal_conductivity/tosi_stackhouse.h b/include/aspect/material_model/thermal_conductivity/tosi_stackhouse.h new file mode 100644 index 00000000000..7942f7a8ab9 --- /dev/null +++ b/include/aspect/material_model/thermal_conductivity/tosi_stackhouse.h @@ -0,0 +1,114 @@ +/* + Copyright (C) 2025 - 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_thermal_conductivity_tosi_stackhouse_h +#define _aspect_material_model_thermal_conductivity_tosi_stackhouse_h + +#include + + +namespace aspect +{ + namespace MaterialModel + { + namespace ThermalConductivity + { + using namespace dealii; + + /** + * A class that implements a pressure- and temperature-dependent thermal conductivity + * following the formulation of + * + * Nicola Tosi, David A. Yuen, Nico de Koker, Renata M. Wentzcovitch, + * Mantle dynamics with pressure- and temperature-dependent thermal expansivity and conductivity, + * Physics of the Earth and Planetary Interiors, Volume 217, + * 2013, Pages 48-58, ISSN 0031-9201, https://doi.org/10.1016/j.pepi.2013.02.004, + * + * and + * + * Stephen Stackhouse, Lars Stixrude, Bijaya B. Karki, + * First-principles calculations of the lattice thermal conductivity of the lower mantle, + * Earth and Planetary Science Letters, Volume 427, 2015, Pages 11-17, ISSN 0012-821X, + * https://doi.org/10.1016/j.epsl.2015.06.050. + * + * The thermal conductivity parameter sets can be chosen in such a + * way that either the Stackhouse or the Tosi relations are used. + * The conductivity description can consist of several layers with + * different sets of parameters. Note that the Stackhouse + * parametrization is only valid for the lower mantle (bridgmanite). + * + * The default parameters of this class use the Tosi parametrization in the upper + * mantle and the Stackhouse parametrization in the lower mantle, which is how + * it was used in the publication + * + * Juliane Dannberg, Rene Gassmöller, Daniele Thallner, Frederick LaCombe, Courtney Sprain, + * Changes in core–mantle boundary heat flux patterns throughout the supercontinent cycle, + * Geophysical Journal International, Volume 237, Issue 3, June 2024, Pages 1251–1274, https://doi.org/10.1093/gji/ggae075, + * + * which introduced this implementation. + * + * @ingroup MaterialModels + */ + template + class TosiStackhouse : public Interface, public aspect::SimulatorAccess + { + public: + /** + * Function to compute the thermal conductivities in @p out given the + * inputs in @p in. + */ + void evaluate (const MaterialModel::MaterialModelInputs &in, + MaterialModel::MaterialModelOutputs &out) const override; + + /** + * Declare the parameters this plugin takes through input files. + */ + static + void + declare_parameters (ParameterHandler &prm); + + /** + * Read the parameters from the parameter file. + */ + void + parse_parameters (ParameterHandler &prm) override; + + private: + /** + * Parameters for the temperature and pressure dependence of the + * thermal conductivity. + */ + std::vector conductivity_transition_depths; + std::vector reference_thermal_conductivities; + std::vector conductivity_pressure_dependencies; + std::vector conductivity_reference_temperatures; + std::vector conductivity_exponents; + std::vector saturation_scaling; + + /** + * The maximum allowed thermal conductivity. + */ + double maximum_conductivity; + }; + } + } +} + +#endif diff --git a/source/material_model/entropy_model.cc b/source/material_model/entropy_model.cc index c6a7daa1a17..db5d9762a9d 100644 --- a/source/material_model/entropy_model.cc +++ b/source/material_model/entropy_model.cc @@ -20,6 +20,8 @@ #include +#include +#include #include #include @@ -116,50 +118,6 @@ namespace aspect - template - double - EntropyModel:: - thermal_conductivity (const double temperature, - const double pressure, - const Point &position) const - { - if (conductivity_formulation == constant) - return thermal_conductivity_value; - - else if (conductivity_formulation == p_T_dependent) - { - // Find the conductivity layer that corresponds to the depth of the evaluation point. - const double depth = this->get_geometry_model().depth(position); - unsigned int layer_index = std::distance(conductivity_transition_depths.begin(), - std::lower_bound(conductivity_transition_depths.begin(),conductivity_transition_depths.end(), depth)); - - const double p_dependence = reference_thermal_conductivities[layer_index] + conductivity_pressure_dependencies[layer_index] * pressure; - - // Make reasonably sure we will not compute any invalid values due to the temperature-dependence. - // Since both the temperature-dependence and the saturation time scale with (Tref/T), we have to - // make sure we can compute the square of this number. If the temperature is small enough to - // be close to yielding NaN values, the conductivity will be set to the maximum value anyway. - const double T = std::max(temperature, std::sqrt(std::numeric_limits::min()) * conductivity_reference_temperatures[layer_index]); - const double T_dependence = std::pow(conductivity_reference_temperatures[layer_index] / T, conductivity_exponents[layer_index]); - - // Function based on the theory of Roufosse and Klemens (1974) that accounts for saturation. - // For the Tosi formulation, the scaling should be zero so that this term is 1. - double saturation_function = 1.0; - if (1./T_dependence > 1.) - saturation_function = (1. - saturation_scaling[layer_index]) - + saturation_scaling[layer_index] * (2./3. * std::sqrt(T_dependence) + 1./3. * 1./T_dependence); - - return std::min(p_dependence * saturation_function * T_dependence, maximum_conductivity); - } - else - { - AssertThrow(false, ExcNotImplemented()); - return numbers::signaling_nan(); - } - } - - - template void EntropyModel::evaluate(const MaterialModel::MaterialModelInputs &in, @@ -179,6 +137,10 @@ namespace aspect std::vector volume_fractions (material_file_names.size()); std::vector mass_fractions (material_file_names.size()); + // We need to make a copy of the material model inputs because we want to replace the + // temperature with the temperature from the lookup table. + MaterialModel::MaterialModelInputs adjusted_inputs(in); + for (unsigned int i=0; i < in.n_evaluation_points(); ++i) { // Use the adiabatic pressure instead of the real one, @@ -188,6 +150,7 @@ namespace aspect // Also convert pressure from Pa to bar, bar is used in the table. const double entropy = in.composition[i][entropy_index]; const double pressure = this->get_adiabatic_conditions().pressure(in.position[i]) / 1.e5; + adjusted_inputs.temperature[i] = entropy_reader[0]->temperature(entropy,pressure); // Loop over all material files, and store the looked-up values for all compositions. for (unsigned int j=0; jtemperature(entropy,pressure); - out.thermal_conductivities[i] = thermal_conductivity(temperature_lookup, in.pressure[i], in.position[i]); - out.entropy_derivative_pressure[i] = 0.; out.entropy_derivative_temperature[i] = 0.; for (unsigned int c=0; c *prescribed_temperature_out = out.template get_additional_output>()) { - prescribed_temperature_out->prescribed_temperature_outputs[i] = temperature_lookup; + prescribed_temperature_out->prescribed_temperature_outputs[i] = adjusted_inputs.temperature[i]; } // Calculate Viscosity @@ -257,14 +216,14 @@ namespace aspect : this->get_parameters().adiabatic_surface_temperature; - const double delta_temperature = temperature_lookup-reference_temperature; + const double delta_temperature = adjusted_inputs.temperature[i] - reference_temperature; // Steinberger & Calderwood viscosity - if (temperature_lookup*reference_temperature == 0) + if (adjusted_inputs.temperature[i]*reference_temperature == 0) out.viscosities[i] = max_eta; else { - double vis_lateral = std::exp(-lateral_viscosity_prefactor_lookup->lateral_viscosity(depth)*delta_temperature/(temperature_lookup*reference_temperature)); + double vis_lateral = std::exp(-lateral_viscosity_prefactor_lookup->lateral_viscosity(depth)*delta_temperature/(adjusted_inputs.temperature[i]*reference_temperature)); // lateral vis variation vis_lateral = std::max(std::min((vis_lateral),max_lateral_eta_variation),1/max_lateral_eta_variation); @@ -313,6 +272,10 @@ namespace aspect seismic_out->vs[i] = MaterialUtilities::average_value (volume_fractions, vs, MaterialUtilities::arithmetic); } } + + // Evaluate thermal conductivity. This has to happen after + // the evaluation of the equation of state and calculation of temperature. + thermal_conductivity->evaluate(adjusted_inputs, out); } @@ -367,10 +330,9 @@ namespace aspect "The value of the cohesion, $C$. The extremely large default" "cohesion value (1e20 Pa) prevents the viscous stress from " "exceeding the yield stress. Units: \\si{\\pascal}."); - prm.declare_entry ("Thermal conductivity", "4.7", - Patterns::Double (0), - "The value of the thermal conductivity $k$. " - "Units: \\si{\\watt\\per\\meter\\per\\kelvin}."); + + // Thermal conductivity parameters + ThermalConductivity::Constant::declare_parameters(prm); prm.declare_entry ("Thermal conductivity formulation", "constant", Patterns::Selection("constant|p-T-dependent"), "Which law should be used to compute the thermal conductivity. " @@ -388,55 +350,8 @@ namespace aspect "The conductivity description can consist of several layers with " "different sets of parameters. Note that the Stackhouse " "parametrization is only valid for the lower mantle (bridgmanite)."); - prm.declare_entry ("Thermal conductivity transition depths", "410000, 520000, 660000", - Patterns::List(Patterns::Double (0.)), - "A list of depth values that indicate where the transitions between " - "the different conductivity parameter sets should occur in the " - "'p-T-dependent' Thermal conductivity formulation (in most cases, " - "this will be the depths of major mantle phase transitions). " - "Units: \\si{\\meter}."); - prm.declare_entry ("Reference thermal conductivities", "2.47, 3.81, 3.52, 4.9", - Patterns::List(Patterns::Double (0.)), - "A list of base values of the thermal conductivity for each of the " - "horizontal layers in the 'p-T-dependent' thermal conductivity " - "formulation. Pressure- and temperature-dependence will be applied" - "on top of this base value, according to the parameters 'Pressure " - "dependencies of thermal conductivity' and 'Reference temperatures " - "for thermal conductivity'. " - "Units: \\si{\\watt\\per\\meter\\per\\kelvin}"); - prm.declare_entry ("Pressure dependencies of thermal conductivity", "3.3e-10, 3.4e-10, 3.6e-10, 1.05e-10", - Patterns::List(Patterns::Double ()), - "A list of values that determine the linear scaling of the " - "thermal conductivity with the pressure in the 'p-T-dependent' " - "thermal conductivity formulation. " - "Units: \\si{\\watt\\per\\meter\\per\\kelvin\\per\\pascal}."); - prm.declare_entry ("Reference temperatures for thermal conductivity", "300, 300, 300, 1200", - Patterns::List(Patterns::Double (0.)), - "A list of values of reference temperatures used to determine " - "the temperature-dependence of the thermal conductivity in the " - "'p-T-dependent' thermal conductivity formulation. " - "Units: \\si{\\kelvin}."); - prm.declare_entry ("Thermal conductivity exponents", "0.48, 0.56, 0.61, 1.0", - Patterns::List(Patterns::Double (0.)), - "A list of exponents in the temperature-dependent term of the " - "'p-T-dependent' thermal conductivity formulation. Note that this " - "exponent is not used (and should have a value of 1) in the " - "formulation of Stackhouse et al. (2015). " - "Units: none."); - prm.declare_entry ("Saturation prefactors", "0, 0, 0, 1", - Patterns::List(Patterns::Double (0., 1.)), - "A list of values that indicate how a given layer in the " - "conductivity formulation should take into account the effects " - "of saturation on the temperature-dependence of the thermal " - "conducitivity. This factor is multiplied with a saturation function " - "based on the theory of Roufosse and Klemens, 1974. A value of 1 " - "reproduces the formulation of Stackhouse et al. (2015), a value of " - "0 reproduces the formulation of Tosi et al., (2013). " - "Units: none."); - prm.declare_entry ("Maximum thermal conductivity", "1000", - Patterns::Double (0.), - "The maximum thermal conductivity that is allowed in the " - "model. Larger values will be cut off."); + ThermalConductivity::TosiStackhouse::declare_parameters(prm); + prm.leave_subsection(); } @@ -464,38 +379,20 @@ namespace aspect min_eta = prm.get_double ("Minimum viscosity"); max_eta = prm.get_double ("Maximum viscosity"); max_lateral_eta_variation = prm.get_double ("Maximum lateral viscosity variation"); - thermal_conductivity_value = prm.get_double ("Thermal conductivity"); + // Thermal conductivity parameters if (prm.get ("Thermal conductivity formulation") == "constant") - conductivity_formulation = constant; + thermal_conductivity = std::make_unique>(); else if (prm.get ("Thermal conductivity formulation") == "p-T-dependent") - conductivity_formulation = p_T_dependent; + { + thermal_conductivity = std::make_unique>(); + if (SimulatorAccess *sim = dynamic_cast*>(thermal_conductivity.get())) + sim->initialize_simulator (this->get_simulator()); + } else AssertThrow(false, ExcMessage("Not a valid thermal conductivity formulation")); - conductivity_transition_depths = Utilities::string_to_double - (Utilities::split_string_list(prm.get ("Thermal conductivity transition depths"))); - const unsigned int n_conductivity_layers = conductivity_transition_depths.size() + 1; - AssertThrow (std::is_sorted(conductivity_transition_depths.begin(), conductivity_transition_depths.end()), - ExcMessage("The list of 'Thermal conductivity transition depths' must " - "be sorted such that the values increase monotonically.")); - - reference_thermal_conductivities = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Reference thermal conductivities"))), - n_conductivity_layers, - "Reference thermal conductivities"); - conductivity_pressure_dependencies = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Pressure dependencies of thermal conductivity"))), - n_conductivity_layers, - "Pressure dependencies of thermal conductivity"); - conductivity_reference_temperatures = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Reference temperatures for thermal conductivity"))), - n_conductivity_layers, - "Reference temperatures for thermal conductivity"); - conductivity_exponents = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Thermal conductivity exponents"))), - n_conductivity_layers, - "Thermal conductivity exponents"); - saturation_scaling = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Saturation prefactors"))), - n_conductivity_layers, - "Saturation prefactors"); - maximum_conductivity = prm.get_double ("Maximum thermal conductivity"); + thermal_conductivity->parse_parameters(prm); angle_of_internal_friction = prm.get_double ("Angle of internal friction") * constants::degree_to_radians; cohesion = prm.get_double("Cohesion"); diff --git a/source/material_model/steinberger.cc b/source/material_model/steinberger.cc index 8cf8ec03f2a..245137ce7e5 100644 --- a/source/material_model/steinberger.cc +++ b/source/material_model/steinberger.cc @@ -21,6 +21,8 @@ #include #include +#include +#include #include #include #include @@ -216,50 +218,6 @@ namespace aspect - template - double - Steinberger:: - thermal_conductivity (const double temperature, - const double pressure, - const Point &position) const - { - if (conductivity_formulation == constant) - return thermal_conductivity_value; - - else if (conductivity_formulation == p_T_dependent) - { - // Find the conductivity layer that corresponds to the depth of the evaluation point. - const double depth = this->get_geometry_model().depth(position); - unsigned int layer_index = std::distance(conductivity_transition_depths.begin(), - std::lower_bound(conductivity_transition_depths.begin(),conductivity_transition_depths.end(), depth)); - - const double p_dependence = reference_thermal_conductivities[layer_index] + conductivity_pressure_dependencies[layer_index] * pressure; - - // Make reasonably sure we will not compute any invalid values due to the temperature-dependence. - // Since both the temperature-dependence and the saturation term scale with (Tref/T), we have to - // make sure we can compute the square of this number. If the temperature is small enough to - // be close to yielding NaN values, the conductivity will be set to the maximum value anyway. - const double T = std::max(temperature, std::sqrt(std::numeric_limits::min()) * conductivity_reference_temperatures[layer_index]); - const double T_dependence = std::pow(conductivity_reference_temperatures[layer_index] / T, conductivity_exponents[layer_index]); - - // Function based on the theory of Roufosse and Klemens (1974) that accounts for saturation. - // For the Tosi formulation, the scaling should be zero so that this term is 1. - double saturation_function = 1.0; - if (1./T_dependence > 1.) - saturation_function = (1. - saturation_scaling[layer_index]) - + saturation_scaling[layer_index] * (2./3. * std::sqrt(T_dependence) + 1./3. * 1./T_dependence); - - return std::min(p_dependence * saturation_function * T_dependence, maximum_conductivity); - } - else - { - AssertThrow(false, ExcNotImplemented()); - return numbers::signaling_nan(); - } - } - - - template void Steinberger::evaluate(const MaterialModel::MaterialModelInputs &in, @@ -276,10 +234,10 @@ namespace aspect // Evaluate the equation of state properties over all evaluation points equation_of_state.evaluate(eos_in, eos_outputs); + thermal_conductivity->evaluate(in, out); for (unsigned int i=0; i < in.n_evaluation_points(); ++i) { - out.thermal_conductivities[i] = thermal_conductivity(in.temperature[i], in.pressure[i], in.position[i]); for (unsigned int c=0; c::declare_parameters (prm); prm.declare_entry ("Thermal conductivity formulation", "constant", Patterns::Selection("constant|p-T-dependent"), "Which law should be used to compute the thermal conductivity. " @@ -416,55 +372,7 @@ namespace aspect "The conductivity description can consist of several layers with " "different sets of parameters. Note that the Stackhouse " "parametrization is only valid for the lower mantle (bridgmanite)."); - prm.declare_entry ("Thermal conductivity transition depths", "410000, 520000, 660000", - Patterns::List(Patterns::Double (0.)), - "A list of depth values that indicate where the transitions between " - "the different conductivity parameter sets should occur in the " - "'p-T-dependent' Thermal conductivity formulation (in most cases, " - "this will be the depths of major mantle phase transitions). " - "Units: \\si{\\meter}."); - prm.declare_entry ("Reference thermal conductivities", "2.47, 3.81, 3.52, 4.9", - Patterns::List(Patterns::Double (0.)), - "A list of base values of the thermal conductivity for each of the " - "horizontal layers in the 'p-T-dependent' Thermal conductivity " - "formulation. Pressure- and temperature-dependence will be applied" - "on top of this base value, according to the parameters 'Pressure " - "dependencies of thermal conductivity' and 'Reference temperatures " - "for thermal conductivity'. " - "Units: \\si{\\watt\\per\\meter\\per\\kelvin}"); - prm.declare_entry ("Pressure dependencies of thermal conductivity", "3.3e-10, 3.4e-10, 3.6e-10, 1.05e-10", - Patterns::List(Patterns::Double ()), - "A list of values that determine the linear scaling of the " - "thermal conductivity with the pressure in the 'p-T-dependent' " - "Thermal conductivity formulation. " - "Units: \\si{\\watt\\per\\meter\\per\\kelvin\\per\\pascal}."); - prm.declare_entry ("Reference temperatures for thermal conductivity", "300, 300, 300, 1200", - Patterns::List(Patterns::Double (0.)), - "A list of values of reference temperatures used to determine " - "the temperature-dependence of the thermal conductivity in the " - "'p-T-dependent' Thermal conductivity formulation. " - "Units: \\si{\\kelvin}."); - prm.declare_entry ("Thermal conductivity exponents", "0.48, 0.56, 0.61, 1.0", - Patterns::List(Patterns::Double (0.)), - "A list of exponents in the temperature-dependent term of the " - "'p-T-dependent' Thermal conductivity formulation. Note that this " - "exponent is not used (and should have a value of 1) in the " - "formulation of Stackhouse et al. (2015). " - "Units: none."); - prm.declare_entry ("Saturation prefactors", "0, 0, 0, 1", - Patterns::List(Patterns::Double (0., 1.)), - "A list of values that indicate how a given layer in the " - "conductivity formulation should take into account the effects " - "of saturation on the temperature-dependence of the thermal " - "conducitivity. This factor is multiplied with a saturation function " - "based on the theory of Roufosse and Klemens, 1974. A value of 1 " - "reproduces the formulation of Stackhouse et al. (2015), a value of " - "0 reproduces the formulation of Tosi et al., (2013). " - "Units: none."); - prm.declare_entry ("Maximum thermal conductivity", "1000", - Patterns::Double (0.), - "The maximum thermal conductivity that is allowed in the " - "model. Larger values will be cut off."); + ThermalConductivity::TosiStackhouse::declare_parameters (prm); // Table lookup parameters EquationOfState::ThermodynamicTableLookup::declare_parameters(prm); @@ -495,40 +403,20 @@ namespace aspect max_lateral_eta_variation = prm.get_double ("Maximum lateral viscosity variation"); viscosity_averaging_scheme = MaterialUtilities::parse_compositional_averaging_operation ("Viscosity averaging scheme", prm); - thermal_conductivity_value = prm.get_double ("Thermal conductivity"); - // Rheological parameters + // Thermal conductivity parameters if (prm.get ("Thermal conductivity formulation") == "constant") - conductivity_formulation = constant; + thermal_conductivity = std::make_unique>(); else if (prm.get ("Thermal conductivity formulation") == "p-T-dependent") - conductivity_formulation = p_T_dependent; + { + thermal_conductivity = std::make_unique>(); + if (SimulatorAccess *sim = dynamic_cast*>(thermal_conductivity.get())) + sim->initialize_simulator (this->get_simulator()); + } else AssertThrow(false, ExcMessage("Not a valid thermal conductivity formulation")); - conductivity_transition_depths = Utilities::string_to_double - (Utilities::split_string_list(prm.get ("Thermal conductivity transition depths"))); - const unsigned int n_conductivity_layers = conductivity_transition_depths.size() + 1; - - AssertThrow (std::is_sorted(conductivity_transition_depths.begin(), conductivity_transition_depths.end()), - ExcMessage("The list of 'Thermal conductivity transition depths' must " - "be sorted such that the values increase monotonically.")); - - reference_thermal_conductivities = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Reference thermal conductivities"))), - n_conductivity_layers, - "Reference thermal conductivities"); - conductivity_pressure_dependencies = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Pressure dependencies of thermal conductivity"))), - n_conductivity_layers, - "Pressure dependencies of thermal conductivity"); - conductivity_reference_temperatures = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Reference temperatures for thermal conductivity"))), - n_conductivity_layers, - "Reference temperatures for thermal conductivity"); - conductivity_exponents = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Thermal conductivity exponents"))), - n_conductivity_layers, - "Thermal conductivity exponents"); - saturation_scaling = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Saturation prefactors"))), - n_conductivity_layers, - "Saturation prefactors"); - maximum_conductivity = prm.get_double ("Maximum thermal conductivity"); + thermal_conductivity->parse_parameters(prm); // Parse the table lookup parameters equation_of_state.initialize_simulator (this->get_simulator()); diff --git a/source/material_model/thermal_conductivity/tosi_stackhouse.cc b/source/material_model/thermal_conductivity/tosi_stackhouse.cc new file mode 100644 index 00000000000..7b7624fa3dd --- /dev/null +++ b/source/material_model/thermal_conductivity/tosi_stackhouse.cc @@ -0,0 +1,166 @@ +/* + Copyright (C) 2025 - 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 + +namespace aspect +{ + namespace MaterialModel + { + namespace ThermalConductivity + { + template + void + TosiStackhouse::evaluate (const MaterialModel::MaterialModelInputs &in, + MaterialModel::MaterialModelOutputs &out) const + { + const unsigned int n_evaluation_points = in.n_evaluation_points(); + for (unsigned int i=0; iget_geometry_model().depth(in.position[i]); + const unsigned int layer_index = std::distance(conductivity_transition_depths.begin(), + std::lower_bound(conductivity_transition_depths.begin(),conductivity_transition_depths.end(), depth)); + + const double p_dependence = reference_thermal_conductivities[layer_index] + conductivity_pressure_dependencies[layer_index] * in.pressure[i]; + + // Make reasonably sure we will not compute any invalid values due to the temperature-dependence. + // Since both the temperature-dependence and the saturation term scale with (Tref/T), we have to + // make sure we can compute the square of this number. If the temperature is small enough to + // be close to yielding NaN values, the conductivity will be set to the maximum value anyway. + const double T = std::max(in.temperature[i], std::sqrt(std::numeric_limits::min()) * conductivity_reference_temperatures[layer_index]); + const double T_dependence = std::pow(conductivity_reference_temperatures[layer_index] / T, conductivity_exponents[layer_index]); + + // Function based on the theory of Roufosse and Klemens (1974) that accounts for saturation. + // For the Tosi formulation, the scaling should be zero so that this term is 1. + double saturation_function = 1.0; + if (1./T_dependence > 1.) + saturation_function = (1. - saturation_scaling[layer_index]) + + saturation_scaling[layer_index] * (2./3. * std::sqrt(T_dependence) + 1./3. * 1./T_dependence); + + out.thermal_conductivities[i] = std::min(p_dependence * saturation_function * T_dependence, maximum_conductivity); + } + } + + + + template + void + TosiStackhouse::declare_parameters (ParameterHandler &prm) + { + prm.declare_entry ("Thermal conductivity transition depths", "410000, 520000, 660000", + Patterns::List(Patterns::Double (0.)), + "A list of depth values that indicate where the transitions between " + "the different conductivity parameter sets should occur (in most cases, " + "these will be the depths of major phase transitions). " + "Units: \\si{\\meter}."); + prm.declare_entry ("Reference thermal conductivities", "2.47, 3.81, 3.52, 4.9", + Patterns::List(Patterns::Double (0.)), + "A list of base values of the thermal conductivity for each of the " + "horizontal layers. Pressure- and temperature-dependence will be applied " + "on top of this base value, according to the parameters 'Pressure " + "dependencies of thermal conductivity' and 'Reference temperatures " + "for thermal conductivity'. " + "Units: \\si{\\watt\\per\\meter\\per\\kelvin}"); + prm.declare_entry ("Pressure dependencies of thermal conductivity", "3.3e-10, 3.4e-10, 3.6e-10, 1.05e-10", + Patterns::List(Patterns::Double ()), + "A list of values that determine the linear scaling of the " + "thermal conductivity with pressure. " + "Units: \\si{\\watt\\per\\meter\\per\\kelvin\\per\\pascal}."); + prm.declare_entry ("Reference temperatures for thermal conductivity", "300, 300, 300, 1200", + Patterns::List(Patterns::Double (0.)), + "A list of values of reference temperatures used to determine " + "the temperature-dependence of the thermal conductivity. " + "Units: \\si{\\kelvin}."); + prm.declare_entry ("Thermal conductivity exponents", "0.48, 0.56, 0.61, 1.0", + Patterns::List(Patterns::Double (0.)), + "A list of exponents in the temperature-dependent term of the " + "conductivity formulation. Note that this " + "exponent is not used (and should have a value of 1) in the " + "formulation of Stackhouse et al. (2015). " + "Units: none."); + prm.declare_entry ("Saturation prefactors", "0, 0, 0, 1", + Patterns::List(Patterns::Double (0., 1.)), + "A list of values that indicate how a given layer " + "should take into account the effects " + "of saturation on the temperature-dependence of the thermal " + "conductivity. This factor is multiplied with a saturation function " + "based on the theory of Roufosse and Klemens, 1974. A value of 1 " + "reproduces the formulation of Stackhouse et al. (2015), a value of " + "0 reproduces the formulation of Tosi et al., (2013). " + "Units: none."); + prm.declare_entry ("Maximum thermal conductivity", "1000", + Patterns::Double (0.), + "The maximum thermal conductivity that is allowed in the " + "model. Larger values will be cut off."); + } + + + + template + void + TosiStackhouse::parse_parameters (ParameterHandler &prm) + { + conductivity_transition_depths = Utilities::string_to_double + (Utilities::split_string_list(prm.get ("Thermal conductivity transition depths"))); + const unsigned int n_conductivity_layers = conductivity_transition_depths.size() + 1; + + AssertThrow (std::is_sorted(conductivity_transition_depths.begin(), conductivity_transition_depths.end()), + ExcMessage("The list of 'Thermal conductivity transition depths' must " + "be sorted such that the values increase monotonically.")); + + reference_thermal_conductivities = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Reference thermal conductivities"))), + n_conductivity_layers, + "Reference thermal conductivities"); + conductivity_pressure_dependencies = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Pressure dependencies of thermal conductivity"))), + n_conductivity_layers, + "Pressure dependencies of thermal conductivity"); + conductivity_reference_temperatures = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Reference temperatures for thermal conductivity"))), + n_conductivity_layers, + "Reference temperatures for thermal conductivity"); + conductivity_exponents = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Thermal conductivity exponents"))), + n_conductivity_layers, + "Thermal conductivity exponents"); + saturation_scaling = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Saturation prefactors"))), + n_conductivity_layers, + "Saturation prefactors"); + maximum_conductivity = prm.get_double ("Maximum thermal conductivity"); + } + } + } +} + +// explicit instantiations +namespace aspect +{ + namespace MaterialModel + { + namespace ThermalConductivity + { +#define INSTANTIATE(dim) \ + template class TosiStackhouse; + + ASPECT_INSTANTIATE(INSTANTIATE) + +#undef INSTANTIATE + } + } +}