-
Notifications
You must be signed in to change notification settings - Fork 241
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Refactor the Tian Approximation in the Reactive Fluid Transport Model #6161
Open
danieldouglas92
wants to merge
3
commits into
geodynamics:main
Choose a base branch
from
danieldouglas92:refactor_reactive_fluid_transport_model
base: main
Could not load branches
Branch not found: {{ refName }}
Loading
Could not load tags
Nothing to show
Loading
Are you sure you want to change the base?
Some commits from the old base branch may be removed from the timeline,
and old review comments may become outdated.
Open
Changes from all commits
Commits
Show all changes
3 commits
Select commit
Hold shift + click to select a range
File filter
Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
There are no files selected for viewing
66 changes: 66 additions & 0 deletions
66
contrib/utilities/update_scripts/prm_files/move_tian_reaction_model.py
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,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)) |
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,5 @@ | ||
| Changed: Separated out the tian2019 solubility | ||
| reaction model into it's own module independent | ||
| of the reactive fluid transport material model. | ||
| <br> | ||
| (Daniel Douglas, 2025/02/06) |
155 changes: 155 additions & 0 deletions
155
include/aspect/material_model/reaction_model/tian2019_solubility.h
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,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 | ||
| <http://www.gnu.org/licenses/>. | ||
| */ | ||
|
|
||
| #ifndef _aspect_material_reaction_model_melt_tian2019_solubility_h | ||
| #define _aspect_material_reaction_model_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 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 <int dim> | ||
| class Tian2019Solubility : public ::aspect::SimulatorAccess<dim> | ||
| { | ||
| 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<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 the quadrature point index | ||
| */ | ||
| std::vector<double> tian_equilibrium_bound_water_content(const MaterialModel::MaterialModelInputs<dim> &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<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}; | ||
|
|
||
| /** | ||
| * 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<double,4 > pressure_cutoffs {{10, 26, 16, 50}}; | ||
|
|
||
| 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.
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
std::vector<double>toconstexpr std::array<double,N>where N is the number of values in this array. Initialization and access to these arrays can remain unchanged.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 had tried to do this with arrays originally but I ultimately ran into issues because all of the arrays are different lengths. Because I ultimately end up storing each of these
poly_coeffsarrays into a vector farther down (Lines 138-148). It seems like there might be a way forward with tuples, when I tried this the looping over eachpoly_coeffsarray seemed a little strange and this type of approach doesn't seem like it's done anywhere else in the code.I think the only way to get around the way that I'm currently storing these constant values would be to alter the way the code within the
tian_equilibrium_bound_water_contentfunction. But because all of these arrays are not necessarily the same size, I'm not sure that there's a good way to do this without adding a bunch of if statements.