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 "GeochemicalDatabaseReader.h"
13 : #include "GeochemistryKineticRateCalculator.h"
14 : #include <unordered_map>
15 : #include "DenseMatrix.h"
16 : #include <libmesh/dense_vector.h>
17 :
18 : /**
19 : * Data structure designed to hold information related to sorption via surface complexation.
20 : * Since this is uncommon in geochemistry models, and this information only needs to be retrieved
21 : * once per node at the start of the Newton process, it uses an std::map, which is slow
22 : * compared with a std::vector, but saves memory compared with storing a lot of "zeroes" for every
23 : * non-sorbing species in the model
24 : */
25 351 : struct SurfaceComplexationInfo
26 : {
27 216 : SurfaceComplexationInfo(){};
28 :
29 : bool operator==(const SurfaceComplexationInfo & rhs) const
30 : {
31 1 : return surface_area == rhs.surface_area && sorption_sites == rhs.sorption_sites;
32 : };
33 :
34 : Real surface_area;
35 : std::map<std::string, Real> sorption_sites;
36 : };
37 :
38 : /**
39 : * A single rate expression for the kinetic species with index kinetic_species_index.
40 : * @param kinetic_species_index index of the kinetic species that is governed by this rate
41 : * @param promoting_indices the kinetic rate is multiplied by the produce over all basis and
42 : * equilibrium species of m^promoting_indices / (m^promoting_indices +
43 : * promoting_half_saturation^promoting_indices)^promoting_monod_indices, where m is the molality or
44 : * activity of the species
45 : * @param promoting_monod_indices the kinetic rate is multiplied by the produce over all basis and
46 : * equilibrium species of m^promoting_indices / (m^promoting_indices +
47 : * promoting_half_saturation^promoting_indices)^promoting_monod_indices, where m is the molality or
48 : * activity of the species
49 : * @param promoting_half_saturation the kinetic rate is multiplied by the produce over all basis and
50 : * equilibrium species of m^promoting_indices / (m^promoting_indices +
51 : * promoting_half_saturation^promoting_indices)^promoting_monod_indices, where m is the molality or
52 : * activity of the species
53 : * @param progeny_index the index of the basis or equilibrium species in the current basis that is
54 : * produced by the kinetic reaction (usually this is 0, and description.progeny_efficiency = 0, so
55 : * there are no progeny effects)
56 : * @param description the KineticRateUserDescription of this rate
57 : */
58 : struct KineticRateDefinition
59 : {
60 108 : KineticRateDefinition(unsigned kinetic_species_index,
61 : const std::vector<Real> & promoting_indices,
62 : const std::vector<Real> & promoting_monod_indices,
63 : const std::vector<Real> & promoting_half_saturation,
64 : unsigned progeny_index,
65 : const KineticRateUserDescription & description)
66 108 : : kinetic_species_index(kinetic_species_index),
67 108 : promoting_indices(promoting_indices),
68 108 : promoting_monod_indices(promoting_monod_indices),
69 108 : promoting_half_saturation(promoting_half_saturation),
70 108 : progeny_index(progeny_index),
71 108 : description(description){};
72 :
73 0 : bool operator==(const KineticRateDefinition & rhs) const
74 : {
75 0 : return (kinetic_species_index == rhs.kinetic_species_index) &&
76 0 : (promoting_indices == rhs.promoting_indices) &&
77 0 : (promoting_monod_indices == rhs.promoting_monod_indices) &&
78 0 : (promoting_half_saturation == rhs.promoting_half_saturation) &&
79 0 : (progeny_index == rhs.progeny_index) && (description == rhs.description);
80 : };
81 :
82 : unsigned kinetic_species_index;
83 : std::vector<Real> promoting_indices;
84 : std::vector<Real> promoting_monod_indices;
85 : std::vector<Real> promoting_half_saturation;
86 : unsigned progeny_index;
87 : KineticRateUserDescription description;
88 : };
89 :
90 : /**
91 : * Data structure to hold all relevant information from the database file.
92 : * Generally, the database file contains information on a lot more species than any numerical model
93 : * considers. The ModelGeochemicalDatabase only holds the minimal information required for the
94 : * numerical model. It also holds the information as std::vector and DenseMatrix data structures
95 : * for numerical efficiency. This makes the ModelGeochemicalDatabase a little more obscure compared
96 : * with the database file, but considering this information is used hundreds, maybe millions, of
97 : * times per node during a single timestep, numerical efficiency is paramount.
98 : */
99 : struct ModelGeochemicalDatabase
100 : {
101 : /**
102 : * Constructor sets original_database. Also initializes swap_to_original_basis to "nothing" in an
103 : * attempt to reduce memory consumption
104 : */
105 1178 : ModelGeochemicalDatabase(const GeochemicalDatabaseReader & db)
106 1178 : : original_database(&db), swap_to_original_basis(DenseMatrix<Real>()){};
107 :
108 2 : bool operator==(const ModelGeochemicalDatabase & rhs) const
109 : {
110 2 : return (original_database == rhs.original_database) &&
111 2 : (basis_species_index == rhs.basis_species_index) &&
112 2 : (basis_species_name == rhs.basis_species_name) &&
113 2 : (basis_species_mineral == rhs.basis_species_mineral) &&
114 2 : (basis_species_gas == rhs.basis_species_gas) &&
115 2 : (basis_species_transported == rhs.basis_species_transported) &&
116 2 : (basis_species_charge == rhs.basis_species_charge) &&
117 2 : (basis_species_radius == rhs.basis_species_radius) &&
118 2 : (basis_species_molecular_weight == rhs.basis_species_molecular_weight) &&
119 2 : (basis_species_molecular_volume == rhs.basis_species_molecular_volume) &&
120 2 : (eqm_species_index == rhs.eqm_species_index) &&
121 2 : (eqm_species_name == rhs.eqm_species_name) &&
122 2 : (eqm_species_mineral == rhs.eqm_species_mineral) &&
123 2 : (eqm_species_gas == rhs.eqm_species_gas) &&
124 2 : (eqm_species_transported == rhs.eqm_species_transported) &&
125 2 : (eqm_species_charge == rhs.eqm_species_charge) &&
126 2 : (eqm_species_radius == rhs.eqm_species_radius) &&
127 2 : (eqm_species_molecular_weight == rhs.eqm_species_molecular_weight) &&
128 2 : (eqm_species_molecular_volume == rhs.eqm_species_molecular_volume) &&
129 2 : (eqm_stoichiometry == rhs.eqm_stoichiometry) && (eqm_log10K == rhs.eqm_log10K) &&
130 2 : (surface_sorption_name == rhs.surface_sorption_name) &&
131 2 : (surface_sorption_area == rhs.surface_sorption_area) &&
132 2 : (surface_sorption_related == rhs.surface_sorption_related) &&
133 2 : (surface_sorption_number == rhs.surface_sorption_number) &&
134 2 : (redox_lhs == rhs.redox_lhs) && (redox_stoichiometry == rhs.redox_stoichiometry) &&
135 2 : (redox_log10K == rhs.redox_log10K) &&
136 2 : (surface_complexation_info == rhs.surface_complexation_info) &&
137 2 : (gas_chi == rhs.gas_chi) && (kin_species_index == rhs.kin_species_index) &&
138 2 : (kin_species_name == rhs.kin_species_name) &&
139 2 : (kin_species_mineral == rhs.kin_species_mineral) &&
140 2 : (kin_species_transported == rhs.kin_species_transported) &&
141 2 : (kin_species_charge == rhs.kin_species_charge) &&
142 2 : (kin_species_molecular_weight == rhs.kin_species_molecular_weight) &&
143 2 : (kin_species_molecular_volume == rhs.kin_species_molecular_volume) &&
144 2 : (kin_log10K == rhs.kin_log10K) && (kin_stoichiometry == rhs.kin_stoichiometry) &&
145 2 : (kin_rate == rhs.kin_rate) &&
146 2 : (have_swapped_out_of_basis == rhs.have_swapped_out_of_basis) &&
147 4 : (have_swapped_into_basis == rhs.have_swapped_into_basis) &&
148 2 : (swap_to_original_basis == rhs.swap_to_original_basis);
149 : };
150 :
151 : /// a pointer to the original database used to build this ModelGeochemicalDatabase
152 : const GeochemicalDatabaseReader * original_database;
153 :
154 : /**
155 : * basis_species_index[name] = index of the basis species, within all ModelGeochemicalDatabase
156 : * internal datastrcutres, with given name
157 : */
158 : std::unordered_map<std::string, unsigned> basis_species_index;
159 :
160 : /// basis_species_name[j] = name of the j^th basis species
161 : std::vector<std::string> basis_species_name;
162 :
163 : /// basis_species_mineral[j] = true iff the j^th basis species is a mineral
164 : std::vector<bool> basis_species_mineral;
165 :
166 : /// basis_species_gas[j] = true iff the j^th basis species is a gas
167 : std::vector<bool> basis_species_gas;
168 :
169 : /// basis_species_transported[j] = true iff the j^th basis species is transported in reactive-transport sims
170 : std::vector<bool> basis_species_transported;
171 :
172 : /// all quantities have a charge (mineral charge = 0, gas charge = 0, oxide charge = 0)
173 : std::vector<Real> basis_species_charge;
174 :
175 : /// all quantities have an ionic radius (Angstrom) for computing activity (mineral radius = 0, gas radius = 0, surface species radius = 0)
176 : std::vector<Real> basis_species_radius;
177 :
178 : /// all quantities have a molecular weight (g)
179 : std::vector<Real> basis_species_molecular_weight;
180 :
181 : /// all quantities have a molecular volume (cm^3) (only nonzero for minerals, however)
182 : std::vector<Real> basis_species_molecular_volume;
183 :
184 : /**
185 : * eqm_species_index[name] = index of the equilibrium species (secondary aqueous species, redox
186 : * couples in equilibrium with the aqueous solution, minerals in equilibrium with the aqueous
187 : * solution, gases in equilibrium with the aqueous solution) within all ModelGeochemicalDatabase
188 : * internal datastrcutres, with given name
189 : */
190 : std::unordered_map<std::string, unsigned> eqm_species_index;
191 :
192 : /// eqm_species_name[i] = name of the i^th eqm species
193 : std::vector<std::string> eqm_species_name;
194 :
195 : /// eqm_species_mineral[i] = true iff the i^th equilibrium species is a mineral
196 : std::vector<bool> eqm_species_mineral;
197 :
198 : /// eqm_species_gas[i] = true iff the i^th equilibrium species is a gas
199 : std::vector<bool> eqm_species_gas;
200 :
201 : /// eqm_species_transported[i] = true iff the i^th eqm species is transported in reactive-transport sims
202 : std::vector<bool> eqm_species_transported;
203 :
204 : /// all quantities have a charge (mineral charge = 0, gas charge = 0, oxide charge = 0)
205 : std::vector<Real> eqm_species_charge;
206 :
207 : /// all quantities have an ionic radius (Angstrom) for computing activity (mineral radius = 0, gas radius = 0, surface species radius = 0)
208 : std::vector<Real> eqm_species_radius;
209 :
210 : /// all quantities have a molecular weight (g)
211 : std::vector<Real> eqm_species_molecular_weight;
212 :
213 : /// all quantities have a molecular volume (cm^3) (only nonzero for minerals, however)
214 : std::vector<Real> eqm_species_molecular_volume;
215 :
216 : /**
217 : * eqm_stoichiometry(i, j) = stoichiometric coefficient for equilibrium species "i" in terms of
218 : * the basis species "j"
219 : */
220 : DenseMatrix<Real> eqm_stoichiometry;
221 :
222 : /**
223 : * eqm_log10K(i, j) = log10(equilibrium constant) for i^th equilibrium species at the j^th
224 : * temperature point
225 : */
226 : DenseMatrix<Real> eqm_log10K;
227 :
228 : /**
229 : * surface_sorption_name[k] = name of the mineral involved in surface sorption. Each of these
230 : * will be associated with a unique surface potential.
231 : */
232 : std::vector<std::string> surface_sorption_name;
233 :
234 : /**
235 : * surface_sorption_area[k] = specific surface area [m^2/g] for the k^th mineral involved in
236 : * surface sorption. Each mineral involved in surface sorption must have a surface_sorption_area
237 : * prescribed it (in the database) and will be associated with a unique surface potential.
238 : */
239 : std::vector<Real> surface_sorption_area;
240 :
241 : /// surface_sorption_related[j] = true iff the j^th equilibrium species is involved in surface sorption
242 : std::vector<bool> surface_sorption_related;
243 :
244 : /**
245 : * surface_sorption_number[j] = the index of the surface potential that should be used to modify
246 : * the equilibrium constant for the j^th equilibrium species. surface_sorption_number is only
247 : * meaningful if surface_sorption_related[j] = true. 0 <= surface_sorption_number[:] <
248 : * surface_sorption_name.size() = number of minerals involved in surface sorption = number of
249 : * surface potentials in the simulation = surface_sorption_area.size().
250 : */
251 : std::vector<unsigned> surface_sorption_number;
252 :
253 : /**
254 : * the name of the species on the left-hand side of the redox equations. Upon creation of the
255 : * model this is e-, but it may change due to swaps
256 : */
257 : std::string redox_lhs;
258 :
259 : /**
260 : * redox_stoichiometry(i, j) = stoichiometric coefficients for i^th redox species that is in
261 : * disequilibrium in terms of the j basis species. These equations are all written with their
262 : * left-hand sides being e-.
263 : * For instance, the database contains the reactions
264 : * Fe+++ = -0.5H20 + Fe++ + H+ + 0.25*O2(aq)
265 : * e- = 0.5H2O - 0.25O2(aq)- H+
266 : * The first is a redox reaction, and assume the user has specified it is in disquilibrium (by
267 : * specifying it is in the basis). If Fe++, H+ and O2(aq) are also in the basis then the two
268 : * equations yield the single redox equation
269 : * e- = -Fe+++ + Fe++
270 : * Then redox_stoichiometry will be -1 for j = Fe+++, and +1 for j = Fe++.
271 : */
272 : DenseMatrix<Real> redox_stoichiometry;
273 :
274 : /**
275 : * redox_log10K(i, j) = log10(equilibrium constant) for i^th redox species at the j^th temperature
276 : * point
277 : */
278 : DenseMatrix<Real> redox_log10K;
279 :
280 : /**
281 : * Holds info on surface complexation, if any, in the model. Note this is slow compared with
282 : * storing information in a std::vector or DenseMatrix, but it also saves memory over storing a
283 : * lot of "zeroes" for all species not involved in surface complexation
284 : */
285 : std::unordered_map<std::string, SurfaceComplexationInfo> surface_complexation_info;
286 :
287 : /**
288 : * Holds info on gas fugacity "chi" parameters. Note this is slow compared with storing info in a
289 : * DenseMatrix, but it also saves memory over storing a lot of "zeroes" for non-gas species. I
290 : * believe the slowness will not present a problem because fugacity parameters are typically only
291 : * queried at the start of the simulation, or at most once per timestep. If the lookup becomes
292 : * burdensome, change from unordered_map to DenseMatrix
293 : */
294 : std::unordered_map<std::string, std::vector<Real>> gas_chi;
295 :
296 : /**
297 : * kin_species_index[name] = index of the kinetic species, within all ModelGeochemicalDatabase
298 : * internal datastrcutres, with given name
299 : */
300 : std::unordered_map<std::string, unsigned> kin_species_index;
301 :
302 : /// kin_species_name[j] = name of the j^th kinetic species
303 : std::vector<std::string> kin_species_name;
304 :
305 : /// kin_species_mineral[j] = true iff the j^th kinetic species is a mineral
306 : std::vector<bool> kin_species_mineral;
307 :
308 : /// kin_species_transported[j] = true iff the j^th kinetic species is transported in reactive-transport sims
309 : std::vector<bool> kin_species_transported;
310 :
311 : /// all kinetic quantities have a charge (mineral charge = 0)
312 : std::vector<Real> kin_species_charge;
313 :
314 : /// all quantities have a molecular weight (g/mol)
315 : std::vector<Real> kin_species_molecular_weight;
316 :
317 : /// all quantities have a molecular volume (cm^3/mol) (only nonzero for minerals, however)
318 : std::vector<Real> kin_species_molecular_volume;
319 :
320 : /**
321 : * kin_log10K(i, j) = log10(equilibrium constant for the i^th kinetic species at the j^th
322 : * temperature point
323 : */
324 : DenseMatrix<Real> kin_log10K;
325 :
326 : /**
327 : * kin_stoichiometry(i, j) = stoichiometric coefficient for kinetic species "i" in terms of
328 : * the basis species "j"
329 : */
330 : DenseMatrix<Real> kin_stoichiometry;
331 :
332 : /**
333 : * rates given to kinetic species. See the method addKineticRate for a detailed description. This
334 : * quantity is organised in such a way that a solver can loop through kin_rate, calculting the
335 : * rates and applying them to the kin_rate[i].kinetic_species_index species
336 : */
337 : std::vector<KineticRateDefinition> kin_rate;
338 :
339 : /**
340 : * Species that have been swapped out of the basis. Every time a swap is performed on the
341 : * ModelGeochemicalDatabase, the basis-index of the species removed from the basis is appended to
342 : * this list
343 : */
344 : std::vector<unsigned> have_swapped_out_of_basis;
345 :
346 : /**
347 : * Species that have been swapped into the basis. Every time a swap is performed on the
348 : * ModelGeochemicalDatabase, the equilibrium-index of the species added to the basis is appended
349 : * to this list
350 : */
351 : std::vector<unsigned> have_swapped_into_basis;
352 :
353 : /**
354 : * Swap matrix that allows expression in terms of the original basis. When a swap is performed
355 : * bulk_new = S^-1^T * bulk_old,
356 : * where S is the swap matrix (and S^-1^T is the transposed inverse of S)
357 : * Hence, upon multiple swaps (S_1 followed by S_2 followed by S_3, etc)
358 : * bulk_new = S_n^-1^T * S_{n-1}^-1^T * ... * S_1^-1^T * bulk_original
359 : * So
360 : * bulk_original = (S_n * S_{n-1} * ... * S_1)^T * bulk_new
361 : * swap_to_original_basis = S_n * S_{n-1} * ... * S_1
362 : * swap_to_original_basis is initialized to DenseMatrix<Real>(), ie, a zero-sized matrix, to
363 : * reduce memory usage. The first time a swap is performed, it is set to the swap matrix
364 : * generated in the GeochemistrySpeciesSwapper, and every time a swap is performed on the
365 : * ModelGeochemicalDatabase this is updated
366 : */
367 : DenseMatrix<Real> swap_to_original_basis;
368 : };
369 :
370 : /**
371 : * Constructs and stores a minimal amount of information that is pertinent to the user-defined
372 : * geochemical system. Most importantly, all basis species, secondary species, mineral species,
373 : * etc, that are defined in the geochemical database but are irrelevant to the user-defined system
374 : * are eliminated from further consideration. This reduces the amount of information considerably.
375 : * The final result is stored in a ModelGeochemicalDatabase structure. This structure is designed
376 : * to be computationally efficient. A "getter" that copies this structure is provided as a public
377 : * method. It is intended that a MOOSE simulation will construct one PertinentGeochemicalSystem
378 : * object, and then copy the information to the nodes during the initial setup. If different nodes
379 : * are allowed different swaps (likely to be the case) each node will need a different copy of the
380 : * ModelGeochemicalDatabase structure.
381 : */
382 : class PertinentGeochemicalSystem
383 : {
384 : public:
385 : /**
386 : * @param db the database reader, which will have parsed the database file
387 : *
388 : * @param basis_species A list of basis components relevant to the aqueous-equilibrium problem.
389 : * "H2O" must appear first in this list. No element must appear more than once in this list.
390 : * These components must be chosen from the "basis species" in the database, the sorbing sites (if
391 : * any) and the decoupled redox states that are in disequilibrium (if any). Any redox pair that
392 : * is not in this list or the kinetic_redox list, will be assumed to be at equilibrium with the
393 : * aqueous solution and will be considered a secondary species. All these species, except H2O,
394 : * may be later swapped out of this list, either by a manual user-prescribed swap (and replaced by
395 : * a mineral or a gas of fixed fugacity, for instance), or during the numerical solve.
396 : *
397 : * @param minerals A list of minerals that are in equilibrium with the aqueous solution. This can
398 : * only include the "minerals" in the database file. No element can appear more than once in this
399 : * list. Their equilibrium reaction must consist of only the basis_species, and
400 : * secondary species and non-kinetically-controlled redox couples that can be expressed in terms
401 : * of the basis_species. If they are also "sorbing minerals" in the database then their sorption
402 : * sites must consist of the basis_species only. During simulation, the user can compute the
403 : * saturation index of these minerals, and these minerals can be "swapped" into the basis if
404 : * desired (or required during the numerical solve). If the user performs a manual "swap" then an
405 : * initial condition must be provided for the mineral. The user choose whether these minerals are
406 : * allowed to precipitate or not - that is, they can be "supressed". This list, along with the
407 : * kinetic_minerals list, comprises the entire list of minerals in the problem: all others are
408 : * eliminated from consideration.
409 : *
410 : * @param gases A list of gases that are in equilibrium with the aqueous solution and can have
411 : * their fugacities fixed, at least at some time and spatial location. All members of this list
412 : * must be a "gas" in the database file. No gas must appear more than once in this list. The
413 : * equilibrium reaction of each gas must involve only the basis_species, or secondary species or
414 : * non-kinetically-controlled redox couples that can be expressed in terms of the basis_species.
415 : *
416 : * @param kinetic_minerals A list of minerals that whose dynamics are governed by a rate law.
417 : * These are not in equilibrium with the aqueous solution. This can only include the "minerals"
418 : * in the database file. No element can appear more than once in this list. Their equilibrium
419 : * reaction must involve only the basis_species, or secondary species or
420 : * non-kinetically-controlled redox couples that can be expressed in terms of the basis_species.
421 : * If they are also "sorbing minerals" in the database then their sorption sites must consist of
422 : * the basis_speices only. No members of this list must be in the minerals list. They can never
423 : * be "swapped" into the basis, nor can they be "supressed".
424 : *
425 : * @param kinetic_redox A list of redox pairs whose dynamics are governed by a rate law.
426 : * These are not in equilibrium with the aqueous solution. Each element of this list must appear
427 : * in the "redox couples" section of the database. No element can appear more than once in this
428 : * list. Their reaction must involve only the basis_species, or secondary species or
429 : * non-kinetically-controlled redox couples that can be expressed in terms of the basis_species.
430 : * No members of this list must be in the basis_species list. They can never be "swapped" into
431 : * the basis.
432 : *
433 : * @param kinetic_surface_species A list of surface sorbing species whose dynamics are governed by
434 : * a rate law. These are not in equilibrium with the aqueous solution. All elements of this list
435 : * must appear as a "surface species" in the database. No member must appear more than twice in
436 : * this list. Their reaction must involve only the basis_species, or secondary species or
437 : * non-kinetically-controlled redox couples that can be expressed in terms of the basis_species.
438 : * They can never be "swapped" into the basis.
439 : *
440 : * @param redox_ox The name of the oxygen species, eg O2(aq), that appears in redox reactions. For
441 : * redox pairs that are in disequilibrium to be correctly recorded, and hence their Nernst
442 : * potentials to be computed easily, redox_ox must be a basis_species and it must appear in the
443 : * reaction for each redox pair.
444 : *
445 : * @param redox_e The name of the free electron, eg e-. For redox pairs that are in
446 : * disequilibrium to be correctly recorded, and hence their Nernst potentials to be computed
447 : * easily, the equilibrium reaction for redox_e must involve redox_ox, and the basis species must
448 : * be chosen so that redox_e is an equilibrium species according to the database reader
449 : */
450 : PertinentGeochemicalSystem(const GeochemicalDatabaseReader & db,
451 : const std::vector<std::string> & basis_species,
452 : const std::vector<std::string> & minerals,
453 : const std::vector<std::string> & gases,
454 : const std::vector<std::string> & kinetic_minerals,
455 : const std::vector<std::string> & kinetic_redox,
456 : const std::vector<std::string> & kinetic_surface_species,
457 : const std::string & redox_ox,
458 : const std::string & redox_e);
459 :
460 : /// Return a reference to the ModelGeochemicalDatabase structure
461 : const ModelGeochemicalDatabase & modelGeochemicalDatabase() const;
462 :
463 : /**
464 : * Adds a rate description for kinetic_species. Note that a single kinetic species can have
465 : * multiple rates prescribed to it (by calling this method multiple times): they are added
466 : * together to give an overall rate.
467 : */
468 : void addKineticRate(const KineticRateUserDescription & description);
469 :
470 : /**
471 : * @param name species name
472 : * @return the index of the species in the original basis
473 : */
474 : unsigned getIndexOfOriginalBasisSpecies(const std::string & name) const;
475 :
476 : /**
477 : * @return a vector of names of the original basis species
478 : */
479 : std::vector<std::string> originalBasisNames() const;
480 :
481 : private:
482 : /// The database
483 : GeochemicalDatabaseReader _db;
484 :
485 : /// given a species name, return its index in the corresponding "info" std::vector
486 : std::unordered_map<std::string, unsigned> _basis_index;
487 :
488 : /// a vector of all relevant species
489 : std::vector<GeochemistryBasisSpecies> _basis_info;
490 :
491 : /// given a species name, return its index in the corresponding "info" std::vector
492 : std::unordered_map<std::string, unsigned> _mineral_index;
493 :
494 : /// a vector of all relevant species
495 : std::vector<GeochemistryMineralSpecies> _mineral_info;
496 :
497 : /// given a species name, return its index in the corresponding "info" std::vector
498 : std::unordered_map<std::string, unsigned> _gas_index;
499 :
500 : /// a vector of all relevant species
501 : std::vector<GeochemistryGasSpecies> _gas_info;
502 :
503 : /// given a species name, return its index in the corresponding "info" std::vector
504 : std::unordered_map<std::string, unsigned> _kinetic_mineral_index;
505 :
506 : /// a vector of all relevant species
507 : std::vector<GeochemistryMineralSpecies> _kinetic_mineral_info;
508 :
509 : /// given a species name, return its index in the corresponding "info" std::vector
510 : std::unordered_map<std::string, unsigned> _kinetic_redox_index;
511 :
512 : /// a vector of all relevant species
513 : std::vector<GeochemistryRedoxSpecies> _kinetic_redox_info;
514 :
515 : /// given a species name, return its index in the corresponding "info" std::vector
516 : std::unordered_map<std::string, unsigned> _kinetic_surface_index;
517 :
518 : /// a vector of all relevant species
519 : std::vector<GeochemistrySurfaceSpecies> _kinetic_surface_info;
520 :
521 : /// given a species name, return its index in the corresponding "info" std::vector
522 : std::unordered_map<std::string, unsigned> _secondary_index;
523 :
524 : /// a vector of all relevant species
525 : std::vector<GeochemistryEquilibriumSpecies> _secondary_info;
526 :
527 : /**
528 : * The name of the oxygen in all disequilibrium-redox equations, eg O2(aq), which must be a basis
529 : * species if disequilibrium redox reactions are to be recorded
530 : */
531 : const std::string _redox_ox;
532 :
533 : /// The name of the free electron involved in redox reactions
534 : const std::string _redox_e;
535 :
536 : /// The important datastructure built by this class
537 : ModelGeochemicalDatabase _model;
538 :
539 : /**
540 : * using the basis_species list, this method builds _basis_index and _basis_info
541 : */
542 : void buildBasis(const std::vector<std::string> & basis_species);
543 :
544 : /**
545 : * using the minerals list, this method builds _mineral_index and _mineral_info, unless minerals =
546 : * {"*"}, in which case buildAllMinerals is used instead
547 : */
548 : void buildMinerals(const std::vector<std::string> & minerals);
549 :
550 : /**
551 : * If minerals = {"*"} then populate _mineral_index and _mineral_info with all relevant minerals
552 : * This is called in the constructor after buildSecondarySpecies and buildKineticMinerals because
553 : * it adds to _mineral_index and _mineral_info all minerals that:
554 : * - are not kinetic minerals (these have been identified by buildKineticMinerals) and
555 : * - whose stoichiometry depends on the basis species specified in the constructor or on secondary
556 : * species (which are redox couples, secondary species and surface species) that depend on these
557 : * basis species. These secondary species have already been identified by buildSecondarySpecies.
558 : */
559 : void buildAllMinerals(const std::vector<std::string> & minerals);
560 :
561 : /**
562 : * using the gas list, this method builds _gas_index and _gas_info
563 : */
564 : void buildGases(const std::vector<std::string> & gases);
565 :
566 : /**
567 : * using the kinetic_minerals list, this method builds _kinetic_mineral_index and
568 : * _kinetic_mineral_info
569 : */
570 : void buildKineticMinerals(const std::vector<std::string> & kinetic_minerals);
571 :
572 : /**
573 : * using the kinetic_redox list, this method builds _kinetic_redox_index and
574 : * _kinetic_redox_info
575 : */
576 : void buildKineticRedox(const std::vector<std::string> & kinetic_redox);
577 :
578 : /**
579 : * using the kinetic_surface list, this method builds _kinetic_surface_index and
580 : * _kinetic_surface_info
581 : */
582 : void buildKineticSurface(const std::vector<std::string> & kinetic_surface);
583 :
584 : /**
585 : * Extract all relevant "redox couples" and "secondary species" and "surface species" from the
586 : * database. These are all species whose reaction involves only the basis_species and are not
587 : * kinetically controlled or already in the basis_species list.
588 : */
589 : void buildSecondarySpecies();
590 :
591 : /**
592 : * @return true if _redox_e appears as a secondary species in the database and that its
593 : * stoichiometry can be expressed in terms of basis species
594 : */
595 : bool checkRedoxe();
596 :
597 : /**
598 : * Check that all minerals in mineral_info have reactions that involve only the
599 : * basis_species or secondary_species. Check that if a mineral in this list is also a
600 : * "sorbing mineral", its sorbing sites are present in the basis_species list
601 : */
602 : void checkMinerals(const std::vector<GeochemistryMineralSpecies> & mineral_info) const;
603 :
604 : /**
605 : * Check that all gases in the "gases" list have reactions that involve only the
606 : * basis_species or secondary_species.
607 : */
608 : void checkGases() const;
609 :
610 : /**
611 : * Check that all kinetic redox species in the _kinetic_redox list have reactions that involve
612 : * only the basis_species or secondary_species.
613 : */
614 : void checkKineticRedox() const;
615 :
616 : /**
617 : * Check that all kinetic surface species in the _kinetic_surface_species list have reactions that
618 : * involve only the basis_species or secondary_species.
619 : */
620 : void checkKineticSurfaceSpecies() const;
621 :
622 : /**
623 : * Fully populate the ModelGeochemicalDatabase
624 : */
625 : void createModel();
626 :
627 : /**
628 : * Extract the stoichiometry and log10K for the _redox_e species. This is called during
629 : * createModel()
630 : * @param redox_e_stoichiometry upon exit will contain the stoichiometry for the _redox_e species
631 : * @param redox_e_log10K upon exit will contain the log10K info for the _redox_e species
632 : */
633 : void buildRedoxeInfo(std::vector<Real> & redox_e_stoichiometry,
634 : std::vector<Real> & redox_e_log10K);
635 : };
|
