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 "GeochemistrySpeciesSwapper.h"
11 :
12 1615 : GeochemistrySpeciesSwapper::GeochemistrySpeciesSwapper(unsigned basis_size, Real stoi_tol)
13 1615 : : _stoi_tol(stoi_tol),
14 1615 : _swap_matrix(basis_size, basis_size),
15 1615 : _true_swap_matrix(basis_size, basis_size),
16 1615 : _inv_swap_matrix(basis_size, basis_size),
17 1615 : _swap_sigma(basis_size),
18 1615 : _swap_U(basis_size, basis_size),
19 1615 : _swap_VT(basis_size, basis_size)
20 : {
21 1615 : }
22 :
23 : void
24 845 : GeochemistrySpeciesSwapper::checkSwap(const ModelGeochemicalDatabase & mgd,
25 : const std::string & replace_this,
26 : const std::string & with_this)
27 : {
28 845 : if (replace_this == "H2O")
29 1 : mooseError("Cannot remove H2O from the basis");
30 : if (mgd.basis_species_index.count(replace_this) == 0)
31 3 : mooseError(replace_this + " is not in the basis, so cannot be removed from the basis");
32 : if (mgd.eqm_species_index.count(with_this) == 0)
33 1 : mooseError(with_this + " is not an equilibrium species, so cannot be "
34 : "removed from the equilibrium species list");
35 :
36 840 : checkSwap(mgd, mgd.basis_species_index.at(replace_this), mgd.eqm_species_index.at(with_this));
37 835 : }
38 :
39 : void
40 1249 : GeochemistrySpeciesSwapper::checkSwap(const ModelGeochemicalDatabase & mgd,
41 : unsigned basis_index_to_replace,
42 : unsigned eqm_index_to_insert)
43 : {
44 1249 : const unsigned num_cols = mgd.basis_species_index.size();
45 1249 : const unsigned num_rows = mgd.eqm_species_index.size();
46 1249 : if (basis_index_to_replace == 0)
47 1 : mooseError("Cannot remove H2O from the basis");
48 1248 : if (basis_index_to_replace >= num_cols)
49 1 : mooseError(basis_index_to_replace, " exceeds the number of basis species in the problem");
50 1247 : if (eqm_index_to_insert >= num_rows)
51 1 : mooseError(eqm_index_to_insert, " exceeds the number of equilibrium species in the problem");
52 1246 : if (mgd.surface_sorption_related[eqm_index_to_insert])
53 1 : mooseError(
54 : "Equilibrium species ",
55 : mgd.eqm_species_name[eqm_index_to_insert],
56 : " is involved in surface sorption so cannot be swapped into the basis. If this is truly "
57 : "necessary, code enhancements will need to be made including: recording whether basis "
58 : "species are involved in surface sorption, including them in the surface-potential "
59 : "calculations, and carefully swapping surface-potential-modified equilibrium constants");
60 :
61 1245 : constructInverseMatrix(mgd, basis_index_to_replace, eqm_index_to_insert);
62 1241 : }
63 :
64 : void
65 1245 : GeochemistrySpeciesSwapper::constructInverseMatrix(const ModelGeochemicalDatabase & mgd,
66 : unsigned basis_index_to_replace,
67 : unsigned eqm_index_to_insert)
68 : {
69 1245 : const unsigned num_cols = mgd.basis_species_index.size();
70 :
71 1245 : if (_swap_matrix.n() != num_cols)
72 1 : mooseError("GeochemistrySpeciesSwapper constructed with incorrect "
73 : "basis_species size");
74 :
75 : // This is a private method, called from checkSwap, so we know that all the inputs have been
76 : // sanity-checked. The only way the swap could be invalid is that the swap matrix is not
77 : // invertible. This could be due to the solve algorithm attempting to perform an invalid swap, so
78 : // only mooseExceptions are thrown below, instead of mooseError, to allow the solve algorithm a
79 : // chance to attempt another swap.
80 :
81 : // construct the swap matrix. new_basis = swap_matrix * old_basis
82 : // We shove the equilibrium species into the "slot" of the basis species it is
83 : // replacing
84 1244 : _swap_matrix.zero();
85 12744 : for (unsigned i = 0; i < num_cols; ++i)
86 11500 : _swap_matrix(i, i) = 1.0;
87 12744 : for (unsigned i = 0; i < num_cols; ++i)
88 11500 : _swap_matrix(basis_index_to_replace, i) = mgd.eqm_stoichiometry(eqm_index_to_insert, i);
89 :
90 : // record the swap matrix since the SVD trashes it
91 1244 : _true_swap_matrix = _swap_matrix;
92 :
93 : // In order to invert _swap_matrix, perform the SVD: A = U * D * VT
94 : // Although every matrix has a SVD, this process might fail due to numerical errors. Therefore
95 : // the following could be wrapped in a try-catch block, but I have never been able to trigger a
96 : // failure since _swap_matrix is always so well-formed due to all the checks on the database while
97 : // parsing it
98 1244 : _swap_matrix.svd(_swap_sigma, _swap_U, _swap_VT);
99 1244 : const Real l1norm = _swap_sigma.l1_norm();
100 12741 : for (unsigned i = 0; i < num_cols; ++i)
101 11500 : if (std::abs(_swap_sigma(i) / l1norm) < _stoi_tol)
102 3 : mooseException("Matrix is not invertible, which signals an invalid basis swap");
103 :
104 : // Find the inverse, which is VT^-1 * D^-1 * U^-1 = V * D^-1 * UT
105 12724 : for (unsigned i = 0; i < num_cols; ++i)
106 11483 : _swap_U.scale_column(i, 1.0 / _swap_sigma(i)); // (scale the columns of U)^T = D^-1 * UT
107 12724 : for (unsigned i = 0; i < num_cols; ++i)
108 140354 : for (unsigned j = 0; j < num_cols; ++j)
109 128871 : _inv_swap_matrix(i, j) = _swap_VT(j, i); // _inv_swap_matrix = V
110 1241 : _inv_swap_matrix.right_multiply_transpose(_swap_U); // V * (scale the columns of U)^T
111 1241 : }
112 :
113 : void
114 839 : GeochemistrySpeciesSwapper::performSwap(ModelGeochemicalDatabase & mgd,
115 : const std::string & replace_this,
116 : const std::string & with_this)
117 : {
118 : // check the swap is valid (if not then mooseError or mooseException)
119 : // and if it's valid, create the inverse swap matrix
120 839 : checkSwap(mgd, replace_this, with_this);
121 :
122 : // perform the swap inside the MGD datastructure
123 835 : alterMGD(mgd, mgd.basis_species_index.at(replace_this), mgd.eqm_species_index.at(with_this));
124 835 : }
125 :
126 : void
127 406 : GeochemistrySpeciesSwapper::performSwap(ModelGeochemicalDatabase & mgd,
128 : unsigned basis_index_to_replace,
129 : unsigned eqm_index_to_insert)
130 : {
131 : // check the swap is valid (if not then mooseError or mooseException)
132 : // and if it's valid, create the inverse swap matrix
133 406 : checkSwap(mgd, basis_index_to_replace, eqm_index_to_insert);
134 :
135 : // perform the swap inside the MGD datastructure
136 406 : alterMGD(mgd, basis_index_to_replace, eqm_index_to_insert);
137 406 : }
138 :
139 : void
140 3 : GeochemistrySpeciesSwapper::performSwap(ModelGeochemicalDatabase & mgd,
141 : DenseVector<Real> & bulk_composition,
142 : const std::string & replace_this,
143 : const std::string & with_this)
144 : {
145 3 : performSwap(mgd, replace_this, with_this);
146 : // compute the bulk composition expressed in the new basis
147 3 : alterBulkComposition(mgd.basis_species_index.size(), bulk_composition);
148 2 : }
149 :
150 : void
151 406 : GeochemistrySpeciesSwapper::performSwap(ModelGeochemicalDatabase & mgd,
152 : DenseVector<Real> & bulk_composition,
153 : unsigned basis_index_to_replace,
154 : unsigned eqm_index_to_insert)
155 : {
156 406 : performSwap(mgd, basis_index_to_replace, eqm_index_to_insert);
157 : // compute the bulk composition expressed in the new basis
158 406 : alterBulkComposition(mgd.basis_species_index.size(), bulk_composition);
159 406 : }
160 :
161 : void
162 1241 : GeochemistrySpeciesSwapper::alterMGD(ModelGeochemicalDatabase & mgd,
163 : unsigned basis_index_to_replace,
164 : unsigned eqm_index_to_insert)
165 : {
166 1241 : const unsigned num_cols = mgd.basis_species_index.size();
167 1241 : const unsigned num_rows = mgd.eqm_species_index.size();
168 : const unsigned num_temperatures = mgd.eqm_log10K.n();
169 : const unsigned num_redox = mgd.redox_stoichiometry.m();
170 : const unsigned kin_rows = mgd.kin_stoichiometry.m();
171 1241 : const unsigned num_rate = mgd.kin_rate.size();
172 1241 : const unsigned pro_ind_eqm =
173 1241 : num_cols + eqm_index_to_insert; // index of the (eqm species that is being swapped into the
174 : // basis) in the kinetic-rate promoting_indices vectors
175 :
176 : // change names
177 1241 : const std::string basis_name = mgd.basis_species_name[basis_index_to_replace];
178 1241 : const std::string eqm_name = mgd.eqm_species_name[eqm_index_to_insert];
179 : mgd.basis_species_name[basis_index_to_replace] = eqm_name;
180 : mgd.basis_species_index.erase(basis_name);
181 1241 : mgd.basis_species_index[eqm_name] = basis_index_to_replace;
182 : mgd.eqm_species_name[eqm_index_to_insert] = basis_name;
183 : mgd.eqm_species_index.erase(eqm_name);
184 1241 : mgd.eqm_species_index[basis_name] = eqm_index_to_insert;
185 :
186 : // flag indicating whether the redox_lhs is the equilibrium species that is being put into the
187 : // basis
188 1241 : const bool redox_lhs_going_to_basis = (eqm_name == mgd.redox_lhs);
189 : mooseAssert(
190 : !redox_lhs_going_to_basis,
191 : "Should not be able to swap the redox LHS into the basis"); // At present, it is impossible to
192 : // swap redox_lhs into the basis,
193 : // since redox_lhs cannot be an
194 : // equilibrium species, because it
195 : // is explicitly excluded in
196 : // PertinentGeochemicalSystem when
197 : // building the secondary species.
198 : /*
199 : In past versions of the code, it was possible to swap the redox LHS into the basis.
200 : If you want to make it possible again then mgd needs to be modified as follows
201 : if (redox_lhs_going_to_basis)
202 : {
203 : // need to make the redox_lhs equal to the species that is being taken from the basis
204 : mgd.redox_lhs = basis_name;
205 : for (unsigned red = 0; red < num_redox; ++red)
206 : {
207 : const Real alpha = mgd.eqm_stoichiometry(
208 : eqm_index_to_insert, basis_index_to_replace); // must be nonzero due to valid swap
209 : const Real alpha_r = mgd.redox_stoichiometry(red, basis_index_to_replace);
210 : mooseAssert(alpha != alpha_r,
211 : "Cannot swap equilibrium species "
212 : << eqm_name << " with basis species " << basis_name
213 : << " because a redox reaction would result in 0 <-> 1");
214 : const Real coef = 1.0 / (alpha - alpha_r);
215 : for (unsigned i = 0; i < num_cols; ++i)
216 : mgd.redox_stoichiometry(red, i) = coef * (mgd.redox_stoichiometry(red, i) -
217 : mgd.eqm_stoichiometry(eqm_index_to_insert, i));
218 : mgd.redox_stoichiometry(red, basis_index_to_replace) = 0.0;
219 : for (unsigned t = 0; t < num_temperatures; ++t)
220 : mgd.redox_log10K(red, t) =
221 : coef * (mgd.redox_log10K(red, t) - mgd.eqm_log10K(eqm_index_to_insert, t));
222 : }
223 : }
224 : */
225 :
226 : // record that the swap has occurred
227 1241 : if (mgd.swap_to_original_basis.n() == 0)
228 542 : mgd.swap_to_original_basis = _true_swap_matrix;
229 : else
230 699 : mgd.swap_to_original_basis.left_multiply(_true_swap_matrix);
231 1241 : mgd.have_swapped_out_of_basis.push_back(basis_index_to_replace);
232 1241 : mgd.have_swapped_into_basis.push_back(eqm_index_to_insert);
233 :
234 : // express stoichiometry in new basis
235 12724 : for (unsigned i = 0; i < num_cols; ++i)
236 11483 : mgd.eqm_stoichiometry(eqm_index_to_insert, i) = 0.0;
237 1241 : mgd.eqm_stoichiometry(eqm_index_to_insert, basis_index_to_replace) =
238 : 1.0; // 1 * replace_this <-> 1 * replace_this
239 1241 : mgd.eqm_stoichiometry.right_multiply(_inv_swap_matrix);
240 64149 : for (unsigned i = 0; i < num_rows; ++i)
241 857623 : for (unsigned j = 0; j < num_cols; ++j)
242 794715 : if (std::abs(mgd.eqm_stoichiometry(i, j)) < _stoi_tol)
243 535232 : mgd.eqm_stoichiometry(i, j) = 0.0;
244 :
245 : // if the redox_lhs is not changed by the swap, alter the redox stoichiometry
246 1241 : if (!redox_lhs_going_to_basis)
247 1241 : mgd.redox_stoichiometry.right_multiply(_inv_swap_matrix);
248 2036 : for (unsigned red = 0; red < num_redox; ++red)
249 11721 : for (unsigned j = 0; j < num_cols; ++j)
250 10926 : if (std::abs(mgd.redox_stoichiometry(red, j)) < _stoi_tol)
251 7655 : mgd.redox_stoichiometry(red, j) = 0.0;
252 :
253 : // express kinetic stoichiometry in new basis
254 1241 : mgd.kin_stoichiometry.right_multiply(_inv_swap_matrix);
255 1392 : for (unsigned i = 0; i < kin_rows; ++i)
256 1614 : for (unsigned j = 0; j < num_cols; ++j)
257 1463 : if (std::abs(mgd.kin_stoichiometry(i, j)) < _stoi_tol)
258 937 : mgd.kin_stoichiometry(i, j) = 0.0;
259 :
260 : // alter equilibrium constants for each temperature point
261 11169 : for (unsigned t = 0; t < num_temperatures; ++t)
262 : {
263 9928 : const Real log10k_eqm_species = mgd.eqm_log10K(eqm_index_to_insert, t);
264 9928 : mgd.eqm_log10K(eqm_index_to_insert, t) =
265 : 0.0; // 1 * replace_this <-> 1 * replace_this with log10K = 0
266 513192 : for (unsigned row = 0; row < num_rows; ++row)
267 503264 : mgd.eqm_log10K(row, t) -=
268 503264 : mgd.eqm_stoichiometry(row, basis_index_to_replace) * log10k_eqm_species;
269 :
270 : // similar for the redox equations
271 9928 : if (!redox_lhs_going_to_basis)
272 16288 : for (unsigned red = 0; red < num_redox; ++red)
273 6360 : mgd.redox_log10K(red, t) -=
274 6360 : mgd.redox_stoichiometry(red, basis_index_to_replace) * log10k_eqm_species;
275 :
276 : // similar for kinetic
277 11136 : for (unsigned kin = 0; kin < kin_rows; ++kin)
278 1208 : mgd.kin_log10K(kin, t) -=
279 1208 : mgd.kin_stoichiometry(kin, basis_index_to_replace) * log10k_eqm_species;
280 : }
281 :
282 : // swap the "is mineral" information
283 : const bool basis_was_mineral = mgd.basis_species_mineral[basis_index_to_replace];
284 1241 : mgd.basis_species_mineral[basis_index_to_replace] = mgd.eqm_species_mineral[eqm_index_to_insert];
285 1241 : mgd.eqm_species_mineral[eqm_index_to_insert] = basis_was_mineral;
286 :
287 : // swap the "is gas" information
288 : const bool basis_was_gas = mgd.basis_species_gas[basis_index_to_replace];
289 1241 : mgd.basis_species_gas[basis_index_to_replace] = mgd.eqm_species_gas[eqm_index_to_insert];
290 1241 : mgd.eqm_species_gas[eqm_index_to_insert] = basis_was_gas;
291 :
292 : // swap the "is transported" information
293 : const bool basis_was_transported = mgd.basis_species_transported[basis_index_to_replace];
294 2482 : mgd.basis_species_transported[basis_index_to_replace] =
295 1241 : mgd.eqm_species_transported[eqm_index_to_insert];
296 1241 : mgd.eqm_species_transported[eqm_index_to_insert] = basis_was_transported;
297 :
298 : // swap the "charge" information
299 1241 : const Real basis_charge = mgd.basis_species_charge[basis_index_to_replace];
300 1241 : mgd.basis_species_charge[basis_index_to_replace] = mgd.eqm_species_charge[eqm_index_to_insert];
301 1241 : mgd.eqm_species_charge[eqm_index_to_insert] = basis_charge;
302 :
303 : // swap the "radius" information
304 1241 : const Real basis_radius = mgd.basis_species_radius[basis_index_to_replace];
305 1241 : mgd.basis_species_radius[basis_index_to_replace] = mgd.eqm_species_radius[eqm_index_to_insert];
306 1241 : mgd.eqm_species_radius[eqm_index_to_insert] = basis_radius;
307 :
308 : // swap the "molecular weight" information
309 1241 : const Real basis_molecular_weight = mgd.basis_species_molecular_weight[basis_index_to_replace];
310 1241 : mgd.basis_species_molecular_weight[basis_index_to_replace] =
311 : mgd.eqm_species_molecular_weight[eqm_index_to_insert];
312 1241 : mgd.eqm_species_molecular_weight[eqm_index_to_insert] = basis_molecular_weight;
313 :
314 : // swap the "molecular volume" information
315 1241 : const Real basis_molecular_volume = mgd.basis_species_molecular_volume[basis_index_to_replace];
316 1241 : mgd.basis_species_molecular_volume[basis_index_to_replace] =
317 : mgd.eqm_species_molecular_volume[eqm_index_to_insert];
318 1241 : mgd.eqm_species_molecular_volume[eqm_index_to_insert] = basis_molecular_volume;
319 :
320 : // No need to swap surface_complexation_info or gas_chi because they are not
321 : // bound to basis or eqm species
322 :
323 : // swap promoting indices in the rates
324 1385 : for (unsigned r = 0; r < num_rate; ++r)
325 : {
326 : const Real promoting_index_of_original_basis =
327 144 : mgd.kin_rate[r].promoting_indices[basis_index_to_replace];
328 144 : mgd.kin_rate[r].promoting_indices[basis_index_to_replace] =
329 144 : mgd.kin_rate[r].promoting_indices[pro_ind_eqm];
330 144 : mgd.kin_rate[r].promoting_indices[pro_ind_eqm] = promoting_index_of_original_basis;
331 : const Real promoting_monod_index_of_original_basis =
332 144 : mgd.kin_rate[r].promoting_monod_indices[basis_index_to_replace];
333 144 : mgd.kin_rate[r].promoting_monod_indices[basis_index_to_replace] =
334 : mgd.kin_rate[r].promoting_monod_indices[pro_ind_eqm];
335 144 : mgd.kin_rate[r].promoting_monod_indices[pro_ind_eqm] = promoting_monod_index_of_original_basis;
336 : const Real promoting_half_saturation_of_original_basis =
337 144 : mgd.kin_rate[r].promoting_half_saturation[basis_index_to_replace];
338 144 : mgd.kin_rate[r].promoting_half_saturation[basis_index_to_replace] =
339 : mgd.kin_rate[r].promoting_half_saturation[pro_ind_eqm];
340 144 : mgd.kin_rate[r].promoting_half_saturation[pro_ind_eqm] =
341 : promoting_half_saturation_of_original_basis;
342 : }
343 :
344 : // swap progeny, if needed
345 1385 : for (unsigned r = 0; r < num_rate; ++r)
346 : {
347 144 : if (mgd.kin_rate[r].progeny_index == pro_ind_eqm)
348 1 : mgd.kin_rate[r].progeny_index = basis_index_to_replace;
349 143 : else if (mgd.kin_rate[r].progeny_index == basis_index_to_replace)
350 1 : mgd.kin_rate[r].progeny_index = pro_ind_eqm;
351 : }
352 1241 : }
353 :
354 : void
355 409 : GeochemistrySpeciesSwapper::alterBulkComposition(unsigned basis_size,
356 : DenseVector<Real> & bulk_composition) const
357 : {
358 409 : if (bulk_composition.size() != basis_size)
359 1 : mooseError("GeochemistrySpeciesSwapper: bulk_composition has size ",
360 1 : bulk_composition.size(),
361 : " which differs from the basis size");
362 408 : DenseVector<Real> result;
363 408 : _inv_swap_matrix.vector_mult_transpose(result, bulk_composition);
364 : bulk_composition = result;
365 408 : }
366 :
367 : bool
368 80 : GeochemistrySpeciesSwapper::findBestEqmSwap(unsigned basis_ind,
369 : const ModelGeochemicalDatabase & mgd,
370 : const std::vector<Real> & eqm_molality,
371 : bool minerals_allowed,
372 : bool gas_allowed,
373 : bool sorption_allowed,
374 : unsigned & best_eqm_species) const
375 : {
376 80 : const unsigned num_eqm = mgd.eqm_species_name.size();
377 80 : if (eqm_molality.size() != num_eqm)
378 1 : mooseError("Size of eqm_molality is ",
379 1 : eqm_molality.size(),
380 : " which is not equal to the number of equilibrium species ",
381 : num_eqm);
382 79 : if (basis_ind >= mgd.basis_species_name.size())
383 1 : mooseError("basis index ", basis_ind, " must be less than ", mgd.basis_species_name.size());
384 78 : best_eqm_species = 0;
385 : bool legitimate_swap_found = false;
386 : Real best_stoi = 0.0;
387 5715 : for (unsigned j = 0; j < num_eqm; ++j)
388 : {
389 5637 : if (mgd.eqm_stoichiometry(j, basis_ind) == 0.0)
390 2810 : continue;
391 2827 : if (!minerals_allowed && mgd.eqm_species_mineral[j])
392 212 : continue;
393 2615 : if (!gas_allowed && mgd.eqm_species_gas[j])
394 15 : continue;
395 2600 : if (!sorption_allowed && mgd.surface_sorption_related[j])
396 2 : continue;
397 2598 : const Real stoi = std::abs(mgd.eqm_stoichiometry(j, basis_ind)) * eqm_molality[j];
398 2598 : if (stoi >= best_stoi)
399 : {
400 : best_stoi = stoi;
401 485 : best_eqm_species = j;
402 : legitimate_swap_found = true;
403 : }
404 : }
405 78 : return legitimate_swap_found;
406 : }
|
