From 873090f6caf96b6efb573bf96d1f327e49726be8 Mon Sep 17 00:00:00 2001 From: Juliane Dannberg Date: Tue, 28 May 2019 16:59:41 -0700 Subject: [PATCH] add multicomponent equation of state --- .../aspect/material_model/dynamic_friction.h | 18 +-- .../multicomponent_incompressible.h | 118 +++++++++++++++ .../aspect/material_model/multicomponent.h | 18 +-- include/aspect/material_model/visco_plastic.h | 6 +- include/aspect/material_model/viscoelastic.h | 16 +-- .../material_model/viscoelastic_plastic.h | 16 +-- source/material_model/dynamic_friction.cc | 86 ++++------- .../multicomponent_incompressible.cc | 136 ++++++++++++++++++ source/material_model/multicomponent.cc | 97 ++++--------- source/material_model/visco_plastic.cc | 66 +++------ source/material_model/viscoelastic.cc | 71 +++------ source/material_model/viscoelastic_plastic.cc | 62 +++----- 12 files changed, 376 insertions(+), 334 deletions(-) create mode 100644 include/aspect/material_model/equation_of_state/multicomponent_incompressible.h create mode 100644 source/material_model/equation_of_state/multicomponent_incompressible.cc diff --git a/include/aspect/material_model/dynamic_friction.h b/include/aspect/material_model/dynamic_friction.h index c54165da63b..8383d68a27b 100644 --- a/include/aspect/material_model/dynamic_friction.h +++ b/include/aspect/material_model/dynamic_friction.h @@ -23,7 +23,7 @@ #include #include -#include +#include namespace aspect { @@ -137,10 +137,7 @@ namespace aspect MaterialUtilities::CompositionalAveragingOperation viscosity_averaging; - /** - * Vector for field densities, read from parameter file . - */ - std::vector densities; + EquationOfState::MulticomponentIncompressible equation_of_state; /** * The dynamic coefficient of friction @@ -174,21 +171,10 @@ namespace aspect double reference_strain_rate; double minimum_strain_rate; - /** - * Vector for field thermal expnsivities, read from parameter file. - */ - std::vector thermal_expansivities; - /** * Vector for field thermal conductivities, read from parameter file. */ std::vector thermal_conductivities; - - /** - * Vector for field specific heats, read from parameter file. - */ - std::vector specific_heats; - }; } diff --git a/include/aspect/material_model/equation_of_state/multicomponent_incompressible.h b/include/aspect/material_model/equation_of_state/multicomponent_incompressible.h new file mode 100644 index 00000000000..19f0b936faf --- /dev/null +++ b/include/aspect/material_model/equation_of_state/multicomponent_incompressible.h @@ -0,0 +1,118 @@ +/* + Copyright (C) 2011 - 2019 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_equation_of_state_multicomponent_incompressible_h +#define _aspect_material_model_equation_of_state_multicomponent_incompressible_h + +#include +#include +#include + + +namespace aspect +{ + namespace MaterialModel + { + namespace EquationOfState + { + using namespace dealii; + + /** + * An incompressible equation of state that is intended for use with multiple compositional + * fields. For each material property, the user supplies a comma delimited list of + * length N+1, where N is the number of compositional fields used in the computation. + * The first entry corresponds to the "background" (which is also why there are N+1 entries). + * + * If a single value is given, then all the compositional fields are given + * that value. Other lengths of lists are not allowed. For a given + * compositional field the material parameters are treated as constant, + * except density, which varies linearly with temperature according to the equation: + * + * $\rho(p,T,\mathfrak c) = \left(1-\alpha_i (T-T_0)\right) \rho_0(\mathfrak c_i).$ + * + * There is no pressure-dependence of the density. + */ + template + class MulticomponentIncompressible : public ::aspect::SimulatorAccess + { + public: + /** + * A function that computes the output of the equation of state @p out + * for all compositions, given the inputs in @p in and an index q that + * determines which entry of the vector of inputs is used. + */ + void evaluate(const MaterialModel::MaterialModelInputs &in, + const unsigned int q, + MaterialModel::EquationOfStateOutputs &out) const; + + /** + * Return whether the model is compressible or not. Incompressibility + * does not necessarily imply that the density is constant; rather, it + * may still depend on temperature or pressure. In the current + * context, compressibility means whether we should solve the continuity + * equation as $\nabla \cdot (\rho \mathbf u)=0$ (compressible Stokes) + * or as $\nabla \cdot \mathbf{u}=0$ (incompressible Stokes). + */ + bool is_compressible () const; + + /** + * Declare the parameters this class takes through input files. + * The optional parameter @p n_compositions determines the maximum + * number of compositions the equation of state is set up with, + * in other words, how many compositional fields influence the + * density. + */ + static + void + declare_parameters (ParameterHandler &prm, + const double default_thermal_expansion = 3.5e-5); + + /** + * Read the parameters this class declares from the parameter file. + * The optional parameter @p n_compositions determines the maximum + * number of compositions the equation of state is set up with, + * and should have the same value as the parameter with the same + * name in the declare_parameters() function. + */ + void + parse_parameters (ParameterHandler &prm); + + + private: + // Vector of reference densities $\rho_0$ with one entry per composition, + // used in the computation of the density. + std::vector densities; + + // The reference temperature $T_0$ used in the computation of the density. + // All component use the same reference temperature. + double reference_T; + + // Vector of constant thermal expansivities $\alpha$ with one entry per composition, + // used in the computation of the density. + std::vector thermal_expansivities; + + // Vector of specific heat capacities with one entry per composition. + std::vector specific_heats; + }; + } + } +} + +#endif diff --git a/include/aspect/material_model/multicomponent.h b/include/aspect/material_model/multicomponent.h index 41e1d4d6ec4..87959adb313 100644 --- a/include/aspect/material_model/multicomponent.h +++ b/include/aspect/material_model/multicomponent.h @@ -23,7 +23,8 @@ #include #include -#include +#include + namespace aspect { @@ -128,30 +129,17 @@ namespace aspect */ MaterialUtilities::CompositionalAveragingOperation viscosity_averaging; - /** - * Vector for field densities, read from parameter file . - */ - std::vector densities; - /** * Vector for field viscosities, read from parameter file. */ std::vector viscosities; - /** - * Vector for field thermal expansivities, read from parameter file. - */ - std::vector thermal_expansivities; - /** * Vector for field thermal conductivities, read from parameter file. */ std::vector thermal_conductivities; - /** - * Vector for field specific heats, read from parameter file. - */ - std::vector specific_heats; + EquationOfState::MulticomponentIncompressible equation_of_state; }; } diff --git a/include/aspect/material_model/visco_plastic.h b/include/aspect/material_model/visco_plastic.h index 19e2293818a..973e4ec514d 100644 --- a/include/aspect/material_model/visco_plastic.h +++ b/include/aspect/material_model/visco_plastic.h @@ -23,6 +23,7 @@ #include #include +#include #include @@ -178,10 +179,9 @@ namespace aspect double grain_size; - std::vector densities; - std::vector thermal_expansivities; std::vector thermal_diffusivities; - std::vector heat_capacities; + + EquationOfState::MulticomponentIncompressible equation_of_state; /** * Enumeration for selecting which viscosity averaging scheme to use. diff --git a/include/aspect/material_model/viscoelastic.h b/include/aspect/material_model/viscoelastic.h index 5eb5cc627c3..70258cc9a3e 100644 --- a/include/aspect/material_model/viscoelastic.h +++ b/include/aspect/material_model/viscoelastic.h @@ -23,6 +23,7 @@ #include #include +#include namespace aspect { @@ -249,31 +250,18 @@ namespace aspect const double dte) const; - /** - * Vector for field densities, read from parameter file. - */ - std::vector densities; + EquationOfState::MulticomponentIncompressible equation_of_state; /** * Vector for field viscosities, read from parameter file. */ std::vector viscosities; - /** - * Vector for field thermal expnsivities, read from parameter file. - */ - std::vector thermal_expansivities; - /** * Vector for field thermal conductivities, read from parameter file. */ std::vector thermal_conductivities; - /** - * Vector for field specific heats, read from parameter file. - */ - std::vector specific_heats; - /** * Vector for field elastic shear moduli, read from parameter file. */ diff --git a/include/aspect/material_model/viscoelastic_plastic.h b/include/aspect/material_model/viscoelastic_plastic.h index f3c349945a0..e8ffa9852d6 100644 --- a/include/aspect/material_model/viscoelastic_plastic.h +++ b/include/aspect/material_model/viscoelastic_plastic.h @@ -24,6 +24,7 @@ #include #include #include +#include namespace aspect { @@ -268,26 +269,13 @@ namespace aspect */ MaterialUtilities::CompositionalAveragingOperation viscosity_averaging; - /** - * Vector for field densities, read from parameter file. - */ - std::vector densities; - - /** - * Vector for field thermal expansivities, read from parameter file. - */ - std::vector thermal_expansivities; + EquationOfState::MulticomponentIncompressible equation_of_state; /** * Vector for field thermal conductivities, read from parameter file. */ std::vector thermal_conductivities; - /** - * Vector for field specific heats, read from parameter file. - */ - std::vector specific_heats; - /** * Vector for field linear (fixed) viscosities, read from parameter file. */ diff --git a/source/material_model/dynamic_friction.cc b/source/material_model/dynamic_friction.cc index 160ab92c067..cbc71ce3be2 100644 --- a/source/material_model/dynamic_friction.cc +++ b/source/material_model/dynamic_friction.cc @@ -78,51 +78,35 @@ namespace aspect evaluate(const MaterialModel::MaterialModelInputs &in, MaterialModel::MaterialModelOutputs &out) const { + EquationOfStateOutputs eos_outputs (this->n_compositional_fields()+1); + for (unsigned int i=0; i < in.position.size(); ++i) { - - const std::vector composition = in.composition[i]; - const std::vector volume_fractions = MaterialUtilities::compute_volume_fractions(composition); + const std::vector volume_fractions = MaterialUtilities::compute_volume_fractions(in.composition[i]); if (in.strain_rate.size() > 0) { const std::vector viscosities = compute_viscosities(in.pressure[i], in.strain_rate[i]); out.viscosities[i] = MaterialUtilities::average_value (volume_fractions, viscosities, viscosity_averaging); } - out.specific_heat[i] = MaterialUtilities::average_value (volume_fractions, specific_heats, MaterialUtilities::arithmetic); + equation_of_state.evaluate(in, i, eos_outputs); + + // The averaging is not strictly correct if thermal expansivities are different, since we are interpreting + // these compositions as volume fractions, but the error introduced should not be too bad. + out.densities[i] = MaterialUtilities::average_value (volume_fractions, eos_outputs.densities, MaterialUtilities::arithmetic); + out.thermal_expansion_coefficients[i] = MaterialUtilities::average_value (volume_fractions, eos_outputs.thermal_expansion_coefficients, MaterialUtilities::arithmetic); + out.specific_heat[i] = MaterialUtilities::average_value (volume_fractions, eos_outputs.specific_heat_capacities, MaterialUtilities::arithmetic); // Arithmetic averaging of thermal conductivities // This may not be strictly the most reasonable thing, but for most Earth materials we hope // that they do not vary so much that it is a big problem. out.thermal_conductivities[i] = MaterialUtilities::average_value (volume_fractions, thermal_conductivities, MaterialUtilities::arithmetic); - double density = 0.0; - for (unsigned int j=0; j < volume_fractions.size(); ++j) - { - // not strictly correct if thermal expansivities are different, since we are interpreting - // these compositions as volume fractions, but the error introduced should not be too bad. - const double temperature_factor = (1.0 - thermal_expansivities[j] * (in.temperature[i] - reference_T)); - density += volume_fractions[j] * densities[j] * temperature_factor; - } - out.densities[i] = density; - - - out.thermal_expansion_coefficients[i] = MaterialUtilities::average_value (volume_fractions, thermal_expansivities, MaterialUtilities::arithmetic); - + out.compressibilities[i] = MaterialUtilities::average_value (volume_fractions, eos_outputs.compressibilities, MaterialUtilities::arithmetic); + out.entropy_derivative_pressure[i] = MaterialUtilities::average_value (volume_fractions, eos_outputs.entropy_derivative_pressure, MaterialUtilities::arithmetic); + out.entropy_derivative_temperature[i] = MaterialUtilities::average_value (volume_fractions, eos_outputs.entropy_derivative_temperature, MaterialUtilities::arithmetic); - // Compressibility at the given positions. - // The compressibility is given as - // $\frac 1\rho \frac{\partial\rho}{\partial p}$. - // (here we use an incompressible medium) - out.compressibilities[i] = 0.0; - // Pressure derivative of entropy at the given positions. - out.entropy_derivative_pressure[i] = 0.0; - // Temperature derivative of entropy at the given positions. - out.entropy_derivative_temperature[i] = 0.0; - // Change in composition due to chemical reactions at the - // given positions. The term reaction_terms[i][c] is the - // change in compositional field c at point i. for (unsigned int c=0; c:: is_compressible () const { - return false; + return equation_of_state.is_compressible(); } template @@ -153,35 +137,22 @@ namespace aspect { prm.enter_subsection("Dynamic Friction"); { + EquationOfState::MulticomponentIncompressible::declare_parameters (prm, 4.e-5); + prm.declare_entry ("Reference temperature", "293", Patterns::Double (0), "The reference temperature $T_0$. Units: $\\text{K}$."); - prm.declare_entry ("Densities", "3300.", - Patterns::List(Patterns::Double(0)), - "List of densities for background mantle and compositional fields," - "for a total of N+1 values, where N is the number of compositional fields." - "If only one value is given, then all use the same value. Units: $kg / m^3$"); - prm.declare_entry ("Thermal expansivities", "4.e-5", - Patterns::List(Patterns::Double(0)), - "List of thermal expansivities for background mantle and compositional fields," - "for a total of N+1 values, where N is the number of compositional fields." - "If only one value is given, then all use the same value. Units: $1/K$"); - prm.declare_entry ("Specific heats", "1250.", - Patterns::List(Patterns::Double(0)), - "List of specific heats $C_p$ for background mantle and compositional fields," - "for a total of N+1 values, where N is the number of compositional fields." - "If only one value is given, then all use the same value. Units: $J /kg /K$"); prm.declare_entry ("Thermal conductivities", "4.7", Patterns::List(Patterns::Double(0)), "List of thermal conductivities for background mantle and compositional fields," "for a total of N+1 values, where N is the number of compositional fields." "If only one value is given, then all use the same value. Units: $W/m/K$."); - prm.declare_entry("Viscosity averaging scheme", "harmonic", - Patterns::Selection("arithmetic|harmonic|geometric|maximum composition"), - "When more than one compositional field is present at a point " - "with different viscosities, we need to come up with an average " - "viscosity at that point. Select a weighted harmonic, arithmetic, " - "geometric, or maximum composition."); + prm.declare_entry ("Viscosity averaging scheme", "harmonic", + Patterns::Selection("arithmetic|harmonic|geometric|maximum composition"), + "When more than one compositional field is present at a point " + "with different viscosities, we need to come up with an average " + "viscosity at that point. Select a weighted harmonic, arithmetic, " + "geometric, or maximum composition."); prm.enter_subsection("Viscosities"); { prm.declare_entry ("Minimum viscosity", "1e19", @@ -237,24 +208,17 @@ namespace aspect { prm.enter_subsection("Dynamic Friction"); { + equation_of_state.initialize_simulator (this->get_simulator()); + equation_of_state.parse_parameters (prm); + reference_T = prm.get_double ("Reference temperature"); viscosity_averaging = MaterialUtilities::parse_compositional_averaging_operation ("Viscosity averaging scheme", prm); - // Parse DynamicFriction properties - densities = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Densities"))), - n_fields, - "Densities"); thermal_conductivities = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Thermal conductivities"))), n_fields, "Thermal conductivities"); - thermal_expansivities = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Thermal expansivities"))), - n_fields, - "Thermal expansivities"); - specific_heats = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Specific heats"))), - n_fields, - "Specific heats"); prm.enter_subsection("Viscosities"); { minimum_viscosity = prm.get_double ("Minimum viscosity"); diff --git a/source/material_model/equation_of_state/multicomponent_incompressible.cc b/source/material_model/equation_of_state/multicomponent_incompressible.cc new file mode 100644 index 00000000000..a999e0254bb --- /dev/null +++ b/source/material_model/equation_of_state/multicomponent_incompressible.cc @@ -0,0 +1,136 @@ +/* + Copyright (C) 2011 - 2019 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 EquationOfState + { + template + void + MulticomponentIncompressible:: + evaluate(const MaterialModel::MaterialModelInputs &in, + const unsigned int q, + MaterialModel::EquationOfStateOutputs &out) const + { + for (unsigned int c=0; c < out.densities.size(); ++c) + { + out.densities[c] = densities[c] * (1 - thermal_expansivities[c] * (in.temperature[q] - reference_T)); + out.thermal_expansion_coefficients[c] = thermal_expansivities[c]; + out.specific_heat_capacities[c] = specific_heats[c]; + out.compressibilities[c] = 0.0; + out.entropy_derivative_pressure[c] = 0.0; + out.entropy_derivative_temperature[c] = 0.0; + } + } + + + + template + bool + MulticomponentIncompressible:: + is_compressible () const + { + return false; + } + + + + template + void + MulticomponentIncompressible::declare_parameters (ParameterHandler &prm, + const double default_thermal_expansion) + { + prm.declare_entry ("Reference temperature", "293", + Patterns::Double (0), + "The reference temperature $T_0$. Units: $\\text{K}$."); + prm.declare_entry ("Densities", "3300.", + Patterns::Anything(), + "List of densities for background mantle and compositional fields," + "for a total of N+1 values, where N is the number of compositional fields." + "If only one value is given, then all use the same value. Units: $kg / m^3$"); + prm.declare_entry ("Thermal expansivities", std::to_string(default_thermal_expansion), + Patterns::Anything(), + "List of thermal expansivities for background mantle and compositional fields," + "for a total of N+1 values, where N is the number of compositional fields." + "If only one value is given, then all use the same value. Units: $1/K$"); + prm.declare_entry ("Heat capacities", "1250.", + Patterns::Anything(), + "List of specific heats $C_p$ for background mantle and compositional fields," + "for a total of N+1 values, where N is the number of compositional fields." + "If only one value is given, then all use the same value. Units: $J /kg /K$"); + prm.declare_alias ("Heat capacities", "Specific heats"); + } + + + + template + void + MulticomponentIncompressible::parse_parameters (ParameterHandler &prm) + { + reference_T = prm.get_double ("Reference temperature"); + + // Establish that a background field is required here + const bool has_background_field = true; + + // Retrieve the list of composition names + const std::vector list_of_composition_names = this->introspection().get_composition_names(); + + // Parse multicomponent properties + densities = Utilities::parse_map_to_double_array (prm.get("Densities"), + list_of_composition_names, + has_background_field, + "Densities"); + + thermal_expansivities = Utilities::parse_map_to_double_array (prm.get("Thermal expansivities"), + list_of_composition_names, + has_background_field, + "Thermal expansivities"); + + specific_heats = Utilities::parse_map_to_double_array (prm.get("Heat capacities"), + list_of_composition_names, + has_background_field, + "Specific heats"); + } + + } + } +} + + +// explicit instantiations +namespace aspect +{ + namespace MaterialModel + { + namespace EquationOfState + { +#define INSTANTIATE(dim) \ + template class MulticomponentIncompressible; + ASPECT_INSTANTIATE(INSTANTIATE) + } + } +} diff --git a/source/material_model/multicomponent.cc b/source/material_model/multicomponent.cc index 553db1c81a9..22b54cf3b50 100644 --- a/source/material_model/multicomponent.cc +++ b/source/material_model/multicomponent.cc @@ -20,7 +20,6 @@ #include -#include #include #include @@ -36,47 +35,30 @@ namespace aspect evaluate(const MaterialModel::MaterialModelInputs &in, MaterialModel::MaterialModelOutputs &out) const { + EquationOfStateOutputs eos_outputs (this->n_compositional_fields()+1); + for (unsigned int i=0; i < in.temperature.size(); ++i) { - const double temperature = in.temperature[i]; - const std::vector composition = in.composition[i]; - const std::vector volume_fractions = MaterialUtilities::compute_volume_fractions(composition); + const std::vector volume_fractions = MaterialUtilities::compute_volume_fractions(in.composition[i]); - out.viscosities[i] = MaterialUtilities::average_value (volume_fractions, viscosities, viscosity_averaging); - out.specific_heat[i] = MaterialUtilities::average_value (volume_fractions, specific_heats, MaterialUtilities::arithmetic); + equation_of_state.evaluate(in, i, eos_outputs); + out.viscosities[i] = MaterialUtilities::average_value (volume_fractions, viscosities, viscosity_averaging); + out.specific_heat[i] = MaterialUtilities::average_value (volume_fractions, eos_outputs.specific_heat_capacities, MaterialUtilities::arithmetic); + out.thermal_expansion_coefficients[i] = MaterialUtilities::average_value (volume_fractions, eos_outputs.thermal_expansion_coefficients, MaterialUtilities::arithmetic); + out.compressibilities[i] = MaterialUtilities::average_value (volume_fractions, eos_outputs.compressibilities, MaterialUtilities::arithmetic); + out.entropy_derivative_pressure[i] = MaterialUtilities::average_value (volume_fractions, eos_outputs.entropy_derivative_pressure, MaterialUtilities::arithmetic); + out.entropy_derivative_temperature[i] = MaterialUtilities::average_value (volume_fractions, eos_outputs.entropy_derivative_temperature, MaterialUtilities::arithmetic); // Arithmetic averaging of thermal conductivities // This may not be strictly the most reasonable thing, but for most Earth materials we hope // that they do not vary so much that it is a big problem. out.thermal_conductivities[i] = MaterialUtilities::average_value (volume_fractions, thermal_conductivities, MaterialUtilities::arithmetic); - double density = 0.0; - for (unsigned int j=0; j < volume_fractions.size(); ++j) - { - // not strictly correct if thermal expansivities are different, since we are interpreting - // these compositions as volume fractions, but the error introduced should not be too bad. - const double temperature_factor= (1.0 - thermal_expansivities[j] * (temperature - reference_T)); - density += volume_fractions[j] * densities[j] * temperature_factor; - } - out.densities[i] = density; - - - out.thermal_expansion_coefficients[i] = MaterialUtilities::average_value (volume_fractions, thermal_expansivities, MaterialUtilities::arithmetic); - - - // Compressibility at the given positions. - // The compressibility is given as - // $\frac 1\rho \frac{\partial\rho}{\partial p}$. - // (here we use an incompressible medium) - out.compressibilities[i] = 0.0; - // Pressure derivative of entropy at the given positions. - out.entropy_derivative_pressure[i] = 0.0; - // Temperature derivative of entropy at the given positions. - out.entropy_derivative_temperature[i] = 0.0; - // Change in composition due to chemical reactions at the - // given positions. The term reaction_terms[i][c] is the - // change in compositional field c at point i. + // not strictly correct if thermal expansivities are different, since we are interpreting + // these compositions as volume fractions, but the error introduced should not be too bad. + out.densities[i] = MaterialUtilities::average_value (volume_fractions, eos_outputs.densities, MaterialUtilities::arithmetic); + for (unsigned int c=0; c:: is_compressible () const { - return false; + return equation_of_state.is_compressible(); } template @@ -107,40 +89,27 @@ namespace aspect { prm.enter_subsection("Multicomponent"); { + EquationOfState::MulticomponentIncompressible::declare_parameters (prm, 4.e-5); + prm.declare_entry ("Reference temperature", "293", Patterns::Double (0), "The reference temperature $T_0$. Units: $\\text{K}$."); - prm.declare_entry ("Densities", "3300.", - Patterns::Anything(), - "List of densities for background mantle and compositional fields," - "for a total of N+1 values, where N is the number of compositional fields." - "If only one value is given, then all use the same value. Units: $kg / m^3$"); prm.declare_entry ("Viscosities", "1.e21", Patterns::Anything(), "List of viscosities for background mantle and compositional fields," "for a total of N+1 values, where N is the number of compositional fields." "If only one value is given, then all use the same value. Units: $Pa \\, s$"); - prm.declare_entry ("Thermal expansivities", "4.e-5", - Patterns::Anything(), - "List of thermal expansivities for background mantle and compositional fields," - "for a total of N+1 values, where N is the number of compositional fields." - "If only one value is given, then all use the same value. Units: $1/K$"); - prm.declare_entry ("Specific heats", "1250.", - Patterns::Anything(), - "List of specific heats $C_p$ for background mantle and compositional fields," - "for a total of N+1 values, where N is the number of compositional fields." - "If only one value is given, then all use the same value. Units: $J /kg /K$"); prm.declare_entry ("Thermal conductivities", "4.7", Patterns::Anything(), "List of thermal conductivities for background mantle and compositional fields," "for a total of N+1 values, where N is the number of compositional fields." "If only one value is given, then all use the same value. Units: $W/m/K$."); - prm.declare_entry("Viscosity averaging scheme", "harmonic", - Patterns::Selection("arithmetic|harmonic|geometric|maximum composition"), - "When more than one compositional field is present at a point " - "with different viscosities, we need to come up with an average " - "viscosity at that point. Select a weighted harmonic, arithmetic, " - "geometric, or maximum composition."); + prm.declare_entry ("Viscosity averaging scheme", "harmonic", + Patterns::Selection("arithmetic|harmonic|geometric|maximum composition"), + "When more than one compositional field is present at a point " + "with different viscosities, we need to come up with an average " + "viscosity at that point. Select a weighted harmonic, arithmetic, " + "geometric, or maximum composition."); } prm.leave_subsection(); } @@ -155,6 +124,9 @@ namespace aspect { prm.enter_subsection("Multicomponent"); { + equation_of_state.initialize_simulator (this->get_simulator()); + equation_of_state.parse_parameters (prm); + reference_T = prm.get_double ("Reference temperature"); viscosity_averaging = MaterialUtilities::parse_compositional_averaging_operation ("Viscosity averaging scheme", @@ -166,12 +138,6 @@ namespace aspect // Retrieve the list of composition names const std::vector list_of_composition_names = this->introspection().get_composition_names(); - // Parse multicomponent properties - densities = Utilities::parse_map_to_double_array (prm.get("Densities"), - list_of_composition_names, - has_background_field, - "Densities"); - viscosities = Utilities::parse_map_to_double_array (prm.get("Viscosities"), list_of_composition_names, has_background_field, @@ -181,17 +147,6 @@ namespace aspect list_of_composition_names, has_background_field, "Thermal conductivities"); - - thermal_expansivities = Utilities::parse_map_to_double_array (prm.get("Thermal expansivities"), - list_of_composition_names, - has_background_field, - "Thermal expansivities"); - - specific_heats = Utilities::parse_map_to_double_array (prm.get("Specific heats"), - list_of_composition_names, - has_background_field, - "Specific heats"); - } prm.leave_subsection(); } diff --git a/source/material_model/visco_plastic.cc b/source/material_model/visco_plastic.cc index 302e5c6cb2e..54995525239 100644 --- a/source/material_model/visco_plastic.cc +++ b/source/material_model/visco_plastic.cc @@ -584,27 +584,25 @@ namespace aspect // Store which components do not represent volumetric compositions (e.g. strain components). const ComponentMask volumetric_compositions = get_volumetric_composition_mask(); + EquationOfStateOutputs eos_outputs (this->n_compositional_fields()+1); + // Loop through all requested points for (unsigned int i=0; i < in.temperature.size(); ++i) { + equation_of_state.evaluate(in, i, eos_outputs); + // First compute the equation of state variables and thermodynamic properties - out.densities[i] = 0.0; - out.thermal_expansion_coefficients[i] = 0.0; - out.specific_heat[i] = 0.0; + const std::vector volume_fractions = MaterialUtilities::compute_volume_fractions(in.composition[i], volumetric_compositions); + + // not strictly correct if thermal expansivities are different, since we are interpreting + // these compositions as volume fractions, but the error introduced should not be too bad. + out.densities[i] = MaterialUtilities::average_value (volume_fractions, eos_outputs.densities, MaterialUtilities::arithmetic); + out.thermal_expansion_coefficients[i] = MaterialUtilities::average_value (volume_fractions, eos_outputs.thermal_expansion_coefficients, MaterialUtilities::arithmetic); + out.specific_heat[i] = MaterialUtilities::average_value (volume_fractions, eos_outputs.specific_heat_capacities, MaterialUtilities::arithmetic); double thermal_diffusivity = 0.0; - const std::vector volume_fractions = MaterialUtilities::compute_volume_fractions(in.composition[i], volumetric_compositions); for (unsigned int j=0; j < volume_fractions.size(); ++j) - { - // not strictly correct if thermal expansivities are different, since we are interpreting - // these compositions as volume fractions, but the error introduced should not be too bad. - const double temperature_factor = (1.0 - thermal_expansivities[j] * (in.temperature[i] - reference_T)); - - out.densities[i] += volume_fractions[j] * densities[j] * temperature_factor; - out.thermal_expansion_coefficients[i] += volume_fractions[j] * thermal_expansivities[j]; - out.specific_heat[i] += volume_fractions[j] * heat_capacities[j]; - thermal_diffusivity += volume_fractions[j] * thermal_diffusivities[j]; - } + thermal_diffusivity += volume_fractions[j] * thermal_diffusivities[j]; // Thermal conductivity at the given positions. If the temperature equation uses // the reference density profile formulation, use the reference density to @@ -618,9 +616,9 @@ namespace aspect else out.thermal_conductivities[i] = thermal_diffusivity * out.specific_heat[i] * out.densities[i]; - out.compressibilities[i] = 0.0; - out.entropy_derivative_pressure[i] = 0.0; - out.entropy_derivative_temperature[i] = 0.0; + out.compressibilities[i] = MaterialUtilities::average_value (volume_fractions, eos_outputs.compressibilities, MaterialUtilities::arithmetic); + out.entropy_derivative_pressure[i] = MaterialUtilities::average_value (volume_fractions, eos_outputs.entropy_derivative_pressure, MaterialUtilities::arithmetic); + out.entropy_derivative_temperature[i] = MaterialUtilities::average_value (volume_fractions, eos_outputs.entropy_derivative_temperature, MaterialUtilities::arithmetic); // Compute the effective viscosity if requested and retrieve whether the material is plastically yielding bool plastic_yielding = false; @@ -703,7 +701,7 @@ namespace aspect ViscoPlastic:: is_compressible () const { - return false; + return equation_of_state.is_compressible(); } template @@ -721,6 +719,8 @@ namespace aspect { prm.enter_subsection ("Visco Plastic"); { + EquationOfState::MulticomponentIncompressible::declare_parameters (prm); + // Reference and minimum/maximum values prm.declare_entry ("Reference temperature", "293", Patterns::Double(0), "For calculating density by thermal expansivity. Units: $\\text{K}$"); @@ -760,21 +760,6 @@ namespace aspect "List of thermal diffusivities, for background material and compositional fields, " "for a total of N+1 values, where N is the number of compositional fields. " "If only one value is given, then all use the same value. Units: $m^2/s$"); - prm.declare_entry ("Heat capacities", "1.25e3", - Patterns::List(Patterns::Double(0)), - "List of heat capacities $C_p$, for background material and compositional fields, " - "for a total of N+1 values, where N is the number of compositional fields. " - "If only one value is given, then all use the same value. Units: $J/kg/K$"); - prm.declare_entry ("Densities", "3300.", - Patterns::List(Patterns::Double(0)), - "List of densities, $\\rho$, for background material and compositional fields, " - "for a total of N+1 values, where N is the number of compositional fields. " - "If only one value is given, then all use the same value. Units: $kg / m^3$"); - prm.declare_entry ("Thermal expansivities", "3.5e-5", - Patterns::List(Patterns::Double(0)), - "List of thermal expansivities for background material and compositional fields, " - "for a total of N+1 values, where N is the number of compositional fields. " - "If only one value is given, then all use the same value. Units: $1 / K$"); // Strain weakening parameters prm.declare_entry ("Strain weakening mechanism", "default", @@ -1042,6 +1027,10 @@ namespace aspect { prm.enter_subsection ("Visco Plastic"); { + // Equation of state parameters + equation_of_state.initialize_simulator (this->get_simulator()); + equation_of_state.parse_parameters (prm); + // Reference and minimum/maximum values reference_T = prm.get_double("Reference temperature"); min_strain_rate = prm.get_double("Minimum strain rate"); @@ -1050,22 +1039,11 @@ namespace aspect max_visc = prm.get_double ("Maximum viscosity"); ref_visc = prm.get_double ("Reference viscosity"); - // Equation of state parameters thermal_diffusivities = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Thermal diffusivities"))), n_fields, "Thermal diffusivities"); - heat_capacities = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Heat capacities"))), - n_fields, - "Heat capacities"); - // ---- Compositional parameters grain_size = prm.get_double("Grain size"); - densities = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Densities"))), - n_fields, - "Densities"); - thermal_expansivities = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Thermal expansivities"))), - n_fields, - "Thermal expansivities"); // Strain weakening parameters if (prm.get ("Strain weakening mechanism") == "none") diff --git a/source/material_model/viscoelastic.cc b/source/material_model/viscoelastic.cc index 338b28e8918..ea21e8d15b4 100644 --- a/source/material_model/viscoelastic.cc +++ b/source/material_model/viscoelastic.cc @@ -109,6 +109,8 @@ namespace aspect MaterialModel::ElasticOutputs *force_out = out.template get_additional_output >(); + EquationOfStateOutputs eos_outputs (this->n_compositional_fields()+1); + // Store which components to exclude during volume fraction computation. ComponentMask composition_mask(this->n_compositional_fields(),true); // assign compositional fields associated with viscoelastic stress a value of 0 @@ -134,41 +136,26 @@ namespace aspect for (unsigned int i=0; i < in.temperature.size(); ++i) { - const double temperature = in.temperature[i]; const std::vector composition = in.composition[i]; const std::vector volume_fractions = MaterialUtilities::compute_volume_fractions(composition, composition_mask); - out.specific_heat[i] = MaterialUtilities::average_value(volume_fractions, specific_heats, MaterialUtilities::arithmetic); + equation_of_state.evaluate(in, i, eos_outputs); // Arithmetic averaging of thermal conductivities // This may not be strictly the most reasonable thing, but for most Earth materials we hope // that they do not vary so much that it is a big problem. out.thermal_conductivities[i] = MaterialUtilities::average_value(volume_fractions, thermal_conductivities, MaterialUtilities::arithmetic); - double density = 0.0; - for (unsigned int j=0; j < volume_fractions.size(); ++j) - { - // not strictly correct if thermal expansivities are different, since we are interpreting - // these compositions as volume fractions, but the error introduced should not be too bad. - const double temperature_factor= (1.0 - thermal_expansivities[j] * (temperature - reference_T)); - density += volume_fractions[j] * densities[j] * temperature_factor; - } - out.densities[i] = density; - - out.thermal_expansion_coefficients[i] = MaterialUtilities::average_value(volume_fractions, thermal_expansivities, MaterialUtilities::arithmetic); - - // Compressibility at the given positions. - // The compressibility is given as - // $\frac 1\rho \frac{\partial\rho}{\partial p}$. - // (here we use an incompressible medium) - out.compressibilities[i] = 0.0; - // Pressure derivative of entropy at the given positions. - out.entropy_derivative_pressure[i] = 0.0; - // Temperature derivative of entropy at the given positions. - out.entropy_derivative_temperature[i] = 0.0; - // Change in composition due to chemical reactions at the - // given positions. The term reaction_terms[i][c] is the - // change in compositional field c at point i. + // not strictly correct if thermal expansivities are different, since we are interpreting + // these compositions as volume fractions, but the error introduced should not be too bad. + out.densities[i] = MaterialUtilities::average_value (volume_fractions, eos_outputs.densities, MaterialUtilities::arithmetic); + out.thermal_expansion_coefficients[i] = MaterialUtilities::average_value (volume_fractions, eos_outputs.thermal_expansion_coefficients, MaterialUtilities::arithmetic); + out.specific_heat[i] = MaterialUtilities::average_value (volume_fractions, eos_outputs.specific_heat_capacities, MaterialUtilities::arithmetic); + + out.compressibilities[i] = MaterialUtilities::average_value (volume_fractions, eos_outputs.compressibilities, MaterialUtilities::arithmetic); + out.entropy_derivative_pressure[i] = MaterialUtilities::average_value (volume_fractions, eos_outputs.entropy_derivative_pressure, MaterialUtilities::arithmetic); + out.entropy_derivative_temperature[i] = MaterialUtilities::average_value (volume_fractions, eos_outputs.entropy_derivative_temperature, MaterialUtilities::arithmetic); + for (unsigned int c=0; c:: is_compressible () const { - return false; + return equation_of_state.is_compressible(); } template @@ -300,29 +287,16 @@ namespace aspect { prm.enter_subsection("Viscoelastic"); { + EquationOfState::MulticomponentIncompressible::declare_parameters (prm); + prm.declare_entry ("Reference temperature", "293", Patterns::Double (0), "The reference temperature $T_0$. Units: $\\text{K}$."); - prm.declare_entry ("Densities", "3300.", - Patterns::List(Patterns::Double(0)), - "List of densities for background mantle and compositional fields, " - "for a total of N+1 values, where N is the number of compositional fields. " - "If only one value is given, then all use the same value. Units: $kg / m^3$"); prm.declare_entry ("Viscosities", "1.e21", Patterns::List(Patterns::Double(0)), "List of viscosities for background mantle and compositional fields, " "for a total of N+1 values, where N is the number of compositional fields. " "If only one value is given, then all use the same value. Units: $Pa s$"); - prm.declare_entry ("Thermal expansivities", "4.e-5", - Patterns::List(Patterns::Double(0)), - "List of thermal expansivities for background mantle and compositional fields, " - "for a total of N+1 values, where N is the number of compositional fields. " - "If only one value is given, then all use the same value. Units: $1/K$"); - prm.declare_entry ("Specific heats", "1250.", - Patterns::List(Patterns::Double(0)), - "List of specific heats $C_p$ for background mantle and compositional fields, " - "for a total of N+1 values, where N is the number of compositional fields. " - "If only one value is given, then all use the same value. Units: $J /kg /K$"); prm.declare_entry ("Thermal conductivities", "4.7", Patterns::List(Patterns::Double(0)), "List of thermal conductivities for background mantle and compositional fields, " @@ -379,27 +353,22 @@ namespace aspect { prm.enter_subsection("Viscoelastic"); { + // Equation of state parameters + equation_of_state.initialize_simulator (this->get_simulator()); + equation_of_state.parse_parameters (prm); + reference_T = prm.get_double ("Reference temperature"); viscosity_averaging = MaterialUtilities::parse_compositional_averaging_operation ("Viscosity averaging scheme", prm); // Parse viscoelastic properties - densities = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Densities"))), - n_fields, - "Densities"); viscosities = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Viscosities"))), n_fields, "Viscosities"); thermal_conductivities = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Thermal conductivities"))), n_fields, "Thermal conductivities"); - thermal_expansivities = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Thermal expansivities"))), - n_fields, - "Thermal expansivities"); - specific_heats = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Specific heats"))), - n_fields, - "Specific heats"); elastic_shear_moduli = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Elastic shear moduli"))), n_fields, "Elastic shear moduli"); diff --git a/source/material_model/viscoelastic_plastic.cc b/source/material_model/viscoelastic_plastic.cc index 65792ae11b1..be2aff96bd2 100644 --- a/source/material_model/viscoelastic_plastic.cc +++ b/source/material_model/viscoelastic_plastic.cc @@ -49,6 +49,8 @@ namespace aspect MaterialModel::ElasticOutputs *force_out = out.template get_additional_output >(); + EquationOfStateOutputs eos_outputs (this->n_compositional_fields()+1); + // The elastic time step (dte) is equal to the numerical time step if the time step number // is greater than 0 and the parameter 'use_fixed_elastic_time_step' is set to false. // On the first (0) time step the elastic time step is always equal to the value @@ -62,37 +64,28 @@ namespace aspect for (unsigned int i=0; i < in.temperature.size(); ++i) { - - const double temperature = in.temperature[i]; const std::vector composition = in.composition[i]; const double pressure = in.pressure[i]; const SymmetricTensor<2,dim> strain_rate = in.strain_rate[i]; const std::vector volume_fractions = MaterialUtilities::compute_volume_fractions(composition, composition_mask); + equation_of_state.evaluate(in, i, eos_outputs); + // Arithmetic averaging of specific heat. // This may not be strictly the most reasonable thing, but for most Earth materials we hope // that they do not vary so much that it is a big problem. This statement also applies to // the arithmetic averaging of density and thermal conductivity below. - out.specific_heat[i] = MaterialUtilities::average_value(volume_fractions, specific_heats, MaterialUtilities::arithmetic); + out.densities[i] = MaterialUtilities::average_value (volume_fractions, eos_outputs.densities, MaterialUtilities::arithmetic); + out.thermal_expansion_coefficients[i] = MaterialUtilities::average_value (volume_fractions, eos_outputs.thermal_expansion_coefficients, MaterialUtilities::arithmetic); + out.specific_heat[i] = MaterialUtilities::average_value (volume_fractions, eos_outputs.specific_heat_capacities, MaterialUtilities::arithmetic); out.thermal_conductivities[i] = MaterialUtilities::average_value(volume_fractions, thermal_conductivities, MaterialUtilities::arithmetic); - double density = 0.0; - for (unsigned int j=0; j < volume_fractions.size(); ++j) - { - // not strictly correct if thermal expansivities are different, since we are interpreting - // these compositions as volume fractions, but the error introduced should not be too bad. - const double temperature_factor = (1.0 - thermal_expansivities[j] * (temperature - reference_temperature)); - density += volume_fractions[j] * densities[j] * temperature_factor; - } - out.densities[i] = density; - - out.thermal_expansion_coefficients[i] = MaterialUtilities::average_value(volume_fractions, thermal_expansivities, MaterialUtilities::arithmetic); - // Set properties that are not relevant for this material model to 0 - out.compressibilities[i] = 0.0; - out.entropy_derivative_pressure[i] = 0.0; - out.entropy_derivative_temperature[i] = 0.0; + out.compressibilities[i] = MaterialUtilities::average_value (volume_fractions, eos_outputs.compressibilities, MaterialUtilities::arithmetic); + out.entropy_derivative_pressure[i] = MaterialUtilities::average_value (volume_fractions, eos_outputs.entropy_derivative_pressure, MaterialUtilities::arithmetic); + out.entropy_derivative_temperature[i] = MaterialUtilities::average_value (volume_fractions, eos_outputs.entropy_derivative_temperature, MaterialUtilities::arithmetic); + for (unsigned int c=0; c:: is_compressible () const { - return false; + return equation_of_state.is_compressible(); } template @@ -269,6 +262,8 @@ namespace aspect { prm.enter_subsection("Viscoelastic Plastic"); { + // Equation of state parameters + EquationOfState::MulticomponentIncompressible::declare_parameters (prm); // Reference and minimum/maximum values prm.declare_entry ("Reference temperature", "293", @@ -304,22 +299,6 @@ namespace aspect "\n\n" "Units: $Pa \\, s$"); - // Equation of state parameters - prm.declare_entry ("Densities", "3300.", - Patterns::List(Patterns::Double(0)), - "List of densities for background material and compositional fields, " - "for a total of N+1 values, where N is the number of compositional fields. " - "If only one value is given, then all use the same value. Units: $kg / m^3$"); - prm.declare_entry ("Thermal expansivities", "4.e-5", - Patterns::List(Patterns::Double(0)), - "List of thermal expansivities for background material and compositional fields, " - "for a total of N+1 values, where N is the number of compositional fields. " - "If only one value is given, then all use the same value. Units: $1/K$"); - prm.declare_entry ("Specific heats", "1250.", - Patterns::List(Patterns::Double(0)), - "List of specific heats $C_p$ for background material and compositional fields, " - "for a total of N+1 values, where N is the number of compositional fields. " - "If only one value is given, then all use the same value. Units: $J /kg /K$"); prm.declare_entry ("Thermal conductivities", "4.7", Patterns::List(Patterns::Double(0)), "List of thermal conductivities for background material and compositional fields, " @@ -397,10 +376,12 @@ namespace aspect { prm.enter_subsection("Viscoelastic Plastic"); { - AssertThrow(this->get_parameters().enable_elasticity == true, ExcMessage ("Material model Viscoelastic only works if 'Enable elasticity' is set to true")); + equation_of_state.initialize_simulator (this->get_simulator()); + equation_of_state.parse_parameters (prm); + reference_temperature = prm.get_double ("Reference temperature"); minimum_strain_rate = prm.get_double("Minimum strain rate"); reference_strain_rate = prm.get_double("Reference strain rate"); @@ -424,21 +405,12 @@ namespace aspect "Cohesions"); // Parse additional material properties - densities = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Densities"))), - n_fields, - "Densities"); linear_viscosities = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Linear viscosities"))), n_fields, "Viscosities"); thermal_conductivities = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Thermal conductivities"))), n_fields, "Thermal conductivities"); - thermal_expansivities = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Thermal expansivities"))), - n_fields, - "Thermal expansivities"); - specific_heats = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Specific heats"))), - n_fields, - "Specific heats"); elastic_shear_moduli = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Elastic shear moduli"))), n_fields, "Elastic shear moduli");