Line data Source code
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 :
10 : #include "GeochemistryReactorBase.h"
11 :
12 : InputParameters
13 2268 : GeochemistryReactorBase::sharedParams()
14 : {
15 2268 : InputParameters params = emptyInputParameters();
16 4536 : 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 4536 : 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 2268 : "free_mineral activity log10activity fugacity log10fugacity");
33 4536 : params.addRequiredParam<MultiMooseEnum>(
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 2268 : "g_per_kg_solvent mg_per_kg_solvent ug_per_kg_solvent cm3");
55 4536 : params.addRequiredParam<MultiMooseEnum>(
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 4536 : 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 4536 : params.addRequiredParam<std::vector<Real>>(
72 : "constraint_value", "Numerical value of the containts on constraint_species");
73 6804 : params.addRangeCheckedParam<Real>(
74 4536 : "max_ionic_strength", 3.0, "max_ionic_strength >= 0.0", "Maximum value of ionic strength");
75 4536 : params.addParam<unsigned>(
76 : "extra_iterations_to_make_consistent",
77 4536 : 0,
78 : "Extra iterations to make the molalities, activities, etc, consistent "
79 : "before commencing the Newton process to find the aqueous configuration");
80 4536 : 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 4536 : 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 4536 : params.addParam<Real>(
92 : "abs_tol",
93 4536 : 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 4536 : params.addParam<Real>("rel_tol",
97 4536 : 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 6804 : params.addRangeCheckedParam<Real>("min_initial_molality",
102 4536 : 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 4536 : params.addParam<unsigned>(
107 : "max_iter",
108 4536 : 100,
109 : "Maximum number of Newton iterations allowed when finding the aqueous configuration");
110 4536 : params.addParam<Real>(
111 : "max_initial_residual",
112 4536 : 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 6804 : params.addRangeCheckedParam<Real>(
116 : "swap_threshold",
117 4536 : 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 4536 : params.addParam<unsigned>(
123 : "max_swaps_allowed",
124 4536 : 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 4536 : params.addParam<unsigned>(
129 : "ramp_max_ionic_strength_initial",
130 4536 : 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 4536 : params.addParam<bool>(
136 : "ionic_str_using_basis_only",
137 4536 : 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 4536 : params.addParam<bool>("stoichiometric_ionic_str_using_Cl_only",
146 4536 : 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 2268 : return params;
151 2268 : }
152 :
153 : InputParameters
154 1593 : GeochemistryReactorBase::validParams()
155 : {
156 1593 : InputParameters params = NodalUserObject::validParams();
157 1593 : params += GeochemistryReactorBase::sharedParams();
158 :
159 3186 : params.addRequiredParam<UserObjectName>(
160 : "model_definition", "The name of the GeochemicalModelDefinition user object.");
161 4779 : params.addRangeCheckedParam<Real>(
162 : "stoichiometry_tolerance",
163 3186 : 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 1593 : params.addClassDescription("Base class for UserObject to solve geochemistry reactions");
170 1593 : return params;
171 0 : }
172 :
173 933 : GeochemistryReactorBase::GeochemistryReactorBase(const InputParameters & parameters)
174 : : NodalUserObject(parameters),
175 933 : _num_my_nodes(_subproblem.mesh().getMesh().n_local_nodes()),
176 933 : _mgd(getUserObject<GeochemicalModelDefinition>("model_definition").getDatabase()),
177 1866 : _pgs(getUserObject<GeochemicalModelDefinition>("model_definition")
178 933 : .getPertinentGeochemicalSystem()),
179 933 : _num_basis(_mgd.basis_species_name.size()),
180 933 : _num_eqm(_mgd.eqm_species_name.size()),
181 933 : _initial_max_ionic_str(getParam<Real>("max_ionic_strength") /
182 1866 : (1.0 + getParam<unsigned>("ramp_max_ionic_strength_initial"))),
183 2799 : _is(_initial_max_ionic_str,
184 933 : _initial_max_ionic_str,
185 1866 : getParam<bool>("ionic_str_using_basis_only"),
186 1866 : getParam<bool>("stoichiometric_ionic_str_using_Cl_only")),
187 933 : _gac(_is,
188 933 : getUserObject<GeochemicalModelDefinition>("model_definition").getOriginalFullDatabase()),
189 1866 : _max_swaps_allowed(getParam<unsigned>("max_swaps_allowed")),
190 1866 : _swapper(_num_basis, getParam<Real>("stoichiometry_tolerance")),
191 2799 : _small_molality(getParam<Real>("swap_threshold") * getParam<Real>("abs_tol")),
192 933 : _solver_output(_num_my_nodes),
193 933 : _tot_iter(_num_my_nodes, 0),
194 1866 : _abs_residual(_num_my_nodes, 0.0)
195 : {
196 933 : }
|
