https://mooseframework.inl.gov
GeochemistryReactorBase.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://mooseframework.inl.gov
3 //*
4 //* All rights reserved, see COPYRIGHT for full restrictions
5 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 //*
7 //* Licensed under LGPL 2.1, please see LICENSE for details
8 //* https://www.gnu.org/licenses/lgpl-2.1.html
9 
11 
14 {
16  params.addParam<std::vector<std::string>>(
17  "swap_out_of_basis",
18  {},
19  "Species that should be removed from the model_definition's basis and be replaced with the "
20  "swap_into_basis species");
21  params.addParam<std::vector<std::string>>(
22  "swap_into_basis",
23  {},
24  "Species that should be removed from the model_definition's equilibrium species list and "
25  "added to the basis. There must be the same number of species in swap_out_of_basis and "
26  "swap_into_basis. These swaps are performed before any other computations during the "
27  "initial problem setup. If this list contains more than one species, the swapping is "
28  "performed one-by-one, starting with the first pair (swap_out_of_basis[0] and "
29  "swap_into_basis[0]), then the next pair, etc");
30  MultiMooseEnum constraint_user_meaning(
31  "kg_solvent_water bulk_composition bulk_composition_with_kinetic free_concentration "
32  "free_mineral activity log10activity fugacity log10fugacity");
34  "constraint_meaning",
35  constraint_user_meaning,
36  "Meanings of the numerical values given in constraint_value. kg_solvent_water: can only be "
37  "applied to H2O and units must be kg. bulk_composition: can be applied to all non-gas "
38  "species, and represents the total amount of the basis species contained as free species as "
39  "well as the amount found in secondary species but not in kinetic species, and units must be "
40  "moles or mass (kg, g, etc). bulk_composition_with_kinetic: can be applied to all non-gas "
41  "species, and represents the total amount of the basis species contained as free species as "
42  "well as the amount found in secondary species and in kinetic species, and units must be "
43  "moles or mass (kg, g, etc). free_concentration: can be applied to all basis species that "
44  "are not gas and not H2O and not mineral, and represents the total amount of the basis "
45  "species existing freely (not as secondary species) within the solution, and units must be "
46  "molal or mass_per_kg_solvent. free_mineral: can be applied to all mineral basis species, "
47  "and represents the total amount of the mineral existing freely (precipitated) within the "
48  "solution, and units must be moles, mass or cm3. activity and log10activity: can be applied "
49  "to basis species that are not gas and not mineral and not sorbing sites, and represents the "
50  "activity of the basis species (recall pH = -log10activity), and units must be "
51  "dimensionless. fugacity and log10fugacity: can be applied to gases, and units must be "
52  "dimensionless");
53  MultiMooseEnum constraint_unit("dimensionless moles molal kg g mg ug kg_per_kg_solvent "
54  "g_per_kg_solvent mg_per_kg_solvent ug_per_kg_solvent cm3");
56  "constraint_unit",
57  constraint_unit,
58  "Units of the numerical values given in constraint_value. Dimensionless: should only be "
59  "used for activity or fugacity constraints. Moles: mole number. Molal: moles per kg "
60  "solvent water. kg: kilograms. g: grams. mg: milligrams. ug: micrograms. "
61  "kg_per_kg_solvent: kilograms per kg solvent water. g_per_kg_solvent: grams per kg solvent "
62  "water. mg_per_kg_solvent: milligrams per kg solvent water. ug_per_kg_solvent: micrograms "
63  "per kg solvent water. cm3: cubic centimeters");
64  params.addRequiredParam<std::vector<std::string>>(
65  "constraint_species",
66  "Names of the species that have their values fixed to constraint_value with meaning "
67  "constraint_meaning. All basis species (after swap_into_basis and swap_out_of_basis) must "
68  "be provided with exactly one constraint. These constraints are used to compute the "
69  "configuration during the initial problem setup, and in time-dependent simulations they may "
70  "be modified as time progresses.");
71  params.addRequiredParam<std::vector<Real>>(
72  "constraint_value", "Numerical value of the containts on constraint_species");
73  params.addRangeCheckedParam<Real>(
74  "max_ionic_strength", 3.0, "max_ionic_strength >= 0.0", "Maximum value of ionic strength");
75  params.addParam<unsigned>(
76  "extra_iterations_to_make_consistent",
77  0,
78  "Extra iterations to make the molalities, activities, etc, consistent "
79  "before commencing the Newton process to find the aqueous configuration");
80  params.addRequiredParam<std::string>(
81  "charge_balance_species",
82  "Charge balance will be enforced on this basis species. This means that its bulk mole "
83  "number may be changed from the initial value you provide in order to ensure charge "
84  "neutrality. After the initial swaps have been performed, this must be in the basis, and it "
85  "must be provided with a bulk_composition constraint_meaning.");
86  params.addParam<std::vector<std::string>>(
87  "prevent_precipitation",
88  {},
89  "Mineral species in this list will be prevented from precipitating, irrespective of their "
90  "saturation index (unless they are initially in the basis)");
91  params.addParam<Real>(
92  "abs_tol",
93  1E-10,
94  "If the residual of the algebraic system (measured in mol) is lower than this value, the "
95  "Newton process (that finds the aqueous configuration) is deemed to have converged");
96  params.addParam<Real>("rel_tol",
97  1E-200,
98  "If the residual of the algebraic system (measured in mol) is lower than "
99  "this value times the initial residual, the Newton process (that finds the "
100  "aqueous configuration) is deemed to have converged");
101  params.addRangeCheckedParam<Real>("min_initial_molality",
102  1E-20,
103  "min_initial_molality > 0.0",
104  "Minimum value of the initial-guess molality used in the "
105  "Newton process to find the aqueous configuration");
106  params.addParam<unsigned>(
107  "max_iter",
108  100,
109  "Maximum number of Newton iterations allowed when finding the aqueous configuration");
110  params.addParam<Real>(
111  "max_initial_residual",
112  1E3,
113  "Attempt to alter the initial-guess molalities so that the initial residual for the Newton "
114  "process (that finds the aqueous configuration) is less than this number of moles");
115  params.addRangeCheckedParam<Real>(
116  "swap_threshold",
117  0.1,
118  "swap_threshold >= 0.0",
119  "If the molality of a basis species in the algebraic system falls below swap_threshold * "
120  "abs_tol then it is swapped out of the basis. The dimensions of swap_threshold are "
121  "1/kg(solvent water)");
122  params.addParam<unsigned>(
123  "max_swaps_allowed",
124  20,
125  "Maximum number of swaps allowed during a single attempt at finding the aqueous "
126  "configuration. Usually only a handful of swaps are used: this parameter prevents endless "
127  "cyclic swapping that prevents the algorithm from progressing");
128  params.addParam<unsigned>(
129  "ramp_max_ionic_strength_initial",
130  20,
131  "The number of iterations over which to progressively increase the maximum ionic strength "
132  "(from zero to max_ionic_strength) during the initial equilibration. Increasing this can "
133  "help in convergence of the Newton process, at the cost of spending more time finding the "
134  "aqueous configuration.");
135  params.addParam<bool>(
136  "ionic_str_using_basis_only",
137  false,
138  "If set to true, ionic strength and stoichiometric ionic strength will be computed using "
139  "only the basis molalities, ignoring molalities of equilibrium species. Since basis "
140  "molality is usually greater than equilibrium molality, and the whole Debye-Huckel concept "
141  "of activity coefficients depending on ionic strength is only approximate in practice, "
142  "setting this parameter true often results in a reasonable approximation. It can aid in "
143  "convergence since it eliminates problems associated with unphysical huge equilibrium "
144  "molalities that can occur during Newton iteration to the solution");
145  params.addParam<bool>("stoichiometric_ionic_str_using_Cl_only",
146  false,
147  "If set to true, the stoichiometric ionic strength will be set equal to "
148  "Cl- molality (or max_ionic_strength if the Cl- molality is too "
149  "big). This flag overrides ionic_str_using_basis_molality_only");
150  return params;
151 }
152 
155 {
158 
159  params.addRequiredParam<UserObjectName>(
160  "model_definition", "The name of the GeochemicalModelDefinition user object.");
161  params.addRangeCheckedParam<Real>(
162  "stoichiometry_tolerance",
163  1E-6,
164  "stoichiometry_tolerance >= 0.0",
165  "Swapping involves inverting matrices via a singular value decomposition. During this "
166  "process: (1) if abs(singular value) < stoi_tol * L1norm(singular values), then the "
167  "matrix is deemed singular (so the basis swap is deemed invalid); (2) if abs(any "
168  "stoichiometric coefficient) < stoi_tol then it is set to zero.");
169  params.addClassDescription("Base class for UserObject to solve geochemistry reactions");
170  return params;
171 }
172 
174  : NodalUserObject(parameters),
175  _num_my_nodes(_subproblem.mesh().getMesh().n_local_nodes()),
176  _mgd(getUserObject<GeochemicalModelDefinition>("model_definition").getDatabase()),
177  _pgs(getUserObject<GeochemicalModelDefinition>("model_definition")
178  .getPertinentGeochemicalSystem()),
179  _num_basis(_mgd.basis_species_name.size()),
180  _num_eqm(_mgd.eqm_species_name.size()),
181  _initial_max_ionic_str(getParam<Real>("max_ionic_strength") /
182  (1.0 + getParam<unsigned>("ramp_max_ionic_strength_initial"))),
183  _is(_initial_max_ionic_str,
184  _initial_max_ionic_str,
185  getParam<bool>("ionic_str_using_basis_only"),
186  getParam<bool>("stoichiometric_ionic_str_using_Cl_only")),
187  _gac(_is,
188  getUserObject<GeochemicalModelDefinition>("model_definition").getOriginalFullDatabase()),
189  _max_swaps_allowed(getParam<unsigned>("max_swaps_allowed")),
190  _swapper(_num_basis, getParam<Real>("stoichiometry_tolerance")),
191  _small_molality(getParam<Real>("swap_threshold") * getParam<Real>("abs_tol")),
192  _solver_output(_num_my_nodes),
193  _tot_iter(_num_my_nodes, 0),
194  _abs_residual(_num_my_nodes, 0.0)
195 {
196 }
static InputParameters validParams()
T & getMesh(MooseMesh &mesh)
function to cast mesh
Definition: SCM.h:35
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
MeshBase & mesh
static InputParameters validParams()
static InputParameters sharedParams()
contains params that are shared with AddGeochemistrySolverAction and its children ...
void addRequiredParam(const std::string &name, const std::string &doc_string)
InputParameters emptyInputParameters()
User object that parses a geochemical database file, and only retains information relevant to the cur...
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
GeochemistryReactorBase(const InputParameters &parameters)
void addClassDescription(const std::string &doc_string)
void addRangeCheckedParam(const std::string &name, const T &value, const std::string &parsed_function, const std::string &doc_string)