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 "GeochemicalSolver.h"
11 : #include "GeochemistrySortedIndices.h"
12 :
13 901 : GeochemicalSolver::GeochemicalSolver(unsigned num_basis,
14 : unsigned num_kin,
15 : GeochemistryIonicStrength & is,
16 : Real abs_tol,
17 : Real rel_tol,
18 : unsigned max_iter,
19 : Real max_initial_residual,
20 : Real swap_threshold,
21 : unsigned max_swaps_allowed,
22 : const std::vector<std::string> & prevent_precipitation,
23 : Real max_ionic_strength,
24 : unsigned ramp_max_ionic_strength,
25 901 : bool evaluate_kin_always)
26 901 : : _is(is),
27 901 : _num_basis(num_basis),
28 901 : _num_kin(num_kin),
29 901 : _num_basis_in_algebraic_system(_num_basis),
30 901 : _num_in_algebraic_system(_num_basis + _num_kin),
31 901 : _residual(_num_in_algebraic_system),
32 901 : _abs_residual(0.0),
33 901 : _jacobian(_num_in_algebraic_system, _num_in_algebraic_system),
34 901 : _new_mol(_num_in_algebraic_system),
35 901 : _abs_tol(abs_tol),
36 901 : _rel_tol(rel_tol),
37 901 : _res0_times_rel(0.0),
38 901 : _max_iter(max_iter),
39 901 : _max_initial_residual(max_initial_residual),
40 901 : _swap_threshold(swap_threshold),
41 901 : _max_swaps_allowed(max_swaps_allowed),
42 901 : _prevent_precipitation(prevent_precipitation),
43 901 : _max_ionic_strength(max_ionic_strength),
44 901 : _ramp_max_ionic_strength(ramp_max_ionic_strength),
45 901 : _evaluate_kin_always(evaluate_kin_always),
46 901 : _input_mole_additions(_num_basis + _num_kin),
47 901 : _input_dmole_additions(_num_basis + _num_kin, _num_basis + _num_kin)
48 : {
49 901 : if (_ramp_max_ionic_strength > _max_iter)
50 1 : mooseError("GeochemicalSolver: ramp_max_ionic_strength must be less than max_iter");
51 900 : if (max_initial_residual <= 0.0)
52 1 : mooseError("GeochemicalSolver: max_initial_residual must be positive");
53 899 : if (max_ionic_strength < 0.0)
54 1 : mooseError("GeochemicalSolver: max_ionic_strength must not be negative");
55 898 : if (abs_tol < 0.0)
56 1 : mooseError("GeochemicalSolver: abs_tol must not be negative");
57 897 : if (rel_tol < 0.0)
58 1 : mooseError("GeochemicalSolver: rel_tol must not be negative");
59 896 : if (rel_tol == 0.0 && abs_tol == 0.0)
60 1 : mooseError("GeochemicalSolver: either rel_tol or abs_tol must be positive");
61 913 : }
62 :
63 : void
64 2 : GeochemicalSolver::setMaxInitialResidual(Real max_initial_residual)
65 : {
66 2 : if (max_initial_residual <= 0.0)
67 1 : mooseError("GeochemicalSolver: max_initial_residual must be positive");
68 1 : _max_initial_residual = max_initial_residual;
69 1 : }
70 :
71 : Real
72 2 : GeochemicalSolver::getMaxInitialResidual() const
73 : {
74 2 : return _max_initial_residual;
75 : }
76 :
77 : Real
78 49173 : GeochemicalSolver::computeResidual(const GeochemicalSystem & egs,
79 : DenseVector<Real> & residual,
80 : const DenseVector<Real> & mole_additions) const
81 : {
82 302541 : for (unsigned a_ind = 0; a_ind < _num_in_algebraic_system; ++a_ind)
83 253368 : residual(a_ind) = egs.getResidualComponent(a_ind, mole_additions);
84 49173 : return residual.l1_norm();
85 : }
86 :
87 : void
88 40579 : GeochemicalSolver::solveAndUnderrelax(const GeochemicalSystem & egs,
89 : DenseMatrix<Real> & jacobian,
90 : DenseVector<Real> & new_mol) const
91 : {
92 40579 : jacobian.lu_solve(_residual, new_mol);
93 :
94 : // at this point we want to do molality = molality - new_mol, but
95 : // Bethke recommends underrelaxation - probably want to do PETSc variational bounds in the
96 : // future
97 40579 : Real one_over_delta = 1.0;
98 40579 : const std::vector<Real> current_molality_and_pot = egs.getAlgebraicVariableValues();
99 233742 : for (unsigned a_ind = 0; a_ind < _num_in_algebraic_system; ++a_ind)
100 193163 : one_over_delta =
101 199402 : std::max(one_over_delta, new_mol(a_ind) * 2.0 / current_molality_and_pot[a_ind]);
102 233742 : for (unsigned a_ind = 0; a_ind < _num_in_algebraic_system; ++a_ind)
103 193163 : new_mol(a_ind) = current_molality_and_pot[a_ind] - new_mol(a_ind) / one_over_delta;
104 40579 : }
105 :
106 : bool
107 4879 : GeochemicalSolver::swapNeeded(const GeochemicalSystem & egs,
108 : unsigned & swap_out_of_basis,
109 : unsigned & swap_into_basis,
110 : std::stringstream & ss) const
111 : {
112 : bool swap_needed = false;
113 :
114 4879 : const ModelGeochemicalDatabase & mgd = egs.getModelGeochemicalDatabase();
115 :
116 : // check if any basis minerals have negative free number of moles
117 4879 : const std::vector<Real> & basis_molality = egs.getSolventMassAndFreeMolalityAndMineralMoles();
118 : const std::vector<unsigned> molality_order =
119 4879 : GeochemistrySortedIndices::sortedIndices(basis_molality, true);
120 :
121 : // if the Newton did not converge then check for small molalities in the non-minerals
122 4879 : if (_abs_residual > _res0_times_rel && _abs_residual > _abs_tol)
123 : {
124 29 : for (const auto & i : molality_order)
125 : {
126 29 : if (basis_molality[i] >= _swap_threshold)
127 : {
128 : // since we're going through basis_molality in ascending order, as soon as we get >=
129 : // _swap_threshold, we're safe
130 : break;
131 : }
132 10 : else if (!mgd.basis_species_mineral[i] && egs.getInAlgebraicSystem()[i] &&
133 5 : egs.getChargeBalanceBasisIndex() != i)
134 : {
135 : // a non-mineral in the algebraic system has super low molality: try to find a legitimate
136 : // swap
137 5 : swap_out_of_basis = i;
138 5 : bool legitimate_swap_found = egs.getSwapper().findBestEqmSwap(swap_out_of_basis,
139 : mgd,
140 : egs.getEquilibriumMolality(),
141 : true,
142 : false,
143 : false,
144 : swap_into_basis);
145 5 : if (legitimate_swap_found)
146 : {
147 5 : ss << "Basis species " << mgd.basis_species_name[swap_out_of_basis]
148 10 : << " has very low molality of " << basis_molality[i] << " compared to "
149 5 : << _swap_threshold << ". Swapping it with equilibrium species "
150 5 : << mgd.eqm_species_name[swap_into_basis] << std::endl;
151 : swap_needed = true;
152 : break;
153 : }
154 : // if no legitimate swap is found, then loop around to the next basis species
155 : }
156 : }
157 : }
158 : if (swap_needed)
159 : return swap_needed;
160 :
161 : // now look through the molalities for minerals that are consumed
162 5186 : for (const auto & i : molality_order)
163 : {
164 5186 : if (basis_molality[i] > 0.0)
165 : {
166 : // since we're going through basis_molality in ascending order, as soon as we get >0, we're
167 : // safe
168 : break;
169 : }
170 373 : else if (mgd.basis_species_mineral[i])
171 : {
172 : swap_needed = true;
173 61 : swap_out_of_basis = i;
174 61 : bool legitimate_swap_found = egs.getSwapper().findBestEqmSwap(swap_out_of_basis,
175 : mgd,
176 : egs.getEquilibriumMolality(),
177 : false,
178 : false,
179 : false,
180 : swap_into_basis);
181 61 : if (!legitimate_swap_found)
182 0 : mooseException("Cannot find a legitimate swap for mineral ",
183 : mgd.basis_species_name[swap_out_of_basis]);
184 61 : ss << "Mineral " << mgd.basis_species_name[swap_out_of_basis]
185 : << " consumed. Swapping it with equilibrium species "
186 122 : << mgd.eqm_species_name[swap_into_basis] << std::endl;
187 : break;
188 : }
189 : }
190 : if (swap_needed)
191 : return swap_needed;
192 :
193 : // check maximum saturation index is not positive
194 4813 : const std::vector<Real> & eqm_SI = egs.getSaturationIndices();
195 : const std::vector<unsigned> mineral_SI_order =
196 4813 : GeochemistrySortedIndices::sortedIndices(eqm_SI, false);
197 :
198 76145 : for (const auto & j : mineral_SI_order)
199 71557 : if (eqm_SI[j] > 0.0)
200 507 : if (mgd.eqm_species_mineral[j])
201 507 : if (std::find(_prevent_precipitation.begin(),
202 : _prevent_precipitation.end(),
203 : mgd.eqm_species_name[j]) == _prevent_precipitation.end())
204 : {
205 : // mineral has positive saturation index and user is not preventing its precipitation.
206 : // determine the basis species to swap out
207 : swap_needed = true;
208 225 : swap_into_basis = j;
209 : bool legitimate_swap_found = false;
210 225 : swap_out_of_basis = 0;
211 : Real best_stoi = 0.0;
212 2744 : for (unsigned i = 1; i < _num_basis; ++i) // never swap water (i=0)
213 : {
214 2519 : if (basis_molality[i] > 0.0 && i != egs.getChargeBalanceBasisIndex() &&
215 4769 : !mgd.basis_species_gas[i] && !egs.getBasisActivityKnown()[i])
216 : {
217 : // don't want to swap out the charge-balance species or any gases of fixed fugacity
218 : // or any species with fixed activity
219 1751 : const Real stoi = std::abs(mgd.eqm_stoichiometry(j, i)) / basis_molality[i];
220 1751 : if (stoi > best_stoi)
221 : {
222 : legitimate_swap_found = true;
223 : best_stoi = stoi;
224 428 : swap_out_of_basis = i;
225 : }
226 : }
227 : }
228 225 : if (!legitimate_swap_found)
229 0 : mooseException(
230 : "Cannot find a legitimate swap for the supersaturated equilibrium species ",
231 : mgd.eqm_species_name[j]);
232 225 : ss << "Mineral " << mgd.eqm_species_name[j]
233 : << " supersaturated. Swapping it with basis species "
234 450 : << mgd.basis_species_name[swap_out_of_basis] << std::endl;
235 : break;
236 : }
237 : return swap_needed;
238 : }
239 :
240 : bool
241 1874 : GeochemicalSolver::reduceInitialResidual(GeochemicalSystem & egs,
242 : Real dt,
243 : DenseVector<Real> & mole_additions,
244 : DenseMatrix<Real> & dmole_additions)
245 : {
246 1874 : const Real initial_r = _abs_residual;
247 :
248 : // to get an indication of whether we should increase or decrease molalities in the algorithm
249 : // below, find the median of the original molalities
250 1874 : const std::vector<Real> & original_molality = egs.getAlgebraicBasisValues();
251 : const std::vector<unsigned> mol_order =
252 1874 : GeochemistrySortedIndices::sortedIndices(original_molality, false);
253 1874 : const Real median_molality = original_molality[mol_order[_num_basis_in_algebraic_system / 2]];
254 :
255 : // get the index order of the residual vector (largest first): ignore residuals for surface
256 : // potentials (we're only using _num_basis_in_algebraic_system, not _num_in_algebraic_system)
257 : unsigned ind = 0;
258 1874 : std::vector<unsigned> res_order(_num_basis_in_algebraic_system);
259 : std::iota(res_order.begin(), res_order.end(), ind++);
260 1874 : std::sort(res_order.begin(),
261 : res_order.end(),
262 58215 : [&](int i, int j) { return std::abs(_residual(i)) > std::abs(_residual(j)); });
263 :
264 1874 : const std::vector<Real> & original_molality_and_pot = egs.getAlgebraicVariableValues();
265 : DenseVector<Real> new_molality_and_pot(original_molality_and_pot);
266 2149 : for (const auto & a : res_order)
267 : {
268 2149 : if (std::abs(_residual(a)) < _max_initial_residual)
269 : return false; // haven't managed to find a suitable new molality, and all remaining residual
270 : // components are less than the cutoff, so cannot appreciably reduce from now
271 : // on
272 :
273 2092 : const Real multiplier = (original_molality_and_pot[a] > median_molality) ? 0.5 : 2.0;
274 : // try using the multiplier
275 2092 : new_molality_and_pot(a) = multiplier * original_molality_and_pot[a];
276 2092 : egs.setAlgebraicVariables(new_molality_and_pot);
277 2092 : if (_evaluate_kin_always)
278 : {
279 : mole_additions = _input_mole_additions;
280 338 : dmole_additions = _input_dmole_additions;
281 338 : egs.addKineticRates(dt, mole_additions, dmole_additions);
282 : }
283 2092 : _abs_residual = computeResidual(egs, _residual, mole_additions);
284 2092 : if (_abs_residual < initial_r)
285 : return true; // success: found a new molality that decreases the initial |R|
286 :
287 : // the above approach did not decrease |R|, so try using 1/multiplier
288 1348 : new_molality_and_pot(a) = original_molality_and_pot[a] / multiplier;
289 1348 : egs.setAlgebraicVariables(new_molality_and_pot);
290 1348 : if (_evaluate_kin_always)
291 : {
292 : mole_additions = _input_mole_additions;
293 2 : dmole_additions = _input_dmole_additions;
294 2 : egs.addKineticRates(dt, mole_additions, dmole_additions);
295 : }
296 1348 : _abs_residual = computeResidual(egs, _residual, mole_additions);
297 1348 : if (_abs_residual < initial_r)
298 : return true; // success: found a new molality that decreases the initial |R|
299 :
300 : // the new molalities did not decrease |R|, so revert to the original molality, and move to
301 : // next-biggest residual component
302 275 : new_molality_and_pot(a) = original_molality_and_pot[a];
303 275 : egs.setAlgebraicVariables(new_molality_and_pot);
304 275 : if (_evaluate_kin_always)
305 : {
306 : mole_additions = _input_mole_additions;
307 2 : dmole_additions = _input_dmole_additions;
308 2 : egs.addKineticRates(dt, mole_additions, dmole_additions);
309 : }
310 275 : _abs_residual = computeResidual(egs, _residual, mole_additions);
311 : }
312 : return false;
313 : }
314 :
315 : void
316 4589 : GeochemicalSolver::solveSystem(GeochemicalSystem & egs,
317 : std::stringstream & ss,
318 : unsigned & tot_iter,
319 : Real & abs_residual,
320 : Real dt,
321 : DenseVector<Real> & mole_additions,
322 : DenseMatrix<Real> & dmole_additions)
323 : {
324 4589 : ss.str("");
325 4589 : tot_iter = 0;
326 : unsigned num_swaps = 0;
327 :
328 : // capture the inputs expressed in the current basis
329 : _input_mole_additions = mole_additions;
330 4589 : _input_dmole_additions = dmole_additions;
331 :
332 4589 : const ModelGeochemicalDatabase & mgd = egs.getModelGeochemicalDatabase();
333 :
334 : bool still_swapping = true;
335 9468 : while (still_swapping && num_swaps <= _max_swaps_allowed)
336 : {
337 4879 : _num_basis_in_algebraic_system = egs.getNumBasisInAlgebraicSystem();
338 4879 : _num_in_algebraic_system = egs.getNumInAlgebraicSystem();
339 4879 : _residual = DenseVector<Real>(_num_in_algebraic_system);
340 4879 : _jacobian = DenseMatrix<Real>(_num_in_algebraic_system, _num_in_algebraic_system);
341 4879 : _new_mol = DenseVector<Real>(_num_in_algebraic_system);
342 :
343 : unsigned iter = 0;
344 : const Real max_is0 =
345 4879 : std::min(1.0, (iter + 1.0) / (_ramp_max_ionic_strength + 1.0)) * _max_ionic_strength;
346 4879 : _is.setMaxIonicStrength(max_is0);
347 4879 : _is.setMaxStoichiometricIonicStrength(max_is0);
348 :
349 : mole_additions = _input_mole_additions;
350 4879 : dmole_additions = _input_dmole_additions;
351 4879 : egs.addKineticRates(dt, mole_additions, dmole_additions);
352 :
353 4879 : _abs_residual = computeResidual(egs, _residual, mole_additions);
354 4879 : bool reducing_initial_molalities = (_abs_residual > _max_initial_residual);
355 : const unsigned max_tries =
356 4879 : _num_basis *
357 4879 : unsigned(
358 4879 : std::log(_abs_residual / _max_initial_residual) /
359 4879 : std::log(1.1)); // assume each successful call to reduceInitialResidual reduces residual
360 : // by about 1.1, but that the correct basis species has to be chosen
361 : unsigned tries = 0;
362 6753 : while (reducing_initial_molalities && ++tries <= max_tries)
363 1874 : reducing_initial_molalities = reduceInitialResidual(egs, dt, mole_additions, dmole_additions);
364 :
365 9758 : ss << "iter = " << iter << " |R| = " << _abs_residual << std::endl;
366 4879 : _res0_times_rel = _abs_residual * _rel_tol;
367 45458 : while ((_abs_residual >= _res0_times_rel && _abs_residual >= _abs_tol && iter < _max_iter) ||
368 20590 : iter < _ramp_max_ionic_strength)
369 : {
370 40579 : iter += 1;
371 40579 : tot_iter += 1;
372 40579 : egs.computeJacobian(_residual, _jacobian, mole_additions, dmole_additions);
373 40579 : solveAndUnderrelax(egs, _jacobian, _new_mol);
374 : const Real max_is =
375 40579 : std::min(1.0, (iter + 1.0) / (_ramp_max_ionic_strength + 1.0)) * _max_ionic_strength;
376 40579 : _is.setMaxIonicStrength(max_is);
377 40579 : _is.setMaxStoichiometricIonicStrength(max_is);
378 40579 : egs.setAlgebraicVariables(_new_mol);
379 40579 : if (egs.alterChargeBalanceSpecies(_swap_threshold))
380 : ss << "Changed change balance species to "
381 7 : << mgd.basis_species_name.at(egs.getChargeBalanceBasisIndex()) << std::endl;
382 40579 : if (_evaluate_kin_always)
383 : {
384 : mole_additions = _input_mole_additions;
385 34807 : dmole_additions = _input_dmole_additions;
386 34807 : egs.addKineticRates(dt, mole_additions, dmole_additions);
387 : }
388 40579 : _abs_residual = computeResidual(egs, _residual, mole_additions);
389 81158 : ss << "iter = " << iter << " |R| = " << _abs_residual << std::endl;
390 : }
391 :
392 4879 : abs_residual = _abs_residual; // record the final residual
393 :
394 4879 : if (iter >= _max_iter)
395 53 : ss << std::endl << "Warning: Number of iterations exceeds " << _max_iter << std::endl;
396 :
397 4879 : unsigned swap_out_of_basis = 0;
398 4879 : unsigned swap_into_basis = 0;
399 : try
400 : {
401 : // to ensure basis molality is correct:
402 33430 : for (unsigned i = 0; i < _num_basis; ++i)
403 28551 : egs.addToBulkMoles(i, mole_additions(i));
404 : // now check for small basis molalities, minerals consumed and minerals precipitated
405 4879 : still_swapping = swapNeeded(egs, swap_out_of_basis, swap_into_basis, ss);
406 : }
407 0 : catch (const MooseException & e)
408 : {
409 0 : mooseException(e.what());
410 0 : }
411 4879 : if (still_swapping)
412 : {
413 : // need to do a swap and re-solve
414 : try
415 : {
416 : // before swapping, remove any basis mole_additions that came from kinetics. The
417 : // following loop over i, combined with the above egs.addToBulkMoles(mole_additions), where
418 : // mole_additions = _input_mole_additions + kinetic_additions, means that the bulk
419 : // moles in egs will have only been incremented by _input_mole_additions. Hence,
420 : // _input_mole_additions can be set to zero in preparation for the next solve in the new
421 : // basis.
422 : // NOTE: if _input_mole_additions depend on molalities, this approach introduces error,
423 : // because the following call to addToBulkMoles and subsequent _input_mole_additions = 0
424 : // means the mole additions are FIXED
425 3850 : for (unsigned i = 0; i < _num_basis; ++i)
426 : {
427 3559 : egs.addToBulkMoles(i, _input_mole_additions(i) - mole_additions(i));
428 3559 : _input_mole_additions(i) = 0.0;
429 51412 : for (unsigned j = 0; j < _num_basis + _num_kin; ++j)
430 47853 : _input_dmole_additions(i, j) = 0.0;
431 : }
432 291 : egs.performSwap(swap_out_of_basis, swap_into_basis);
433 291 : num_swaps += 1;
434 : }
435 0 : catch (const MooseException & e)
436 : {
437 0 : mooseException(e.what());
438 0 : }
439 : }
440 : }
441 :
442 4589 : if (num_swaps > _max_swaps_allowed)
443 : {
444 1 : mooseException("Maximum number of swaps performed during solve");
445 : }
446 4588 : else if (_abs_residual >= _res0_times_rel && _abs_residual >= _abs_tol)
447 : {
448 19 : mooseException("Failed to converge");
449 : }
450 : else
451 : {
452 : // the basis species additions have been added above with addtoBulkMoles: now add the kinetic
453 : // additions:
454 29504 : for (unsigned i = 0; i < _num_basis; ++i)
455 24935 : mole_additions(i) = 0.0;
456 4569 : egs.updateOldWithCurrent(mole_additions);
457 4569 : egs.enforceChargeBalance();
458 : }
459 4569 : }
460 :
461 : void
462 3258 : GeochemicalSolver::setRampMaxIonicStrength(unsigned ramp_max_ionic_strength)
463 : {
464 3258 : if (ramp_max_ionic_strength > _max_iter)
465 1 : mooseError("GeochemicalSolver: ramp_max_ionic_strength must be less than max_iter");
466 3257 : _ramp_max_ionic_strength = ramp_max_ionic_strength;
467 3257 : }
468 :
469 : unsigned
470 2 : GeochemicalSolver::getRampMaxIonicStrength() const
471 : {
472 2 : return _ramp_max_ionic_strength;
473 : }
|
