diff --git a/contrib/utilities/update_scripts/prm_files/move_tian_reaction_model.py b/contrib/utilities/update_scripts/prm_files/move_tian_reaction_model.py new file mode 100644 index 00000000000..d05ce3c07eb --- /dev/null +++ b/contrib/utilities/update_scripts/prm_files/move_tian_reaction_model.py @@ -0,0 +1,66 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +""" This script reformats the given .prm files to move the Tian 2019 reaction +mode parameters to the correct subsection. +""" + +import sys +import os +import re +import argparse + +__author__ = 'The authors of the ASPECT code' +__copyright__ = 'Copyright 2024, ASPECT' +__license__ = 'GNU GPL 2 or later' + +# Add the ASPECT root directory to the path so we can import from the aspect_data module +sys.path.append(os.path.join(os.path.dirname(__file__), '..', '..', '..')) +import python.scripts.aspect_input as aspect + + + +def create_tian2019_reaction_subsection(parameters): + """ Move the Tian 2019 reaction parameters to their own subsection. """ + + + parameters_to_move = ["Maximum weight percent water in sediment", \ + "Maximum weight percent water in MORB", \ + "Maximum weight percent water in gabbro", \ + "Maximum weight percent water in peridotite"] + + # Collect existing parameters and delete old entries + reactive_fluid_params = dict({}) + if "Reactive Fluid Transport Model" in parameters["Material model"]["value"]: + reactive_fluid_params = parameters["Material model"]["value"]["Reactive Fluid Transport Model"] + for param in parameters_to_move: + if "Tian 2019 model" not in reactive_fluid_params["value"]: + reactive_fluid_params["value"]["Tian 2019 model"] = {"comment": "", "value" : dict({}), "type": "subsection"} + + if param in reactive_fluid_params["value"]: + reactive_fluid_params["value"]["Tian 2019 model"]["value"][param] = reactive_fluid_params["value"][param] + del reactive_fluid_params["value"][param] + + return parameters + + + +def main(input_file, output_file): + """Echo the input arguments to standard output""" + parameters = aspect.read_parameter_file(input_file) + + parameters = create_tian2019_reaction_subsection(parameters) + + aspect.write_parameter_file(parameters, output_file) + + + +if __name__ == '__main__': + parser = argparse.ArgumentParser( + prog='ASPECT .prm file reformatter', + description='Reformats ASPECT .prm files to follow our general formatting guidelines. See the documentation of this script for details.') + parser.add_argument('input_file', type=str, help='The .prm file to reformat') + parser.add_argument('output_file', type=str, help='The .prm file to write the reformatted file to') + args = parser.parse_args() + + sys.exit(main(args.input_file, args.output_file)) diff --git a/cookbooks/tian_parameterization_kinematic_slab/coupled-two-phase-tian-parameterization-kinematic-slab.prm b/cookbooks/tian_parameterization_kinematic_slab/coupled-two-phase-tian-parameterization-kinematic-slab.prm index 2ff209c9701..bcb5c4a746d 100644 --- a/cookbooks/tian_parameterization_kinematic_slab/coupled-two-phase-tian-parameterization-kinematic-slab.prm +++ b/cookbooks/tian_parameterization_kinematic_slab/coupled-two-phase-tian-parameterization-kinematic-slab.prm @@ -229,10 +229,12 @@ subsection Material model # values to encourage water to hydrate the overlying mantle. The polynomials defined # in Tian et al., 2019 also reach very large values at low P-T conditions, and so limiting # the weight percent to reasonable values is recommended. - set Maximum weight percent water in peridotite = 2 - set Maximum weight percent water in gabbro = 1 - set Maximum weight percent water in MORB = 2 - set Maximum weight percent water in sediment = 3 + subsection Tian 2019 model + set Maximum weight percent water in peridotite = 2 + set Maximum weight percent water in gabbro = 1 + set Maximum weight percent water in MORB = 2 + set Maximum weight percent water in sediment = 3 + end end subsection Visco Plastic diff --git a/doc/modules/changes/20250206_danieldouglas92 b/doc/modules/changes/20250206_danieldouglas92 new file mode 100644 index 00000000000..1949fdf1d9b --- /dev/null +++ b/doc/modules/changes/20250206_danieldouglas92 @@ -0,0 +1,5 @@ +Changed: Separated out the tian2019 solubility +reaction model into it's own module independent +of the reactive fluid transport material model. +
+(Daniel Douglas, 2025/02/06) diff --git a/include/aspect/material_model/reaction_model/tian2019_solubility.h b/include/aspect/material_model/reaction_model/tian2019_solubility.h new file mode 100644 index 00000000000..c033a84c611 --- /dev/null +++ b/include/aspect/material_model/reaction_model/tian2019_solubility.h @@ -0,0 +1,155 @@ +/* + Copyright (C) 2024 by the authors of the ASPECT code. + + This file is part of ASPECT. + + ASPECT is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2, or (at your option) + any later version. + + ASPECT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ASPECT; see the file LICENSE. If not see + . +*/ + +#ifndef _aspect_material_reaction_model_melt_tian2019_solubility_h +#define _aspect_material_reaction_model_melt_tian2019_solubility_h + +#include +#include +#include +#include + + +namespace aspect +{ + namespace MaterialModel + { + using namespace dealii; + + namespace ReactionModel + { + + /** + * A melt model that calculates the solubility of water according to + * parameterized phase diagrams for four lithologies: + * 1) sediment + * 2) mid-ocean ridge basalt (MORB) + * 3) gabbro + * 4) peridotite + * from Tian, 2019 https://doi.org/10.1029/2019GC008488. + * + * These functions can be used in the calculation of reactive fluid transport + * of water. + * + * @ingroup ReactionModel + */ + template + class Tian2019Solubility : public ::aspect::SimulatorAccess + { + public: + + /** + * Compute the free fluid fraction that is present in the material based on the + * fluid content of the material and the fluid solubility for the given input conditions. + * @p in and @p melt_fraction need to have the same size. + * + * @param in Object that contains the current conditions. + * @param porosity_idx the index of the "porosity" composition + * @param q the quadrature point index + */ + double + melt_fraction (const MaterialModel::MaterialModelInputs &in, + const unsigned int porosity_idx, + unsigned int q) const; + + /** + * Compute the maximum allowed bound water content at the input + * pressure and temperature conditions. This is used to determine + * how free water interacts with the solid phase. + * @param in Object that contains the current conditions. + * @param q the quadrature point index + */ + std::vector tian_equilibrium_bound_water_content(const MaterialModel::MaterialModelInputs &in, + unsigned int q) const; + + /** + * Declare the parameters this function takes through input files. + */ + static + void + declare_parameters (ParameterHandler &prm); + + /** + * Read the parameters from the parameter file. + */ + void + parse_parameters (ParameterHandler &prm); + + private: + + /** + * The maximum water content for each of the 4 rock types in the tian approximation + * method. These are important for keeping the polynomial bounded within reasonable + * values. + */ + double tian_max_peridotite_water; + double tian_max_gabbro_water; + double tian_max_MORB_water; + double tian_max_sediment_water; + + /** + * + * The following coefficients are taken from a publication from Tian et al., 2019, and can be found + * in Table 3 (Gabbro), Table B1 (MORB), Table B2 (Sediments) and Table B3 (peridotite). + * LR refers to the effective enthalpy change for devolatilization reactions, + * csat is the saturated mass fraction of water in the solid, and Td is the + * onset temperature of devolatilization for water. + */ + std::vector LR_peridotite_poly_coeffs {-19.0609, 168.983, -630.032, 1281.84, -1543.14, 1111.88, -459.142, 95.4143, 1.97246}; + std::vector csat_peridotite_poly_coeffs {0.00115628, 2.42179}; + std::vector Td_peridotite_poly_coeffs {-15.4627, 94.9716, 636.603}; + + std::vector LR_gabbro_poly_coeffs {-1.81745, 7.67198, -10.8507, 5.09329, 8.14519}; + std::vector csat_gabbro_poly_coeffs {-0.0176673, 0.0893044, 1.52732}; + std::vector Td_gabbro_poly_coeffs {-1.72277, 20.5898, 637.517}; + + std::vector LR_MORB_poly_coeffs {-1.78177, 7.50871, -10.4840, 5.19725, 7.96365}; + std::vector csat_MORB_poly_coeffs {0.0102725, -0.115390, 0.324452, 1.41588}; + std::vector Td_MORB_poly_coeffs {-3.81280, 22.7809, 638.049}; + + std::vector LR_sediment_poly_coeffs {-2.03283, 10.8186, -21.2119, 18.3351, -6.48711, 8.32459}; + std::vector csat_sediment_poly_coeffs {-0.150662, 0.301807, 1.01867}; + std::vector Td_sediment_poly_coeffs {2.83277, -24.7593, 85.9090, 524.898}; + + /** + * The polynomials breakdown above certain pressures, 10 GPa for peridotite, 26 GPa for gabbro, 16 GPa for MORB, + * and 50 GPa for sediment. These cutoff pressures were determined by extending the pressure range in Tian et al. (2019) + * and observing where the maximum allowed water contents jump towards infinite values. + */ + const std::array pressure_cutoffs {{10, 26, 16, 50}}; + + std::vector> devolatilization_enthalpy_changes {LR_peridotite_poly_coeffs, LR_gabbro_poly_coeffs, \ + LR_MORB_poly_coeffs, LR_sediment_poly_coeffs + }; + + std::vector> water_mass_fractions {csat_peridotite_poly_coeffs, csat_gabbro_poly_coeffs, \ + csat_MORB_poly_coeffs, csat_sediment_poly_coeffs + }; + + std::vector> devolatilization_onset_temperatures {Td_peridotite_poly_coeffs, Td_gabbro_poly_coeffs, \ + Td_MORB_poly_coeffs, Td_sediment_poly_coeffs + }; + }; + } + + } +} + +#endif diff --git a/include/aspect/material_model/reactive_fluid_transport.h b/include/aspect/material_model/reactive_fluid_transport.h index 34a62af4acb..093c3e2b265 100644 --- a/include/aspect/material_model/reactive_fluid_transport.h +++ b/include/aspect/material_model/reactive_fluid_transport.h @@ -28,11 +28,8 @@ #include #include -#include -#include #include - - +#include namespace aspect { @@ -64,18 +61,6 @@ namespace aspect * @} */ - - /** - * Compute the maximum allowed bound water content at the input - * pressure and temperature conditions. This is used to determine - * how free water interacts with the solid phase. - * @param in Object that contains the current conditions. - * @param q unsigned int from 0-3 indexing which rock phase the equilbrium - * bound water content is being calculated for - */ - std::vector tian_equilibrium_bound_water_content(const MaterialModel::MaterialModelInputs &in, - unsigned int q) const; - /** * Compute the free fluid fraction that can be present in the material based on the * fluid content of the material and the fluid solubility for the given input conditions. @@ -160,64 +145,16 @@ namespace aspect */ double fluid_reaction_time_scale; - /** - * The maximum water content for each of the 4 rock types in the tian approximation - * method. These are important for keeping the polynomial bounded within reasonable - * values. - */ - double tian_max_peridotite_water; - double tian_max_gabbro_water; - double tian_max_MORB_water; - double tian_max_sediment_water; - - /** - * - * The following coefficients are taken from a publication from Tian et al., 2019, and can be found - * in Table 3 (Gabbro), Table B1 (MORB), Table B2 (Sediments) and Table B3 (peridotite). - * LR refers to the effective enthalpy change for devolatilization reactions, - * csat is the saturated mass fraction of water in the solid, and Td is the - * onset temperature of devolatilization for water. - */ - std::vector LR_peridotite_poly_coeffs {-19.0609, 168.983, -630.032, 1281.84, -1543.14, 1111.88, -459.142, 95.4143, 1.97246}; - std::vector csat_peridotite_poly_coeffs {0.00115628, 2.42179}; - std::vector Td_peridotite_poly_coeffs {-15.4627, 94.9716, 636.603}; - - std::vector LR_gabbro_poly_coeffs {-1.81745, 7.67198, -10.8507, 5.09329, 8.14519}; - std::vector csat_gabbro_poly_coeffs {-0.0176673, 0.0893044, 1.52732}; - std::vector Td_gabbro_poly_coeffs {-1.72277, 20.5898, 637.517}; - - std::vector LR_MORB_poly_coeffs {-1.78177, 7.50871, -10.4840, 5.19725, 7.96365}; - std::vector csat_MORB_poly_coeffs {0.0102725, -0.115390, 0.324452, 1.41588}; - std::vector Td_MORB_poly_coeffs {-3.81280, 22.7809, 638.049}; - - std::vector LR_sediment_poly_coeffs {-2.03283, 10.8186, -21.2119, 18.3351, -6.48711, 8.32459}; - std::vector csat_sediment_poly_coeffs {-0.150662, 0.301807, 1.01867}; - std::vector Td_sediment_poly_coeffs {2.83277, -24.7593, 85.9090, 524.898}; - - /** - * The polynomials breakdown above certain pressures, 10 GPa for peridotite, 26 GPa for gabbro, 16 GPa for MORB, - * and 50 GPa for sediment. These cutoff pressures were determined by extending the pressure range in Tian et al. (2019) - * and observing where the maximum allowed water contents jump towards infinite values. - */ - const std::vector pressure_cutoffs {10, 26, 16, 50}; - - std::vector> devolatilization_enthalpy_changes {LR_peridotite_poly_coeffs, LR_gabbro_poly_coeffs, \ - LR_MORB_poly_coeffs, LR_sediment_poly_coeffs - }; - - std::vector> water_mass_fractions {csat_peridotite_poly_coeffs, csat_gabbro_poly_coeffs, \ - csat_MORB_poly_coeffs, csat_sediment_poly_coeffs - }; - - std::vector> devolatilization_onset_temperatures {Td_peridotite_poly_coeffs, Td_gabbro_poly_coeffs, \ - Td_MORB_poly_coeffs, Td_sediment_poly_coeffs - }; - /* * Object for computing Katz 2003 melt parameters */ ReactionModel::Katz2003MantleMelting katz2003_model; + /* + * Object for computing Tian 2019 parameterized solubility parameters + */ + ReactionModel::Tian2019Solubility tian2019_model; + /** * Enumeration for selecting which type of scheme to use for * reactions between fluids and solids. The available diff --git a/source/material_model/reaction_model/tian2019_solubility.cc b/source/material_model/reaction_model/tian2019_solubility.cc new file mode 100644 index 00000000000..30f5d3cd955 --- /dev/null +++ b/source/material_model/reaction_model/tian2019_solubility.cc @@ -0,0 +1,207 @@ +/* + Copyright (C) 2024 by the authors of the ASPECT code. + + This file is part of ASPECT. + + ASPECT is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2, or (at your option) + any later version. + + ASPECT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ASPECT; see the file LICENSE. If not see + . +*/ + + +#include +#include + + +namespace aspect +{ + namespace MaterialModel + { + namespace ReactionModel + { + + + template + double + Tian2019Solubility:: + melt_fraction (const MaterialModel::MaterialModelInputs &in, + const unsigned int porosity_idx, + unsigned int q) const + { + double melt_fractions; + // The bound fluid content is calculated using parametrized phase + // diagrams for four different rock types: sediment, MORB, gabbro, and + // peridotite. + const unsigned int bound_fluid_idx = this->introspection().compositional_index_for_name("bound_fluid"); + const unsigned int sediment_idx = this->introspection().compositional_index_for_name("sediment"); + const unsigned int MORB_idx = this->introspection().compositional_index_for_name("MORB"); + const unsigned int gabbro_idx = this->introspection().compositional_index_for_name("gabbro"); + const unsigned int peridotite_idx = this->introspection().compositional_index_for_name("peridotite"); + + // Initialize a vector that stores the compositions (mass fractions) for + // the four different rock compositions, + std::vector tracked_rock_mass_fractions(4); + tracked_rock_mass_fractions[0] = (in.composition[q][peridotite_idx]); + tracked_rock_mass_fractions[1] = (in.composition[q][gabbro_idx]); + tracked_rock_mass_fractions[2] = (in.composition[q][MORB_idx]); + tracked_rock_mass_fractions[3] = (in.composition[q][sediment_idx]); + + // The bound water content (water within the solid phase) for the four different rock types + std::vector tian_eq_bound_water_content = tian_equilibrium_bound_water_content(in, q); + + // average the water content between the four different rock types + double average_eq_bound_water_content = MaterialUtilities::average_value (tracked_rock_mass_fractions, tian_eq_bound_water_content, MaterialUtilities::arithmetic); + + // The fluid volume fraction in equilibrium with the solid (stored in the melt_fractions vector) + // is equal to the sum of the porosity and the change in bound fluid content + // (current bound fluid - updated average bound fluid). + melt_fractions = std::max(in.composition[q][bound_fluid_idx] + in.composition[q][porosity_idx] - average_eq_bound_water_content, 0.0); + + return melt_fractions; + } + + + + template + std::vector + Tian2019Solubility:: + tian_equilibrium_bound_water_content (const MaterialModel::MaterialModelInputs &in, + unsigned int q) const + { + // Create arrays that will store the values of the polynomials at the current pressure + std::vector LR_values(4); + std::vector csat_values(4); + std::vector Td_values(4); + + // Loop over the four rock types (peridotite, gabbro, MORB, sediment) and the polynomial + // coefficients to fill the vectors defined above. The polynomials for LR are defined in + // equations 13, B2, B10, and B18. csat polynomials are defined in equations 14, B1, B9, and B17. + // Td polynomials are defined in equations 15, B3, B11, and B19. + for (unsigned int i = 0; i eq_bound_water_content(4); + + // Define the maximum bound water content allowed for the four different rock compositions + std::vector max_bound_water_content = {tian_max_peridotite_water, tian_max_gabbro_water, tian_max_MORB_water, tian_max_sediment_water}; + + // Loop over all rock compositions and fill the equilibrium bound water content, divide by 100 to convert + // from percentage to fraction (equation 1) + for (unsigned int k = 0; k + void + Tian2019Solubility::declare_parameters (ParameterHandler &prm) + { + prm.declare_entry ("Maximum weight percent water in sediment", "3", + Patterns::Double (0), + "The maximum allowed weight percent that the sediment composition can hold."); + prm.declare_entry ("Maximum weight percent water in MORB", "2", + Patterns::Double (0), + "The maximum allowed weight percent that the sediment composition can hold."); + prm.declare_entry ("Maximum weight percent water in gabbro", "1", + Patterns::Double (0), + "The maximum allowed weight percent that the sediment composition can hold."); + prm.declare_entry ("Maximum weight percent water in peridotite", "8", + Patterns::Double (0), + "The maximum allowed weight percent that the sediment composition can hold."); + } + + + template + void + Tian2019Solubility::parse_parameters (ParameterHandler &prm) + { + AssertThrow(this->introspection().compositional_name_exists("sediment"), + ExcMessage("The Tian approximation only works " + "if there is a compositional field called sediment.")); + AssertThrow(this->introspection().compositional_name_exists("MORB"), + ExcMessage("The Tian approximation only works " + "if there is a compositional field called MORB.")); + AssertThrow(this->introspection().compositional_name_exists("gabbro"), + ExcMessage("The Tian approximation only works " + "if there is a compositional field called gabbro.")); + AssertThrow(this->introspection().compositional_name_exists("peridotite"), + ExcMessage("The Tian approximation only works " + "if there is a compositional field called peridotite.")); + tian_max_peridotite_water = prm.get_double ("Maximum weight percent water in peridotite"); + tian_max_gabbro_water = prm.get_double ("Maximum weight percent water in gabbro"); + tian_max_MORB_water = prm.get_double ("Maximum weight percent water in MORB"); + tian_max_sediment_water = prm.get_double ("Maximum weight percent water in sediment"); + } + } + } +} + + +// explicit instantiations +namespace aspect +{ + namespace MaterialModel + { +#define INSTANTIATE(dim) \ + namespace ReactionModel \ + { \ + template class Tian2019Solubility; \ + } + + ASPECT_INSTANTIATE(INSTANTIATE) + +#undef INSTANTIATE + } +} diff --git a/source/material_model/reactive_fluid_transport.cc b/source/material_model/reactive_fluid_transport.cc index da6d3b00eea..12de82700a0 100644 --- a/source/material_model/reactive_fluid_transport.cc +++ b/source/material_model/reactive_fluid_transport.cc @@ -19,10 +19,7 @@ */ #include -#include #include -#include -#include #include #include @@ -57,79 +54,6 @@ namespace aspect } - - template - std::vector - ReactiveFluidTransport:: - tian_equilibrium_bound_water_content (const MaterialModel::MaterialModelInputs &in, - unsigned int q) const - { - // Create arrays that will store the values of the polynomials at the current pressure - std::vector LR_values(4); - std::vector csat_values(4); - std::vector Td_values(4); - - // Loop over the four rock types (peridotite, gabbro, MORB, sediment) and the polynomial - // coefficients to fill the vectors defined above. The polynomials for LR are defined in - // equations 13, B2, B10, and B18. csat polynomials are defined in equations 14, B1, B9, and B17. - // Td polynomials are defined in equations 15, B3, B11, and B19. - for (unsigned int i = 0; i eq_bound_water_content(4); - - // Define the maximum bound water content allowed for the four different rock compositions - std::vector max_bound_water_content = {tian_max_peridotite_water, tian_max_gabbro_water, tian_max_MORB_water, tian_max_sediment_water}; - - // Loop over all rock compositions and fill the equilibrium bound water content, divide by 100 to convert - // from percentage to fraction (equation 1) - for (unsigned int k = 0; k void ReactiveFluidTransport:: @@ -160,33 +84,7 @@ namespace aspect } case tian_approximation: { - // The bound fluid content is calculated using parametrized phase - // diagrams for four different rock types: sediment, MORB, gabbro, and - // peridotite. - const unsigned int bound_fluid_idx = this->introspection().compositional_index_for_name("bound_fluid"); - const unsigned int sediment_idx = this->introspection().compositional_index_for_name("sediment"); - const unsigned int MORB_idx = this->introspection().compositional_index_for_name("MORB"); - const unsigned int gabbro_idx = this->introspection().compositional_index_for_name("gabbro"); - const unsigned int peridotite_idx = this->introspection().compositional_index_for_name("peridotite"); - - // Initialize a vector that stores the compositions (mass fractions) for - // the four different rock compositions, - std::vector tracked_rock_mass_fractions(4); - tracked_rock_mass_fractions[0] = (in.composition[q][peridotite_idx]); - tracked_rock_mass_fractions[1] = (in.composition[q][gabbro_idx]); - tracked_rock_mass_fractions[2] = (in.composition[q][MORB_idx]); - tracked_rock_mass_fractions[3] = (in.composition[q][sediment_idx]); - - // The bound water content (water within the solid phase) for the four different rock types - std::vector tian_eq_bound_water_content = tian_equilibrium_bound_water_content(in, q); - - // average the water content between the four different rock types - double average_eq_bound_water_content = MaterialUtilities::average_value (tracked_rock_mass_fractions, tian_eq_bound_water_content, MaterialUtilities::arithmetic); - - // The fluid volume fraction in equilibrium with the solid (stored in the melt_fractions vector) - // is equal to the sum of the porosity and the change in bound fluid content - // (current bound fluid - updated average bound fluid). - melt_fractions[q] = std::max(in.composition[q][bound_fluid_idx] + in.composition[q][porosity_idx] - average_eq_bound_water_content, 0.0); + melt_fractions[q] = tian2019_model.melt_fraction(in, porosity_idx, q); break; } case katz2003: @@ -335,6 +233,12 @@ namespace aspect } prm.leave_subsection(); + prm.enter_subsection("Tian 2019 model"); + { + ReactionModel::Tian2019Solubility::declare_parameters(prm); + } + prm.leave_subsection(); + prm.declare_entry("Base model","visco plastic", Patterns::Selection(MaterialModel::get_valid_model_names_pattern()), "The name of a material model incorporating the " @@ -391,18 +295,6 @@ namespace aspect "computed. If the model does not use operator splitting, this parameter is not used. " "Units: yr or s, depending on the ``Use years " "in output instead of seconds'' parameter."); - prm.declare_entry ("Maximum weight percent water in sediment", "3", - Patterns::Double (0), - "The maximum allowed weight percent that the sediment composition can hold."); - prm.declare_entry ("Maximum weight percent water in MORB", "2", - Patterns::Double (0), - "The maximum allowed weight percent that the sediment composition can hold."); - prm.declare_entry ("Maximum weight percent water in gabbro", "1", - Patterns::Double (0), - "The maximum allowed weight percent that the sediment composition can hold."); - prm.declare_entry ("Maximum weight percent water in peridotite", "8", - Patterns::Double (0), - "The maximum allowed weight percent that the sediment composition can hold."); prm.declare_entry ("Fluid-solid reaction scheme", "no reaction", Patterns::Selection("no reaction|zero solubility|tian approximation|katz2003"), "Select what type of scheme to use for reactions between fluid and solid phases. " @@ -449,11 +341,6 @@ namespace aspect fluid_reaction_time_scale = prm.get_double ("Fluid reaction time scale for operator splitting"); reference_T = prm.get_double ("Reference temperature"); - tian_max_peridotite_water = prm.get_double ("Maximum weight percent water in peridotite"); - tian_max_gabbro_water = prm.get_double ("Maximum weight percent water in gabbro"); - tian_max_MORB_water = prm.get_double ("Maximum weight percent water in MORB"); - tian_max_sediment_water = prm.get_double ("Maximum weight percent water in sediment"); - // Create the base model and initialize its SimulatorAccess base // class; it will get a chance to read its parameters below after we // leave the current section. @@ -475,19 +362,14 @@ namespace aspect } else if (prm.get ("Fluid-solid reaction scheme") == "tian approximation") { - AssertThrow(this->introspection().compositional_name_exists("sediment"), - ExcMessage("The Tian approximation only works " - "if there is a compositional field called sediment.")); - AssertThrow(this->introspection().compositional_name_exists("MORB"), - ExcMessage("The Tian approximation only works " - "if there is a compositional field called MORB.")); - AssertThrow(this->introspection().compositional_name_exists("gabbro"), - ExcMessage("The Tian approximation only works " - "if there is a compositional field called gabbro.")); - AssertThrow(this->introspection().compositional_name_exists("peridotite"), - ExcMessage("The Tian approximation only works " - "if there is a compositional field called peridotite.")); fluid_solid_reaction_scheme = tian_approximation; + + prm.enter_subsection("Tian 2019 model"); + { + tian2019_model.initialize_simulator (this->get_simulator()); + tian2019_model.parse_parameters(prm); + } + prm.leave_subsection(); } else if (prm.get ("Fluid-solid reaction scheme") == "katz2003") { diff --git a/tests/tian_MORB.prm b/tests/tian_MORB.prm index 4b0706a9f3e..c59c3a6ac46 100644 --- a/tests/tian_MORB.prm +++ b/tests/tian_MORB.prm @@ -22,6 +22,8 @@ end subsection Material model subsection Reactive Fluid Transport Model - set Maximum weight percent water in MORB = 5.3 + subsection Tian 2019 model + set Maximum weight percent water in MORB = 5.3 + end end end diff --git a/tests/tian_gabbro.prm b/tests/tian_gabbro.prm index b96ed026bf8..90f44caddd9 100644 --- a/tests/tian_gabbro.prm +++ b/tests/tian_gabbro.prm @@ -22,6 +22,8 @@ end subsection Material model subsection Reactive Fluid Transport Model - set Maximum weight percent water in gabbro = 5.1 + subsection Tian 2019 model + set Maximum weight percent water in gabbro = 5.1 + end end end diff --git a/tests/tian_peridotite.prm b/tests/tian_peridotite.prm index df0f6571a49..5a2c309ca55 100644 --- a/tests/tian_peridotite.prm +++ b/tests/tian_peridotite.prm @@ -99,7 +99,9 @@ subsection Material model set Reference fluid density = 3000 set Fluid reaction time scale for operator splitting = 1e2 set Fluid-solid reaction scheme = tian approximation - set Maximum weight percent water in peridotite = 11 + subsection Tian 2019 model + set Maximum weight percent water in peridotite = 11 + end end subsection Simpler model diff --git a/tests/tian_sediment.prm b/tests/tian_sediment.prm index 9e2f6973ba7..56811dac761 100644 --- a/tests/tian_sediment.prm +++ b/tests/tian_sediment.prm @@ -22,6 +22,8 @@ end subsection Material model subsection Reactive Fluid Transport Model - set Maximum weight percent water in sediment = 3.2 + subsection Tian 2019 model + set Maximum weight percent water in sediment = 3.2 + end end end