diff --git a/include/aspect/material_model/entropy_model.h b/include/aspect/material_model/entropy_model.h index 0f5e7d52dfe..bc876011a45 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 @@ -125,31 +126,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 8ffd2ce9e31..c54bd21314c 100644 --- a/include/aspect/material_model/steinberger.h +++ b/include/aspect/material_model/steinberger.h @@ -23,6 +23,7 @@ #include #include +#include #include #include @@ -118,7 +119,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 following a + * pressure-temperature dependent approximation and the other * parameters are provided via lookup tables from the software PERPLEX. * * @ingroup MaterialModels @@ -204,19 +206,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 @@ -250,31 +239,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/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());