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 1376 : GeochemistrySpeciesSwapper::GeochemistrySpeciesSwapper(unsigned basis_size, Real stoi_tol)
13 1376 : : _stoi_tol(stoi_tol),
14 1376 : _swap_matrix(basis_size, basis_size),
15 1376 : _true_swap_matrix(basis_size, basis_size),
16 1376 : _inv_swap_matrix(basis_size, basis_size),
17 1376 : _swap_sigma(basis_size),
18 1376 : _swap_U(basis_size, basis_size),
19 1376 : _swap_VT(basis_size, basis_size)
20 : {
21 1376 : }
22 :
23 : void
24 649 : GeochemistrySpeciesSwapper::checkSwap(const ModelGeochemicalDatabase & mgd,
25 : const std::string & replace_this,
26 : const std::string & with_this)
27 : {
28 649 : 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 644 : checkSwap(mgd, mgd.basis_species_index.at(replace_this), mgd.eqm_species_index.at(with_this));
37 639 : }
38 :
39 : void
40 1001 : GeochemistrySpeciesSwapper::checkSwap(const ModelGeochemicalDatabase & mgd,
41 : unsigned basis_index_to_replace,
42 : unsigned eqm_index_to_insert)
43 : {
44 1001 : const unsigned num_cols = mgd.basis_species_index.size();
45 1001 : const unsigned num_rows = mgd.eqm_species_index.size();
46 1001 : if (basis_index_to_replace == 0)
47 1 : mooseError("Cannot remove H2O from the basis");
48 1000 : if (basis_index_to_replace >= num_cols)
49 1 : mooseError(basis_index_to_replace, " exceeds the number of basis species in the problem");
50 999 : if (eqm_index_to_insert >= num_rows)
51 1 : mooseError(eqm_index_to_insert, " exceeds the number of equilibrium species in the problem");
52 998 : if (mgd.surface_sorption_related[eqm_index_to_insert])
53 1 : mooseError(
54 : "Equilibrium species ",
55 1 : 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 997 : constructInverseMatrix(mgd, basis_index_to_replace, eqm_index_to_insert);
62 993 : }
63 :
64 : void
65 997 : GeochemistrySpeciesSwapper::constructInverseMatrix(const ModelGeochemicalDatabase & mgd,
66 : unsigned basis_index_to_replace,
67 : unsigned eqm_index_to_insert)
68 : {
69 997 : const unsigned num_cols = mgd.basis_species_index.size();
70 :
71 997 : 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 : _swap_matrix.zero();
85 10060 : for (unsigned i = 0; i < num_cols; ++i)
86 9064 : _swap_matrix(i, i) = 1.0;
87 10060 : for (unsigned i = 0; i < num_cols; ++i)
88 9064 : _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 996 : _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 996 : _swap_matrix.svd(_swap_sigma, _swap_U, _swap_VT);
99 996 : const Real l1norm = _swap_sigma.l1_norm();
100 10057 : for (unsigned i = 0; i < num_cols; ++i)
101 9064 : 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 10040 : for (unsigned i = 0; i < num_cols; ++i)
106 9047 : _swap_U.scale_column(i, 1.0 / _swap_sigma(i)); // (scale the columns of U)^T = D^-1 * UT
107 10040 : for (unsigned i = 0; i < num_cols; ++i)
108 109452 : for (unsigned j = 0; j < num_cols; ++j)
109 100405 : _inv_swap_matrix(i, j) = _swap_VT(j, i); // _inv_swap_matrix = V
110 993 : _inv_swap_matrix.right_multiply_transpose(_swap_U); // V * (scale the columns of U)^T
111 993 : }
112 :
113 : void
114 643 : 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 643 : checkSwap(mgd, replace_this, with_this);
121 :
122 : // perform the swap inside the MGD datastructure
123 639 : alterMGD(mgd, mgd.basis_species_index.at(replace_this), mgd.eqm_species_index.at(with_this));
124 639 : }
125 :
126 : void
127 354 : 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 354 : checkSwap(mgd, basis_index_to_replace, eqm_index_to_insert);
134 :
135 : // perform the swap inside the MGD datastructure
136 354 : alterMGD(mgd, basis_index_to_replace, eqm_index_to_insert);
137 354 : }
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 354 : GeochemistrySpeciesSwapper::performSwap(ModelGeochemicalDatabase & mgd,
152 : DenseVector<Real> & bulk_composition,
153 : unsigned basis_index_to_replace,
154 : unsigned eqm_index_to_insert)
155 : {
156 354 : performSwap(mgd, basis_index_to_replace, eqm_index_to_insert);
157 : // compute the bulk composition expressed in the new basis
158 354 : alterBulkComposition(mgd.basis_species_index.size(), bulk_composition);
159 354 : }
160 :
161 : void
162 993 : GeochemistrySpeciesSwapper::alterMGD(ModelGeochemicalDatabase & mgd,
163 : unsigned basis_index_to_replace,
164 : unsigned eqm_index_to_insert)
165 : {
166 993 : const unsigned num_cols = mgd.basis_species_index.size();
167 993 : 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 993 : const unsigned num_rate = mgd.kin_rate.size();
172 993 : const unsigned pro_ind_eqm =
173 993 : 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 993 : const std::string basis_name = mgd.basis_species_name[basis_index_to_replace];
178 993 : 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 993 : 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 993 : 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 993 : 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 993 : if (mgd.swap_to_original_basis.n() == 0)
228 430 : mgd.swap_to_original_basis = _true_swap_matrix;
229 : else
230 563 : mgd.swap_to_original_basis.left_multiply(_true_swap_matrix);
231 993 : mgd.have_swapped_out_of_basis.push_back(basis_index_to_replace);
232 993 : mgd.have_swapped_into_basis.push_back(eqm_index_to_insert);
233 :
234 : // express stoichiometry in new basis
235 10040 : for (unsigned i = 0; i < num_cols; ++i)
236 9047 : mgd.eqm_stoichiometry(eqm_index_to_insert, i) = 0.0;
237 993 : mgd.eqm_stoichiometry(eqm_index_to_insert, basis_index_to_replace) =
238 : 1.0; // 1 * replace_this <-> 1 * replace_this
239 993 : mgd.eqm_stoichiometry.right_multiply(_inv_swap_matrix);
240 50031 : for (unsigned i = 0; i < num_rows; ++i)
241 665160 : for (unsigned j = 0; j < num_cols; ++j)
242 616122 : if (std::abs(mgd.eqm_stoichiometry(i, j)) < _stoi_tol)
243 410993 : mgd.eqm_stoichiometry(i, j) = 0.0;
244 :
245 : // if the redox_lhs is not changed by the swap, alter the redox stoichiometry
246 993 : if (!redox_lhs_going_to_basis)
247 993 : mgd.redox_stoichiometry.right_multiply(_inv_swap_matrix);
248 1586 : for (unsigned red = 0; red < num_redox; ++red)
249 8622 : for (unsigned j = 0; j < num_cols; ++j)
250 8029 : if (std::abs(mgd.redox_stoichiometry(red, j)) < _stoi_tol)
251 5576 : mgd.redox_stoichiometry(red, j) = 0.0;
252 :
253 : // express kinetic stoichiometry in new basis
254 993 : mgd.kin_stoichiometry.right_multiply(_inv_swap_matrix);
255 1121 : for (unsigned i = 0; i < kin_rows; ++i)
256 1364 : for (unsigned j = 0; j < num_cols; ++j)
257 1236 : if (std::abs(mgd.kin_stoichiometry(i, j)) < _stoi_tol)
258 775 : mgd.kin_stoichiometry(i, j) = 0.0;
259 :
260 : // alter equilibrium constants for each temperature point
261 8937 : for (unsigned t = 0; t < num_temperatures; ++t)
262 : {
263 7944 : const Real log10k_eqm_species = mgd.eqm_log10K(eqm_index_to_insert, t);
264 7944 : mgd.eqm_log10K(eqm_index_to_insert, t) =
265 : 0.0; // 1 * replace_this <-> 1 * replace_this with log10K = 0
266 400248 : for (unsigned row = 0; row < num_rows; ++row)
267 392304 : mgd.eqm_log10K(row, t) -=
268 392304 : mgd.eqm_stoichiometry(row, basis_index_to_replace) * log10k_eqm_species;
269 :
270 : // similar for the redox equations
271 7944 : if (!redox_lhs_going_to_basis)
272 12688 : for (unsigned red = 0; red < num_redox; ++red)
273 4744 : mgd.redox_log10K(red, t) -=
274 4744 : mgd.redox_stoichiometry(red, basis_index_to_replace) * log10k_eqm_species;
275 :
276 : // similar for kinetic
277 8968 : for (unsigned kin = 0; kin < kin_rows; ++kin)
278 1024 : mgd.kin_log10K(kin, t) -=
279 1024 : 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 993 : mgd.basis_species_mineral[basis_index_to_replace] = mgd.eqm_species_mineral[eqm_index_to_insert];
285 993 : 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 993 : mgd.basis_species_gas[basis_index_to_replace] = mgd.eqm_species_gas[eqm_index_to_insert];
290 993 : 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 1986 : mgd.basis_species_transported[basis_index_to_replace] =
295 993 : mgd.eqm_species_transported[eqm_index_to_insert];
296 993 : mgd.eqm_species_transported[eqm_index_to_insert] = basis_was_transported;
297 :
298 : // swap the "charge" information
299 993 : const Real basis_charge = mgd.basis_species_charge[basis_index_to_replace];
300 993 : mgd.basis_species_charge[basis_index_to_replace] = mgd.eqm_species_charge[eqm_index_to_insert];
301 993 : mgd.eqm_species_charge[eqm_index_to_insert] = basis_charge;
302 :
303 : // swap the "radius" information
304 993 : const Real basis_radius = mgd.basis_species_radius[basis_index_to_replace];
305 993 : mgd.basis_species_radius[basis_index_to_replace] = mgd.eqm_species_radius[eqm_index_to_insert];
306 993 : mgd.eqm_species_radius[eqm_index_to_insert] = basis_radius;
307 :
308 : // swap the "molecular weight" information
309 993 : const Real basis_molecular_weight = mgd.basis_species_molecular_weight[basis_index_to_replace];
310 993 : mgd.basis_species_molecular_weight[basis_index_to_replace] =
311 : mgd.eqm_species_molecular_weight[eqm_index_to_insert];
312 993 : mgd.eqm_species_molecular_weight[eqm_index_to_insert] = basis_molecular_weight;
313 :
314 : // swap the "molecular volume" information
315 993 : const Real basis_molecular_volume = mgd.basis_species_molecular_volume[basis_index_to_replace];
316 993 : mgd.basis_species_molecular_volume[basis_index_to_replace] =
317 : mgd.eqm_species_molecular_volume[eqm_index_to_insert];
318 993 : 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 1114 : for (unsigned r = 0; r < num_rate; ++r)
325 : {
326 : const Real promoting_index_of_original_basis =
327 121 : mgd.kin_rate[r].promoting_indices[basis_index_to_replace];
328 121 : mgd.kin_rate[r].promoting_indices[basis_index_to_replace] =
329 121 : mgd.kin_rate[r].promoting_indices[pro_ind_eqm];
330 121 : mgd.kin_rate[r].promoting_indices[pro_ind_eqm] = promoting_index_of_original_basis;
331 : const Real promoting_monod_index_of_original_basis =
332 121 : mgd.kin_rate[r].promoting_monod_indices[basis_index_to_replace];
333 121 : mgd.kin_rate[r].promoting_monod_indices[basis_index_to_replace] =
334 : mgd.kin_rate[r].promoting_monod_indices[pro_ind_eqm];
335 121 : 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 121 : mgd.kin_rate[r].promoting_half_saturation[basis_index_to_replace];
338 121 : mgd.kin_rate[r].promoting_half_saturation[basis_index_to_replace] =
339 : mgd.kin_rate[r].promoting_half_saturation[pro_ind_eqm];
340 121 : 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 1114 : for (unsigned r = 0; r < num_rate; ++r)
346 : {
347 121 : if (mgd.kin_rate[r].progeny_index == pro_ind_eqm)
348 1 : mgd.kin_rate[r].progeny_index = basis_index_to_replace;
349 120 : 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 993 : }
353 :
354 : void
355 357 : GeochemistrySpeciesSwapper::alterBulkComposition(unsigned basis_size,
356 : DenseVector<Real> & bulk_composition) const
357 : {
358 357 : 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 356 : DenseVector<Real> result;
363 356 : _inv_swap_matrix.vector_mult_transpose(result, bulk_composition);
364 : bulk_composition = result;
365 356 : }
366 :
367 : bool
368 73 : 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 73 : const unsigned num_eqm = mgd.eqm_species_name.size();
377 73 : 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 72 : 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 71 : best_eqm_species = 0;
385 : bool legitimate_swap_found = false;
386 : Real best_stoi = 0.0;
387 4911 : for (unsigned j = 0; j < num_eqm; ++j)
388 : {
389 4840 : if (mgd.eqm_stoichiometry(j, basis_ind) == 0.0)
390 2357 : continue;
391 2483 : if (!minerals_allowed && mgd.eqm_species_mineral[j])
392 178 : continue;
393 2305 : if (!gas_allowed && mgd.eqm_species_gas[j])
394 15 : continue;
395 2290 : if (!sorption_allowed && mgd.surface_sorption_related[j])
396 2 : continue;
397 2288 : const Real stoi = std::abs(mgd.eqm_stoichiometry(j, basis_ind)) * eqm_molality[j];
398 2288 : if (stoi >= best_stoi)
399 : {
400 : best_stoi = stoi;
401 439 : best_eqm_species = j;
402 : legitimate_swap_found = true;
403 : }
404 : }
405 71 : return legitimate_swap_found;
406 : }
|
