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 : #pragma once
11 :
12 : #include "MultiMooseEnum.h"
13 : #include "PertinentGeochemicalSystem.h"
14 : #include "GeochemistrySpeciesSwapper.h"
15 : #include "GeochemistryActivityCoefficientsDebyeHuckel.h"
16 : #include "GeochemistryConstants.h"
17 : #include "GeochemistryUnitConverter.h"
18 :
19 : /**
20 : * This class holds information about bulk composition, molalities, activities, activity
21 : * coefficients, etc of the user-defined geochemical system in PertinentGeochemicalSystem. Hence,
22 : * it extends the generic information in PertinentGeochemicalSystem to a system that is specific to
23 : * one spatio-temporal location. It offers methods to manipulate these molalities, bulk
24 : * compositions, etc, as well as performing basis swaps.
25 : *
26 : * Related to this, this class also builds the so-called "algebraic system" which is the nonlinear
27 : * system of algebraic equations for the unknown basis molalities, surface potentials and mole
28 : * number of kinetic species, and computes the residual vector and jacobian of this algebraic
29 : * system. This class initialises the algebraic variables (the unknown molalities, surface
30 : * potentials and kinetic moles) to reaponable values, but it does not solve the algebraic system:
31 : * instead, it has a "setter" method, setAlgebraicValues that can be used to set the unknowns, hence
32 : * another class can use whatever method it desires to solve the system.
33 : *
34 : * WARNING: because GeochemicalSystem performs swaps, it will change the
35 : * ModelGeochemicalDatabase object.
36 : * WARNING: because GeochemicalSystem sets internal parameters in the
37 : * activity-coefficient object, it will change the GeochemistryActivityCoefficient object.
38 : */
39 : class GeochemicalSystem
40 : {
41 : public:
42 : /**
43 : * Each basis species has one of the following constraints. During the process of units
44 : * conversion (from user-prescribed units to mole-based units) each ConstraintMeaningUserEnum
45 : * supplied by the user is translated to the appropriate ConstraintMeaningEnum
46 : */
47 : enum class ConstraintMeaningEnum
48 : {
49 : MOLES_BULK_WATER,
50 : KG_SOLVENT_WATER,
51 : MOLES_BULK_SPECIES,
52 : FREE_MOLALITY,
53 : FREE_MOLES_MINERAL_SPECIES,
54 : FUGACITY,
55 : ACTIVITY
56 : };
57 :
58 : /// Each basis species must be provided with a constraint, that is chosen by the user from the following enum
59 : enum class ConstraintUserMeaningEnum
60 : {
61 : KG_SOLVENT_WATER,
62 : BULK_COMPOSITION,
63 : BULK_COMPOSITION_WITH_KINETIC,
64 : FREE_CONCENTRATION,
65 : FREE_MINERAL,
66 : ACTIVITY,
67 : LOG10ACTIVITY,
68 : FUGACITY,
69 : LOG10FUGACITY
70 : };
71 :
72 : /**
73 : * Construct the geochemical system, check consistency of the constructor's arguments,
74 : and initialize all internal variables (molalities, bulk compositions, equilibrium constants,
75 : activities, activity coefficients, surface potentials and kinetic mole numbers) and set up the
76 : algebraic system
77 : * @param mgd the model's geochemical database, which is a model-specific version of the full
78 : geochemical database. WARNING: because GeochemicalSystem performs swaps, it will
79 : change mgd
80 : * @param gac the object that computes activity coefficients. WARNING: because
81 : GeochemicalSystem sets internal parameters in gac, it will change it.
82 : * @param is object to compute ionic strengths
83 : * @param swapper object to perform swaps on mgd
84 : * @param swap_out_of_basis A list of basis species that should be swapped out of the basis and
85 : replaced with swap_into_basis, prior to any other computations
86 : * @param swap_into_basis A list of equilibrium species that should be swapped into the basis in
87 : place of swap_out_of_basis
88 : * @param charge_balance_species The charge-balance species, which must be a basis species after
89 : the swaps are performed
90 : * @param constrained_species A list of the basis species after the swaps have been performed
91 : * @param constraint_values Numerical values of the constraints placed on the basis species. Each
92 : basis species must have exactly one constraint
93 : * @param constraint_unit Units of the constraint_values. Each constraint_value must have
94 : exactly one constraint_unit. This is only used during construction, to convert the
95 : constraint_values to mole units
96 : * @param constraint_user_meaning A list the provides physical meaning to the constraint_values
97 : * @param initial_temperature Initial temperature
98 : * @param iters_to_make_consistent The initial equilibrium molalities depend on the activity
99 : coefficients, which depend on the basis and equilibrium molalities. This circular dependence
100 : means it is usually impossible to define an exactly consistent initial configuration for
101 : molalities. An iterative process is used to approach the consistent initial configuration, using
102 : iters_to_make_consistent iterations. Usually iters_to_make_consistent=0 is reasonable, because
103 : there are so many approximations used in the solution process that a slightly inconsistent initial
104 : configuration is fine
105 : * @param min_initial_molality Minimum value of equilibrium molality used in the initial condition
106 : * @param kin_name Names of the kinetic species that are provided with initial conditions in
107 : kin_initial_moles. All kinetic species must be provided with an initial condition
108 : * @param kin_initial Values of the initial mole number or mass or volume (depending on
109 : kin_unit) for the kinetic species
110 : * @param kin_unit The units of the numbers given in kin_initial
111 : */
112 : GeochemicalSystem(ModelGeochemicalDatabase & mgd,
113 : GeochemistryActivityCoefficients & gac,
114 : GeochemistryIonicStrength & is,
115 : GeochemistrySpeciesSwapper & swapper,
116 : const std::vector<std::string> & swap_out_of_basis,
117 : const std::vector<std::string> & swap_into_basis,
118 : const std::string & charge_balance_species,
119 : const std::vector<std::string> & constrained_species,
120 : const std::vector<Real> & constraint_value,
121 : const MultiMooseEnum & constraint_unit,
122 : const MultiMooseEnum & constraint_user_meaning,
123 : Real initial_temperature,
124 : unsigned iters_to_make_consistent,
125 : Real min_initial_molality,
126 : const std::vector<std::string> & kin_name,
127 : const std::vector<Real> & kin_initial_moles,
128 : const MultiMooseEnum & kin_unit);
129 :
130 : GeochemicalSystem(
131 : ModelGeochemicalDatabase & mgd,
132 : GeochemistryActivityCoefficients & gac,
133 : GeochemistryIonicStrength & is,
134 : GeochemistrySpeciesSwapper & swapper,
135 : const std::vector<std::string> & swap_out_of_basis,
136 : const std::vector<std::string> & swap_into_basis,
137 : const std::string & charge_balance_species,
138 : const std::vector<std::string> & constrained_species,
139 : const std::vector<Real> & constraint_value,
140 : const std::vector<GeochemistryUnitConverter::GeochemistryUnit> & constraint_unit,
141 : const std::vector<ConstraintUserMeaningEnum> & constraint_user_meaning,
142 : Real initial_temperature,
143 : unsigned iters_to_make_consistent,
144 : Real min_initial_molality,
145 : const std::vector<std::string> & kin_name,
146 : const std::vector<Real> & kin_initial,
147 : const std::vector<GeochemistryUnitConverter::GeochemistryUnit> & kin_unit);
148 :
149 : /**
150 : * Copy assignment operator. Almost all of the following is trivial. The most important
151 : * non-trivial feature is copying src._mgd into our _mgd. Note this method gets called when
152 : * dest = src
153 : * and dest is already constructed. (Code such as GeochemicalSystem dest = src uses the copy
154 : * constructor which simply sets the _mgd reference in dest equal to the _mgd reference in src,
155 : * and does not make a copy of the data within src._mgd)
156 : */
157 863 : GeochemicalSystem & operator=(const GeochemicalSystem & src)
158 : {
159 863 : if (this == &src) // trivial a=a situation
160 : return *this;
161 : // check for bad assignment situations. Other "const" things that
162 : // needn't be the same (but probably actually are) include: _swap_out and _swap_in
163 862 : if (_num_basis != src._num_basis || _num_eqm != src._num_eqm || _num_redox != src._num_redox ||
164 1725 : _num_surface_pot != src._num_surface_pot || _num_kin != src._num_kin ||
165 861 : _original_charge_balance_species != src._original_charge_balance_species)
166 3 : mooseError("GeochemicalSystem: copy assigment operator called with inconsistent fundamental "
167 : "properties");
168 : // actually do the copying
169 860 : _mgd = src._mgd;
170 860 : _swapper = src._swapper;
171 : _gac = src._gac;
172 860 : _is = src._is;
173 860 : _charge_balance_species = src._charge_balance_species;
174 860 : _charge_balance_basis_index = src._charge_balance_basis_index;
175 860 : _constrained_species = src._constrained_species;
176 860 : _constraint_value = src._constraint_value;
177 860 : _original_constraint_value = src._original_constraint_value;
178 860 : _constraint_unit = src._constraint_unit;
179 860 : _constraint_user_meaning = src._constraint_user_meaning;
180 860 : _constraint_meaning = src._constraint_meaning;
181 860 : _eqm_log10K = src._eqm_log10K;
182 860 : _redox_log10K = src._redox_log10K;
183 860 : _kin_log10K = src._kin_log10K;
184 860 : _num_basis_in_algebraic_system = src._num_basis_in_algebraic_system;
185 860 : _num_in_algebraic_system = src._num_in_algebraic_system;
186 860 : _in_algebraic_system = src._in_algebraic_system;
187 860 : _algebraic_index = src._algebraic_index;
188 860 : _basis_index = src._basis_index;
189 860 : _bulk_moles_old = src._bulk_moles_old;
190 860 : _basis_molality = src._basis_molality;
191 860 : _basis_activity_known = src._basis_activity_known;
192 860 : _basis_activity = src._basis_activity;
193 860 : _eqm_molality = src._eqm_molality;
194 860 : _basis_activity_coef = src._basis_activity_coef;
195 860 : _eqm_activity_coef = src._eqm_activity_coef;
196 860 : _eqm_activity = src._eqm_activity;
197 860 : _surface_pot_expr = src._surface_pot_expr;
198 860 : _sorbing_surface_area = src._sorbing_surface_area;
199 860 : _kin_moles = src._kin_moles;
200 860 : _kin_moles_old = src._kin_moles_old;
201 860 : _iters_to_make_consistent = src._iters_to_make_consistent;
202 860 : _temperature = src._temperature;
203 860 : _min_initial_molality = src._min_initial_molality;
204 860 : _original_redox_lhs = src._original_redox_lhs;
205 860 : return *this;
206 : };
207 :
208 : /**
209 : * Copy constructor.
210 : */
211 2104 : GeochemicalSystem(const GeochemicalSystem & src) = default;
212 :
213 2 : bool operator==(const GeochemicalSystem & rhs) const
214 : {
215 2 : return (_mgd == rhs._mgd) && (_num_basis == rhs._num_basis) && (_num_eqm == rhs._num_eqm) &&
216 2 : (_num_redox == rhs._num_redox) && (_num_surface_pot == rhs._num_surface_pot) &&
217 2 : (_num_kin == rhs._num_kin) && (_swapper == rhs._swapper) &&
218 2 : (_swap_out == rhs._swap_out) && (_swap_in == rhs._swap_in) && (_gac == rhs._gac) &&
219 4 : (_is == rhs._is) && (_charge_balance_species == rhs._charge_balance_species) &&
220 2 : (_original_charge_balance_species == rhs._original_charge_balance_species) &&
221 2 : (_charge_balance_basis_index == rhs._charge_balance_basis_index) &&
222 2 : (_constrained_species == rhs._constrained_species) &&
223 2 : (_constraint_value == rhs._constraint_value) &&
224 2 : (_original_constraint_value == rhs._original_constraint_value) &&
225 2 : (_constraint_unit == rhs._constraint_unit) &&
226 2 : (_constraint_user_meaning == rhs._constraint_user_meaning) &&
227 2 : (_constraint_meaning == rhs._constraint_meaning) && (_eqm_log10K == rhs._eqm_log10K) &&
228 2 : (_redox_log10K == rhs._redox_log10K) && (_kin_log10K == rhs._kin_log10K) &&
229 2 : (_num_basis_in_algebraic_system == rhs._num_basis_in_algebraic_system) &&
230 2 : (_num_in_algebraic_system == rhs._num_in_algebraic_system) &&
231 2 : (_in_algebraic_system == rhs._in_algebraic_system) &&
232 2 : (_algebraic_index == rhs._algebraic_index) && (_basis_index == rhs._basis_index) &&
233 2 : (_bulk_moles_old == rhs._bulk_moles_old) && (_basis_molality == rhs._basis_molality) &&
234 2 : (_basis_activity_known == rhs._basis_activity_known) &&
235 2 : (_basis_activity == rhs._basis_activity) && (_eqm_molality == rhs._eqm_molality) &&
236 2 : (_basis_activity_coef == rhs._basis_activity_coef) &&
237 2 : (_eqm_activity_coef == rhs._eqm_activity_coef) && (_eqm_activity == rhs._eqm_activity) &&
238 2 : (_surface_pot_expr == rhs._surface_pot_expr) &&
239 2 : (_sorbing_surface_area == rhs._sorbing_surface_area) && (_kin_moles == rhs._kin_moles) &&
240 2 : (_kin_moles_old == rhs._kin_moles_old) &&
241 2 : (_iters_to_make_consistent == rhs._iters_to_make_consistent) &&
242 2 : (_temperature == rhs._temperature) &&
243 4 : (_min_initial_molality == rhs._min_initial_molality) &&
244 2 : (_original_redox_lhs == rhs._original_redox_lhs);
245 : }
246 :
247 : /// returns the number of species in the basis
248 : unsigned getNumInBasis() const;
249 :
250 : /// returns the number of species in equilibrium with the basis components
251 : unsigned getNumInEquilibrium() const;
252 :
253 : /// returns the number of kinetic species
254 : unsigned getNumKinetic() const;
255 :
256 : /// initialize all variables, ready for a Newton solve of the algebraic system
257 : void initialize();
258 :
259 : /// return the index of the charge-balance species in the basis list
260 : unsigned getChargeBalanceBasisIndex() const;
261 :
262 : /**
263 : * If the molality of the charge-balance basis species is less than threshold molality, then
264 : * attempt to change the charge-balance species, preferably to another component with opposite
265 : * charge and big molality.
266 : * @param threshold_molality if the charge-balance basis species has molality greater than this,
267 : * this routine does nothing
268 : * @return true if the charge-balance species was changed
269 : */
270 : bool alterChargeBalanceSpecies(Real threshold_molality);
271 :
272 : /**
273 : * @param j the index of the equilibrium species
274 : * @return the value of log10(equilibrium constant) for the j^th equilibrium species
275 : */
276 : Real getLog10K(unsigned j) const;
277 :
278 : /**
279 : * @return the number of redox species in disequibrium
280 : */
281 : unsigned getNumRedox() const;
282 :
283 : /**
284 : * @param red the index of the redox species in disequilibrium (red < _num_redox)
285 : * @return the value of log10(equilibrium constant) for the given disequilibrium-redox species
286 : */
287 : Real getRedoxLog10K(unsigned red) const;
288 :
289 : /**
290 : * @param red the index of the redox species in disequilibrium (red < _num_redox)
291 : * @return the value of log10(activity product) for the given disequilibrium-redox species
292 : */
293 : Real log10RedoxActivityProduct(unsigned red) const;
294 :
295 : /**
296 : * @param kin the index of the kinetic species (must be < _num_kin)
297 : * @return the value of log10(equilibrium constant) for the given kinetic species
298 : */
299 : Real getKineticLog10K(unsigned kin) const;
300 :
301 : /**
302 : * @return the value of log10(equilibrium constant) for the each kinetic species
303 : */
304 : const std::vector<Real> & getKineticLog10K() const;
305 :
306 : /**
307 : * @param kin the index of the kinetic species (must be < _num_kin)
308 : * @return the value of log10(activity product) for the given kinetic species
309 : */
310 : Real log10KineticActivityProduct(unsigned kin) const;
311 :
312 : /**
313 : * @param kin the index of the kinetic species (must be < _num_kin)
314 : * @return the mole number for the given kinetic species
315 : */
316 : Real getKineticMoles(unsigned kin) const;
317 :
318 : /**
319 : * Sets the current AND old mole number for a kinetic species. Note: this does not compute a
320 : * consistent configuration (viz, the bulk mole composition is not updated): use
321 : * computeConsistentConfiguration if you need to.
322 : * @param kin the index of the kinetic species (must be < _num_kin)
323 : * @param moles the number of moles
324 : */
325 : void setKineticMoles(unsigned kin, Real moles);
326 :
327 : /**
328 : * @return the mole number for all kinetic species
329 : */
330 : const std::vector<Real> & getKineticMoles() const;
331 :
332 : /**
333 : * Usually used at the end of a solve, to provide correct initial conditions to the next
334 : * time-step, this method:
335 : * - sets kin_moles_old = kin_moles
336 : * - updates the constraint_values and bulk_moles_old with mole_additions, for BULK-type species
337 : * - ensures charge balance holds if all species are BULK-type species
338 : * - computes basis_molality for the BULK-type mineral species
339 : * - computes bulk_moles_old from the molalities for non-BULK-type species
340 : * Note that it is possible that you would like to enforceChargeBalance() immediately after
341 : * calling this method: that is fine, there will be no negative consequences if the system has
342 : * been solved
343 : * @param mole_additions the increment of mole number of each basis species and kinetic species
344 : * since the last timestep. This must have size _num_basis + _num_kin. Only the first _num_basis
345 : * of these are used.
346 : */
347 : void updateOldWithCurrent(const DenseVector<Real> & mole_additions);
348 :
349 : /**
350 : * @return the number in the algebraic system (number of basis species in algebraic system +
351 : * number of surface potentials + number of kinetic species)
352 : */
353 : unsigned getNumInAlgebraicSystem() const;
354 :
355 : /// return the number of basis species in the algebraic system
356 : unsigned getNumBasisInAlgebraicSystem() const;
357 :
358 : /// return the number of surface potentials
359 : unsigned getNumSurfacePotentials() const;
360 :
361 : /**
362 : * @return a vector of length _num_basis whose entries determine whether the basis species is in
363 : * the algebraic system
364 : */
365 : const std::vector<bool> & getInAlgebraicSystem() const;
366 :
367 : /**
368 : * @return v, a vector of length _num_basis, where v[i] = the basis index of the i^th species in
369 : the algebraic system. Note that for i > _num_basis_in_algebraic_system, the value of v is
370 : undefined.
371 : */
372 : const std::vector<unsigned> & getBasisIndexOfAlgebraicSystem() const;
373 :
374 : /**
375 : * @return v, a vector of length _num_basis, where v[i] = the algebraic index of the i^th basis
376 : * species. Note that unless the i^th basis species is in the algebraic system, this is
377 : * meaningless
378 : */
379 : const std::vector<unsigned> & getAlgebraicIndexOfBasisSystem() const;
380 :
381 : /**
382 : * @return v, a vector of length getNumInAlgebraicSystem, the elements of which are the current
383 : * values of the algebraic variables
384 : */
385 : std::vector<Real> getAlgebraicVariableValues() const;
386 :
387 : /**
388 : * @return v, a vector of length _num_basis_in_algebraic_system, the elements of which are the
389 : * current values of the algebraic variables: molalities only (not surface potentials or kinetic
390 : * species)
391 : */
392 : std::vector<Real> getAlgebraicBasisValues() const;
393 :
394 : /**
395 : * @return v, a DenseVector of length getNumInAlgebraicSystem, the elements of which are the
396 : * current values of the algebraic variables
397 : */
398 : DenseVector<Real> getAlgebraicVariableDenseValues() const;
399 :
400 : /// Returns the mass of solvent water
401 : Real getSolventWaterMass() const;
402 :
403 : /**
404 : * @return the number of bulk-composition moles. Note this will typically be the "old" number of
405 : * bulk moles (from the previous time-step) unless addtoBulkMoles or updateOldWithCurrent or other
406 : * similar methods have just been called. Note that this contains contributions from kinetic
407 : * species, in contrast to the Bethke approach.
408 : */
409 : const std::vector<Real> & getBulkMolesOld() const;
410 :
411 : /**
412 : * @return the number of bulk-composition moles in the original basis. Note this will typically
413 : * be the "old" number of bulk moles (from the previous time-step) unless addtoBulkMoles or
414 : * updateOldWithCurrent or other similar methods have just been called. Note that this contains
415 : * contributions from kinetic species, in contrast to the Bethke approach.
416 : */
417 : DenseVector<Real> getBulkOldInOriginalBasis() const;
418 :
419 : /**
420 : * @return the number of bulk-composition moles that are transported in reactive-transport,
421 : * expressed in the original basis. Note this is computed using the existing molalities, so the
422 : * result might be junk if the system is inconsistent, but will be OK if, for instance, a solve
423 : * has just converged. Also note that this does not include contributions from kinetic species
424 : */
425 : DenseVector<Real> getTransportedBulkInOriginalBasis() const;
426 :
427 : /**
428 : * Computes the value of transported bulk moles for all basis species using the existing
429 : * molalities. Note that this is probably gives rubbish results unless the system is consistent
430 : * (eg, the solve has converged). Also note that this does not include contributions from kinetic
431 : * species
432 : */
433 : void computeTransportedBulkFromMolalities(std::vector<Real> & transported_bulk) const;
434 :
435 : /**
436 : * @return vector v, where v[i] = mass of solvent water (i=0), or v[i] = molality of the basis
437 : * aqueous species i, or v[i] = free moles of basis mineral i, whichever is appropriate
438 : */
439 : const std::vector<Real> & getSolventMassAndFreeMolalityAndMineralMoles() const;
440 :
441 : /**
442 : * @return vector v, where v[i] = true if the activity of basis species i is fixed, either by the
443 : * user providing an activity value, or a fugacity value, or because the i^th basis species is a
444 : * mineral
445 : */
446 : const std::vector<bool> & getBasisActivityKnown() const;
447 :
448 : /**
449 : * @return the activity for the i^th basis species
450 : */
451 : Real getBasisActivity(unsigned i) const;
452 :
453 : /**
454 : * @return the activity for the basis species
455 : */
456 : const std::vector<Real> & getBasisActivity() const;
457 :
458 : /**
459 : * @return the molality of the j^th equilibrium species. This is not defined for minerals or
460 : * gases
461 : */
462 : Real getEquilibriumMolality(unsigned j) const;
463 :
464 : /**
465 : * @return the molalities of the equilibrium species. These are not defined for minerals or
466 : * gases
467 : */
468 : const std::vector<Real> & getEquilibriumMolality() const;
469 :
470 : /**
471 : * @return the activity coefficient for the i^th basis species. This is not defined for water,
472 : * minerals, gases, or aqueous species that have been provided an activity by the user
473 : */
474 : Real getBasisActivityCoefficient(unsigned i) const;
475 :
476 : /**
477 : * @return the activity coefficients for the basis species. These are not defined for water,
478 : * minerals, gases, or aqueous species that have been provided an activity by the user
479 : */
480 : const std::vector<Real> & getBasisActivityCoefficient() const;
481 :
482 : /**
483 : * @return the activity coefficient for the j^th equilibrium species. This is not defined for
484 : * minerals
485 : */
486 : Real getEquilibriumActivityCoefficient(unsigned j) const;
487 :
488 : /**
489 : * @return the activity coefficients for the equilibrium species. These are not defined for
490 : * minerals
491 : */
492 : const std::vector<Real> & getEquilibriumActivityCoefficient() const;
493 :
494 : /**
495 : * @return the total charge in the system. Note this is based on the old values of bulk moles and
496 : * kinetic moles (ie, from the previous time-step) so to get the current values methods like
497 : * addToBulkMoles should be called, or updateOldWithCurrent
498 : */
499 : Real getTotalChargeOld() const;
500 :
501 : /**
502 : * Return the residual of the algebraic system for the given algebraic index
503 : * @param algebraic_ind the algebraic index
504 : * @param mole_additions the increment of mole number of each basis species and kinetic species
505 : * since the last timestep. This must have size _num_basis + _num_kin. For the basis species,
506 : * this is the amount of each species being injected into the system over the timestep. For the
507 : * kinetic species, this is -dt*reaction_rate. Please note: do not decompose the kinetic-species
508 : * additions into basis components and add them to the first slots of mole_additions! This method
509 : * does that decomposition automatically. The first _num_basis slots of mole_additions contain
510 : * kinetic-independent additions, while the last _num_kin slots contain kinetic-rate
511 : * contributions.
512 : */
513 : Real getResidualComponent(unsigned algebraic_ind, const DenseVector<Real> & mole_additions) const;
514 :
515 : /**
516 : * @return reference to the underlying ModelGeochemicalDatabase
517 : */
518 : const ModelGeochemicalDatabase & getModelGeochemicalDatabase() const;
519 :
520 : /**
521 : * Copies a ModelGeochemicalDatabase into our _mgd structure
522 : * @param mgd reference to the ModelGeochemicalDatabase that will be copied into _mgd
523 : */
524 : void setModelGeochemicalDatabase(const ModelGeochemicalDatabase & mgd);
525 :
526 : /**
527 : * Compute the Jacobian for the algebraic system and put it in jac
528 : * @param res The residual of the algebraic system. This is used to speed up computations of the
529 : * jacobian
530 : * @param jac The jacobian entries are placed here
531 : * @param mole_additions The molar additions to the basis species and the kinetic species
532 : * @param dmole_additions d(mole_additions)/d(molality or kinetic_moles)
533 : */
534 : void computeJacobian(const DenseVector<Real> & res,
535 : DenseMatrix<Real> & jac,
536 : const DenseVector<Real> & mole_additions,
537 : const DenseMatrix<Real> & dmole_additions) const;
538 :
539 : /**
540 : * Set the variables in the algebraic system (molalities and potentially the mass of solvent
541 : * water, and surface potentials, and kinetic mole numbers if any) to algebraic_var. All elements
542 : * of this vector must be postivie. This function also calls computeConsistentConfiguration.
543 : * @param algebraic_var the values to set the algebraic variables to
544 : */
545 : void setAlgebraicVariables(const DenseVector<Real> & algebraic_var);
546 :
547 : /**
548 : * Enforces charge balance by altering the constraint_value and bulk_moles_old of the
549 : * charge-balance species. Use with caution, since this overwrites the constraint values provided
550 : * in the constructor, and also changes bulk_moles_old which is usually the value of bulk moles
551 : * from the previous time-step
552 : */
553 : void enforceChargeBalance();
554 :
555 : /**
556 : * Changes the charge-balance species to the original that is specified in the constructor (this
557 : * might not be possible since the original charge-balance species might have been swapped out)
558 : * @return true if the charge-balance species is changed
559 : */
560 : bool revertToOriginalChargeBalanceSpecies();
561 :
562 : /**
563 : * @return vector v, where v[j] = saturation index of equilibrium species j = log10(activity
564 : * product) - log10K. If the equilibrium species is not a mineral then v[j] = 0
565 : */
566 : std::vector<Real> getSaturationIndices() const;
567 :
568 : /**
569 : * Perform the basis swap, and ensure that the resulting system is consistent. Note that this
570 : * alters the bulk_moles_old of species. If you are confident that your configuration is
571 : * reasonably correct, you should therefore ensure bulk_moles_old is set to the actual bulk_moles
572 : * in your system prior to calling performSwap. Alternatively, if you are unsure of the
573 : * correctness, just let performSwap use bulk_moles_old to set the bulk moles of your new basis.
574 : * @param swap_out_of_basis index of basis species to remove from basis
575 : * @param swap_into_basis index of equilibrium species to add to basis
576 : */
577 : void performSwap(unsigned swap_out_of_basis, unsigned swap_into_basis);
578 :
579 : /// Get the ionic strength
580 : Real getIonicStrength() const;
581 :
582 : /// Get the stoichiometric ionic strength
583 : Real getStoichiometricIonicStrength() const;
584 :
585 : /**
586 : * @param sp surface number of interest. sp < _num_surface_pot
587 : * @return the surface potential (units Volts) for the surface number
588 : */
589 : Real getSurfacePotential(unsigned sp) const;
590 :
591 : /**
592 : * @param sp surface number of interest. sp < _num_surface_pot
593 : * @return the specific charge of the surface (units: Coulombs/m^2) for the surface number
594 : */
595 : Real getSurfaceCharge(unsigned sp) const;
596 :
597 : /**
598 : * @return the sorbing surface area (units: m^2) for each sorbing surface
599 : */
600 : const std::vector<Real> & getSorbingSurfaceArea() const;
601 :
602 : /**
603 : * @return the temperature in degC
604 : */
605 : Real getTemperature() const;
606 :
607 : /**
608 : * Sets the temperature to the given quantity. This also:
609 : * - calculates new values of equilibrium constants for equilibrium, redox and kinetic species
610 : * - instructs the activity-coefficient calculator to set its internal temperature and calculate
611 : * new values of the Debye-Huckel parameters.
612 : * However, it does NOT update activity coefficients, basis molalities for species of known
613 : * activities, basis activities, equilibrium molalities, bulk mole number of species with fixed
614 : * free molality, free mineral mole numbers, or surface sorbing area. Hence, the state of the
615 : * equilibrium geochemical system will be left inconsistent after a call to setTemperature, and
616 : * you should computeConsistentConfiguration() to rectify this, if needed.
617 : * @param temperature the temperature in degC
618 : */
619 : void setTemperature(Real temperature);
620 :
621 : /**
622 : * Sets:
623 : * - solvent water mass
624 : * - free molality of all basis species and equilibrium species that are not (water or gas or
625 : * mineral)
626 : * - free number of moles of all minerals (if any)
627 : * - surface potential expressions for all surface-sorption sites (if any)
628 : * - mole number and old mole number of the kinetic species (if any)
629 : * Then computes bulk mole composition, activities, activity coefficients and sorbing surface
630 : * areas so that the GeochemicalSystem is kept self-consistent.
631 : *
632 : * This method is designed to be run at the start of a simulation, or during a "recover" operator.
633 : *
634 : * @param names the names of all the basis species, equilibrium species, surface-sorption
635 : * minerals and kinetic species that are in the system. Note, there must be exactly num_basis +
636 : * num_eqm + num_surface_pot + num_kin of these, ie, all species must be provided with a value
637 : * @param values the values of the molalities (or solvent water mass for water, or free number of
638 : * moles for minerals/kinetics, or surface_potential_expr for surface-sorption sites). There must
639 : * be an equal number of these as "names".
640 : * @param contraints_from_molalities whether the constraints initial provided to the
641 : * GeochemicalSystem should be updated by the values provided. This must have size
642 : * num_basis.
643 : *
644 : * This method is reasonably nontrivial:
645 : * - it assumes that temperature-dependent quantities (equilibrium constants, Debye-Huckel
646 : * parameters, etc) have been set prior to calling this method, eg in the constructor or
647 : * setTemperature()
648 : * - it assumes the algebraic info has been built during instantiation, or during swaps, etc
649 : * - it does not assume the value provided are "good", although since they most usually come from
650 : * a previous solve, they are typically pretty "good". It does check for non-negativity: molality
651 : * for basis gases must be zero; molality for basis minerals must be non-negative; molality for
652 : * every other basis species must be positive; molality for equilibrium gases and equilibrium
653 : * minerals are set to zero irrespective of values provided; molality of other equilibrium species
654 : * must be non-negative; surface_potential_expressions must be positive; mole number of kinetic
655 : * species must be positive. (In the previous sentence, "molality" is molality, kg solvent water,
656 : * or free mineral moles, whichever is appropriate.)
657 : * - if constraints_from_molalities is false then: if the original constraint was
658 : * KG_SOLVENT_WATER, FREE_MOLALITY or FREE_MOLES_MINERAL_SPECIES then the constraint take
659 : * precedence over the molality provided to this method, so the molality is ignored and the
660 : * constraint value used instead.
661 : * - if constraints_from_molalities is true then: if the original constraint was KG_SOLVENT_WATER,
662 : * FREE_MOLALITY or FREE_MOLES_MINERAL_SPECIES then the constraint value is set to the value
663 : * provided by to this method; if the original constraint was MOLES_BULK_WATER or
664 : * MOLES_BULK_SPECIES then the constraint value is set to the value computed from the molalities
665 : * provided to this method; if the original constraint was ACTIVITY, then the constraint value is
666 : * set to activity_coefficient * molality_provided_to_this_method
667 : * - Note the possibilities for ignoring the values provided to this method mentioned in the
668 : * preceeding paragraphs (setting to zero for equilibrium minerals and gases, and constraints
669 : * overriding the values provided)!
670 : */
671 : void setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(
672 : const std::vector<std::string> & names,
673 : const std::vector<Real> & values,
674 : const std::vector<bool> & constraints_from_molalities);
675 :
676 : /**
677 : * @return a vector of the constraint meanings for the basis species
678 : */
679 : const std::vector<GeochemicalSystem::ConstraintMeaningEnum> & getConstraintMeaning() const;
680 :
681 : /**
682 : * Changes a KG_SOLVENT_WATER constraint to MOLES_BULK_WATER (if there is such a constraint) and
683 : * all FREE_MOLALITY and FREE_MOLES_MINERAL_SPECIES to MOLES_BULK_SPECIES using the current values
684 : * of _bulk_moles_old. If, after performing these changes, all constraints are MOLES_BULK_WATER
685 : * or MOLES_BULK_SPECIES then charge-balance is performed by altering the bulk_moles_old for the
686 : * charge-balance species.
687 : */
688 : void closeSystem();
689 :
690 : /**
691 : * Changes the constraint to MOLES_BULK_SPECIES (or MOLES_BULK_WATER if basis_ind = 0) for the
692 : * basis index. The constraint value (the bulk number of moles of the species) is computed from
693 : * the current molality values. Note that if basis_ind corresponds to a gas, then changing the
694 : * constraint involves performing a basis swap with a non-gaseous equilibrium species (eg O2(g) is
695 : * removed from the basis and O2(aq) enters in its place). If, after performing the change to
696 : * MOLES_BULK_SPECIES, all constraints are MOLES_BULK_WATER or MOLES_BULK_SPECIES then
697 : * charge-balance is performed by altering the bulk_moles_old for the charge-balance species.
698 : * @param basis_ind the index of the basis species
699 : */
700 : void changeConstraintToBulk(unsigned basis_ind);
701 :
702 : /**
703 : * Changes the constraint to MOLES_BULK_SPECIES (or MOLES_BULK_WATER if basis_ind = 0) for the
704 : * basis index. The constraint value (the bulk number of moles of the species) is set to value.
705 : * It is an error to call this function when basis_ind corresponds to a gas, because changing the
706 : * constraint involves performing a basis swap with a non-gaseous equilibrium species (eg O2(g) is
707 : * removed from the basis and O2(aq) enters in its place), so the bulk number of moles will be
708 : * computed by the swapper: use changeConstraintToBulk(basis_ind) instead. If, after performing
709 : * the change to MOLES_BULK_SPECIES, all constraints are MOLES_BULK_WATER or MOLES_BULK_SPECIES
710 : * then charge-balance is performed by altering the bulk_moles_old for the charge-balance species.
711 : * @param basis_ind the index of the basis species
712 : * @param value the value to set the bulk number of moles to
713 : */
714 : void changeConstraintToBulk(unsigned basis_ind, Real value);
715 :
716 : /**
717 : * Add to the MOLES_BULK_SPECIES (or MOLES_BULK_WATER if basis_ind = 0) for the basis species.
718 : * Note that if the constraint on the basis species is not MOLES_BULK_SPECIES (or
719 : * MOLES_BULK_WATER) then this will have no impact. Note that charge-balance is performed if all
720 : * constraints are BULK-type.
721 : * @param basis_ind the index of the basis species
722 : * @param value the value to add to the bulk number of moles
723 : */
724 : void addToBulkMoles(unsigned basis_ind, Real value);
725 :
726 : /**
727 : * Set the constraint value for the basis species. If molalities or activities are changed, this
728 : * method uses computeConsistentConfiguration to result in a consistent configuration, while if
729 : * bulk composition is altered it uses alterSystemBecauseBulkChanged to result in a consistent
730 : * configuration (charge-balance is performed if all constraints are BULK-type)
731 : * @param basis_ind the index of the basis species
732 : * @param value the value to set the constraint to
733 : */
734 : void setConstraintValue(unsigned basis_ind, Real value);
735 :
736 : /**
737 : * Compute a reasonably consistent configuration that can be used as an initial condition in a
738 : * Newton process to solve the algebraic system. Note that unless _iters_to_make_consistent is
739 : * very large, the resulting configuration will not be completely consistent, but this is
740 : * conventional in geochemical modelling: there are so many unknowns and approximations used in
741 : * the Newton process, that starting from a completely consistent configuration is unimportant.
742 : *
743 : * Before entering this method, it is assumed that::
744 : * - basis molality has been set for all basis species. For species where this is unknown (viz,
745 : * it will be solved for in the Newton process) a reasonable initial condition should be set. Use
746 : * the initBulkAndFree method.
747 : * - basis_activity_known has been set to true for basis species whose activities are specified by
748 : * the user, for basis gases and for basis minerals. Use the buildKnownBasisActivities method.
749 : * - basis_activity has been set for basis species with basis_activity_known = true. Use the
750 : * buildKnownBasisActivities method.
751 : * - equilibrium molality has been pre-initialised (eg, to zero) for all equilibrium species prior
752 : * to entering this method.
753 : *
754 : * This method computes:
755 : * - ionic strength and stoichiometric ionic strength
756 : * - basis activity coefficients (not for water, minerals or gases)
757 : * - equilibrium activity coefficients (not for minerals or gases)
758 : * - basis molality for basis species where the user has specified the activity
759 : * - basis molality for mineral basis species
760 : * - basis activity for all basis species
761 : * - equilibrium molality for all equilibrium species (molality for minerals and gases are set to
762 : * zero)
763 : * - bulk_moles_old for all basis species (set to the constraint_value for BULK-type constraints,
764 : * otherwise computed from the current values of molality)
765 : * - sorbing surface area
766 : *
767 : * It does not compute activity for equilibrium species: use computeAndGetEquilibriumActivity for
768 : * that
769 : */
770 : void computeConsistentConfiguration();
771 :
772 : /**
773 : * @return the original redox left-hand-side after the initial basis swaps
774 : */
775 : const std::string & getOriginalRedoxLHS() const;
776 :
777 : /**
778 : * Perform the basis swap, and ensure that the resulting system is consistent. Note: no
779 : * sanity-checks regarding swapping out the charge-balance species, etc, are made.
780 : * Note that this
781 : * alters the bulk_moles_old of species. If you are confident that your configuration is
782 : * reasonably correct, you should therefore ensure bulk_moles_old is set to the actual bulk_moles
783 : * in your system prior to calling performSwap. Alternatively, if you are unsure of the
784 : * correctness, just let performSwap use bulk_moles_old to set the bulk moles of your new basis.
785 : * @param swap_out_of_basis index of basis species to remove from basis
786 : * @param swap_into_basis index of equilibrium species to add to basis
787 : */
788 : void performSwapNoCheck(unsigned swap_out_of_basis, unsigned swap_into_basis);
789 :
790 : /**
791 : * @return a reference to the swapper used by this object
792 : */
793 : const GeochemistrySpeciesSwapper & getSwapper() const;
794 :
795 : /**
796 : * Set the free mole number of mineral-related species to the value provided. This sets the mole
797 : * number of basis minerals, the molality of unoccupied sorption sites (whether in the basis or
798 : * not), and the molality of sorbed species to the value provided. It does not set the molality
799 : * of equilibrium minerals, which is always zero.
800 : * NOTE: calling this method can easily result in a
801 : * completely inconsistent GeochemicalSystem because no
802 : * computeConsistentConfiguration() is called. If you need a consistent configuration, please
803 : * call that method after calling this one
804 : * @param value the mole number (or molality, for sorption-related species)
805 : */
806 : void setMineralRelatedFreeMoles(Real value);
807 :
808 : /**
809 : * Computes activity for all equilibrium species (_eqm_activity) and returns a reference to the
810 : * vector. The vector _eqm_activity is only used for output purposes (such as recording the
811 : * fugacity of equilibrium gases) so it is only computed by this method and not internally during
812 : * computeConsistentConfiguration.
813 : * @return equilibrium activities
814 : */
815 : const std::vector<Real> & computeAndGetEquilibriumActivity();
816 :
817 : /**
818 : * Returns the value of activity for the equilibrium species with index eqm_index
819 : * @param eqm_index the index in mgd for the equilibrium species of interest
820 : */
821 : Real getEquilibriumActivity(unsigned eqm_ind) const;
822 :
823 : /**
824 : * Computes the kinetic rates and their derivatives based on the current values of molality, mole
825 : * number, etc, and then, using dt, calculates the appropriate mole_additions and dmole_additions.
826 : * NOTE: this *adds* the kinetic contributes to mole_additions and dmole_additions.
827 : * This method is not const because it modifies _eqm_activity for
828 : * equilibrium gases and eqm species H+ and OH- (if there are any).
829 : * @param dt time-step size
830 : * @param mole_additions The kinetic rates multiplied by dt get placed in the last num_kin slots
831 : * @param dmole_additions d(mole_additions)/d(molality or kinetic_moles)
832 : */
833 : void
834 : addKineticRates(Real dt, DenseVector<Real> & mole_additions, DenseMatrix<Real> & dmole_additions);
835 :
836 : private:
837 : /// The minimal geochemical database
838 : ModelGeochemicalDatabase & _mgd;
839 : /// number of basis species
840 : const unsigned _num_basis;
841 : /// number of equilibrium species
842 : const unsigned _num_eqm;
843 : /// number of redox species
844 : const unsigned _num_redox;
845 : /// number of surface potentials
846 : const unsigned _num_surface_pot;
847 : /// number of kinetic species
848 : const unsigned _num_kin;
849 : /// swapper that can swap species into and out of the basis
850 : GeochemistrySpeciesSwapper & _swapper;
851 : /// Species to immediately remove from the basis in favor of _swap_in
852 : const std::vector<std::string> _swap_out;
853 : /// Species to immediately add to the basis in favor of _swap_out
854 : const std::vector<std::string> _swap_in;
855 : /// Object to compute the activity coefficients and activity of water
856 : GeochemistryActivityCoefficients & _gac;
857 : /// Object that provides the ionic strengths
858 : GeochemistryIonicStrength & _is;
859 : /// The species used to enforce charge balance
860 : std::string _charge_balance_species;
861 : /// The species used to enforce charge balance, as provided in the constructor
862 : const std::string _original_charge_balance_species;
863 : /// The index in the list of basis species corresponding to the charge-balance species. This gets altered as the simulation progresses, due to swaps, and alterChargeBalanceSpecies
864 : unsigned _charge_balance_basis_index;
865 : /// Names of the species in that have their values fixed to _constraint_value with meanings _constraint_meaning. In the constructor, this is ordered to have the same ordering as the basis species.
866 : std::vector<std::string> _constrained_species;
867 : /// Numerical values of the constraints on _constraint_species. In the constructor, this is ordered to have the same ordering as the basis species.
868 : std::vector<Real> _constraint_value;
869 : /// Numerical values of the constraints on _constraint_species. In the constructor, this is ordered to have the same ordering as the basis species. Since values can change due to charge-balance, this holds the original values set by the user.
870 : std::vector<Real> _original_constraint_value;
871 : /// Units of the _constraint_value when the GeochemicalSystem is constructed. This is used during the constructor to convert the constraint_values into mole units
872 : std::vector<GeochemistryUnitConverter::GeochemistryUnit> _constraint_unit;
873 : /// The user-defined meaning of the values in _constraint_value. In the constructor, this is ordered to have the same ordering as the basis species. During the process of unit conversion, from the user-supplied _constraint_unit, to mole-based units, _constraint_user_meaning is used to populate _constraint_meaning, and henceforth usually only _constraint_meaning is used in the code
874 : std::vector<ConstraintUserMeaningEnum> _constraint_user_meaning;
875 : /// The meaning of the values in _constraint_value. In the constructor, this is ordered to have the same ordering as the basis species.
876 : std::vector<ConstraintMeaningEnum> _constraint_meaning;
877 : /// equilibrium constant of the equilibrium species
878 : std::vector<Real> _eqm_log10K;
879 : /// equilibrium constant of the redox species
880 : std::vector<Real> _redox_log10K;
881 : /// equilibrium constant of the kinetic species
882 : std::vector<Real> _kin_log10K;
883 : /// number of unknown molalities in the algebraic system. Note: surface potentials (if any) are extra quantities in the algebraic system
884 : unsigned _num_basis_in_algebraic_system;
885 : /// number of unknowns in the algebraic system (includes molalities and surface potentials)
886 : unsigned _num_in_algebraic_system;
887 : /// if _in_algebraic_system[i] == true then we need to solve for the i^th basis-species molality
888 : std::vector<bool> _in_algebraic_system;
889 : /// _algebraic_index[i] = index in the algebraic system of the basis species i. _basis_index[_algebraic_index[i]] = i
890 : std::vector<unsigned> _algebraic_index;
891 : /// _basis_index[i] = index in the basis of the algebraic system species i, for i<num_basis_in_algebraic_system. _basis_index[_algebraic_index[i]] == i
892 : std::vector<unsigned> _basis_index;
893 : /// Number of bulk moles of basis species (this includes contributions from kinetic species, in contrast to Bethke)
894 : std::vector<Real> _bulk_moles_old;
895 : /**
896 : * IMPORTANT: this is
897 : * - number of kg of solvent water as the first element
898 : * - molality of basis aqueous species, or free moles of basis mineral, whichever is appropriate
899 : * - always zero for gases
900 : */
901 : std::vector<Real> _basis_molality;
902 : /// whether basis_activity[i] is fixed by the user
903 : std::vector<bool> _basis_activity_known;
904 : /// values of activity (for water, minerals and aqueous basis species) or fugacity (for gases)
905 : std::vector<Real> _basis_activity;
906 : /// molality of equilibrium species
907 : std::vector<Real> _eqm_molality;
908 : /// basis activity coefficients
909 : std::vector<Real> _basis_activity_coef;
910 : /// equilibrium activity coefficients
911 : std::vector<Real> _eqm_activity_coef;
912 : /**
913 : * equilibrium activities. NOTE: for computational efficiency, these are not computed until
914 : * computeAndGetEqmActivity is called, and they are currently only ever used for output purposes
915 : * and computing kinetic rates
916 : */
917 : std::vector<Real> _eqm_activity;
918 : /**
919 : * surface potential expressions. These are not the surface potentials themselves. Instead they
920 : * are exp(-surface_potential * Faraday / 2 / R / T_k). Hence _surface_pot_expr >= 0 (like
921 : * molalities) and the surface-potential residual is close to linear in _surface_pot_expr if
922 : * equilibrium molalities are large
923 : */
924 : std::vector<Real> _surface_pot_expr;
925 : /// surface areas of the sorbing minerals
926 : std::vector<Real> _sorbing_surface_area;
927 : /// mole number of kinetic species
928 : std::vector<Real> _kin_moles;
929 : /// old mole number of kinetic species
930 : std::vector<Real> _kin_moles_old;
931 : /**
932 : * Iterations to make the initial configuration consistent. Note that because equilibrium
933 : * molality depends on activity and activity coefficients, and water activity and activity
934 : * coefficients depend on molality, it is a nontrivial task to compute all these so that the
935 : * algorithm starts with a consistent initial condition from which to solve the algebraic system.
936 : * Usually the algorithm doesn't even attempt to make a consistent initial condition (a suitable
937 : * default is iters_to_make_consistent=0), because solving the algebraic system includes so many
938 : * approximations anyway.
939 : */
940 : unsigned _iters_to_make_consistent;
941 : /// The temperature in degC
942 : Real _temperature;
943 : /// Minimum molality ever used in an initial guess
944 : Real _min_initial_molality;
945 : /// The left-hand-side of the redox equations in _mgd after initial swaps
946 : std::string _original_redox_lhs;
947 :
948 : /**
949 : * Using the provided value of temperature, build _eqm_log10K for each eqm species and redox
950 : * species
951 : */
952 : void buildTemperatureDependentQuantities(Real temperature);
953 :
954 : /// Builds in_algebraic_system, algebraic_index and basis_index, and sets num_basis_in_algebraic_system appropriately
955 : void buildAlgebraicInfo(std::vector<bool> & in_algebraic_system,
956 : unsigned & num_basis_in_algebraic_system,
957 : unsigned & num_in_algebraic_system,
958 : std::vector<unsigned> & algebraic_index,
959 : std::vector<unsigned> & basis_index) const;
960 :
961 : /**
962 : * based on _constrained_value and _constraint_meaning, populate nw, bulk_moles_old and
963 : * basis_molality with reasonable initial conditions that may be used during the Newton solve of
964 : * the algebraic system
965 : * @param bulk_moles_old bulk composition number of moles of the basis species
966 : * @param basis_molality zeroth component is mass of solvent water, other components are either
967 : * molality of the basis aqueous species or number of moles of basis mineral, whichever is
968 : * appropriate
969 : */
970 : void initBulkAndFree(std::vector<Real> & bulk_moles_old,
971 : std::vector<Real> & basis_molality) const;
972 :
973 : /**
974 : * Fully populates basis_activity_known, which is true if activity has been set by the user, or
975 : * the fugacity has been set by the user, or the basis species is a mineral. Populates only those
976 : * slots in basis_activity for which basis_activity_known == true
977 : */
978 : void buildKnownBasisActivities(std::vector<bool> & basis_activity_known,
979 : std::vector<Real> & basis_activity) const;
980 :
981 : /**
982 : * Compute the activity of water and put in basis_activity[0], and use the _basis_activity_coef
983 : * and _basis_molality to compute the remaining basis activities (for those that are not minerals,
984 : * gases or aqueous species provided with an activity)
985 : */
986 : void computeRemainingBasisActivities(std::vector<Real> & basis_activity) const;
987 :
988 : /**
989 : * For basis aqueous species (not water, minerals or gases) with activity set by the user, set
990 : * basis_molality = activity / activity_coefficient
991 : */
992 : void updateBasisMolalityForKnownActivity(std::vector<Real> & basis_molality) const;
993 :
994 : /**
995 : * @return log10(activity product) for equilibrium species eqm_j
996 : */
997 : Real log10ActivityProduct(unsigned eqm_j) const;
998 :
999 : /**
1000 : * compute the equilibrium molalities (not for minerals or gases)
1001 : */
1002 : void computeEqmMolalities(std::vector<Real> & eqm_molality) const;
1003 :
1004 : /**
1005 : * If all non-charged basis species are provided with a bulk number of moles, alter the bulk
1006 : * number of moles (old) of the charge_balance_species so that charges are balanced
1007 : * @param constraint_value the values of the constraints provided by the user
1008 : * @param bulk_moles_old the old values of bulk moles (from previous time-step)
1009 : */
1010 : void enforceChargeBalanceIfSimple(std::vector<Real> & constraint_value,
1011 : std::vector<Real> & bulk_moles_old) const;
1012 :
1013 : /**
1014 : * Computes the value of bulk_moles_old for all basis species. For species with BULK-type
1015 : * constraints bulk_moles_old gets set to the constraint value, while for other species
1016 : * bulk_moles_old is computed from the molalities, and, in contrast to Bethke, the kinetic species
1017 : */
1018 : void computeBulk(std::vector<Real> & bulk_moles_old) const;
1019 :
1020 : /**
1021 : * Enforces charge balance by altering the constraint_value and bulk_moles_old of the
1022 : * charge-balance species. Note that this changes the important constraint_value and the old
1023 : * bulk_moles_old (usually from the previous time-step) so should be used with caution.
1024 : */
1025 : void enforceChargeBalance(std::vector<Real> & constraint_value,
1026 : std::vector<Real> & bulk_moles_old) const;
1027 :
1028 : /**
1029 : * Compute the free mineral moles (ie, basis_molality for basis species that are minerals), using
1030 : * the bulk_moles_old. Note that usually computeBulk should have been called immediately
1031 : * preceding this in order to correctly set bulk_moles_old.
1032 : */
1033 : void computeFreeMineralMoles(std::vector<Real> & basis_molality) const;
1034 :
1035 : /**
1036 : * Used during construction: checks for sane inputs and initializes molalities, etc, using the
1037 : * initialize() method
1038 : * @param kin_name names of kinetic species
1039 : * @param kin_initial initial mole numbers (or mass or volume, depending on kin_unit) of the
1040 : * species named in kin_name
1041 : * @param kin_unit units of the numbers provided in kin_initial
1042 : */
1043 : void
1044 : checkAndInitialize(const std::vector<std::string> & kin_name,
1045 : const std::vector<Real> & kin_initial,
1046 : const std::vector<GeochemistryUnitConverter::GeochemistryUnit> & kin_unit);
1047 :
1048 : /**
1049 : * Set the charge-balance species to the basis index provided. No checks are made on the
1050 : * sanity of the desired change. Before setting, the _constraint_value mole number of the old
1051 : * charge-balance species is set to the value provided in the constructor. Then
1052 : * _charge_balance_basis_index and _charge_balance_species is set appropriately
1053 : */
1054 : void setChargeBalanceSpecies(unsigned new_charge_balance_index);
1055 :
1056 : /**
1057 : * @param eqm_j the index of the equilibrium species
1058 : * @return the modifier to the mass-balance equation for equilibrium species j, which is
1059 : * exp(-z * psi * F / R / TK), where z is the charge of the equilibrium species j, psi is the
1060 : * surface potential relevant to the equilibrium species j, F the Faraday constant, R the gas
1061 : * constant, and TK the temperature in Kelvin. If equilibrium species j has no relevant surface
1062 : * potential, unity is returned. Note that exp(-z * psi * F / R / TK) = (_surface_pot_expr)^(2z)
1063 : */
1064 : Real surfaceSorptionModifier(unsigned eqm_j) const;
1065 :
1066 : /**
1067 : * @param sp Surface potential number. sp < _num_surface_pot
1068 : * @return (1/2) * A / F * sqrt(R * T_k * eps * eps_0 * rho * I). This appears in the
1069 : * surface-potential equation.
1070 : */
1071 : Real surfacePotPrefactor(unsigned sp) const;
1072 :
1073 : /// Compute sorbing_surface_area depending on the current molality of the sorbing minerals
1074 : void computeSorbingSurfaceArea(std::vector<Real> & sorbing_surface_area) const;
1075 :
1076 : /**
1077 : * @param basis_ind the index of the basis species in mgd
1078 : * @return the number of bulk moles of the species basis_ind. This depends on the current
1079 : * molality of the basis_ind species, as well as the current molalities of the equilibrium species
1080 : * that depend on this basis species, and, in contrast to the Bethke approach, the current mole
1081 : * number of kinetic species that depend on this basis species.
1082 : */
1083 : Real computeBulkFromMolalities(unsigned basis_ind) const;
1084 :
1085 : /**
1086 : * Alter the GeochemicalSystem to reflect changes in bulk composition constraints that
1087 : * occur through, for instance, setConstraintValue or addToBulkMoles or changeConstraintToBulk
1088 : * (the latter because charge neutrality might be easily enforced, changing constraints and hence
1089 : * bulk moles). This method enforces charge neutrality if simple (ie, all constraints are
1090 : * BULK-type), then sets _bulk_moles_old, and free mineral moles appropriately
1091 : */
1092 : void alterSystemBecauseBulkChanged();
1093 : };
|
