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 "GeochemicalModelInterrogator.h"
11 : #include "GeochemistryFormattedOutput.h"
12 : #include "EquilibriumConstantInterpolator.h"
13 : #include <limits>
14 :
15 : registerMooseObject("GeochemistryApp", GeochemicalModelInterrogator);
16 :
17 : InputParameters
18 222 : GeochemicalModelInterrogator::sharedParams()
19 : {
20 222 : InputParameters params = emptyInputParameters();
21 444 : params.addRequiredParam<UserObjectName>("model_definition",
22 : "The name of the GeochemicalModelDefinition user object");
23 444 : params.addParam<std::vector<std::string>>(
24 : "swap_out_of_basis",
25 : {},
26 : "Species that should be removed from the model_definition's basis and be replaced with the "
27 : "swap_into_basis species. There must be the same number of species in swap_out_of_basis and "
28 : "swap_into_basis. If this list contains more than one species, the swapping is performed "
29 : "one-by-one, starting with the first pair (swap_out_of_basis[0] and swap_into_basis[0]), "
30 : "then the next pair, etc");
31 444 : params.addParam<std::vector<std::string>>(
32 : "swap_into_basis",
33 : {},
34 : "Species that should be removed from the model_definition's "
35 : "equilibrium species list and added to the basis");
36 444 : params.addParam<std::vector<std::string>>(
37 : "activity_species",
38 : {},
39 : "Species that are provided numerical values of activity (or fugacity for gases) in the "
40 : "activity_value input");
41 444 : params.addParam<std::vector<Real>>(
42 : "activity_values",
43 : {},
44 : "Numerical values for the activity (or fugacity) for the "
45 : "species in the activity_species list. These are activity values, not log10(activity).");
46 444 : params.addParam<unsigned int>(
47 : "precision",
48 444 : 4,
49 : "Precision for printing values. Also, if the absolute value of a stoichiometric coefficient "
50 : "is less than 10^(-precision) then it is set to zero. Also, if equilibrium temperatures are "
51 : "desired, they will be computed to a relative error of 10^(-precision)");
52 444 : params.addParam<std::string>("equilibrium_species",
53 : "",
54 : "Only output results for this equilibrium species. If not "
55 : "provided, results for all equilibrium species will be outputted");
56 444 : MooseEnum interrogation_choice("reaction activity eqm_temperature", "reaction");
57 444 : params.addParam<MooseEnum>(
58 : "interrogation",
59 : interrogation_choice,
60 : "Type of interrogation to perform. reaction: Output equilibrium species reactions and "
61 : "log10K. activity: determine activity products at equilibrium. eqm_temperature: determine "
62 : "temperature to ensure equilibrium");
63 444 : params.addParam<Real>(
64 : "temperature",
65 444 : 25,
66 : "Equilibrium constants will be computed at this temperature [degC]. This is "
67 : "ignored if interrogation=eqm_temperature.");
68 666 : params.addRangeCheckedParam<Real>(
69 : "stoichiometry_tolerance",
70 444 : 1E-6,
71 : "stoichiometry_tolerance >= 0.0",
72 : "Swapping involves inverting matrices via a singular value decomposition. During this "
73 : "process: (1) if abs(singular value) < stoi_tol * L1norm(singular values), then the "
74 : "matrix is deemed singular (so the basis swap is deemed invalid); (2) if abs(any "
75 : "stoichiometric coefficient) < stoi_tol then it is set to zero.");
76 222 : return params;
77 222 : }
78 :
79 : InputParameters
80 148 : GeochemicalModelInterrogator::validParams()
81 : {
82 148 : InputParameters params = Output::validParams();
83 148 : params += GeochemicalModelInterrogator::sharedParams();
84 148 : params.addClassDescription("Performing simple manipulations of and querying a "
85 : "geochemical model");
86 :
87 444 : params.set<ExecFlagEnum>("execute_on") = {EXEC_FINAL};
88 148 : return params;
89 148 : }
90 :
91 74 : GeochemicalModelInterrogator::GeochemicalModelInterrogator(const InputParameters & parameters)
92 : : Output(parameters),
93 : UserObjectInterface(this),
94 74 : _mgd(getUserObject<GeochemicalModelDefinition>("model_definition").getDatabase()),
95 148 : _swapper(_mgd.basis_species_index.size(), getParam<Real>("stoichiometry_tolerance")),
96 148 : _swap_out(getParam<std::vector<std::string>>("swap_out_of_basis")),
97 148 : _swap_in(getParam<std::vector<std::string>>("swap_into_basis")),
98 148 : _precision(getParam<unsigned int>("precision")),
99 148 : _interrogation(getParam<MooseEnum>("interrogation").getEnum<InterrogationChoiceEnum>()),
100 148 : _temperature(getParam<Real>("temperature")),
101 148 : _activity_species(getParam<std::vector<std::string>>("activity_species")),
102 148 : _activity_values(getParam<std::vector<Real>>("activity_values")),
103 296 : _equilibrium_species(getParam<std::string>("equilibrium_species"))
104 : {
105 74 : if (_swap_out.size() != _swap_in.size())
106 0 : paramError("swap_out_of_basis must have same length as swap_into_basis");
107 74 : if (_activity_species.size() != _activity_values.size())
108 0 : paramError("activity_species must have same length as activity_values");
109 74 : }
110 :
111 : void
112 74 : GeochemicalModelInterrogator::output()
113 : {
114 518 : for (const auto & sp : eqmSpeciesOfInterest())
115 : {
116 444 : switch (_interrogation)
117 : {
118 404 : case InterrogationChoiceEnum::REACTION:
119 404 : outputReaction(sp);
120 : break;
121 32 : case InterrogationChoiceEnum::ACTIVITY:
122 32 : outputActivity(sp);
123 : break;
124 8 : case InterrogationChoiceEnum::EQM_TEMPERATURE:
125 8 : outputTemperature(sp);
126 : break;
127 : }
128 74 : }
129 :
130 250 : for (unsigned i = 0; i < _swap_out.size(); ++i)
131 : {
132 : // any exception here is an error
133 : try
134 : {
135 178 : _swapper.performSwap(_mgd, _swap_out[i], _swap_in[i]);
136 : }
137 2 : catch (const MooseException & e)
138 : {
139 2 : mooseError(e.what());
140 0 : }
141 1200 : for (const auto & sp : eqmSpeciesOfInterest())
142 : {
143 1024 : switch (_interrogation)
144 : {
145 944 : case InterrogationChoiceEnum::REACTION:
146 944 : outputReaction(sp);
147 : break;
148 72 : case InterrogationChoiceEnum::ACTIVITY:
149 72 : outputActivity(sp);
150 : break;
151 8 : case InterrogationChoiceEnum::EQM_TEMPERATURE:
152 8 : outputTemperature(sp);
153 : break;
154 : }
155 176 : }
156 : }
157 72 : }
158 :
159 : std::vector<std::string>
160 250 : GeochemicalModelInterrogator::eqmSpeciesOfInterest() const
161 : {
162 250 : if (_equilibrium_species == "")
163 58 : return _mgd.eqm_species_name;
164 : else if (_mgd.eqm_species_index.count(_equilibrium_species) == 1)
165 384 : return {_equilibrium_species};
166 0 : return {};
167 192 : }
168 :
169 : void
170 1356 : GeochemicalModelInterrogator::outputReaction(const std::string & eqm_species) const
171 : {
172 : if (_mgd.eqm_species_index.count(eqm_species) == 0)
173 0 : return;
174 1356 : std::stringstream ss;
175 1356 : const unsigned row = _mgd.eqm_species_index.at(eqm_species);
176 1356 : const Real cutoff = std::pow(10.0, -1.0 * _precision);
177 1356 : const std::vector<Real> temps = _mgd.original_database->getTemperatures();
178 1356 : const unsigned numT = temps.size();
179 1356 : const std::string model_type = _mgd.original_database->getLogKModel();
180 : EquilibriumConstantInterpolator log10K(
181 1356 : temps, _mgd.eqm_log10K.sub_matrix(row, 1, 0, numT).get_values(), model_type);
182 1356 : log10K.generate();
183 1356 : const Real log10_eqm_const = log10K.sample(_temperature);
184 1356 : ss << std::setprecision(_precision);
185 1356 : ss << eqm_species << " = ";
186 2712 : ss << GeochemistryFormattedOutput::reaction(
187 1356 : _mgd.eqm_stoichiometry, row, _mgd.basis_species_name, cutoff, _precision);
188 1356 : ss << " . log10(K) = " << log10_eqm_const;
189 : ss << std::endl;
190 1356 : _console << ss.str() << std::flush;
191 2712 : }
192 :
193 : bool
194 616 : GeochemicalModelInterrogator::knownActivity(const std::string & species) const
195 : {
196 896 : for (const auto & sp : _activity_species)
197 352 : if (sp == species)
198 : return true;
199 : if (_mgd.basis_species_index.count(species))
200 424 : return _mgd.basis_species_mineral[_mgd.basis_species_index.at(species)];
201 : if (_mgd.eqm_species_index.count(species))
202 120 : return _mgd.eqm_species_mineral[_mgd.eqm_species_index.at(species)];
203 : return false;
204 : }
205 :
206 : Real
207 304 : GeochemicalModelInterrogator::getActivity(const std::string & species) const
208 : {
209 : unsigned ind = 0;
210 424 : for (const auto & sp : _activity_species)
211 : {
212 192 : if (sp == species)
213 72 : return _activity_values[ind];
214 120 : ind += 1;
215 : }
216 : return 1.0; // must be a mineral
217 : }
218 :
219 : void
220 104 : GeochemicalModelInterrogator::outputActivity(const std::string & eqm_species) const
221 : {
222 : if (_mgd.eqm_species_index.count(eqm_species) == 0)
223 0 : return;
224 :
225 104 : const unsigned row = _mgd.eqm_species_index.at(eqm_species);
226 104 : const unsigned num_cols = _mgd.basis_species_index.size();
227 104 : const Real cutoff = std::pow(10.0, -1.0 * _precision);
228 104 : const std::vector<Real> temps = _mgd.original_database->getTemperatures();
229 104 : const unsigned numT = temps.size();
230 104 : const std::string model_type = _mgd.original_database->getLogKModel();
231 : EquilibriumConstantInterpolator log10K(
232 104 : temps, _mgd.eqm_log10K.sub_matrix(row, 1, 0, numT).get_values(), model_type);
233 104 : log10K.generate();
234 104 : Real rhs = log10K.sample(_temperature);
235 104 : std::stringstream lhs;
236 :
237 104 : lhs << std::setprecision(_precision);
238 104 : if (knownActivity(eqm_species))
239 104 : rhs += std::log10(getActivity(eqm_species));
240 : else
241 0 : lhs << "(A_" << eqm_species << ")^-1 ";
242 :
243 640 : for (unsigned i = 0; i < num_cols; ++i)
244 536 : if (std::abs(_mgd.eqm_stoichiometry(row, i)) > cutoff)
245 : {
246 464 : const std::string sp = _mgd.basis_species_name[i];
247 464 : if (knownActivity(sp))
248 160 : rhs -= _mgd.eqm_stoichiometry(row, i) * std::log10(getActivity(sp));
249 : else
250 912 : lhs << "(A_" << sp << ")^" << _mgd.eqm_stoichiometry(row, i) << " ";
251 : }
252 :
253 104 : lhs << "= 10^" << rhs << std::endl;
254 104 : _console << lhs.str() << std::flush;
255 104 : }
256 :
257 : void
258 16 : GeochemicalModelInterrogator::outputTemperature(const std::string & eqm_species) const
259 : {
260 : if (_mgd.eqm_species_index.count(eqm_species) == 0)
261 8 : return;
262 :
263 16 : const unsigned row = _mgd.eqm_species_index.at(eqm_species);
264 16 : const unsigned num_cols = _mgd.basis_species_index.size();
265 16 : const Real cutoff = std::pow(10.0, -1.0 * _precision);
266 : Real rhs = 1.0;
267 :
268 : // check that we know all the required activities, otherwise compute the RHS
269 16 : if (knownActivity(eqm_species))
270 16 : rhs = -std::log10(getActivity(eqm_species));
271 : else
272 : {
273 0 : _console << "Not enough activites known to compute equilibrium temperature for reaction\n "
274 0 : << std::flush;
275 0 : outputReaction(eqm_species);
276 0 : return;
277 : }
278 :
279 48 : for (unsigned i = 0; i < num_cols; ++i)
280 40 : if (std::abs(_mgd.eqm_stoichiometry(row, i)) > cutoff)
281 : {
282 32 : const std::string sp = _mgd.basis_species_name[i];
283 32 : if (knownActivity(sp))
284 24 : rhs += _mgd.eqm_stoichiometry(row, i) * std::log10(getActivity(sp));
285 : else
286 : {
287 8 : _console << "Not enough activites known to compute equilibrium temperature for reaction\n "
288 8 : << std::flush;
289 8 : outputReaction(eqm_species);
290 : return;
291 : }
292 : }
293 :
294 8 : const Real tsoln = solveForT(
295 16 : _mgd.eqm_log10K.sub_matrix(row, 1, 0, _mgd.original_database->getTemperatures().size()), rhs);
296 8 : _console << eqm_species << ". T = " << tsoln << "degC" << std::endl;
297 : }
298 :
299 : Real
300 8 : GeochemicalModelInterrogator::solveForT(const DenseMatrix<Real> & reference_log10K, Real rhs) const
301 : {
302 8 : const std::vector<Real> temps = _mgd.original_database->getTemperatures();
303 8 : const unsigned numT = temps.size();
304 8 : const std::string model_type = _mgd.original_database->getLogKModel();
305 :
306 : // find the bracket that contains the rhs
307 : unsigned bracket = 0;
308 16 : for (bracket = 0; bracket < numT - 1; ++bracket)
309 : {
310 16 : if (reference_log10K(0, bracket) <= rhs && reference_log10K(0, bracket + 1) > rhs)
311 : break;
312 8 : else if (reference_log10K(0, bracket) >= rhs && reference_log10K(0, bracket + 1) < rhs)
313 : break;
314 : }
315 :
316 8 : if (bracket == numT - 1)
317 : return std::numeric_limits<double>::quiet_NaN();
318 :
319 8 : EquilibriumConstantInterpolator log10K(temps, reference_log10K.get_values(), model_type);
320 8 : log10K.generate();
321 : // now do a Newton-Raphson to find T for which log10K.sample(T) = rhs
322 8 : Real temp = _mgd.original_database->getTemperatures()[bracket + 1];
323 : const Real small_delT =
324 8 : (temp + GeochemistryConstants::CELSIUS_TO_KELVIN) * std::pow(10.0, -1.0 * _precision);
325 : Real del_temp = small_delT;
326 : unsigned iter = 0;
327 32 : while (std::abs(del_temp) >= small_delT && iter++ < 100)
328 : {
329 24 : Real residual = log10K.sample(temp) - rhs;
330 24 : del_temp = -residual / log10K.sampleDerivative(temp);
331 24 : temp += del_temp;
332 : }
333 : return temp;
334 8 : }
|
