Refactor the Tian Approximation in the Reactive Fluid Transport Model #6161
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,159 @@ | ||
| /* | ||
| 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 | ||
| <http://www.gnu.org/licenses/>. | ||
| */ | ||
|
|
||
| #ifndef _aspect_material_reaction_melt_tian2019_solubility_h | ||
| #define _aspect_material_reaction_melt_tian2019_solubility_h | ||
|
|
||
| #include <aspect/material_model/interface.h> | ||
| #include <aspect/simulator_access.h> | ||
| #include <aspect/melt.h> | ||
| #include <aspect/utilities.h> | ||
|
|
||
|
|
||
| namespace aspect | ||
| { | ||
| namespace MaterialModel | ||
| { | ||
| using namespace dealii; | ||
|
|
||
| namespace ReactionModel | ||
| { | ||
|
|
||
| /** | ||
| * A melt model that calculates melt fraction and entropy change | ||
| * according to the melting model for dry peridotite of Katz, 2003. | ||
| * This also includes a computation of the latent heat of melting (if the latent heat | ||
| * heating model is active). | ||
| * | ||
| * These functions can be used in the calculation of melting and melt transport | ||
| * in the melt_simple material model and can be extended to other material models | ||
| * | ||
| * @ingroup ReactionModel | ||
| */ | ||
| template <int dim> | ||
| class Tian2019Solubility : public ::aspect::SimulatorAccess<dim> | ||
| { | ||
| public: | ||
| // constructor | ||
| Tian2019Solubility(); | ||
|
|
||
| /** | ||
| * 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); | ||
|
Comment on lines
+57
to
+68
There was a problem hiding this comment. Please put these functions at the end of the list of public functions here and in the source file. This is just a convention we (try to) follow, but it makes navigating the source files easier. |
||
|
|
||
|
|
||
| /** | ||
| * Compute the free fluid fraction that can be present in the material based on the | ||
|
There was a problem hiding this comment.
|
||
| * fluid content of the material and the fluid solubility for the given input conditions. | ||
| * @p in and @p melt_fractions need to have the same size. | ||
|
There was a problem hiding this comment. there is no |
||
| * | ||
| * @param in Object that contains the current conditions. | ||
| * @param melt_fractions Vector of doubles that is filled with the | ||
| * allowable free fluid fraction for each given input conditions. | ||
|
Comment on lines
+77
to
+78
There was a problem hiding this comment. Update this to the new interface, i.e. explain what |
||
| */ | ||
| double | ||
| melt_fraction (const MaterialModel::MaterialModelInputs<dim> &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 unsigned int from 0-3 indexing which rock phase the equilbrium | ||
| * bound water content is being calculated for | ||
|
Comment on lines
+90
to
+91
There was a problem hiding this comment. Oh, is that really the meaning of |
||
| */ | ||
| std::vector<double> tian_equilibrium_bound_water_content(const MaterialModel::MaterialModelInputs<dim> &in, | ||
| unsigned int q) const; | ||
|
|
||
| private: | ||
| /** | ||
| * Parameters for the solubility model of Tian et al., 2019. | ||
| */ | ||
|
Comment on lines
+97
to
+99
There was a problem hiding this comment. I think this comment here is redundant now. All of these parameters are from Tian et al since this is the Tian et al plugin. |
||
|
|
||
| /** | ||
| * 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<double> LR_peridotite_poly_coeffs {-19.0609, 168.983, -630.032, 1281.84, -1543.14, 1111.88, -459.142, 95.4143, 1.97246}; | ||
| std::vector<double> csat_peridotite_poly_coeffs {0.00115628, 2.42179}; | ||
| std::vector<double> Td_peridotite_poly_coeffs {-15.4627, 94.9716, 636.603}; | ||
|
|
||
| std::vector<double> LR_gabbro_poly_coeffs {-1.81745, 7.67198, -10.8507, 5.09329, 8.14519}; | ||
| std::vector<double> csat_gabbro_poly_coeffs {-0.0176673, 0.0893044, 1.52732}; | ||
| std::vector<double> Td_gabbro_poly_coeffs {-1.72277, 20.5898, 637.517}; | ||
|
|
||
| std::vector<double> LR_MORB_poly_coeffs {-1.78177, 7.50871, -10.4840, 5.19725, 7.96365}; | ||
| std::vector<double> csat_MORB_poly_coeffs {0.0102725, -0.115390, 0.324452, 1.41588}; | ||
| std::vector<double> Td_MORB_poly_coeffs {-3.81280, 22.7809, 638.049}; | ||
|
|
||
| std::vector<double> LR_sediment_poly_coeffs {-2.03283, 10.8186, -21.2119, 18.3351, -6.48711, 8.32459}; | ||
| std::vector<double> csat_sediment_poly_coeffs {-0.150662, 0.301807, 1.01867}; | ||
| std::vector<double> Td_sediment_poly_coeffs {2.83277, -24.7593, 85.9090, 524.898}; | ||
|
Comment on lines
+119
to
+133
There was a problem hiding this comment. All of these numbers are fixed and will never change during the program runtime. You can tell this to the compiler, which will then optimize the code further and also prevent you from accidentally changing the values by changing their type from |
||
|
|
||
| /** | ||
| * 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<double> pressure_cutoffs {10, 26, 16, 50}; | ||
|
There was a problem hiding this comment. same here, make this a |
||
|
|
||
| std::vector<std::vector<double>> devolatilization_enthalpy_changes {LR_peridotite_poly_coeffs, LR_gabbro_poly_coeffs, \ | ||
| LR_MORB_poly_coeffs, LR_sediment_poly_coeffs | ||
| }; | ||
|
|
||
| std::vector<std::vector<double>> water_mass_fractions {csat_peridotite_poly_coeffs, csat_gabbro_poly_coeffs, \ | ||
| csat_MORB_poly_coeffs, csat_sediment_poly_coeffs | ||
| }; | ||
|
|
||
| std::vector<std::vector<double>> devolatilization_onset_temperatures {Td_peridotite_poly_coeffs, Td_gabbro_poly_coeffs, \ | ||
| Td_MORB_poly_coeffs, Td_sediment_poly_coeffs | ||
| }; | ||
| }; | ||
| } | ||
|
|
||
| } | ||
| } | ||
|
|
||
| #endif | ||
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I know you took this from the other reaction model, but that one sets a bad example here. These lines should correspond as closely as possible to the folder structure and file name. This is not technically necessary, but a convention we follow to make them easy to remember: