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 "GeochemicalSystem.h"
11 : #include "GeochemistryActivityCalculators.h"
12 : #include "GeochemistryKineticRateCalculator.h"
13 : #include "EquilibriumConstantInterpolator.h"
14 :
15 1792 : GeochemicalSystem::GeochemicalSystem(ModelGeochemicalDatabase & mgd,
16 : GeochemistryActivityCoefficients & gac,
17 : GeochemistryIonicStrength & is,
18 : GeochemistrySpeciesSwapper & swapper,
19 : const std::vector<std::string> & swap_out_of_basis,
20 : const std::vector<std::string> & swap_into_basis,
21 : const std::string & charge_balance_species,
22 : const std::vector<std::string> & constrained_species,
23 : const std::vector<Real> & constraint_value,
24 : const MultiMooseEnum & constraint_unit,
25 : const MultiMooseEnum & constraint_user_meaning,
26 : Real initial_temperature,
27 : unsigned iters_to_make_consistent,
28 : Real min_initial_molality,
29 : const std::vector<std::string> & kin_name,
30 : const std::vector<Real> & kin_initial,
31 1792 : const MultiMooseEnum & kin_unit)
32 1792 : : _mgd(mgd),
33 1792 : _num_basis(mgd.basis_species_index.size()),
34 1792 : _num_eqm(mgd.eqm_species_index.size()),
35 1792 : _num_redox(mgd.redox_stoichiometry.m()),
36 1792 : _num_surface_pot(mgd.surface_sorption_name.size()),
37 1792 : _num_kin(mgd.kin_species_index.size()),
38 1792 : _swapper(swapper),
39 1792 : _swap_out(swap_out_of_basis),
40 1792 : _swap_in(swap_into_basis),
41 1792 : _gac(gac),
42 1792 : _is(is),
43 1792 : _charge_balance_species(charge_balance_species),
44 1792 : _original_charge_balance_species(charge_balance_species),
45 1792 : _charge_balance_basis_index(0),
46 1792 : _constrained_species(constrained_species),
47 1792 : _constraint_value(constraint_value),
48 1792 : _original_constraint_value(constraint_value),
49 1792 : _constraint_unit(constraint_unit.size()),
50 1792 : _constraint_user_meaning(constraint_user_meaning.size()),
51 1792 : _constraint_meaning(constraint_user_meaning.size()),
52 1792 : _eqm_log10K(_num_eqm),
53 1792 : _redox_log10K(_num_redox),
54 1792 : _kin_log10K(_num_kin),
55 1792 : _num_basis_in_algebraic_system(0),
56 1792 : _num_in_algebraic_system(0),
57 1792 : _in_algebraic_system(_num_basis),
58 1792 : _algebraic_index(_num_basis),
59 1792 : _basis_index(_num_basis),
60 1792 : _bulk_moles_old(_num_basis),
61 1792 : _basis_molality(_num_basis),
62 1792 : _basis_activity_known(_num_basis),
63 1792 : _basis_activity(_num_basis),
64 1792 : _eqm_molality(_num_eqm),
65 1792 : _basis_activity_coef(_num_basis),
66 1792 : _eqm_activity_coef(_num_eqm),
67 1792 : _eqm_activity(_num_eqm),
68 1792 : _surface_pot_expr(_num_surface_pot),
69 1792 : _sorbing_surface_area(_num_surface_pot),
70 1792 : _kin_moles(_num_kin),
71 1792 : _kin_moles_old(_num_kin),
72 1792 : _iters_to_make_consistent(iters_to_make_consistent),
73 1792 : _temperature(initial_temperature),
74 1792 : _min_initial_molality(min_initial_molality),
75 1792 : _original_redox_lhs()
76 : {
77 10205 : for (unsigned i = 0; i < constraint_user_meaning.size(); ++i)
78 8413 : _constraint_user_meaning[i] =
79 8413 : static_cast<ConstraintUserMeaningEnum>(constraint_user_meaning.get(i));
80 10205 : for (unsigned i = 0; i < constraint_unit.size(); ++i)
81 8413 : _constraint_unit[i] =
82 8413 : static_cast<GeochemistryUnitConverter::GeochemistryUnit>(constraint_unit.get(i));
83 1792 : const unsigned ku_size = kin_unit.size();
84 1792 : std::vector<GeochemistryUnitConverter::GeochemistryUnit> k_unit(ku_size);
85 1950 : for (unsigned i = 0; i < ku_size; ++i)
86 158 : k_unit[i] = static_cast<GeochemistryUnitConverter::GeochemistryUnit>(kin_unit.get(i));
87 1792 : checkAndInitialize(kin_name, kin_initial, k_unit);
88 1792 : }
89 :
90 434 : GeochemicalSystem::GeochemicalSystem(
91 : ModelGeochemicalDatabase & mgd,
92 : GeochemistryActivityCoefficients & gac,
93 : GeochemistryIonicStrength & is,
94 : GeochemistrySpeciesSwapper & swapper,
95 : const std::vector<std::string> & swap_out_of_basis,
96 : const std::vector<std::string> & swap_into_basis,
97 : const std::string & charge_balance_species,
98 : const std::vector<std::string> & constrained_species,
99 : const std::vector<Real> & constraint_value,
100 : const std::vector<GeochemistryUnitConverter::GeochemistryUnit> & constraint_unit,
101 : const std::vector<ConstraintUserMeaningEnum> & constraint_user_meaning,
102 : Real initial_temperature,
103 : unsigned iters_to_make_consistent,
104 : Real min_initial_molality,
105 : const std::vector<std::string> & kin_name,
106 : const std::vector<Real> & kin_initial,
107 434 : const std::vector<GeochemistryUnitConverter::GeochemistryUnit> & kin_unit)
108 434 : : _mgd(mgd),
109 434 : _num_basis(mgd.basis_species_index.size()),
110 434 : _num_eqm(mgd.eqm_species_index.size()),
111 434 : _num_redox(mgd.redox_stoichiometry.m()),
112 434 : _num_surface_pot(mgd.surface_sorption_name.size()),
113 434 : _num_kin(mgd.kin_species_index.size()),
114 434 : _swapper(swapper),
115 434 : _swap_out(swap_out_of_basis),
116 434 : _swap_in(swap_into_basis),
117 434 : _gac(gac),
118 434 : _is(is),
119 434 : _charge_balance_species(charge_balance_species),
120 434 : _original_charge_balance_species(charge_balance_species),
121 434 : _charge_balance_basis_index(0),
122 434 : _constrained_species(constrained_species),
123 434 : _constraint_value(constraint_value),
124 434 : _original_constraint_value(constraint_value),
125 434 : _constraint_unit(constraint_unit),
126 434 : _constraint_user_meaning(constraint_user_meaning),
127 472 : _constraint_meaning(constraint_user_meaning.size()),
128 472 : _eqm_log10K(_num_eqm),
129 472 : _redox_log10K(_num_redox),
130 472 : _kin_log10K(_num_kin),
131 434 : _num_basis_in_algebraic_system(0),
132 434 : _num_in_algebraic_system(0),
133 472 : _in_algebraic_system(_num_basis),
134 472 : _algebraic_index(_num_basis),
135 472 : _basis_index(_num_basis),
136 472 : _bulk_moles_old(_num_basis),
137 472 : _basis_molality(_num_basis),
138 472 : _basis_activity_known(_num_basis),
139 472 : _basis_activity(_num_basis),
140 472 : _eqm_molality(_num_eqm),
141 472 : _basis_activity_coef(_num_basis),
142 472 : _eqm_activity_coef(_num_eqm),
143 472 : _eqm_activity(_num_eqm),
144 472 : _surface_pot_expr(_num_surface_pot),
145 472 : _sorbing_surface_area(_num_surface_pot),
146 472 : _kin_moles(_num_kin),
147 472 : _kin_moles_old(_num_kin),
148 434 : _iters_to_make_consistent(iters_to_make_consistent),
149 434 : _temperature(initial_temperature),
150 434 : _min_initial_molality(min_initial_molality),
151 434 : _original_redox_lhs()
152 : {
153 434 : checkAndInitialize(kin_name, kin_initial, kin_unit);
154 510 : }
155 :
156 : void
157 2226 : GeochemicalSystem::checkAndInitialize(
158 : const std::vector<std::string> & kin_name,
159 : const std::vector<Real> & kin_initial,
160 : const std::vector<GeochemistryUnitConverter::GeochemistryUnit> & kin_unit)
161 : {
162 : // initialize every the kinetic species
163 2226 : const unsigned num_kin_name = kin_name.size();
164 2226 : if (!(_num_kin == num_kin_name && _num_kin == kin_initial.size() && _num_kin == kin_unit.size()))
165 5 : mooseError("Initial mole number (or mass or volume) and a unit must be provided for each "
166 : "kinetic species ",
167 5 : _num_kin,
168 : " ",
169 : num_kin_name,
170 : " ",
171 5 : kin_initial.size(),
172 : " ",
173 5 : kin_unit.size());
174 2512 : for (const auto & name_index : _mgd.kin_species_index)
175 : {
176 292 : const unsigned ind = std::distance(
177 292 : kin_name.begin(), std::find(kin_name.begin(), kin_name.end(), name_index.first));
178 292 : if (ind < num_kin_name)
179 : {
180 294 : if (!(kin_unit[ind] == GeochemistryUnitConverter::GeochemistryUnit::MOLES ||
181 46 : kin_unit[ind] == GeochemistryUnitConverter::GeochemistryUnit::KG ||
182 24 : kin_unit[ind] == GeochemistryUnitConverter::GeochemistryUnit::G ||
183 2 : kin_unit[ind] == GeochemistryUnitConverter::GeochemistryUnit::MG ||
184 : kin_unit[ind] == GeochemistryUnitConverter::GeochemistryUnit::UG ||
185 : kin_unit[ind] == GeochemistryUnitConverter::GeochemistryUnit::CM3))
186 1 : mooseError("Kinetic species ",
187 : name_index.first,
188 : ": units must be moles or mass, or volume in the case of minerals");
189 291 : const Real moles = GeochemistryUnitConverter::toMoles(
190 291 : kin_initial[ind], kin_unit[ind], name_index.first, _mgd);
191 291 : setKineticMoles(name_index.second, moles);
192 : }
193 : else
194 0 : mooseError("Initial mole number or mass or volume for kinetic species ",
195 : name_index.first,
196 : " must be provided");
197 : }
198 :
199 : // check sanity of swaps
200 2220 : if (_swap_out.size() != _swap_in.size())
201 1 : mooseError("swap_out_of_basis must have same length as swap_into_basis");
202 2872 : for (unsigned i = 0; i < _swap_out.size(); ++i)
203 654 : if (_swap_out[i] == _charge_balance_species)
204 1 : mooseError("Cannot swap out ",
205 : _charge_balance_species,
206 : " because it is the charge-balance species\n");
207 :
208 : // do swaps desired by user. any exception here is an error
209 2870 : for (unsigned i = 0; i < _swap_out.size(); ++i)
210 : try
211 : {
212 653 : _swapper.performSwap(_mgd, _swap_out[i], _swap_in[i]);
213 : }
214 0 : catch (const MooseException & e)
215 : {
216 0 : mooseError(e.what());
217 0 : }
218 :
219 : // check charge-balance species is in the basis and has a charge
220 2217 : if (_mgd.basis_species_index.count(_charge_balance_species) == 0)
221 1 : mooseError("Cannot enforce charge balance using ",
222 : _charge_balance_species,
223 : " because it is not in the basis");
224 2216 : _charge_balance_basis_index = _mgd.basis_species_index.at(_charge_balance_species);
225 2216 : if (_mgd.basis_species_charge[_charge_balance_basis_index] == 0.0)
226 1 : mooseError("Cannot enforce charge balance using ",
227 : _charge_balance_species,
228 : " because it has zero charge");
229 :
230 : // check that constraint vectors have appropriate sizes
231 2215 : if (_constrained_species.size() != _constraint_value.size())
232 1 : mooseError("Constrained species names must have same length as constraint values");
233 2214 : if (_constrained_species.size() != _constraint_unit.size())
234 1 : mooseError("Constrained species names must have same length as constraint units");
235 2213 : if (_constrained_species.size() != _constraint_user_meaning.size())
236 1 : mooseError("Constrained species names must have same length as constraint meanings");
237 2212 : if (_constrained_species.size() != _num_basis)
238 1 : mooseError("Constrained species names must have same length as the number of species in the "
239 : "basis (each component must be provided with a single constraint");
240 :
241 : // check that each _constrained_species name appears in the basis
242 12798 : for (const auto & name : _mgd.basis_species_name)
243 10589 : if (std::find(_constrained_species.begin(), _constrained_species.end(), name) ==
244 : _constrained_species.end())
245 2 : mooseError("The basis species ", name, " must appear in the constrained species list");
246 :
247 : // order the constraints in the same way as the basis species. This makes the remainder of the
248 : // code much cleaner.
249 2209 : std::vector<std::string> c_s(_num_basis);
250 2209 : std::vector<Real> c_v(_num_basis);
251 2230 : std::vector<GeochemistryUnitConverter::GeochemistryUnit> c_u(_num_basis);
252 2230 : std::vector<ConstraintUserMeaningEnum> c_m(_num_basis);
253 12792 : for (unsigned i = 0; i < _num_basis; ++i)
254 : {
255 10583 : const unsigned basis_ind = _mgd.basis_species_index.at(_constrained_species[i]);
256 10583 : c_s[basis_ind] = _constrained_species[i];
257 10583 : c_v[basis_ind] = _constraint_value[i];
258 10583 : c_u[basis_ind] = _constraint_unit[i];
259 10583 : c_m[basis_ind] = _constraint_user_meaning[i];
260 : }
261 2209 : _constrained_species = c_s;
262 2209 : _constraint_value = c_v;
263 2209 : _constraint_unit = c_u;
264 2209 : _original_constraint_value = c_v;
265 2209 : _constraint_user_meaning = c_m;
266 :
267 : // run through the constraints, checking physical and chemical consistency, converting to mole
268 : // units, and building constraint_meaning
269 12755 : for (unsigned i = 0; i < _constrained_species.size(); ++i)
270 : {
271 10567 : const std::string name = _constrained_species[i];
272 :
273 10567 : switch (_constraint_user_meaning[i])
274 : {
275 749 : case ConstraintUserMeaningEnum::KG_SOLVENT_WATER:
276 : {
277 : // if the mass of solvent water is provided, check it is positive
278 749 : if (_constraint_value[i] <= 0.0)
279 1 : mooseError("Specified mass of solvent water must be positive: you entered ",
280 : _constraint_value[i]);
281 748 : if (_constraint_unit[i] != GeochemistryUnitConverter::GeochemistryUnit::KG)
282 1 : mooseError("Units for kg_solvent_water must be kg");
283 747 : _constraint_meaning[i] = ConstraintMeaningEnum::KG_SOLVENT_WATER;
284 747 : break;
285 : }
286 6789 : case ConstraintUserMeaningEnum::BULK_COMPOSITION:
287 : case ConstraintUserMeaningEnum::BULK_COMPOSITION_WITH_KINETIC:
288 : {
289 : // convert to mole units and specify correct constraint_meaning
290 6789 : _constraint_value[i] = GeochemistryUnitConverter::toMoles(
291 6789 : _constraint_value[i], _constraint_unit[i], name, _mgd);
292 : // add contributions from kinetic moles, if necessary
293 6789 : if (_constraint_user_meaning[i] == ConstraintUserMeaningEnum::BULK_COMPOSITION)
294 7625 : for (unsigned k = 0; k < _mgd.kin_species_name.size(); ++k)
295 849 : _constraint_value[i] += _kin_moles[k] * _mgd.kin_stoichiometry(k, i);
296 6802 : if (!(_constraint_unit[i] == GeochemistryUnitConverter::GeochemistryUnit::MOLES ||
297 234 : _constraint_unit[i] == GeochemistryUnitConverter::GeochemistryUnit::KG ||
298 233 : _constraint_unit[i] == GeochemistryUnitConverter::GeochemistryUnit::G ||
299 : _constraint_unit[i] == GeochemistryUnitConverter::GeochemistryUnit::MG ||
300 : _constraint_unit[i] == GeochemistryUnitConverter::GeochemistryUnit::UG))
301 1 : mooseError("Species ", name, ": units for bulk composition must be moles or mass");
302 6788 : if (name == "H2O")
303 1112 : _constraint_meaning[i] = ConstraintMeaningEnum::MOLES_BULK_WATER;
304 : else
305 5676 : _constraint_meaning[i] = ConstraintMeaningEnum::MOLES_BULK_SPECIES;
306 : break;
307 : }
308 711 : case ConstraintUserMeaningEnum::FREE_CONCENTRATION:
309 : {
310 : // if free concentration, check it is positive and perform the translation to mole units
311 711 : if (_constraint_value[i] <= 0.0)
312 1 : mooseError("Specified free concentration values must be positive: you entered ",
313 : _constraint_value[i]);
314 711 : if (!(_constraint_unit[i] == GeochemistryUnitConverter::GeochemistryUnit::MOLAL ||
315 : _constraint_unit[i] ==
316 2 : GeochemistryUnitConverter::GeochemistryUnit::KG_PER_KG_SOLVENT ||
317 : _constraint_unit[i] ==
318 1 : GeochemistryUnitConverter::GeochemistryUnit::G_PER_KG_SOLVENT ||
319 : _constraint_unit[i] ==
320 : GeochemistryUnitConverter::GeochemistryUnit::MG_PER_KG_SOLVENT ||
321 : _constraint_unit[i] ==
322 : GeochemistryUnitConverter::GeochemistryUnit::UG_PER_KG_SOLVENT))
323 1 : mooseError(
324 : "Species ",
325 : name,
326 : ": units for free concentration quantities must be molal or mass_per_kg_solvent");
327 1418 : _constraint_value[i] = GeochemistryUnitConverter::toMoles(
328 709 : _constraint_value[i], _constraint_unit[i], name, _mgd);
329 709 : _constraint_meaning[i] = ConstraintMeaningEnum::FREE_MOLALITY;
330 709 : break;
331 : }
332 362 : case ConstraintUserMeaningEnum::FREE_MINERAL:
333 : {
334 : // if free mineral, check it is positive and perform the translation to mole units
335 362 : if (_constraint_value[i] <= 0.0)
336 1 : mooseError("Specified free mineral values must be positive: you entered ",
337 : _constraint_value[i]);
338 429 : if (!(_constraint_unit[i] == GeochemistryUnitConverter::GeochemistryUnit::MOLES ||
339 90 : _constraint_unit[i] == GeochemistryUnitConverter::GeochemistryUnit::KG ||
340 90 : _constraint_unit[i] == GeochemistryUnitConverter::GeochemistryUnit::G ||
341 68 : _constraint_unit[i] == GeochemistryUnitConverter::GeochemistryUnit::MG ||
342 : _constraint_unit[i] == GeochemistryUnitConverter::GeochemistryUnit::UG ||
343 : _constraint_unit[i] == GeochemistryUnitConverter::GeochemistryUnit::CM3))
344 1 : mooseError("Species ",
345 : name,
346 : ": units for free mineral quantities must be moles, mass or volume");
347 720 : _constraint_value[i] = GeochemistryUnitConverter::toMoles(
348 360 : _constraint_value[i], _constraint_unit[i], name, _mgd);
349 360 : _constraint_meaning[i] = ConstraintMeaningEnum::FREE_MOLES_MINERAL_SPECIES;
350 360 : break;
351 : }
352 1288 : case ConstraintUserMeaningEnum::ACTIVITY:
353 : {
354 1288 : if (_constraint_value[i] <= 0.0)
355 1 : mooseError("Specified activity values must be positive: you entered ",
356 : _constraint_value[i]);
357 1287 : if (_constraint_unit[i] != GeochemistryUnitConverter::GeochemistryUnit::DIMENSIONLESS)
358 1 : mooseError(
359 : "Species ", name, ": dimensionless units must be used when specifying activity");
360 1286 : _constraint_meaning[i] = ConstraintMeaningEnum::ACTIVITY;
361 1286 : break;
362 : }
363 561 : case ConstraintUserMeaningEnum::LOG10ACTIVITY:
364 : {
365 : // if log10activity is provided, convert it to activity
366 561 : if (_constraint_unit[i] != GeochemistryUnitConverter::GeochemistryUnit::DIMENSIONLESS)
367 1 : mooseError("Species ",
368 : name,
369 : ": dimensionless units must be used when specifying log10activity\n");
370 560 : _constraint_value[i] = std::pow(10.0, _constraint_value[i]);
371 560 : _constraint_meaning[i] = ConstraintMeaningEnum::ACTIVITY;
372 560 : break;
373 : }
374 72 : case ConstraintUserMeaningEnum::FUGACITY:
375 : {
376 : // if fugacity is provided, check it is positive
377 72 : if (_constraint_value[i] <= 0.0)
378 1 : mooseError("Specified fugacity values must be positive: you entered ",
379 : _constraint_value[i]);
380 71 : if (_constraint_unit[i] != GeochemistryUnitConverter::GeochemistryUnit::DIMENSIONLESS)
381 1 : mooseError(
382 : "Species ", name, ": dimensionless units must be used when specifying fugacity\n");
383 70 : _constraint_meaning[i] = ConstraintMeaningEnum::FUGACITY;
384 70 : break;
385 : }
386 35 : case ConstraintUserMeaningEnum::LOG10FUGACITY:
387 : {
388 : // if log10fugacity is provided, convert it to fugacity
389 35 : if (_constraint_unit[i] != GeochemistryUnitConverter::GeochemistryUnit::DIMENSIONLESS)
390 1 : mooseError("Species ",
391 : name,
392 : ": dimensionless units must be used when specifying log10fugacity\n");
393 34 : _constraint_value[i] = std::pow(10.0, _constraint_value[i]);
394 34 : _constraint_meaning[i] = ConstraintMeaningEnum::FUGACITY;
395 34 : break;
396 : }
397 : }
398 :
399 : // check that water is provided with correct meaning
400 10554 : if (name == "H2O")
401 2206 : if (!(_constraint_meaning[i] == ConstraintMeaningEnum::MOLES_BULK_WATER ||
402 : _constraint_meaning[i] == ConstraintMeaningEnum::KG_SOLVENT_WATER ||
403 : _constraint_meaning[i] == ConstraintMeaningEnum::ACTIVITY))
404 1 : mooseError("H2O must be provided with either a mass of solvent water, a bulk composition "
405 : "in moles or mass, or an activity");
406 :
407 : // check that gases are provided with the correct meaning
408 10553 : if (_mgd.basis_species_gas[i])
409 104 : if (_constraint_meaning[i] != ConstraintMeaningEnum::FUGACITY)
410 1 : mooseError("The gas ", name, " must be provided with a fugacity");
411 :
412 : // check that minerals are provided with the correct meaning
413 10552 : if (_mgd.basis_species_mineral[i])
414 503 : if (!(_constraint_meaning[i] == ConstraintMeaningEnum::FREE_MOLES_MINERAL_SPECIES ||
415 : _constraint_meaning[i] == ConstraintMeaningEnum::MOLES_BULK_SPECIES))
416 2 : mooseError("The mineral ",
417 : name,
418 : " must be provided with either: a free number of moles, a free mass or a free "
419 : "volume; or a bulk composition of moles or mass");
420 :
421 : // check that non-water, non-minerals, non-gases are provided with the correct meaning
422 10550 : if (name != "H2O" && !_mgd.basis_species_gas[i] && !_mgd.basis_species_mineral[i])
423 7741 : if (!(_constraint_meaning[i] == ConstraintMeaningEnum::FREE_MOLALITY ||
424 : _constraint_meaning[i] == ConstraintMeaningEnum::ACTIVITY ||
425 : _constraint_meaning[i] == ConstraintMeaningEnum::MOLES_BULK_SPECIES))
426 2 : mooseError("The basis species ",
427 : name,
428 : " must be provided with a free concentration, bulk composition or an activity");
429 :
430 : // check that the charge-balance species has been provided MOLES_BULK_SPECIES
431 10548 : if (name == _charge_balance_species)
432 2203 : if (_constraint_meaning[i] != ConstraintMeaningEnum::MOLES_BULK_SPECIES)
433 2 : mooseError("For code consistency, the species ",
434 : name,
435 : " must be provided with a bulk composition because it is the charge-balance "
436 : "species. The value provided should be a reasonable estimate of the mole "
437 : "number, but will be overridden as the solve progresses");
438 : }
439 2188 : _original_redox_lhs = _mgd.redox_lhs;
440 :
441 2188 : initialize();
442 2209 : }
443 :
444 : void
445 2188 : GeochemicalSystem::initialize()
446 : {
447 2188 : buildTemperatureDependentQuantities(_temperature);
448 2188 : enforceChargeBalanceIfSimple(_constraint_value, _bulk_moles_old);
449 2188 : buildAlgebraicInfo(_in_algebraic_system,
450 2188 : _num_basis_in_algebraic_system,
451 2188 : _num_in_algebraic_system,
452 2188 : _algebraic_index,
453 2188 : _basis_index);
454 2188 : initBulkAndFree(_bulk_moles_old, _basis_molality);
455 2188 : buildKnownBasisActivities(_basis_activity_known, _basis_activity);
456 :
457 2188 : _eqm_molality.assign(_num_eqm, 0.0);
458 2188 : _surface_pot_expr.assign(_num_surface_pot, 1.0);
459 :
460 2188 : computeConsistentConfiguration();
461 2188 : }
462 :
463 : void
464 47936 : GeochemicalSystem::computeConsistentConfiguration()
465 : {
466 : // the steps 1 and 2 below could be iterated for a long time (or a Newton process could even be
467 : // followed) to provide better estimates of activities and molalities, but this is not done in the
468 : // conventional geochemistry approach: there are just too many unknowns and approximations
469 : // employed during the algebraic-system solve to justify iterating towards the perfectly
470 : // consistent initial condition
471 95908 : for (unsigned picard = 0; picard < _iters_to_make_consistent + 1; ++picard)
472 : {
473 : // Step 1: compute ionic strengths and activities using the eqm molalities
474 47972 : _gac.setInternalParameters(_temperature, _mgd, _basis_molality, _eqm_molality, _kin_moles);
475 47972 : _gac.buildActivityCoefficients(_mgd, _basis_activity_coef, _eqm_activity_coef);
476 47972 : updateBasisMolalityForKnownActivity(_basis_molality);
477 47972 : computeRemainingBasisActivities(_basis_activity);
478 :
479 : // Step 2: compute equilibrium molality based on the activities just computed
480 47972 : computeEqmMolalities(_eqm_molality);
481 : }
482 :
483 47936 : computeBulk(_bulk_moles_old);
484 47936 : computeFreeMineralMoles(_basis_molality);
485 47936 : computeSorbingSurfaceArea(_sorbing_surface_area);
486 47936 : }
487 :
488 : unsigned
489 2814 : GeochemicalSystem::getChargeBalanceBasisIndex() const
490 : {
491 2814 : return _charge_balance_basis_index;
492 : }
493 :
494 : Real
495 2516 : GeochemicalSystem::getLog10K(unsigned j) const
496 : {
497 2516 : if (j >= _num_eqm)
498 1 : mooseError("Cannot retrieve log10K for equilibrium species ",
499 : j,
500 : " since there are only ",
501 1 : _num_eqm,
502 : " equilibrium species");
503 2515 : return _eqm_log10K[j];
504 : }
505 :
506 : unsigned
507 8 : GeochemicalSystem::getNumRedox() const
508 : {
509 8 : return _num_redox;
510 : }
511 :
512 : Real
513 143 : GeochemicalSystem::getRedoxLog10K(unsigned red) const
514 : {
515 143 : if (red >= _num_redox)
516 1 : mooseError("Cannot retrieve log10K for redox species ",
517 : red,
518 : " since there are only ",
519 1 : _num_redox,
520 : " redox species");
521 142 : return _redox_log10K[red];
522 : }
523 :
524 : Real
525 143 : GeochemicalSystem::log10RedoxActivityProduct(unsigned red) const
526 : {
527 143 : if (red >= _num_redox)
528 1 : mooseError("Cannot retrieve activity product for redox species ",
529 : red,
530 : " since there are only ",
531 1 : _num_redox,
532 : " redox species");
533 : Real log10ap = 0.0;
534 2268 : for (unsigned basis_i = 0; basis_i < _num_basis; ++basis_i)
535 2126 : log10ap += _mgd.redox_stoichiometry(red, basis_i) * std::log10(_basis_activity[basis_i]);
536 142 : return log10ap;
537 : }
538 :
539 : unsigned
540 681 : GeochemicalSystem::getNumKinetic() const
541 : {
542 681 : return _num_kin;
543 : }
544 :
545 : Real
546 55 : GeochemicalSystem::getKineticLog10K(unsigned kin) const
547 : {
548 55 : if (kin >= _num_kin)
549 1 : mooseError("Cannot retrieve log10K for kinetic species ",
550 : kin,
551 : " since there are only ",
552 1 : _num_kin,
553 : " kinetic species");
554 54 : return _kin_log10K[kin];
555 : }
556 :
557 : const std::vector<Real> &
558 2 : GeochemicalSystem::getKineticLog10K() const
559 : {
560 2 : return _kin_log10K;
561 : }
562 :
563 : Real
564 5963 : GeochemicalSystem::log10KineticActivityProduct(unsigned kin) const
565 : {
566 5963 : if (kin >= _num_kin)
567 1 : mooseError("Cannot retrieve activity product for kinetic species ",
568 : kin,
569 : " since there are only ",
570 1 : _num_kin,
571 : " kinetic species");
572 : Real log10ap = 0.0;
573 52501 : for (unsigned basis_i = 0; basis_i < _num_basis; ++basis_i)
574 46539 : log10ap += _mgd.kin_stoichiometry(kin, basis_i) * std::log10(_basis_activity[basis_i]);
575 5962 : return log10ap;
576 : }
577 :
578 : void
579 3064 : GeochemicalSystem::buildTemperatureDependentQuantities(Real temperature)
580 : {
581 3064 : const std::vector<Real> temps = _mgd.original_database->getTemperatures();
582 3064 : const unsigned numT = temps.size();
583 3064 : const std::string model_type = _mgd.original_database->getLogKModel();
584 :
585 65874 : for (unsigned eqm_j = 0; eqm_j < _num_eqm; ++eqm_j)
586 : {
587 : EquilibriumConstantInterpolator interp(
588 62810 : temps, _mgd.eqm_log10K.sub_matrix(eqm_j, 1, 0, numT).get_values(), model_type);
589 62810 : interp.generate();
590 62810 : _eqm_log10K[eqm_j] = interp.sample(temperature);
591 62810 : }
592 4021 : for (unsigned red = 0; red < _num_redox; ++red)
593 : {
594 : EquilibriumConstantInterpolator interp(
595 957 : temps, _mgd.redox_log10K.sub_matrix(red, 1, 0, numT).get_values(), model_type);
596 957 : interp.generate();
597 957 : _redox_log10K[red] = interp.sample(temperature);
598 957 : }
599 3671 : for (unsigned kin = 0; kin < _num_kin; ++kin)
600 : {
601 : EquilibriumConstantInterpolator interp(
602 607 : temps, _mgd.kin_log10K.sub_matrix(kin, 1, 0, numT).get_values(), model_type);
603 607 : interp.generate();
604 607 : _kin_log10K[kin] = interp.sample(temperature);
605 607 : }
606 3064 : }
607 :
608 : void
609 3000 : GeochemicalSystem::buildAlgebraicInfo(std::vector<bool> & in_algebraic_system,
610 : unsigned & num_basis_in_algebraic_system,
611 : unsigned & num_in_algebraic_system,
612 : std::vector<unsigned> & algebraic_index,
613 : std::vector<unsigned> & basis_index) const
614 : {
615 : // build in_algebraic_system
616 21142 : for (const auto & name_index : _mgd.basis_species_index)
617 : {
618 18142 : const std::string name = name_index.first;
619 18142 : const unsigned basis_ind = name_index.second;
620 18142 : const ConstraintMeaningEnum meaning = _constraint_meaning[basis_ind];
621 18142 : if (name == "H2O")
622 3000 : in_algebraic_system[basis_ind] = (meaning == ConstraintMeaningEnum::MOLES_BULK_WATER);
623 15142 : else if (_mgd.basis_species_gas[basis_ind])
624 : in_algebraic_system[basis_ind] = false;
625 14930 : else if (_mgd.basis_species_mineral[basis_ind])
626 : in_algebraic_system[basis_ind] = false;
627 : else
628 13071 : in_algebraic_system[basis_ind] = (meaning == ConstraintMeaningEnum::MOLES_BULK_SPECIES);
629 : }
630 :
631 : // build algebraic_index and basis_index
632 3000 : num_basis_in_algebraic_system = 0;
633 3000 : algebraic_index.resize(_num_basis, 0);
634 3000 : basis_index.resize(_num_basis, 0);
635 21142 : for (unsigned basis_ind = 0; basis_ind < _num_basis; ++basis_ind)
636 18142 : if (in_algebraic_system[basis_ind])
637 : {
638 11987 : algebraic_index[basis_ind] = _num_basis_in_algebraic_system;
639 11987 : basis_index[_num_basis_in_algebraic_system] = basis_ind;
640 11987 : num_basis_in_algebraic_system += 1;
641 : }
642 :
643 3000 : num_in_algebraic_system = num_basis_in_algebraic_system + _num_surface_pot + _num_kin;
644 3000 : }
645 :
646 : unsigned
647 4946 : GeochemicalSystem::getNumInAlgebraicSystem() const
648 : {
649 4946 : return _num_in_algebraic_system;
650 : }
651 :
652 : unsigned
653 4890 : GeochemicalSystem::getNumBasisInAlgebraicSystem() const
654 : {
655 4890 : return _num_basis_in_algebraic_system;
656 : }
657 :
658 : unsigned
659 342 : GeochemicalSystem::getNumSurfacePotentials() const
660 : {
661 342 : return _num_surface_pot;
662 : }
663 :
664 : const std::vector<bool> &
665 18 : GeochemicalSystem::getInAlgebraicSystem() const
666 : {
667 18 : return _in_algebraic_system;
668 : }
669 :
670 : const std::vector<unsigned> &
671 5 : GeochemicalSystem::getBasisIndexOfAlgebraicSystem() const
672 : {
673 5 : return _basis_index;
674 : }
675 :
676 : const std::vector<unsigned> &
677 5 : GeochemicalSystem::getAlgebraicIndexOfBasisSystem() const
678 : {
679 5 : return _algebraic_index;
680 : }
681 :
682 : std::vector<Real>
683 42470 : GeochemicalSystem::getAlgebraicVariableValues() const
684 : {
685 42470 : std::vector<Real> var(_num_basis_in_algebraic_system + _num_surface_pot + _num_kin);
686 244991 : for (unsigned a = 0; a < _num_basis_in_algebraic_system; ++a)
687 202521 : var[a] = _basis_molality[_basis_index[a]];
688 45629 : for (unsigned s = 0; s < _num_surface_pot; ++s)
689 3159 : var[s + _num_basis_in_algebraic_system] = _surface_pot_expr[s];
690 47731 : for (unsigned k = 0; k < _num_kin; ++k)
691 5261 : var[k + _num_basis_in_algebraic_system + _num_surface_pot] = _kin_moles[k];
692 42470 : return var;
693 : }
694 :
695 : std::vector<Real>
696 1881 : GeochemicalSystem::getAlgebraicBasisValues() const
697 : {
698 1881 : std::vector<Real> var(_num_basis_in_algebraic_system);
699 19530 : for (unsigned a = 0; a < _num_basis_in_algebraic_system; ++a)
700 17649 : var[a] = _basis_molality[_basis_index[a]];
701 1881 : return var;
702 : }
703 :
704 : DenseVector<Real>
705 6 : GeochemicalSystem::getAlgebraicVariableDenseValues() const
706 : {
707 6 : DenseVector<Real> var(_num_in_algebraic_system);
708 16 : for (unsigned a = 0; a < _num_basis_in_algebraic_system; ++a)
709 10 : var(a) = _basis_molality[_basis_index[a]];
710 7 : for (unsigned s = 0; s < _num_surface_pot; ++s)
711 1 : var(s + _num_basis_in_algebraic_system) = _surface_pot_expr[s];
712 9 : for (unsigned k = 0; k < _num_kin; ++k)
713 3 : var(k + _num_basis_in_algebraic_system + _num_surface_pot) = _kin_moles[k];
714 6 : return var;
715 : }
716 :
717 : void
718 2188 : GeochemicalSystem::initBulkAndFree(std::vector<Real> & bulk_moles_old,
719 : std::vector<Real> & basis_molality) const
720 : {
721 12696 : for (unsigned i = 0; i < _num_basis; ++i)
722 : {
723 : // water is done first, so dividing by basis_molality[0] is OK
724 10508 : const Real value = _constraint_value[i];
725 10508 : const ConstraintMeaningEnum meaning = _constraint_meaning[i];
726 10508 : switch (meaning)
727 : {
728 1111 : case ConstraintMeaningEnum::MOLES_BULK_WATER:
729 : {
730 1111 : bulk_moles_old[i] = value;
731 1111 : basis_molality[i] = std::max(
732 1111 : _min_initial_molality,
733 1111 : 0.999 * value /
734 : GeochemistryConstants::MOLES_PER_KG_WATER); // mass of solvent water (water is an
735 : // algebraic variable). Guess used to
736 : // initialize the Newton process
737 1111 : break;
738 : }
739 743 : case ConstraintMeaningEnum::KG_SOLVENT_WATER:
740 : {
741 743 : bulk_moles_old[i] = value * GeochemistryConstants::MOLES_PER_KG_WATER /
742 : 0.999; // initial guess (water is not an algebraic variable). Will be
743 : // determined exactly during the solve
744 743 : basis_molality[i] = value; // mass of solvent water
745 743 : break;
746 : }
747 5660 : case ConstraintMeaningEnum::MOLES_BULK_SPECIES:
748 : {
749 5660 : bulk_moles_old[i] = value;
750 5660 : basis_molality[i] = std::max(
751 5660 : _min_initial_molality,
752 5660 : 0.9 * value / basis_molality[0]); // initial guess (i is an algebraic variable). This
753 : // is what we solve for in the Newton process
754 5660 : break;
755 : }
756 : case ConstraintMeaningEnum::FREE_MOLALITY:
757 : {
758 701 : bulk_moles_old[i] =
759 701 : value * basis_molality[0] / 0.9; // initial guess (i is not an algebraic variable).
760 : // Will be determined exactly during the solve
761 701 : basis_molality[i] = value;
762 701 : break;
763 : }
764 359 : case ConstraintMeaningEnum::FREE_MOLES_MINERAL_SPECIES:
765 : {
766 359 : bulk_moles_old[i] = value / 0.9; // initial guess (i is not an algebraic variable). Will
767 : // be determined exactly during the solve
768 359 : basis_molality[i] = value; // note, this is *moles*, not molality
769 359 : break;
770 : }
771 103 : case ConstraintMeaningEnum::FUGACITY:
772 : {
773 103 : bulk_moles_old[i] = 0.0; // initial guess (i is not an algebraic variable). will be
774 : // determined exactly after the solve
775 103 : basis_molality[i] =
776 : 0.0; // never used in any solve process, but since this is a gas should be zero, and
777 : // setting this explicitly elimiates if(species=gas) checks in various loops
778 103 : break;
779 : }
780 1831 : case ConstraintMeaningEnum::ACTIVITY:
781 : {
782 1831 : bulk_moles_old[i] = value / 0.9; // initial guess (i is not an algebraic variable). Will
783 : // be determined exactly during the solve
784 1831 : if (i == 0)
785 334 : basis_molality[i] = 1.0; // assumption
786 : else
787 1497 : basis_molality[i] =
788 : value /
789 : 0.9; // assume activity_coefficient = 0.9. this will be updated during the solve
790 : break;
791 : }
792 : }
793 : }
794 2188 : }
795 :
796 : Real
797 19 : GeochemicalSystem::getSolventWaterMass() const
798 : {
799 19 : return _basis_molality[0];
800 : }
801 :
802 : const std::vector<Real> &
803 38728 : GeochemicalSystem::getBulkMolesOld() const
804 : {
805 38728 : return _bulk_moles_old;
806 : }
807 :
808 : DenseVector<Real>
809 326 : GeochemicalSystem::getBulkOldInOriginalBasis() const
810 : {
811 326 : DenseVector<Real> result(_bulk_moles_old);
812 326 : if (_mgd.swap_to_original_basis.n() == 0)
813 : return result; // no swaps have been performed
814 156 : _mgd.swap_to_original_basis.vector_mult_transpose(result, _bulk_moles_old);
815 78 : return result;
816 : }
817 :
818 : DenseVector<Real>
819 476 : GeochemicalSystem::getTransportedBulkInOriginalBasis() const
820 : {
821 : std::vector<Real> trans_bulk;
822 476 : computeTransportedBulkFromMolalities(trans_bulk);
823 : DenseVector<Real> result(trans_bulk);
824 476 : if (_mgd.swap_to_original_basis.n() == 0)
825 : return result; // no swaps have been performed
826 364 : _mgd.swap_to_original_basis.vector_mult_transpose(result, trans_bulk);
827 182 : return result;
828 : }
829 :
830 : const std::vector<Real> &
831 75773 : GeochemicalSystem::getSolventMassAndFreeMolalityAndMineralMoles() const
832 : {
833 75773 : return _basis_molality;
834 : }
835 :
836 : void
837 3010 : GeochemicalSystem::buildKnownBasisActivities(std::vector<bool> & basis_activity_known,
838 : std::vector<Real> & basis_activity) const
839 : {
840 3010 : basis_activity_known = std::vector<bool>(_num_basis, false);
841 3010 : basis_activity.resize(_num_basis);
842 :
843 : // all aqueous species with provided activity, and all gases with fugacity have known activity
844 21248 : for (unsigned i = 0; i < _num_basis; ++i)
845 : {
846 18238 : const ConstraintMeaningEnum meaning = _constraint_meaning[i];
847 18238 : if (meaning == ConstraintMeaningEnum::ACTIVITY || meaning == ConstraintMeaningEnum::FUGACITY)
848 : {
849 : basis_activity_known[i] = true;
850 2307 : basis_activity[i] = _constraint_value[i];
851 : }
852 : }
853 :
854 : // all minerals have activity = 1.0
855 21248 : for (unsigned basis_ind = 0; basis_ind < _num_basis; ++basis_ind)
856 18238 : if (_mgd.basis_species_mineral[basis_ind])
857 : {
858 : basis_activity_known[basis_ind] = true;
859 1867 : basis_activity[basis_ind] = 1.0;
860 : }
861 3010 : }
862 :
863 : const std::vector<bool> &
864 2364 : GeochemicalSystem::getBasisActivityKnown() const
865 : {
866 2364 : return _basis_activity_known;
867 : }
868 :
869 : Real
870 44863 : GeochemicalSystem::getBasisActivity(unsigned i) const
871 : {
872 44863 : if (i >= _num_basis)
873 1 : mooseError("Cannot retrieve basis activity for species ",
874 : i,
875 : " since there are only ",
876 1 : _num_basis,
877 : " basis species");
878 44862 : return _basis_activity[i];
879 : }
880 :
881 : const std::vector<Real> &
882 318 : GeochemicalSystem::getBasisActivity() const
883 : {
884 318 : return _basis_activity;
885 : }
886 :
887 : Real
888 247986 : GeochemicalSystem::getEquilibriumMolality(unsigned j) const
889 : {
890 247986 : if (j >= _num_eqm)
891 1 : mooseError("Cannot retrieve molality for equilibrium species ",
892 : j,
893 : " since there are only ",
894 1 : _num_eqm,
895 : " equilibrium species");
896 247985 : return _eqm_molality[j];
897 : }
898 :
899 : const std::vector<Real> &
900 441 : GeochemicalSystem::getEquilibriumMolality() const
901 : {
902 441 : return _eqm_molality;
903 : }
904 :
905 : Real
906 5261 : GeochemicalSystem::getKineticMoles(unsigned k) const
907 : {
908 5261 : if (k >= _num_kin)
909 1 : mooseError("Cannot retrieve moles for kinetic species ",
910 : k,
911 : " since there are only ",
912 1 : _num_kin,
913 : " kinetic species");
914 5260 : return _kin_moles[k];
915 : }
916 :
917 : void
918 299 : GeochemicalSystem::setKineticMoles(unsigned k, Real moles)
919 : {
920 299 : if (k >= _num_kin)
921 1 : mooseError("Cannot set moles for kinetic species ",
922 : k,
923 : " since there are only ",
924 1 : _num_kin,
925 : " kinetic species");
926 298 : if (moles <= 0.0)
927 2 : mooseError("Mole number for kinetic species must be positive, not ", moles);
928 296 : _kin_moles[k] = moles;
929 296 : _kin_moles_old[k] = moles;
930 296 : }
931 :
932 : const std::vector<Real> &
933 480 : GeochemicalSystem::getKineticMoles() const
934 : {
935 480 : return _kin_moles;
936 : }
937 :
938 : void
939 47982 : GeochemicalSystem::computeRemainingBasisActivities(std::vector<Real> & basis_activity) const
940 : {
941 47982 : if (!_basis_activity_known[0])
942 47626 : basis_activity[0] = _gac.waterActivity();
943 323200 : for (unsigned basis_ind = 1; basis_ind < _num_basis; ++basis_ind) // don't loop over water
944 275218 : if (!_basis_activity_known[basis_ind]) // basis_activity_known = true for minerals, gases and
945 : // species with activities provided by the user
946 209979 : basis_activity[basis_ind] = _basis_activity_coef[basis_ind] * _basis_molality[basis_ind];
947 47982 : }
948 :
949 : Real
950 249 : GeochemicalSystem::getEquilibriumActivityCoefficient(unsigned j) const
951 : {
952 249 : if (j >= _num_eqm)
953 1 : mooseError("Cannot retrieve activity coefficient for equilibrium species ",
954 : j,
955 : " since there are only ",
956 1 : _num_eqm,
957 : " equilibrium species");
958 248 : return _eqm_activity_coef[j];
959 : }
960 :
961 : const std::vector<Real> &
962 321 : GeochemicalSystem::getEquilibriumActivityCoefficient() const
963 : {
964 321 : return _eqm_activity_coef;
965 : }
966 :
967 : Real
968 44 : GeochemicalSystem::getBasisActivityCoefficient(unsigned i) const
969 : {
970 44 : if (i >= _num_basis)
971 1 : mooseError("Cannot retrieve basis activity coefficient for species ",
972 : i,
973 : " since there are only ",
974 1 : _num_basis,
975 : " basis species");
976 43 : return _basis_activity_coef[i];
977 : }
978 :
979 : const std::vector<Real> &
980 316 : GeochemicalSystem::getBasisActivityCoefficient() const
981 : {
982 316 : return _basis_activity_coef;
983 : }
984 :
985 : void
986 47982 : GeochemicalSystem::updateBasisMolalityForKnownActivity(std::vector<Real> & basis_molality) const
987 : {
988 323200 : for (unsigned i = 1; i < _num_basis; ++i) // don't loop over water
989 : {
990 275218 : if (_basis_activity_known[i] && !_mgd.basis_species_mineral[i] && !_mgd.basis_species_gas[i])
991 19522 : basis_molality[i] =
992 19522 : _basis_activity[i] / _basis_activity_coef[i]; // just those for
993 : // which activity is provided by the user
994 : }
995 47982 : }
996 :
997 : Real
998 1373827 : GeochemicalSystem::log10ActivityProduct(unsigned eqm_j) const
999 : {
1000 : Real log10ap = 0.0;
1001 18962491 : for (unsigned basis_i = 0; basis_i < _num_basis; ++basis_i)
1002 17588664 : log10ap += _mgd.eqm_stoichiometry(eqm_j, basis_i) * std::log10(_basis_activity[basis_i]);
1003 1373827 : return log10ap;
1004 : }
1005 :
1006 : void
1007 47972 : GeochemicalSystem::computeEqmMolalities(std::vector<Real> & eqm_molality) const
1008 : {
1009 1488049 : for (unsigned eqm_j = 0; eqm_j < _num_eqm; ++eqm_j)
1010 : {
1011 1440077 : if (_mgd.eqm_species_mineral[eqm_j] || _mgd.eqm_species_gas[eqm_j])
1012 70945 : eqm_molality[eqm_j] = 0.0;
1013 : else
1014 : {
1015 : // compute log10 version first, in an attempt to eliminate overflow and underflow problems
1016 : // such as 10^(1000)
1017 1369132 : const Real log10m = log10ActivityProduct(eqm_j) - _eqm_log10K[eqm_j];
1018 1369132 : eqm_molality[eqm_j] =
1019 1369132 : std::pow(10.0, log10m) / _eqm_activity_coef[eqm_j] * surfaceSorptionModifier(eqm_j);
1020 : }
1021 : }
1022 47972 : }
1023 :
1024 : Real
1025 1369132 : GeochemicalSystem::surfaceSorptionModifier(unsigned eqm_j) const
1026 : {
1027 1369132 : if (eqm_j >= _num_eqm)
1028 : return 1.0;
1029 1369132 : if (!_mgd.surface_sorption_related[eqm_j])
1030 : return 1.0;
1031 18714 : return std::pow(_surface_pot_expr[_mgd.surface_sorption_number[eqm_j]],
1032 18714 : 2.0 * _mgd.eqm_species_charge[eqm_j]);
1033 : }
1034 :
1035 : void
1036 54999 : GeochemicalSystem::enforceChargeBalanceIfSimple(std::vector<Real> & constraint_value,
1037 : std::vector<Real> & bulk_moles_old) const
1038 : {
1039 : Real tot_charge = 0.0;
1040 418126 : for (unsigned basis_i = 0; basis_i < _num_basis; ++basis_i)
1041 375298 : if (_mgd.basis_species_charge[basis_i] != 0.0)
1042 : {
1043 225401 : if (_constraint_meaning[basis_i] == ConstraintMeaningEnum::MOLES_BULK_SPECIES)
1044 213230 : tot_charge += _mgd.basis_species_charge[basis_i] * constraint_value[basis_i];
1045 : else
1046 : return;
1047 : }
1048 : // kinetic species are counted in the bulk
1049 : // all charged basis species must have been provided with a MOLES_BULK_SPECIES value, so we can
1050 : // easily enforce charge neutrality
1051 42828 : tot_charge -= _mgd.basis_species_charge[_charge_balance_basis_index] *
1052 : constraint_value[_charge_balance_basis_index];
1053 42828 : constraint_value[_charge_balance_basis_index] =
1054 42828 : -tot_charge / _mgd.basis_species_charge[_charge_balance_basis_index];
1055 42828 : bulk_moles_old[_charge_balance_basis_index] = constraint_value[_charge_balance_basis_index];
1056 : }
1057 :
1058 : Real
1059 4895 : GeochemicalSystem::getTotalChargeOld() const
1060 : {
1061 : Real tot_charge = 0.0;
1062 31663 : for (unsigned basis_i = 0; basis_i < _num_basis; ++basis_i)
1063 26768 : tot_charge += _mgd.basis_species_charge[basis_i] * _bulk_moles_old[basis_i];
1064 : // kinetic species already counted in bulk_moles
1065 4895 : return tot_charge;
1066 : }
1067 :
1068 : void
1069 4570 : GeochemicalSystem::enforceChargeBalance()
1070 : {
1071 4570 : enforceChargeBalance(_constraint_value, _bulk_moles_old);
1072 4570 : }
1073 :
1074 : void
1075 4570 : GeochemicalSystem::enforceChargeBalance(std::vector<Real> & constraint_value,
1076 : std::vector<Real> & bulk_moles_old) const
1077 : {
1078 4570 : const Real tot_charge = getTotalChargeOld();
1079 4570 : constraint_value[_charge_balance_basis_index] -=
1080 4570 : tot_charge / _mgd.basis_species_charge[_charge_balance_basis_index];
1081 4570 : bulk_moles_old[_charge_balance_basis_index] = constraint_value[_charge_balance_basis_index];
1082 4570 : }
1083 :
1084 : void
1085 44338 : GeochemicalSystem::setAlgebraicVariables(const DenseVector<Real> & algebraic_var)
1086 : {
1087 44338 : if (algebraic_var.size() != _num_in_algebraic_system)
1088 1 : mooseError("Incorrect size in setAlgebraicVariables");
1089 274766 : for (unsigned a = 0; a < _num_in_algebraic_system; ++a)
1090 230431 : if (algebraic_var(a) <= 0.0)
1091 2 : mooseError("Cannot set algebraic variables to non-positive values such as ",
1092 : algebraic_var(a));
1093 :
1094 266241 : for (unsigned a = 0; a < _num_basis_in_algebraic_system; ++a)
1095 221906 : _basis_molality[_basis_index[a]] = algebraic_var(a);
1096 47568 : for (unsigned s = 0; s < _num_surface_pot; ++s)
1097 3233 : _surface_pot_expr[s] = algebraic_var(s + _num_basis_in_algebraic_system);
1098 49623 : for (unsigned k = 0; k < _num_kin; ++k)
1099 5288 : _kin_moles[k] = algebraic_var(k + _num_basis_in_algebraic_system + _num_surface_pot);
1100 :
1101 44335 : computeConsistentConfiguration();
1102 44335 : }
1103 :
1104 : void
1105 52515 : GeochemicalSystem::computeBulk(std::vector<Real> & bulk_moles_old) const
1106 : {
1107 400506 : for (unsigned i = 0; i < _num_basis; ++i)
1108 : {
1109 347991 : const Real value = _constraint_value[i];
1110 347991 : const ConstraintMeaningEnum meaning = _constraint_meaning[i];
1111 347991 : switch (meaning)
1112 : {
1113 296370 : case ConstraintMeaningEnum::MOLES_BULK_SPECIES:
1114 : case ConstraintMeaningEnum::MOLES_BULK_WATER:
1115 : {
1116 296370 : bulk_moles_old[i] = value;
1117 296370 : break;
1118 : }
1119 51621 : case ConstraintMeaningEnum::KG_SOLVENT_WATER:
1120 : case ConstraintMeaningEnum::FREE_MOLALITY:
1121 : case ConstraintMeaningEnum::FREE_MOLES_MINERAL_SPECIES:
1122 : case ConstraintMeaningEnum::FUGACITY:
1123 : case ConstraintMeaningEnum::ACTIVITY:
1124 : {
1125 51621 : bulk_moles_old[i] = computeBulkFromMolalities(i);
1126 51621 : break;
1127 : }
1128 : }
1129 : }
1130 52515 : }
1131 :
1132 : Real
1133 72277 : GeochemicalSystem::computeBulkFromMolalities(unsigned basis_ind) const
1134 : {
1135 72277 : const Real nw = _basis_molality[0];
1136 : Real bulk = 0.0;
1137 72277 : if (basis_ind == 0)
1138 : bulk = GeochemistryConstants::MOLES_PER_KG_WATER;
1139 58615 : else if (_mgd.basis_species_mineral[basis_ind])
1140 4596 : bulk = _basis_molality[basis_ind] / nw; // because of multiplication by nw, below
1141 : else
1142 54019 : bulk = _basis_molality[basis_ind];
1143 3182012 : for (unsigned j = 0; j < _num_eqm; ++j)
1144 3109735 : bulk += _mgd.eqm_stoichiometry(j, basis_ind) * _eqm_molality[j];
1145 72277 : bulk *= nw;
1146 79749 : for (unsigned kin = 0; kin < _num_kin; ++kin)
1147 7472 : bulk += _mgd.kin_stoichiometry(kin, basis_ind) * _kin_moles[kin];
1148 72277 : return bulk;
1149 : }
1150 :
1151 : void
1152 480 : GeochemicalSystem::computeTransportedBulkFromMolalities(std::vector<Real> & transported_bulk) const
1153 : {
1154 480 : transported_bulk.resize(_num_basis);
1155 480 : const Real nw = _basis_molality[0];
1156 3521 : for (unsigned i = 0; i < _num_basis; ++i)
1157 : {
1158 3041 : transported_bulk[i] = 0.0;
1159 3041 : if (i == 0)
1160 480 : transported_bulk[i] = GeochemistryConstants::MOLES_PER_KG_WATER;
1161 2561 : else if (_mgd.basis_species_transported[i])
1162 : {
1163 2337 : if (_mgd.basis_species_mineral[i])
1164 0 : transported_bulk[i] = _basis_molality[i] / nw; // because of multiplication by nw, below
1165 : else
1166 2337 : transported_bulk[i] = _basis_molality[i];
1167 : }
1168 : else
1169 : transported_bulk[i] = 0.0;
1170 103789 : for (unsigned j = 0; j < _num_eqm; ++j)
1171 100748 : if (_mgd.eqm_species_transported[j])
1172 95002 : transported_bulk[i] += _mgd.eqm_stoichiometry(j, i) * _eqm_molality[j];
1173 3041 : transported_bulk[i] *= nw;
1174 : }
1175 480 : }
1176 :
1177 : Real
1178 253617 : GeochemicalSystem::getResidualComponent(unsigned algebraic_ind,
1179 : const DenseVector<Real> & mole_additions) const
1180 : {
1181 253617 : if (algebraic_ind >= _num_in_algebraic_system)
1182 2 : mooseError("Cannot retrieve residual for algebraic index ",
1183 : algebraic_ind,
1184 : " because there are only ",
1185 2 : _num_basis_in_algebraic_system,
1186 : " molalities in the algebraic system and ",
1187 2 : _num_surface_pot,
1188 : " surface potentials and ",
1189 2 : _num_kin,
1190 : " kinetic species");
1191 253615 : if (mole_additions.size() != _num_basis + _num_kin)
1192 1 : mooseError("The increment in mole numbers (mole_additions) needs to be of size ",
1193 1 : _num_basis,
1194 : " + ",
1195 1 : _num_kin,
1196 : " but it is of size ",
1197 1 : mole_additions.size());
1198 :
1199 253614 : if (algebraic_ind < _num_basis_in_algebraic_system) // residual for basis molality
1200 : {
1201 : /*
1202 : * Without the special things for water or the charge-balance species, we're trying to solve
1203 : * bulk_new = bulk_old + mole_addition
1204 : * where bulk_old is known (from previous time-step or, for a steady-state problem, the
1205 : * constraint)
1206 : * and mole_addition is given to this function
1207 : * and bulk_new = nw * (m + sum_eqm[stoi * mol_eqm]) + sum_kin[stoi * mole_kin]
1208 : * Hence, the residual is
1209 : * R = -(bulk_old + mole_addition) + nw*(m + sum_eqm[stoi*mol_eqm]) + sum_kin[stoi*mole_kin]
1210 : * This is seemingly different from Bethke Eqns(16.7)-(16.9), because Bethke bulk ("M") do not
1211 : contain the sum_kin[stoi * mol_kin] terms. Converting the above residual to Bethke's
1212 : convention:
1213 : * R = -((nw*(m + sum_eqm[stoi*mol_eqm]) + sum_kin[stoi*mole_kin])_old + mole_addition) + nw*(m
1214 : + sum_eqm[stoi*mol_eqm]) + sum_kin[stoi*mole_kin]
1215 : * = -(M_old + sum_kin[stoi*mole_kin_old] + mole_addition) + M + sum_kin[stoi*mole_kin]
1216 : * = -(M_old + mole_addition) + M + sum_kin[stoi*kin_mole_addition]
1217 : * which is exactly Eqns(16.7)-(16.9).
1218 : * One problem with Bethke's approach is that if kin_mole_addition depends on the mole number
1219 : of the kinetic species, there is a "hidden" variable (moles of kinetic species) which is kept
1220 : constant during the solve. Here, including the moles of kinetic species as additional
1221 : variables allows greater accuracy.
1222 : */
1223 243801 : const unsigned basis_i = _basis_index[algebraic_ind];
1224 : Real res = 0.0;
1225 243801 : if (basis_i == 0)
1226 36884 : res += -(_bulk_moles_old[basis_i] + mole_additions(basis_i)) +
1227 36884 : _basis_molality[0] * GeochemistryConstants::MOLES_PER_KG_WATER;
1228 206917 : else if (basis_i == _charge_balance_basis_index)
1229 : {
1230 49229 : res += _basis_molality[0] * _basis_molality[basis_i];
1231 380244 : for (unsigned i = 0; i < _num_basis; ++i)
1232 : {
1233 331015 : if (i == _charge_balance_basis_index)
1234 49229 : continue;
1235 281786 : else if (_mgd.basis_species_charge[i] ==
1236 : 0.0) // certainly includes water, minerals and gases
1237 119189 : continue;
1238 162597 : else if (_constraint_meaning[i] == ConstraintMeaningEnum::MOLES_BULK_SPECIES)
1239 : {
1240 : // the molalities might not yet have converged so that the bulk moles of this species
1241 : // currently equals _bulk_moles_old + mole_additions, but we know
1242 : // that when the solve has converged it'll have to, so:
1243 142386 : res += _mgd.basis_species_charge[i] * (_bulk_moles_old[i] + mole_additions(i)) /
1244 : _mgd.basis_species_charge[basis_i];
1245 : }
1246 : else
1247 : {
1248 : // do not know the bulk moles of this species: use the current value for molality and
1249 : // kin_moles. Note, there is no mole_additions here since physically any mole_additions
1250 : // get instantly get soaked up by the fixed activity (or fixed molality, etc), and
1251 : // mathematically when we eventually come to add mole_additions to the bulk_moles_old
1252 : // (via addtoBulkMoles, for instance) we immediately return without adding stuff. End
1253 : // result: charge-neutrality should not depend on mole_additions for fixed-activity
1254 : // (molality, etc) species.
1255 20211 : res += _mgd.basis_species_charge[i] * computeBulkFromMolalities(i) /
1256 20211 : _mgd.basis_species_charge[basis_i];
1257 : }
1258 : }
1259 : }
1260 : else
1261 157688 : res += -(_bulk_moles_old[basis_i] + mole_additions(basis_i)) +
1262 157688 : _basis_molality[0] * _basis_molality[basis_i];
1263 14338926 : for (unsigned eqm_j = 0; eqm_j < _num_eqm; ++eqm_j)
1264 14095125 : res += _basis_molality[0] * _mgd.eqm_stoichiometry(eqm_j, basis_i) * _eqm_molality[eqm_j];
1265 269226 : for (unsigned kin = 0; kin < _num_kin; ++kin)
1266 25425 : res += _mgd.kin_stoichiometry(kin, basis_i) * _kin_moles[kin];
1267 243801 : return res;
1268 : }
1269 9813 : else if (algebraic_ind <
1270 9813 : _num_basis_in_algebraic_system + _num_surface_pot) // residual for surface potential
1271 : {
1272 3482 : const unsigned sp = algebraic_ind - _num_basis_in_algebraic_system;
1273 3482 : Real res = surfacePotPrefactor(sp) * (_surface_pot_expr[sp] - 1.0 / _surface_pot_expr[sp]);
1274 44145 : for (unsigned j = 0; j < _num_eqm; ++j)
1275 40663 : if (_mgd.surface_sorption_related[j] && _mgd.surface_sorption_number[j] == sp)
1276 18569 : res += _basis_molality[0] * _mgd.eqm_species_charge[j] * _eqm_molality[j];
1277 3482 : return res;
1278 : }
1279 :
1280 : // else: residual for kinetic
1281 6331 : const unsigned kin = algebraic_ind - _num_basis_in_algebraic_system - _num_surface_pot;
1282 6331 : Real res = _kin_moles[kin] - (_kin_moles_old[kin] + mole_additions(_num_basis + kin));
1283 6331 : return res;
1284 : }
1285 :
1286 : void
1287 40594 : GeochemicalSystem::computeJacobian(const DenseVector<Real> & res,
1288 : DenseMatrix<Real> & jac,
1289 : const DenseVector<Real> & mole_additions,
1290 : const DenseMatrix<Real> & dmole_additions) const
1291 : {
1292 : /* To the reader: yes, this is awfully complicated. The best way of understanding it is to very
1293 : * slowly go through the residual calculation and compute the derivatives by hand. It is quite
1294 : * well tested, but it'd be great if you could add more tests!
1295 : */
1296 40594 : if (res.size() != _num_in_algebraic_system)
1297 1 : mooseError(
1298 1 : "Jacobian: residual size must be ", _num_in_algebraic_system, " but it is ", res.size());
1299 40593 : if (mole_additions.size() != _num_basis + _num_kin)
1300 1 : mooseError("Jacobian: the increment in mole numbers (mole_additions) needs to be of size ",
1301 1 : _num_basis,
1302 : " + ",
1303 1 : _num_kin,
1304 : " but it is of size ",
1305 1 : mole_additions.size());
1306 40592 : if (!(dmole_additions.n() == _num_basis + _num_kin &&
1307 : dmole_additions.m() == _num_basis + _num_kin))
1308 2 : mooseError("Jacobian: the derivative of mole additions (dmole_additions) needs to be of size ",
1309 2 : _num_basis + _num_kin,
1310 : "x",
1311 2 : _num_basis + _num_kin,
1312 : " but it is of size ",
1313 2 : dmole_additions.m(),
1314 : "x",
1315 2 : dmole_additions.n());
1316 : /*
1317 : Note that in this function we make use of the fact that species can only be in the algebraic
1318 : system if their molality is unknown (or in the case of water, the mass of solvent mass is
1319 : unknown). This means that the molality of the equilibrium species depends on
1320 : (activity_coefficient * basis molality)^stoi, so that the derivative is quite simple. Eg, we
1321 : don't have to worry about derivatives with respect to fixed-activity things, or fixed fugacity.
1322 : Also, note that the constructor and the various "set" methods of this class enforce molality > 0,
1323 : so there are no division-by-zero problems. Also, note that we never compute derivatives of the
1324 : activity coefficients, or the activity of water with respect to the molalities, as they are
1325 : assumed to be quite small. Finally, the surface_pot_expr will also always be positive.
1326 : */
1327 : // correctly size and zero
1328 40590 : jac.resize(_num_in_algebraic_system, _num_in_algebraic_system);
1329 40590 : const Real nw = _basis_molality[0];
1330 :
1331 : // jac(a, b) = d(R_a) / d(m_b), where a corresponds to a molality
1332 225465 : for (unsigned a = 0; a < _num_basis_in_algebraic_system; ++a)
1333 : {
1334 184875 : const unsigned basis_of_a = _basis_index[a];
1335 1476670 : for (unsigned b = 0; b < _num_basis_in_algebraic_system; ++b)
1336 : {
1337 1291795 : const unsigned basis_of_b = _basis_index[b];
1338 :
1339 : // contribution from mole_additions
1340 1291795 : if (basis_of_a != _charge_balance_basis_index)
1341 1106920 : jac(a, b) -= dmole_additions(basis_of_a, basis_of_b);
1342 : else
1343 : {
1344 1859238 : for (unsigned i = 0; i < _num_basis; ++i)
1345 : {
1346 1674363 : if (i == _charge_balance_basis_index)
1347 184875 : continue;
1348 1489488 : else if (_mgd.basis_species_charge[i] ==
1349 : 0.0) // certainly includes water, minerals and gases
1350 540313 : continue;
1351 949175 : else if (_constraint_meaning[i] == ConstraintMeaningEnum::MOLES_BULK_SPECIES)
1352 887845 : jac(a, b) += _mgd.basis_species_charge[i] * dmole_additions(i, basis_of_b) /
1353 887845 : _mgd.basis_species_charge[basis_of_a];
1354 : }
1355 : }
1356 :
1357 : // contribution from explicit nw, in case where mass of solvent water is an unknown
1358 1291795 : if (basis_of_b == 0)
1359 : {
1360 122647 : if (basis_of_a != _charge_balance_basis_index)
1361 : {
1362 : // use a short-cut here: dR/dnw = m + sum_eqm[stoi*mol] = (R + bulk + addition - kin)/nw
1363 90098 : Real numerator = res(a) + _bulk_moles_old[basis_of_a] + mole_additions(basis_of_a);
1364 102685 : for (unsigned kin = 0; kin < _num_kin; ++kin)
1365 12587 : numerator -= _mgd.kin_stoichiometry(kin, basis_of_a) * _kin_moles[kin];
1366 90098 : jac(a, 0) += numerator / nw;
1367 : }
1368 : else
1369 : {
1370 198478 : for (unsigned i = 0; i < _num_basis; ++i)
1371 : {
1372 165929 : if (i == _charge_balance_basis_index)
1373 32549 : continue;
1374 133380 : else if (_mgd.basis_species_charge[i] ==
1375 : 0.0) // certainly includes water, minerals and gases
1376 70422 : continue;
1377 62958 : else if (_constraint_meaning[i] == ConstraintMeaningEnum::MOLES_BULK_SPECIES)
1378 52336 : continue;
1379 : else
1380 : {
1381 10622 : Real molal = _basis_molality[i];
1382 49951 : for (unsigned j = 0; j < _num_eqm; ++j)
1383 39329 : molal += _mgd.eqm_stoichiometry(j, i) * _eqm_molality[j];
1384 10622 : jac(a, 0) +=
1385 10622 : _mgd.basis_species_charge[i] * molal / _mgd.basis_species_charge[basis_of_a];
1386 : }
1387 : }
1388 32549 : jac(a, 0) += _basis_molality[basis_of_a];
1389 511886 : for (unsigned j = 0; j < _num_eqm; ++j)
1390 479337 : jac(a, 0) += _mgd.eqm_stoichiometry(j, basis_of_a) * _eqm_molality[j];
1391 : }
1392 : }
1393 :
1394 : // contribution from molality of basis_of_b
1395 : if (basis_of_b != 0)
1396 : {
1397 1169148 : if (a == b)
1398 152326 : jac(a, b) += nw; // remember b != 0
1399 104175398 : for (unsigned eqm_j = 0; eqm_j < _num_eqm; ++eqm_j)
1400 : {
1401 103006250 : jac(a, b) += nw * _mgd.eqm_stoichiometry(eqm_j, basis_of_a) * _eqm_molality[eqm_j] *
1402 103006250 : _mgd.eqm_stoichiometry(eqm_j, basis_of_b) / _basis_molality[basis_of_b];
1403 : }
1404 1169148 : if (basis_of_a == _charge_balance_basis_index)
1405 : {
1406 : // additional terms from the charge-balance additions
1407 1660760 : for (unsigned i = 0; i < _num_basis; ++i)
1408 : {
1409 1508434 : if (i == _charge_balance_basis_index)
1410 152326 : continue;
1411 1356108 : else if (_mgd.basis_species_charge[i] ==
1412 : 0.0) // certainly includes water, minerals and gases
1413 469891 : continue;
1414 886217 : else if (_constraint_meaning[i] == ConstraintMeaningEnum::MOLES_BULK_SPECIES)
1415 835509 : continue;
1416 : else
1417 : {
1418 : const Real prefactor =
1419 50708 : _mgd.basis_species_charge[i] * nw / _mgd.basis_species_charge[basis_of_a];
1420 50708 : if (i == basis_of_b)
1421 0 : jac(a, b) += prefactor;
1422 3068814 : for (unsigned eqm_j = 0; eqm_j < _num_eqm; ++eqm_j)
1423 : {
1424 3018106 : jac(a, b) += prefactor * _mgd.eqm_stoichiometry(eqm_j, i) * _eqm_molality[eqm_j] *
1425 3018106 : _mgd.eqm_stoichiometry(eqm_j, basis_of_b) /
1426 3018106 : _basis_molality[basis_of_b];
1427 : }
1428 : }
1429 : }
1430 : }
1431 : }
1432 : }
1433 : }
1434 :
1435 : // jac(a, b) = d(R_a) / d(surface_pot), where a corresponds to a molality
1436 225465 : for (unsigned a = 0; a < _num_basis_in_algebraic_system; ++a)
1437 : {
1438 184875 : const unsigned basis_of_a = _basis_index[a];
1439 198812 : for (unsigned s = 0; s < _num_surface_pot; ++s)
1440 : {
1441 13937 : const unsigned b = s + _num_basis_in_algebraic_system;
1442 : // derivative of nw * _mgd.eqm_stoichiometry(eqm_j, basis_i) * _eqm_molality[eqm_j],
1443 : // where _eqm_molality = (_surface_pot_expr)^(2 * charge) * stuff
1444 200875 : for (unsigned eqm_j = 0; eqm_j < _num_eqm; ++eqm_j)
1445 186938 : if (_mgd.surface_sorption_related[eqm_j] && _mgd.surface_sorption_number[eqm_j] == s)
1446 82130 : jac(a, b) += nw * _mgd.eqm_stoichiometry(eqm_j, basis_of_a) * 2.0 *
1447 82130 : _mgd.eqm_species_charge[eqm_j] * _eqm_molality[eqm_j] /
1448 82130 : _surface_pot_expr[_mgd.surface_sorption_number[eqm_j]];
1449 : }
1450 : }
1451 :
1452 : // jac(a, b) = d(R_a) / d(_kin_moles) where a corresponds to a molality
1453 225465 : for (unsigned a = 0; a < _num_basis_in_algebraic_system; ++a)
1454 : {
1455 184875 : const unsigned basis_of_a = _basis_index[a];
1456 :
1457 : // contribution from mole_additions
1458 205764 : for (unsigned kin = 0; kin < _num_kin; ++kin)
1459 : {
1460 20889 : const unsigned ind = _num_basis + kin;
1461 20889 : const unsigned b = _num_basis_in_algebraic_system + _num_surface_pot + kin;
1462 20889 : if (basis_of_a != _charge_balance_basis_index)
1463 15638 : jac(a, b) -= dmole_additions(basis_of_a, ind);
1464 : else
1465 : {
1466 47935 : for (unsigned i = 0; i < _num_basis; ++i)
1467 : {
1468 42684 : if (i == _charge_balance_basis_index)
1469 5251 : continue;
1470 37433 : else if (_mgd.basis_species_charge[i] ==
1471 : 0.0) // certainly includes water, minerals and gases
1472 26399 : continue;
1473 11034 : else if (_constraint_meaning[i] == ConstraintMeaningEnum::MOLES_BULK_SPECIES)
1474 9466 : jac(a, b) += _mgd.basis_species_charge[i] * dmole_additions(i, ind) /
1475 9466 : _mgd.basis_species_charge[basis_of_a];
1476 : }
1477 : }
1478 : }
1479 :
1480 : // contribution from sum_kin[stoi*kin_moles]
1481 205764 : for (unsigned kin = 0; kin < _num_kin; ++kin)
1482 : {
1483 20889 : const unsigned b = _num_basis_in_algebraic_system + _num_surface_pot + kin;
1484 20889 : jac(a, b) += _mgd.kin_stoichiometry(kin, basis_of_a);
1485 : }
1486 :
1487 : // special additional contribution for charge-balance from kinetic stuff
1488 184875 : if (basis_of_a == _charge_balance_basis_index)
1489 : {
1490 295338 : for (unsigned i = 0; i < _num_basis; ++i)
1491 : {
1492 254748 : if (i == _charge_balance_basis_index)
1493 40590 : continue;
1494 214158 : else if (_mgd.basis_species_charge[i] ==
1495 : 0.0) // certainly includes water, minerals and gases
1496 98082 : continue;
1497 116076 : else if (_constraint_meaning[i] == ConstraintMeaningEnum::MOLES_BULK_SPECIES)
1498 100370 : continue;
1499 : else
1500 : {
1501 : const Real prefactor =
1502 15706 : _mgd.basis_species_charge[i] / _mgd.basis_species_charge[basis_of_a];
1503 17274 : for (unsigned kin = 0; kin < _num_kin; ++kin)
1504 : {
1505 1568 : const unsigned b = _num_basis_in_algebraic_system + _num_surface_pot + kin;
1506 1568 : jac(a, b) += prefactor * _mgd.kin_stoichiometry(kin, i);
1507 : }
1508 : }
1509 : }
1510 : }
1511 : }
1512 :
1513 : // jac(a, b) = d(R_a) / d(variable_b) where a corresponds to a surface potential
1514 43662 : for (unsigned s = 0; s < _num_surface_pot; ++s)
1515 : {
1516 3072 : const unsigned a = s + _num_basis_in_algebraic_system;
1517 :
1518 17009 : for (unsigned b = 0; b < _num_basis_in_algebraic_system; ++b)
1519 : {
1520 : // derivative of nw * _mgd.eqm_species_charge[j] * _eqm_molality[j];
1521 13937 : const unsigned basis_of_b = _basis_index[b];
1522 13937 : if (basis_of_b == 0) // special case: mass of solvent water is an unknown
1523 : {
1524 14588 : for (unsigned j = 0; j < _num_eqm; ++j)
1525 13125 : if (_mgd.surface_sorption_related[j] && _mgd.surface_sorption_number[j] == s)
1526 5843 : jac(a, b) += _mgd.eqm_species_charge[j] * _eqm_molality[j];
1527 : }
1528 : else
1529 : {
1530 186287 : for (unsigned j = 0; j < _num_eqm; ++j)
1531 173813 : if (_mgd.surface_sorption_related[j] && _mgd.surface_sorption_number[j] == s)
1532 76287 : jac(a, b) += nw * _mgd.eqm_species_charge[j] * _eqm_molality[j] *
1533 76287 : _mgd.eqm_stoichiometry(j, basis_of_b) / _basis_molality[basis_of_b];
1534 : }
1535 : }
1536 :
1537 3072 : const Real coef = surfacePotPrefactor(s);
1538 : // derivative of coef * (x - 1/x) wrt x, where x = _surface_pot_expr
1539 3072 : jac(a, a) += coef * (1.0 + 1.0 / std::pow(_surface_pot_expr[s], 2.0));
1540 : // derivative of nw * _mgd.eqm_species_charge[j] * _eqm_molality[j];
1541 : // where _eqm_molality = (_surface_pot_expr)^(2 * charge) * stuff
1542 37957 : for (unsigned j = 0; j < _num_eqm; ++j)
1543 34885 : if (_mgd.surface_sorption_related[j] && _mgd.surface_sorption_number[j] == s)
1544 16053 : jac(a, a) += nw * std::pow(_mgd.eqm_species_charge[j], 2.0) * 2.0 * _eqm_molality[j] /
1545 : _surface_pot_expr[s];
1546 : }
1547 :
1548 : // jac(a, b) = d(R_a) / d(variable_b) where a corresponds to a kinetic species
1549 45841 : for (unsigned kin = 0; kin < _num_kin; ++kin)
1550 : {
1551 5251 : const unsigned a = _num_basis_in_algebraic_system + _num_surface_pot + kin;
1552 5251 : jac(a, a) += 1.0; // deriv wrt _kin_moles[kin]
1553 :
1554 : // contribution from mole_additions
1555 5251 : const unsigned ind_of_addition = _num_basis + kin;
1556 26140 : for (unsigned b = 0; b < _num_basis_in_algebraic_system; ++b)
1557 : {
1558 20889 : const unsigned basis_of_b = _basis_index[b];
1559 20889 : jac(a, b) -= dmole_additions(ind_of_addition, basis_of_b);
1560 : }
1561 10522 : for (unsigned kinp = 0; kinp < _num_kin; ++kinp)
1562 : {
1563 5271 : const unsigned b = _num_basis_in_algebraic_system + _num_surface_pot + kinp;
1564 5271 : const unsigned indp = _num_basis + kinp;
1565 5271 : jac(a, b) -= dmole_additions(ind_of_addition, indp);
1566 : }
1567 : }
1568 40590 : }
1569 :
1570 : const ModelGeochemicalDatabase &
1571 664963 : GeochemicalSystem::getModelGeochemicalDatabase() const
1572 : {
1573 664963 : return _mgd;
1574 : }
1575 :
1576 : void
1577 100322 : GeochemicalSystem::computeFreeMineralMoles(std::vector<Real> & basis_molality) const
1578 : {
1579 :
1580 100322 : const Real nw = _basis_molality[0];
1581 846125 : for (unsigned i = 0; i < _num_basis; ++i)
1582 745803 : if (_mgd.basis_species_mineral[i])
1583 : {
1584 122448 : if (_constraint_meaning[i] == ConstraintMeaningEnum::FREE_MOLES_MINERAL_SPECIES)
1585 6279 : basis_molality[i] = _constraint_value[i];
1586 : else
1587 : {
1588 116169 : basis_molality[i] = _bulk_moles_old[i];
1589 6611004 : for (unsigned j = 0; j < _num_eqm; ++j)
1590 6494835 : basis_molality[i] -= nw * _mgd.eqm_stoichiometry(j, i) * _eqm_molality[j];
1591 153625 : for (unsigned kin = 0; kin < _num_kin; ++kin)
1592 37456 : basis_molality[i] -= _mgd.kin_stoichiometry(kin, i) * _kin_moles[kin];
1593 : }
1594 : }
1595 100322 : }
1596 :
1597 : std::vector<Real>
1598 5128 : GeochemicalSystem::getSaturationIndices() const
1599 : {
1600 5128 : std::vector<Real> si(_num_eqm);
1601 100330 : for (unsigned j = 0; j < _num_eqm; ++j)
1602 95202 : if (_mgd.eqm_species_mineral[j])
1603 4695 : si[j] = log10ActivityProduct(j) - _eqm_log10K[j];
1604 : else
1605 90507 : si[j] = 0.0;
1606 5128 : return si;
1607 : }
1608 :
1609 : void
1610 304 : GeochemicalSystem::performSwap(unsigned swap_out_of_basis, unsigned swap_into_basis)
1611 : {
1612 304 : if (swap_out_of_basis == 0)
1613 1 : mooseException("GeochemicalSystem: attempting to swap out water and replace it by ",
1614 : _mgd.eqm_species_name[swap_into_basis],
1615 : ". This could be because the algorithm would like to "
1616 : "swap out the charge-balance species, in which case you should choose a "
1617 : "different charge-balance species");
1618 303 : if (swap_out_of_basis == _charge_balance_basis_index)
1619 1 : mooseException("GeochemicalSystem: attempting to swap the charge-balance species out of "
1620 : "the basis");
1621 302 : if (_mgd.basis_species_gas[swap_out_of_basis])
1622 1 : mooseException("GeochemicalSystem: attempting to swap a gas out of the basis");
1623 301 : if (_mgd.eqm_species_gas[swap_into_basis])
1624 1 : mooseException("GeochemicalSystem: attempting to swap a gas into the basis");
1625 : // perform the swap
1626 300 : performSwapNoCheck(swap_out_of_basis, swap_into_basis);
1627 300 : }
1628 :
1629 : void
1630 406 : GeochemicalSystem::performSwapNoCheck(unsigned swap_out_of_basis, unsigned swap_into_basis)
1631 : {
1632 406 : DenseVector<Real> bm = _bulk_moles_old;
1633 406 : _swapper.performSwap(_mgd, bm, swap_out_of_basis, swap_into_basis);
1634 :
1635 : // the swap_into_basis species now has fixed bulk moles irrespective of what the
1636 : // swap_out_of_basis species had fixed
1637 406 : _constraint_meaning[swap_out_of_basis] = ConstraintMeaningEnum::MOLES_BULK_SPECIES;
1638 :
1639 : // the bulk moles will have changed for many components. The _bulk_moles_old values will be
1640 : // set in computeConsistentConfiguration, below
1641 5472 : for (unsigned basis_i = 0; basis_i < _num_basis; ++basis_i)
1642 5066 : if (_constraint_meaning[basis_i] == ConstraintMeaningEnum::MOLES_BULK_SPECIES ||
1643 : _constraint_meaning[basis_i] == ConstraintMeaningEnum::MOLES_BULK_WATER)
1644 : {
1645 4515 : _constraint_value[basis_i] = bm(basis_i);
1646 4515 : _original_constraint_value[basis_i] = bm(basis_i);
1647 : }
1648 :
1649 : // In the following, the molalities are carefully swapped so that the configuration stays
1650 : // reasonably close to the configuration prior to the swap. This may be advantageous as it
1651 : // initializes the Newton process with a better starting guess than what might have been used
1652 : Real molality_of_species_just_swapped_in =
1653 406 : _eqm_molality[swap_into_basis]; // is positive, or zero iff (mineral or gas)
1654 406 : if (_mgd.basis_species_mineral[swap_out_of_basis] || _mgd.basis_species_gas[swap_out_of_basis] ||
1655 : _eqm_molality[swap_into_basis] == 0.0)
1656 : {
1657 : // the species just swapped in is a mineral or a gas, so its equilibium molality is
1658 : // undefined: make a guess for its molality
1659 247 : molality_of_species_just_swapped_in =
1660 431 : std::max(_min_initial_molality, 0.9 * bm(swap_out_of_basis));
1661 : }
1662 : Real molality_of_species_just_swapped_out =
1663 406 : _basis_molality[swap_out_of_basis]; // can be negative if a consumed mineral
1664 406 : _basis_molality[swap_out_of_basis] = molality_of_species_just_swapped_in;
1665 406 : _eqm_molality[swap_into_basis] =
1666 406 : std::max(0.0, molality_of_species_just_swapped_out); // this gets recalculated in
1667 : // computeConsistentConfiguration, below
1668 :
1669 : // depending if minerals were swapped in or out of the basis, the known activity may have
1670 : // changed
1671 406 : buildKnownBasisActivities(_basis_activity_known, _basis_activity);
1672 :
1673 : // the equilibrium constants will have changed due to the swap
1674 406 : buildTemperatureDependentQuantities(_temperature);
1675 :
1676 : // due to re-orderings in mgd, the charge-balance basis index might have changed
1677 406 : _charge_balance_basis_index = _mgd.basis_species_index.at(_charge_balance_species);
1678 : // charge balance might be able to be performed easily
1679 406 : enforceChargeBalanceIfSimple(_constraint_value, _bulk_moles_old);
1680 :
1681 : // the algebraic system has probably changed
1682 406 : buildAlgebraicInfo(_in_algebraic_system,
1683 406 : _num_basis_in_algebraic_system,
1684 406 : _num_in_algebraic_system,
1685 406 : _algebraic_index,
1686 406 : _basis_index);
1687 :
1688 : // finally compute a consistent configuration, based on the basis molalities, etc, above
1689 406 : computeConsistentConfiguration();
1690 406 : }
1691 :
1692 : unsigned
1693 2389 : GeochemicalSystem::getNumInBasis() const
1694 : {
1695 2389 : return _num_basis;
1696 : }
1697 :
1698 : unsigned
1699 528 : GeochemicalSystem::getNumInEquilibrium() const
1700 : {
1701 528 : return _num_eqm;
1702 : }
1703 :
1704 : void
1705 9 : GeochemicalSystem::setChargeBalanceSpecies(unsigned new_charge_balance_index)
1706 : {
1707 : // because the original mole number of the charge balance species may have been changed due to
1708 : // enforcing charge balance:
1709 9 : _constraint_value[_charge_balance_basis_index] =
1710 9 : _original_constraint_value[_charge_balance_basis_index];
1711 9 : _bulk_moles_old[_charge_balance_basis_index] = _constraint_value[_charge_balance_basis_index];
1712 : // now change the charge balance info
1713 9 : _charge_balance_basis_index = new_charge_balance_index;
1714 9 : _charge_balance_species = _mgd.basis_species_name[_charge_balance_basis_index];
1715 : // enforce charge-balance if easily possible
1716 9 : enforceChargeBalanceIfSimple(_constraint_value, _bulk_moles_old);
1717 9 : }
1718 :
1719 : bool
1720 40580 : GeochemicalSystem::alterChargeBalanceSpecies(Real threshold_molality)
1721 : {
1722 40580 : if (_basis_molality[_charge_balance_basis_index] > threshold_molality)
1723 : return false;
1724 : unsigned best_species_opp_sign = 0;
1725 : unsigned best_species_same_sign = 0;
1726 : Real best_molality_opp_sign = 0.0;
1727 : Real best_molality_same_sign = 0.0;
1728 1965 : for (unsigned basis_i = 0; basis_i < _num_basis; ++basis_i)
1729 : {
1730 1761 : if (basis_i == _charge_balance_basis_index)
1731 204 : continue;
1732 1557 : else if (_constraint_meaning[basis_i] != ConstraintMeaningEnum::MOLES_BULK_SPECIES)
1733 490 : continue;
1734 1067 : else if (_mgd.basis_species_charge[basis_i] == 0.0)
1735 203 : continue;
1736 864 : else if (_basis_molality[basis_i] <= threshold_molality)
1737 772 : continue;
1738 : // we know basis_i is a charged species with set bulk moles and high molality
1739 92 : if (_mgd.basis_species_charge[basis_i] *
1740 : _mgd.basis_species_charge[_charge_balance_basis_index] <
1741 : 0.0)
1742 : {
1743 : // charge of opposite sign
1744 78 : if (_basis_molality[basis_i] > best_molality_opp_sign)
1745 : {
1746 : best_species_opp_sign = basis_i;
1747 : best_molality_opp_sign = _basis_molality[basis_i];
1748 : }
1749 : }
1750 : else
1751 : {
1752 : // charge of same sign
1753 14 : if (_basis_molality[basis_i] > best_molality_same_sign)
1754 : {
1755 : best_species_same_sign = basis_i;
1756 : best_molality_same_sign = _basis_molality[basis_i];
1757 : }
1758 : }
1759 : }
1760 204 : if (best_species_opp_sign != 0)
1761 : {
1762 : // this is preferred over the same-sign version
1763 8 : setChargeBalanceSpecies(best_species_opp_sign);
1764 8 : return true;
1765 : }
1766 196 : else if (best_species_same_sign != 0)
1767 : {
1768 : // this is not preferred, but is better than no charge-balance species change
1769 0 : setChargeBalanceSpecies(best_species_same_sign);
1770 0 : return true;
1771 : }
1772 : return false;
1773 : }
1774 :
1775 : bool
1776 1 : GeochemicalSystem::revertToOriginalChargeBalanceSpecies()
1777 : {
1778 1 : if (_mgd.basis_species_index.find(_original_charge_balance_species) ==
1779 1 : _mgd.basis_species_index.end())
1780 : return false; // original charge-balance species no longer in basis
1781 1 : const unsigned original_index = _mgd.basis_species_index.at(_original_charge_balance_species);
1782 1 : if (original_index == _charge_balance_basis_index)
1783 : return false; // current charge-balance species is the original
1784 1 : setChargeBalanceSpecies(original_index);
1785 1 : return true;
1786 : }
1787 :
1788 : Real
1789 321 : GeochemicalSystem::getIonicStrength() const
1790 : {
1791 321 : return _is.ionicStrength(_mgd, _basis_molality, _eqm_molality, _kin_moles);
1792 : }
1793 :
1794 : Real
1795 316 : GeochemicalSystem::getStoichiometricIonicStrength() const
1796 : {
1797 316 : return _is.stoichiometricIonicStrength(_mgd, _basis_molality, _eqm_molality, _kin_moles);
1798 : }
1799 :
1800 : Real
1801 7022 : GeochemicalSystem::surfacePotPrefactor(unsigned sp) const
1802 : {
1803 7022 : return 0.5 * _sorbing_surface_area[sp] / GeochemistryConstants::FARADAY *
1804 14044 : std::sqrt(8.0 * GeochemistryConstants::GAS_CONSTANT *
1805 7022 : (_temperature + GeochemistryConstants::CELSIUS_TO_KELVIN) *
1806 7022 : GeochemistryConstants::PERMITTIVITY_FREE_SPACE *
1807 7022 : GeochemistryConstants::DIELECTRIC_CONSTANT_WATER *
1808 : GeochemistryConstants::DENSITY_WATER *
1809 7022 : _is.ionicStrength(_mgd, _basis_molality, _eqm_molality, _kin_moles));
1810 : }
1811 :
1812 : Real
1813 482 : GeochemicalSystem::getSurfacePotential(unsigned sp) const
1814 : {
1815 482 : if (sp >= _num_surface_pot)
1816 2 : mooseError("Cannot retrieve the surface potential for surface ",
1817 : sp,
1818 : " since there are only ",
1819 2 : _num_surface_pot,
1820 : " surfaces involved in surface complexation");
1821 480 : return -2.0 * GeochemistryConstants::GAS_CONSTANT *
1822 480 : (_temperature + GeochemistryConstants::CELSIUS_TO_KELVIN) /
1823 480 : GeochemistryConstants::FARADAY * std::log(_surface_pot_expr[sp]);
1824 : }
1825 :
1826 : Real
1827 470 : GeochemicalSystem::getSurfaceCharge(unsigned sp) const
1828 : {
1829 470 : if (sp >= _num_surface_pot)
1830 2 : mooseError("Cannot retrieve the surface charge for surface ",
1831 : sp,
1832 : " since there are only ",
1833 2 : _num_surface_pot,
1834 : " surfaces involved in surface complexation");
1835 : // pre = mol of charge per square metre. To convert to Coulombs/m^2:
1836 468 : const Real pre = surfacePotPrefactor(sp) / _sorbing_surface_area[sp] *
1837 468 : (-_surface_pot_expr[sp] + 1.0 / _surface_pot_expr[sp]);
1838 468 : return pre * GeochemistryConstants::FARADAY;
1839 : }
1840 :
1841 : void
1842 47946 : GeochemicalSystem::computeSorbingSurfaceArea(std::vector<Real> & sorbing_surface_area) const
1843 : {
1844 51472 : for (unsigned sp = 0; sp < _num_surface_pot; ++sp)
1845 : {
1846 3526 : sorbing_surface_area[sp] = _mgd.surface_sorption_area[sp];
1847 : if (_mgd.basis_species_index.count(_mgd.surface_sorption_name[sp]) == 1)
1848 : {
1849 2971 : const unsigned basis_ind = _mgd.basis_species_index.at(_mgd.surface_sorption_name[sp]);
1850 : const Real grams =
1851 2971 : _mgd.basis_species_molecular_weight[basis_ind] * _basis_molality[basis_ind];
1852 2971 : sorbing_surface_area[sp] *= grams;
1853 : }
1854 : }
1855 47946 : }
1856 :
1857 : const std::vector<Real> &
1858 60 : GeochemicalSystem::getSorbingSurfaceArea() const
1859 : {
1860 60 : return _sorbing_surface_area;
1861 : }
1862 :
1863 : Real
1864 8271 : GeochemicalSystem::getTemperature() const
1865 : {
1866 8271 : return _temperature;
1867 : }
1868 :
1869 : void
1870 470 : GeochemicalSystem::setTemperature(Real temperature)
1871 : {
1872 470 : _temperature = temperature;
1873 470 : buildTemperatureDependentQuantities(_temperature);
1874 470 : _gac.setInternalParameters(_temperature, _mgd, _basis_molality, _eqm_molality, _kin_moles);
1875 470 : }
1876 :
1877 : void
1878 28 : GeochemicalSystem::setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(
1879 : const std::vector<std::string> & names,
1880 : const std::vector<Real> & values,
1881 : const std::vector<bool> & constraints_from_molalities)
1882 : {
1883 : // assume temperature-dependent quantities have been built during instantiation
1884 : // assume algebraic info has been built during instantiation
1885 :
1886 : /*
1887 : * STEP 0: Check sizes are correct
1888 : */
1889 28 : const unsigned num_names = names.size();
1890 28 : if (num_names != values.size())
1891 1 : mooseError("When setting all molalities, names and values must have same size");
1892 27 : if (num_names != _num_basis + _num_eqm + _num_surface_pot + _num_kin)
1893 2 : mooseError("When setting all molalities, values must be provided for every species and surface "
1894 : "potentials");
1895 25 : if (constraints_from_molalities.size() != _num_basis)
1896 1 : mooseError("constraints_from_molalities has size ",
1897 1 : constraints_from_molalities.size(),
1898 : " while number of basis species is ",
1899 1 : _num_basis);
1900 :
1901 : /*
1902 : * STEP 1: Read values into _basis_molality, _eqm_molality and _surface_pot_expr and
1903 : * _kin_moles
1904 : */
1905 : // The is no guarantee is made that the user-supplied molalities are "good", except they must
1906 : // be non-negative. There are lots of "finds" in the loops below, which is slow, but it
1907 : // ensures all species are given molalities. This method is designed to be called only every
1908 : // once-in-a-while (eg, at the start of a simulation)
1909 176 : for (const auto & name_index : _mgd.basis_species_index)
1910 : {
1911 : const unsigned ind =
1912 156 : std::distance(names.begin(), std::find(names.begin(), names.end(), name_index.first));
1913 156 : if (ind < num_names)
1914 : {
1915 155 : if (_mgd.basis_species_gas[name_index.second])
1916 : {
1917 8 : if (values[ind] != 0.0)
1918 1 : mooseError("Molality for gas ",
1919 : name_index.first,
1920 : " cannot be ",
1921 : values[ind],
1922 : ": it must be zero");
1923 : else
1924 7 : _basis_molality[name_index.second] = values[ind];
1925 : }
1926 147 : else if (_mgd.basis_species_mineral[name_index.second])
1927 : {
1928 15 : if (values[ind] < 0.0)
1929 1 : mooseError("Molality for mineral ",
1930 : name_index.first,
1931 : " cannot be ",
1932 : values[ind],
1933 : ": it must be non-negative");
1934 : else
1935 14 : _basis_molality[name_index.second] = values[ind];
1936 : }
1937 132 : else if (values[ind] <= 0.0)
1938 1 : mooseError("Molality for species ",
1939 : name_index.first,
1940 : " cannot be ",
1941 : values[ind],
1942 : ": it must be positive");
1943 : else
1944 131 : _basis_molality[name_index.second] = values[ind];
1945 : }
1946 : else
1947 1 : mooseError("Molality (or free mineral moles, etc - whatever is appropriate) for species ",
1948 : name_index.first,
1949 : " needs to be provided when setting all molalities");
1950 : }
1951 486 : for (const auto & name_index : _mgd.eqm_species_index)
1952 : {
1953 : const unsigned ind =
1954 470 : std::distance(names.begin(), std::find(names.begin(), names.end(), name_index.first));
1955 470 : if (ind < num_names)
1956 : {
1957 469 : if (_mgd.eqm_species_gas[name_index.second] || _mgd.eqm_species_mineral[name_index.second])
1958 30 : _eqm_molality[name_index.second] =
1959 : 0.0; // Note: explicitly setting to zero, irrespective of user input. The reason
1960 : // for doing this is that during a restore, a basis species with positive
1961 : // molality could become a secondary species, which should have zero molality
1962 439 : else if (values[ind] < 0.0)
1963 3 : mooseError("Molality for species ",
1964 : name_index.first,
1965 : " cannot be ",
1966 : values[ind],
1967 : ": it must be non-negative");
1968 : else
1969 436 : _eqm_molality[name_index.second] = values[ind];
1970 : }
1971 : else
1972 1 : mooseError("Molality for species ",
1973 : name_index.first,
1974 : " needs to be provided when setting all molalities");
1975 : }
1976 22 : for (unsigned sp = 0; sp < _num_surface_pot; ++sp)
1977 : {
1978 : const unsigned ind =
1979 8 : std::distance(names.begin(),
1980 : std::find(names.begin(),
1981 : names.end(),
1982 8 : _mgd.surface_sorption_name[sp] + "_surface_potential_expr"));
1983 8 : if (ind < num_names)
1984 : {
1985 7 : if (values[ind] <= 0.0)
1986 1 : mooseError("Surface-potential expression for mineral ",
1987 1 : _mgd.surface_sorption_name[sp],
1988 : " cannot be ",
1989 : values[ind],
1990 : ": it must be positive");
1991 6 : _surface_pot_expr[sp] = values[ind];
1992 : }
1993 : else
1994 1 : mooseError("Surface potential for mineral ",
1995 1 : _mgd.surface_sorption_name[sp],
1996 : " needs to be provided when setting all molalities");
1997 : }
1998 16 : for (const auto & name_index : _mgd.kin_species_index)
1999 : {
2000 : const unsigned ind =
2001 4 : std::distance(names.begin(), std::find(names.begin(), names.end(), name_index.first));
2002 4 : if (ind < num_names)
2003 3 : setKineticMoles(name_index.second, values[ind]);
2004 : else
2005 1 : mooseError("Moles for species ",
2006 : name_index.first,
2007 : " needs to be provided when setting all molalities");
2008 : }
2009 :
2010 : /*
2011 : * STEP 2: Alter _constraint_values if necessary
2012 : */
2013 : // If some of the constraints_from_molalities are false, then the molalities provided to this
2014 : // method may have to be modified to satisfy the contraints: this alters _basis_molality so
2015 : // must occur first
2016 119 : for (unsigned i = 0; i < _num_basis; ++i)
2017 : {
2018 107 : const ConstraintMeaningEnum meaning = _constraint_meaning[i];
2019 107 : if (meaning == ConstraintMeaningEnum::KG_SOLVENT_WATER ||
2020 107 : meaning == ConstraintMeaningEnum::FREE_MOLALITY ||
2021 : meaning == ConstraintMeaningEnum::FREE_MOLES_MINERAL_SPECIES)
2022 : {
2023 25 : if (constraints_from_molalities[i])
2024 : {
2025 : // molalities provided to this method dictate the contraints:
2026 16 : _constraint_value[i] = _basis_molality[i];
2027 16 : _original_constraint_value[i] = _constraint_value[i];
2028 : }
2029 : else
2030 : // contraints take precedence over the molalities provided to this method:
2031 9 : _basis_molality[i] = _constraint_value[i];
2032 : }
2033 : }
2034 :
2035 : // Potentially alter _constraint_value for the BULK contraints:
2036 119 : for (unsigned i = 0; i < _num_basis; ++i)
2037 : {
2038 107 : const ConstraintMeaningEnum meaning = _constraint_meaning[i];
2039 107 : if (meaning == ConstraintMeaningEnum::MOLES_BULK_WATER ||
2040 107 : meaning == ConstraintMeaningEnum::MOLES_BULK_SPECIES)
2041 70 : if (constraints_from_molalities[i])
2042 : {
2043 : // the constraint value should be overridden by the molality-computed bulk mole number
2044 40 : _constraint_value[i] = computeBulkFromMolalities(i);
2045 40 : _original_constraint_value[i] = _constraint_value[i];
2046 : }
2047 : }
2048 :
2049 : // Potentially alter _constraint_value for ACTIVITY contraints:
2050 111 : for (unsigned i = 0; i < _num_basis; ++i)
2051 101 : if (constraints_from_molalities[i])
2052 : {
2053 56 : const ConstraintMeaningEnum meaning = _constraint_meaning[i];
2054 56 : if (meaning == ConstraintMeaningEnum::FUGACITY)
2055 1 : mooseError("Gas fugacity cannot be determined from molality: constraints_from_molalities "
2056 : "must be set false for all gases");
2057 55 : else if (meaning == ConstraintMeaningEnum::ACTIVITY)
2058 : {
2059 5 : if (i == 0)
2060 1 : mooseError("Water activity cannot be determined from molalities: "
2061 : "constraints_from_molalities "
2062 : "must be set to false for water if activity of water is fixed");
2063 : // the constraint value should be overidden by the molality provided to this method
2064 4 : _constraint_value[i] = _basis_activity_coef[i] * _basis_molality[i];
2065 4 : _original_constraint_value[i] = _constraint_value[i];
2066 : }
2067 : }
2068 :
2069 : /*
2070 : * STEP 3: Follow the initialize() and computeConsistentConfiguration() methods
2071 : */
2072 : // assume done already: buildTemperatureDependentQuantities
2073 10 : enforceChargeBalanceIfSimple(_constraint_value, _bulk_moles_old);
2074 : // assume done already: buildAlgebraicInfo
2075 : // should not be done, as basis_molality is set by this method instead: initBulkAndFree
2076 10 : buildKnownBasisActivities(_basis_activity_known, _basis_activity);
2077 : // should not be done, as these are set by this method: _eqm_molality.assign
2078 : // should not be done, as these are set by this method: _surface_pot_expr.assign
2079 10 : _gac.setInternalParameters(_temperature, _mgd, _basis_molality, _eqm_molality, _kin_moles);
2080 10 : _gac.buildActivityCoefficients(_mgd, _basis_activity_coef, _eqm_activity_coef);
2081 : // for constraints_from_molalities = false then the following: (1) overrides the basis
2082 : // molality provided to this method; (2) produces a slightly inconsistent equilibrium
2083 : // geochemical system because basis_activity_coef was computed on the basis of the basis
2084 : // molalities provided to this method
2085 10 : updateBasisMolalityForKnownActivity(_basis_molality);
2086 10 : computeRemainingBasisActivities(_basis_activity);
2087 : // should not be done, as these are set by this method: computeEqmMolalities
2088 10 : computeBulk(_bulk_moles_old);
2089 : // should not be done, as these are set by this method: computeFreeMineralMoles
2090 10 : computeSorbingSurfaceArea(_sorbing_surface_area);
2091 10 : }
2092 :
2093 : const std::vector<GeochemicalSystem::ConstraintMeaningEnum> &
2094 745 : GeochemicalSystem::getConstraintMeaning() const
2095 : {
2096 745 : return _constraint_meaning;
2097 : }
2098 :
2099 : void
2100 711 : GeochemicalSystem::closeSystem()
2101 : {
2102 3360 : for (unsigned basis_ind = 0; basis_ind < _num_basis; ++basis_ind)
2103 2649 : if (_constraint_meaning[basis_ind] == ConstraintMeaningEnum::KG_SOLVENT_WATER ||
2104 5068 : _constraint_meaning[basis_ind] == ConstraintMeaningEnum::FREE_MOLALITY ||
2105 : _constraint_meaning[basis_ind] == ConstraintMeaningEnum::FREE_MOLES_MINERAL_SPECIES)
2106 306 : changeConstraintToBulk(basis_ind);
2107 711 : }
2108 :
2109 : void
2110 412 : GeochemicalSystem::changeConstraintToBulk(unsigned basis_ind)
2111 : {
2112 412 : if (basis_ind >= _num_basis)
2113 1 : mooseError("Cannot changeConstraintToBulk for species ",
2114 : basis_ind,
2115 : " because there are only ",
2116 1 : _num_basis,
2117 : " basis species");
2118 411 : if (_mgd.basis_species_gas[basis_ind])
2119 : {
2120 : // this is a special case where we have to swap out the gas in favour of an equilibrium
2121 : // species
2122 : unsigned swap_into_basis = 0;
2123 : bool legitimate_swap_found = false;
2124 : Real best_stoi = 0.0;
2125 392 : for (unsigned j = 0; j < _num_eqm; ++j)
2126 : {
2127 386 : if (_mgd.eqm_species_gas[j] || _mgd.eqm_stoichiometry(j, basis_ind) == 0.0 ||
2128 : _mgd.surface_sorption_related[j])
2129 198 : continue;
2130 : const Real stoi = std::abs(_mgd.eqm_stoichiometry(j, basis_ind));
2131 188 : if (stoi > best_stoi)
2132 : {
2133 : best_stoi = stoi;
2134 : swap_into_basis = j;
2135 : legitimate_swap_found = true;
2136 : }
2137 : }
2138 6 : if (legitimate_swap_found)
2139 6 : performSwapNoCheck(basis_ind, swap_into_basis);
2140 : else
2141 0 : mooseError("Attempting to change constraint of gas ",
2142 : _mgd.basis_species_name[basis_ind],
2143 : " to MOLES_BULK_SPECIES, which involves a search for a suitable non-gas species "
2144 : "to swap with. No such species was found");
2145 : }
2146 : else
2147 405 : changeConstraintToBulk(basis_ind, computeBulkFromMolalities(basis_ind));
2148 411 : }
2149 :
2150 : void
2151 408 : GeochemicalSystem::changeConstraintToBulk(unsigned basis_ind, Real value)
2152 : {
2153 408 : if (basis_ind >= _num_basis)
2154 1 : mooseError("Cannot changeConstraintToBulk for species ",
2155 : basis_ind,
2156 : " because there are only ",
2157 1 : _num_basis,
2158 : " basis species");
2159 407 : if (_mgd.basis_species_gas[basis_ind])
2160 1 : mooseError("Attempting to changeConstraintToBulk for a gas species. Since a swap is involved, "
2161 : "you cannot specify a value for the bulk number of moles. You must use "
2162 : "changeConstraintToBulk(basis_ind) method instead of "
2163 : "changeConstraintToBulk(basis_ind, value)");
2164 406 : if (basis_ind == 0)
2165 135 : _constraint_meaning[basis_ind] = ConstraintMeaningEnum::MOLES_BULK_WATER;
2166 : else
2167 271 : _constraint_meaning[basis_ind] = ConstraintMeaningEnum::MOLES_BULK_SPECIES;
2168 406 : setConstraintValue(basis_ind, value);
2169 :
2170 : // it is possible that FIXED_ACTIVITY just became MOLES_BULK_SPECIES
2171 406 : buildKnownBasisActivities(_basis_activity_known, _basis_activity);
2172 :
2173 : // it is likely the algebraic system has changed
2174 406 : buildAlgebraicInfo(_in_algebraic_system,
2175 406 : _num_basis_in_algebraic_system,
2176 406 : _num_in_algebraic_system,
2177 406 : _algebraic_index,
2178 406 : _basis_index);
2179 406 : }
2180 :
2181 : void
2182 57419 : GeochemicalSystem::addToBulkMoles(unsigned basis_ind, Real value)
2183 : {
2184 57419 : if (basis_ind >= _num_basis)
2185 1 : mooseError("Cannot addToBulkMoles for species ",
2186 : basis_ind,
2187 : " because there are only ",
2188 1 : _num_basis,
2189 : " basis species");
2190 57418 : if (!(_constraint_meaning[basis_ind] ==
2191 : GeochemicalSystem::ConstraintMeaningEnum::MOLES_BULK_WATER ||
2192 : _constraint_meaning[basis_ind] ==
2193 : GeochemicalSystem::ConstraintMeaningEnum::MOLES_BULK_SPECIES))
2194 : return;
2195 51976 : setConstraintValue(basis_ind, _constraint_value[basis_ind] + value);
2196 : }
2197 :
2198 : void
2199 52924 : GeochemicalSystem::setConstraintValue(unsigned basis_ind, Real value)
2200 : {
2201 52924 : if (basis_ind >= _num_basis)
2202 1 : mooseError("Cannot setConstraintValue for species ",
2203 : basis_ind,
2204 : " because there are only ",
2205 1 : _num_basis,
2206 : " basis species");
2207 52923 : _constraint_value[basis_ind] = value;
2208 52923 : _original_constraint_value[basis_ind] = value;
2209 52923 : switch (_constraint_meaning[basis_ind])
2210 : {
2211 52386 : case ConstraintMeaningEnum::MOLES_BULK_SPECIES:
2212 : case ConstraintMeaningEnum::MOLES_BULK_WATER:
2213 : {
2214 52386 : alterSystemBecauseBulkChanged();
2215 52386 : break;
2216 : }
2217 1 : case ConstraintMeaningEnum::KG_SOLVENT_WATER:
2218 : case ConstraintMeaningEnum::FREE_MOLALITY:
2219 : case ConstraintMeaningEnum::FREE_MOLES_MINERAL_SPECIES:
2220 : {
2221 : // the changes resulting from this change are very similar to setAlgebraicVariables
2222 1 : _basis_molality[basis_ind] = value;
2223 1 : computeConsistentConfiguration();
2224 1 : break;
2225 : }
2226 50 : case ConstraintMeaningEnum::FUGACITY:
2227 : {
2228 : // the changes resulting from this change are very similar to setAlgebraicVariables
2229 50 : _basis_activity[basis_ind] = value;
2230 50 : _basis_molality[basis_ind] = 0.0;
2231 50 : computeConsistentConfiguration();
2232 50 : break;
2233 : }
2234 486 : case ConstraintMeaningEnum::ACTIVITY:
2235 : {
2236 : // the changes resulting from this change are very similar to setAlgebraicVariables
2237 486 : _basis_activity[basis_ind] = value;
2238 486 : _basis_molality[basis_ind] = _basis_activity[basis_ind] / _basis_activity_coef[basis_ind];
2239 486 : computeConsistentConfiguration();
2240 486 : break;
2241 : }
2242 : }
2243 52923 : }
2244 :
2245 : void
2246 52386 : GeochemicalSystem::alterSystemBecauseBulkChanged()
2247 : {
2248 : // Altering the constraints on bulk number of moles impacts various other things.
2249 : // Here, we follow the initialize() and computeConsistentConfiguration() methods, picking out
2250 : // what might have changed.
2251 : // Because the constraint meanings have changed, enforcing charge-neutrality might be easy
2252 52386 : enforceChargeBalanceIfSimple(_constraint_value, _bulk_moles_old);
2253 : // Because the constraint values have changed, either through charge-neutrality or directly
2254 : // through changing the constraint values, the bulk moles must be updated:
2255 475229 : for (unsigned i = 0; i < _num_basis; ++i)
2256 422843 : if (_constraint_meaning[i] == ConstraintMeaningEnum::MOLES_BULK_SPECIES ||
2257 : _constraint_meaning[i] == ConstraintMeaningEnum::MOLES_BULK_WATER)
2258 394573 : _bulk_moles_old[i] = _constraint_value[i];
2259 : // Because the bulk mineral moles might have changed, the free mineral moles might have
2260 : // changed
2261 52386 : computeFreeMineralMoles(_basis_molality); // free mineral moles might be changed
2262 52386 : }
2263 :
2264 : const std::string &
2265 51 : GeochemicalSystem::getOriginalRedoxLHS() const
2266 : {
2267 51 : return _original_redox_lhs;
2268 : }
2269 :
2270 : void
2271 62 : GeochemicalSystem::setModelGeochemicalDatabase(const ModelGeochemicalDatabase & mgd)
2272 : {
2273 62 : _mgd = mgd;
2274 62 : }
2275 :
2276 : const GeochemistrySpeciesSwapper &
2277 71 : GeochemicalSystem::getSwapper() const
2278 : {
2279 71 : return _swapper;
2280 : }
2281 :
2282 : void
2283 107 : GeochemicalSystem::setMineralRelatedFreeMoles(Real value)
2284 : {
2285 : // mole numbers of basis minerals
2286 916 : for (unsigned i = 0; i < _num_basis; ++i)
2287 809 : if (_mgd.basis_species_mineral[i])
2288 77 : _basis_molality[i] = value;
2289 :
2290 : // mole numbers of sorption sites
2291 107 : for (const auto & name_info :
2292 109 : _mgd.surface_complexation_info) // all minerals involved in surface complexation
2293 2 : for (const auto & name_frac :
2294 6 : name_info.second.sorption_sites) // all sorption sites on the given mineral
2295 : {
2296 4 : if (_mgd.basis_species_index.count(name_frac.first))
2297 4 : _basis_molality[_mgd.basis_species_index.at(name_frac.first)] = value;
2298 : else
2299 0 : _eqm_molality[_mgd.eqm_species_index.at(name_frac.first)] = value;
2300 : }
2301 :
2302 : // mole numbers of sorbed equilibrium species
2303 2749 : for (unsigned j = 0; j < _num_eqm; ++j)
2304 : {
2305 2642 : if (_mgd.eqm_species_mineral[j])
2306 361 : _eqm_molality[j] = 0.0; // Note: explicitly setting to zero, irrespective of user input,
2307 : // to ensure consistency with the rest of the code
2308 2642 : if (_mgd.surface_sorption_related[j])
2309 3 : _eqm_molality[j] = value;
2310 : }
2311 :
2312 : // mole numbers of kinetic species
2313 110 : for (unsigned k = 0; k < _num_kin; ++k)
2314 3 : if (_mgd.kin_species_mineral[k])
2315 3 : _kin_moles[k] = value;
2316 107 : }
2317 :
2318 : Real
2319 132029 : GeochemicalSystem::getEquilibriumActivity(unsigned eqm_ind) const
2320 : {
2321 132029 : if (eqm_ind >= _num_eqm)
2322 1 : mooseError("Cannot retrieve activity for equilibrium species ",
2323 : eqm_ind,
2324 : " since there are only ",
2325 1 : _num_eqm,
2326 : " equilibrium species");
2327 132028 : if (_mgd.eqm_species_mineral[eqm_ind])
2328 : return 1.0;
2329 126776 : else if (_mgd.eqm_species_gas[eqm_ind])
2330 : {
2331 : Real log10f = 0.0;
2332 2536 : for (unsigned basis_i = 0; basis_i < _num_basis; ++basis_i)
2333 2232 : log10f += _mgd.eqm_stoichiometry(eqm_ind, basis_i) * std::log10(_basis_activity[basis_i]);
2334 304 : log10f -= getLog10K(eqm_ind);
2335 304 : return std::pow(10.0, log10f);
2336 : }
2337 : else
2338 126472 : return _eqm_activity_coef[eqm_ind] * _eqm_molality[eqm_ind];
2339 : }
2340 :
2341 : const std::vector<Real> &
2342 1 : GeochemicalSystem::computeAndGetEquilibriumActivity()
2343 : {
2344 7 : for (unsigned j = 0; j < _num_eqm; ++j)
2345 6 : _eqm_activity[j] = getEquilibriumActivity(j);
2346 1 : return _eqm_activity;
2347 : }
2348 :
2349 : void
2350 4570 : GeochemicalSystem::updateOldWithCurrent(const DenseVector<Real> & mole_additions)
2351 : {
2352 4570 : if (mole_additions.size() != _num_basis + _num_kin)
2353 1 : mooseError("The increment in mole numbers (mole_additions) needs to be of size ",
2354 1 : _num_basis,
2355 : " + ",
2356 1 : _num_kin,
2357 : " but it is of size ",
2358 1 : mole_additions.size());
2359 :
2360 : // Copy the current kin_moles to the old values
2361 5540 : for (unsigned k = 0; k < _num_kin; ++k)
2362 971 : _kin_moles_old[k] = _kin_moles[k];
2363 :
2364 : // The following:
2365 : // - just returns if constraint is not BULK type.
2366 : // - ensures charge balance if simple
2367 : // - sets the constraint value and bulk_moles_old (for BULK-type species)
2368 : // - computes basis_molality for the BULK-type mineral species
2369 29504 : for (unsigned i = 0; i < _num_basis; ++i)
2370 24935 : addToBulkMoles(i, mole_additions(i));
2371 :
2372 : // The following:
2373 : // - sets bulk_moles_old to the Constraint value for BULK-type species (duplicating the above
2374 : // function)
2375 : // - for the other species, computes bulk_moles_old from the molalities
2376 4569 : computeBulk(_bulk_moles_old);
2377 :
2378 : // It is possible that the user would also like to enforceChargeBalance() now, and that is fine
2379 : // - there will be no negative consequences
2380 4569 : }
2381 :
2382 : void
2383 40045 : GeochemicalSystem::addKineticRates(Real dt,
2384 : DenseVector<Real> & mole_additions,
2385 : DenseMatrix<Real> & dmole_additions)
2386 : {
2387 40045 : if (_num_kin == 0)
2388 33750 : return;
2389 :
2390 : // check sizes
2391 6295 : const unsigned tot = mole_additions.size();
2392 6295 : if (!(tot == _num_kin + _num_basis && dmole_additions.m() == tot && dmole_additions.n() == tot))
2393 3 : mooseError("addKineticRates: incorrectly sized additions: ",
2394 : tot,
2395 : " ",
2396 3 : dmole_additions.m(),
2397 : " ",
2398 3 : dmole_additions.n());
2399 :
2400 : // construct eqm_activity for species that we need
2401 187143 : for (unsigned j = 0; j < _num_eqm; ++j)
2402 539843 : if (_mgd.eqm_species_gas[j] || _mgd.eqm_species_name[j] == "H+" ||
2403 178141 : _mgd.eqm_species_name[j] == "OH-")
2404 8460 : _eqm_activity[j] = getEquilibriumActivity(j);
2405 :
2406 : // calculate the rates and put into appropriate slots
2407 : Real rate;
2408 : Real drate_dkin;
2409 6292 : std::vector<Real> drate_dmol(_num_basis);
2410 12200 : for (const auto & krd : _mgd.kin_rate)
2411 : {
2412 5908 : const unsigned kin = krd.kinetic_species_index;
2413 11816 : GeochemistryKineticRateCalculator::calculateRate(krd.promoting_indices,
2414 5908 : krd.promoting_monod_indices,
2415 5908 : krd.promoting_half_saturation,
2416 5908 : krd.description,
2417 5908 : _mgd.basis_species_name,
2418 5908 : _mgd.basis_species_gas,
2419 5908 : _basis_molality,
2420 5908 : _basis_activity,
2421 5908 : _basis_activity_known,
2422 5908 : _mgd.eqm_species_name,
2423 5908 : _mgd.eqm_species_gas,
2424 5908 : _eqm_molality,
2425 5908 : _eqm_activity,
2426 5908 : _mgd.eqm_stoichiometry,
2427 : _kin_moles[kin],
2428 5908 : _mgd.kin_species_molecular_weight[kin],
2429 5908 : _kin_log10K[kin],
2430 : log10KineticActivityProduct(kin),
2431 5908 : _mgd.kin_stoichiometry,
2432 : kin,
2433 : _temperature,
2434 : rate,
2435 : drate_dkin,
2436 : drate_dmol);
2437 :
2438 : /* the following block of code may be confusing at first sight.
2439 : * (1) The usual case is that the kinetic reaction is written
2440 : * kinetic -> basis_species, and direction = BOTH, and kinetic_bio_efficiency = -1
2441 : * In this case, the following does mole_additions(kinetic_species) += -1 * rate * dt
2442 : * If rate > 0 then the kinetic species decreases in mass.
2443 : * The solver will then detect that the kinetic_species has changed, and adjust basis_species
2444 : * accordingly, viz it will add rate * dt of basis species to the aqueous solution.
2445 : * (2) The common biologically-catalysed case is the reaction is written
2446 : * reactants -> products, catalysed by microbe. In the database file this is written as
2447 : * microbe -> -reactants + products
2448 : * The input file will set direction = BOTH, and kinetic_bio_efficiency != -1
2449 : * Suppose that rate > 0, corresponding to reactants being converted to products.
2450 : * With kinetic_bio_efficiency = -1) this would correspond to microbe decreasing in mass,
2451 : * however with kinetic_bio_efficiency > 0, this is not true.
2452 : * The following block of code sets
2453 : * mole_additions(microbe) = kinetic_bio_efficiency * rate * dt > 0
2454 : * that is, the microbe mass is increasing.
2455 : * However, this means the solver will decrease products and increase reactants,
2456 : * by kinetic_bio_efficiency * rate * dt because it sees the reaction microbe -> -reactants +
2457 : * products in the database file.
2458 : * This is clearly wrong, because we actually want the end result to be that the products have
2459 : * increased by rate * dt, and the reactants to have decreased by rate * dt.
2460 : * So, to counter the solver, we add (-1 - kinetic_bio_efficiency) * rate * dt of reactants and
2461 : * add (1 + kinetic_bio_efficienty) * rate * dt of products.
2462 : * (3) The other situation is the biological death, where direction = DEATH.
2463 : * In this case the database file contains microbe -> -reactants + products
2464 : * However, independent of rate, we don't want to change the molality of reactants or products.
2465 : * Following the logic of (2), to counter the solver, we add -kinetic_bio_efficiency * rate * dt
2466 : * of reactants and kinetic_bio_efficiency * rate * dt of products
2467 : */
2468 5908 : const unsigned ind = kin + _num_basis;
2469 5908 : mole_additions(ind) += krd.description.kinetic_bio_efficiency * rate * dt;
2470 5908 : dmole_additions(ind, ind) += krd.description.kinetic_bio_efficiency * drate_dkin * dt;
2471 5908 : const Real extra_additions = (krd.description.direction == DirectionChoiceEnum::DEATH)
2472 5908 : ? krd.description.kinetic_bio_efficiency
2473 5675 : : krd.description.kinetic_bio_efficiency + 1.0;
2474 52151 : for (unsigned i = 0; i < _num_basis; ++i)
2475 : {
2476 46243 : dmole_additions(ind, i) += krd.description.kinetic_bio_efficiency * drate_dmol[i] * dt;
2477 46243 : const Real stoi_fac = _mgd.kin_stoichiometry(kin, i) * extra_additions * dt;
2478 46243 : mole_additions(i) += stoi_fac * rate;
2479 46243 : dmole_additions(i, ind) += stoi_fac * drate_dkin;
2480 451184 : for (unsigned j = 0; j < _num_basis; ++j)
2481 404941 : dmole_additions(i, j) += stoi_fac * drate_dmol[j];
2482 : }
2483 :
2484 5908 : const Real eff = krd.description.progeny_efficiency;
2485 5908 : if (eff != 0.0)
2486 : {
2487 166 : const unsigned bio_i = krd.progeny_index;
2488 166 : if (bio_i < _num_basis)
2489 : {
2490 165 : mole_additions(bio_i) += eff * rate * dt;
2491 165 : dmole_additions(bio_i, ind) += eff * drate_dkin * dt;
2492 1805 : for (unsigned i = 0; i < _num_basis; ++i)
2493 1640 : dmole_additions(bio_i, i) += eff * drate_dmol[i] * dt;
2494 : }
2495 : else
2496 : {
2497 6 : for (unsigned i = 0; i < _num_basis; ++i)
2498 : {
2499 5 : mole_additions(i) += _mgd.eqm_stoichiometry(bio_i - _num_basis, i) * eff * rate * dt;
2500 5 : dmole_additions(i, ind) +=
2501 5 : _mgd.eqm_stoichiometry(bio_i - _num_basis, i) * eff * drate_dkin * dt;
2502 30 : for (unsigned j = 0; j < _num_basis; ++j)
2503 25 : dmole_additions(i, j) +=
2504 25 : _mgd.eqm_stoichiometry(bio_i - _num_basis, i) * eff * drate_dmol[j] * dt;
2505 : }
2506 : }
2507 : }
2508 : }
2509 : }
|
