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 "nlohmann/json.h"
13 : #include "MooseTypes.h"
14 :
15 : /**
16 : * Data structure for basis (primary) species.
17 : * Members are:
18 : * * species name
19 : * * elemental composition (elements and weights)
20 : * * ionic radius (Angstrom)
21 : * * charge
22 : * * molecular weight (g)
23 : */
24 : struct GeochemistryBasisSpecies
25 : {
26 13252 : GeochemistryBasisSpecies(){};
27 :
28 : std::string name;
29 : std::map<std::string, Real> elements;
30 : Real radius;
31 : Real charge;
32 : Real molecular_weight;
33 : };
34 :
35 : /**
36 : * Data structure for secondary equilibrium species.
37 : * Members are:
38 : * * Species name
39 : * * Basis species and stoichiometric coefficients present
40 : * * Equilibrium constant at given temperatures
41 : * * ionic radius (Angstrom)
42 : * * charge
43 : * * molecular weight (g)
44 : */
45 : struct GeochemistryEquilibriumSpecies
46 : {
47 302720 : GeochemistryEquilibriumSpecies(){};
48 :
49 : std::string name;
50 : std::map<std::string, Real> basis_species;
51 : std::vector<Real> equilibrium_const;
52 : Real radius;
53 : Real charge;
54 : Real molecular_weight;
55 : };
56 :
57 : /**
58 : * Data structure for mineral species.
59 : * Members are:
60 : * * Species name
61 : * * Molecular volume (units?)
62 : * * Basis species and stoichiometric coefficients present
63 : * * Equilibrium constant at given temperatures ()
64 : * * molecular weight (g)
65 : * * specific surface area (m2/g) (only nonzero for sorbing minerals)
66 : * * basis species (sorbing sites) and their site density (only populated for sorbing minerals)
67 : */
68 : struct GeochemistryMineralSpecies
69 : {
70 1966 : GeochemistryMineralSpecies(){};
71 :
72 : std::string name;
73 : Real molecular_volume;
74 : std::map<std::string, Real> basis_species;
75 : std::vector<Real> equilibrium_const;
76 : Real molecular_weight;
77 : Real surface_area;
78 : std::map<std::string, Real> sorption_sites;
79 : };
80 :
81 : /**
82 : * Data structure for mineral species.
83 : * Members are:
84 : * * Species name
85 : * * Basis species and stoichiometric coefficients present
86 : * * Equilibrium constant at given temperatures ()
87 : * * molecular weight (g)
88 : * * Spycher-Reed fugacity coefficients
89 : * * Pcrit, Tcrit and omega - Tsonopoulos fugacity coefficients
90 : */
91 : struct GeochemistryGasSpecies
92 : {
93 343 : GeochemistryGasSpecies(){};
94 :
95 : std::string name;
96 : std::map<std::string, Real> basis_species;
97 : std::vector<Real> equilibrium_const;
98 : Real molecular_weight;
99 : std::vector<Real> chi;
100 : Real Pcrit;
101 : Real Tcrit;
102 : Real omega;
103 : };
104 :
105 : /**
106 : * Data structure for redox species.
107 : * Members are:
108 : * * Species name
109 : * * Basis species and stoichiometric coefficients present
110 : * * Equilibrium constant at given temperatures
111 : * * ionic radius (Angstrom)
112 : * * charge
113 : * * molecular weight (g)
114 : */
115 : struct GeochemistryRedoxSpecies
116 : {
117 29728 : GeochemistryRedoxSpecies(){};
118 :
119 : std::string name;
120 : std::map<std::string, Real> basis_species;
121 : std::vector<Real> equilibrium_const;
122 : Real radius;
123 : Real charge;
124 : Real molecular_weight;
125 : };
126 :
127 : /**
128 : * Data structure for oxide species.
129 : * Members are:
130 : * * Species name
131 : * * Basis species and stoichiometric coefficients present
132 : * * molecular weight (g)
133 : */
134 : struct GeochemistryOxideSpecies
135 : {
136 2 : GeochemistryOxideSpecies(){};
137 :
138 : std::string name;
139 : std::map<std::string, Real> basis_species;
140 : Real molecular_weight;
141 : };
142 :
143 : /**
144 : * Data structure for sorbing surface species.
145 : * * Species name
146 : * * Basis species and stoichiometric coefficients present
147 : * * log10(Equilibrium constant) and dlog10K/dT
148 : * * charge
149 : * * molecular weight (g)
150 : */
151 : struct GeochemistrySurfaceSpecies
152 : {
153 21850 : GeochemistrySurfaceSpecies(){};
154 :
155 : std::string name;
156 : std::map<std::string, Real> basis_species;
157 : Real charge;
158 : Real molecular_weight;
159 : Real log10K;
160 : Real dlog10KdT;
161 : };
162 :
163 : /**
164 : * Data structure for Debye-Huckel activity coefficients.
165 : * Members are:
166 : * * adh: Debye-Huckel a parameter
167 : * * bdh: Debye-Huckel b parameter
168 : * * bdot: Debye-Huckel bdot parameter
169 : */
170 : struct GeochemistryDebyeHuckel
171 : {
172 949 : GeochemistryDebyeHuckel(){};
173 :
174 : std::vector<Real> adh;
175 : std::vector<Real> bdh;
176 : std::vector<Real> bdot;
177 : };
178 :
179 : /**
180 : * Data structure for elements.
181 : * Members are:
182 : * * name: Element name
183 : * * molecular_weight: Element molecular weight
184 : */
185 9 : struct GeochemistryElements
186 : {
187 2 : GeochemistryElements(){};
188 :
189 : std::string name;
190 : Real molecular_weight;
191 : };
192 :
193 : /**
194 : * Data structure for neutral species activity coefficients.
195 : * Members are a, b, c and d, which are temperature dependent coefficients in the Debye-Huckel
196 : * activity model.
197 : */
198 : struct GeochemistryNeutralSpeciesActivity
199 : {
200 2 : GeochemistryNeutralSpeciesActivity(){};
201 1893 : GeochemistryNeutralSpeciesActivity(std::vector<std::vector<Real>> coeffs)
202 1893 : : a(coeffs[0]), b(coeffs[1]), c(coeffs[2]), d(coeffs[3]){};
203 :
204 : std::vector<Real> a;
205 : std::vector<Real> b;
206 : std::vector<Real> c;
207 : std::vector<Real> d;
208 : };
209 :
210 : /**
211 : * Class for reading geochemical reactions from a MOOSE geochemical database
212 : */
213 : class GeochemicalDatabaseReader
214 : {
215 : public:
216 : /**
217 : * Parse the file
218 : * @param filename Moose geochemical database file
219 : * @param reexpress_free_electron If true, and if the free electron in the database file has an
220 : * equilibrium reaction expressed in terms of O2(g), and O2(g) exists as a gas in the database
221 : * file, and O2(g)'s equilibrium reaction is O2(g)=O2(eq), and O2(aq) exists as a basis species in
222 : * the database file, then reexpress the free electron's equilibrium reaction in terms of O2(aq)
223 : * @param use_piecewise_interpolation If true then set the "logk model" to "piecewise-linear"
224 : * regardless of the value found in the filename. This is designed to make testing easy (because
225 : * logK and Debye-Huckel parameters will be exactly as set in the filename instead of from a 4-th
226 : * order least-squares fit) but should rarely be used for real geochemical simulations
227 : */
228 : GeochemicalDatabaseReader(const FileName filename,
229 : const bool reexpress_free_electron = true,
230 : const bool use_piecewise_interpolation = false,
231 : const bool remove_all_extrapolated_secondary_species = false);
232 :
233 : /**
234 : * Parse the thermodynamic database
235 : * @param filename Name of thermodynamic database file
236 : */
237 : void read(const FileName filename);
238 :
239 : /**
240 : * Validate the thermodynamic database
241 : * @param filename Name of thermodynamic database file
242 : * @param db JSON database read from filename
243 : */
244 : void validate(const FileName filename, const nlohmann::json & db);
245 :
246 : /**
247 : * Sometimes the free electron's equilibrium reaction is defined in terms of O2(g) which is not a
248 : * basis species. If this is the case, re-express it in terms of O2(aq), if O2(g) is a gas and
249 : * O2(aq) is a basis species.
250 : */
251 : void reexpressFreeElectron();
252 :
253 : /**
254 : * Get the activity model type
255 : * @return activity model
256 : */
257 : std::string getActivityModel() const;
258 :
259 : /**
260 : * Get the fugacity model type
261 : * @return fugacity model
262 : */
263 : std::string getFugacityModel() const;
264 :
265 : /**
266 : * Get the equilibrium constant model type
267 : * @return equilibrium constant model
268 : */
269 : std::string getLogKModel() const;
270 :
271 : /**
272 : * Get the list of basis (primary) species read from database
273 : * @return list of primary species names
274 : */
275 : std::vector<std::string> getBasisSpeciesNames() const { return _bs_names; };
276 :
277 : /**
278 : * Get the list of secondary equilibrium species read from database
279 : * @return list of equilibrium species names
280 : */
281 : std::vector<std::string> getEquilibriumSpeciesNames() const { return _es_names; };
282 :
283 : /**
284 : * Get the list of secondary mineral species read from database
285 : * @return list of mineral species names
286 : */
287 : std::vector<std::string> getMineralSpeciesNames() const { return _ms_names; };
288 :
289 : /**
290 : * Get the temperature points that the equilibrium constant is defined at.
291 : * @return vector of temperature points (C)
292 : */
293 : const std::vector<Real> & getTemperatures() const;
294 :
295 : /**
296 : * Get the pressure points that the equilibrium constant is defined at.
297 : * @return vector of pressure points (C)
298 : */
299 : std::vector<Real> getPressures();
300 :
301 : /**
302 : * Get the Debye-Huckel activity coefficients
303 : * @return vectors of adh, bdh and bdot
304 : */
305 : const GeochemistryDebyeHuckel & getDebyeHuckel() const;
306 :
307 : /**
308 : * Get the basis (primary) species information
309 : * @param names list of basis species
310 : * @return basis species structure
311 : */
312 : std::map<std::string, GeochemistryBasisSpecies>
313 : getBasisSpecies(const std::vector<std::string> & names);
314 :
315 : /**
316 : * Get the secondary equilibrium species information
317 : * @param names list of equilibrium species
318 : * @return secondary species structure
319 : */
320 : std::map<std::string, GeochemistryEquilibriumSpecies>
321 : getEquilibriumSpecies(const std::vector<std::string> & names);
322 :
323 : /**
324 : * Get the mineral species information
325 : * @param names list of mineral species
326 : * @return mineral species structure
327 : */
328 : std::map<std::string, GeochemistryMineralSpecies>
329 : getMineralSpecies(const std::vector<std::string> & names);
330 :
331 : /**
332 : * Get all the elements
333 : * @return elements species structure
334 : */
335 : std::map<std::string, GeochemistryElements> getElements();
336 :
337 : /**
338 : * Get the gas species information
339 : * @param names list of gs species
340 : * @return gas species structure
341 : */
342 : std::map<std::string, GeochemistryGasSpecies>
343 : getGasSpecies(const std::vector<std::string> & names);
344 :
345 : /**
346 : * Get the redox species (couples) information
347 : * @param names list of gs species
348 : * @return redox species structure
349 : */
350 : std::map<std::string, GeochemistryRedoxSpecies>
351 : getRedoxSpecies(const std::vector<std::string> & names);
352 :
353 : /**
354 : * Get the oxide species information
355 : * @param names list of gs species
356 : * @return oxide species structure
357 : */
358 : std::map<std::string, GeochemistryOxideSpecies>
359 : getOxideSpecies(const std::vector<std::string> & names);
360 :
361 : /**
362 : * Get the surface sorbing species information
363 : * @param names list of surface sorbing species
364 : * @return surface sorbing species structure
365 : */
366 : std::map<std::string, GeochemistrySurfaceSpecies>
367 : getSurfaceSpecies(const std::vector<std::string> & names);
368 :
369 : /**
370 : * Get the neutral species activity coefficients
371 : * @return neutral species activity coefficients
372 : */
373 : const std::map<std::string, GeochemistryNeutralSpeciesActivity> &
374 : getNeutralSpeciesActivity() const;
375 :
376 : /**
377 : * Generates a formatted vector of strings representing all aqueous equilibrium
378 : * reactions
379 : * @param names list of equilibrium species
380 : * return formatted equilibrium reactions
381 : */
382 : std::vector<std::string> equilibriumReactions(const std::vector<std::string> & names) const;
383 :
384 : /**
385 : * Generates a formatted vector of strings representing all mineral reactions
386 : * @param names list of mineral species
387 : * @preturn formatted mineral reactions
388 : */
389 : std::vector<std::string> mineralReactions(const std::vector<std::string> & names) const;
390 :
391 : /**
392 : * Generates a formatted vector of strings representing all gas reactions
393 : * @param names list of gas species
394 : * @preturn formatted gas reactions
395 : */
396 : std::vector<std::string> gasReactions(const std::vector<std::string> & names) const;
397 :
398 : /**
399 : * Generates a formatted vector of strings representing all redox reactions
400 : * @param names list of redox species
401 : * @preturn formatted redox reactions
402 : */
403 : std::vector<std::string> redoxReactions(const std::vector<std::string> & names) const;
404 :
405 : /**
406 : * Generates a formatted vector of strings representing all oxide reactions
407 : * @param names list of oxide species
408 : * @preturn formatted oxide reactions
409 : */
410 : std::vector<std::string> oxideReactions(const std::vector<std::string> & names) const;
411 :
412 : /**
413 : * String representation of JSON species object contents
414 : * @param name name of species
415 : * @return styled string of species information
416 : */
417 : std::string getSpeciesData(const std::string name) const;
418 :
419 : /**
420 : * Filename of database
421 : * @return filename
422 : */
423 : const FileName & filename() const;
424 :
425 : /// Returns true if name is a "secondary species" or "free electron" in the database
426 : bool isSecondarySpecies(const std::string & name) const;
427 :
428 : /**
429 : * Checks if species is of given type
430 : * @param name species name
431 : * @return true iff species is of given type
432 : */
433 : bool isBasisSpecies(const std::string & name) const;
434 : bool isRedoxSpecies(const std::string & name) const;
435 : bool isGasSpecies(const std::string & name) const;
436 : bool isMineralSpecies(const std::string & name) const;
437 : bool isOxideSpecies(const std::string & name) const;
438 : bool isSurfaceSpecies(const std::string & name) const;
439 :
440 : /// returns True iff name is the name of a sorbing mineral
441 : bool isSorbingMineral(const std::string & name) const;
442 :
443 : /// Returns a list of all the names of the "mineral species" in the database
444 : std::vector<std::string> mineralSpeciesNames() const;
445 :
446 : /// Returns a list of all the names of the "secondary species" and "free electron" in the database
447 : std::vector<std::string> secondarySpeciesNames() const;
448 :
449 : /// Returns a list of all the names of the "redox couples" in the database
450 : std::vector<std::string> redoxCoupleNames() const;
451 :
452 : /// Returns a list of all the names of the "surface species" in the database
453 : std::vector<std::string> surfaceSpeciesNames() const;
454 :
455 : protected:
456 : /**
457 : * After parsing the database file, remove any secondary species that have extrapolated
458 : * equilibrium constants. This is called in the constructor if the
459 : * remove_all_extrapolated_secondary_species flag is true
460 : */
461 : void removeExtrapolatedSecondarySpecies();
462 :
463 : /**
464 : * Copy the temperature points (if any) found in the database into _temperature_points. This
465 : * method is called in the constructor
466 : */
467 : void setTemperatures();
468 :
469 : /**
470 : * Copy the Debye-Huckel parameters (if any) found in the database into _debye_huckel. This
471 : * method is called in the constructor
472 : */
473 : void setDebyeHuckel();
474 :
475 : /**
476 : * Copy the Debye-Huckel parameters for computing neutral species activity (if any) found in the
477 : * database into _neutral_species_activity. This method is called in the constructor
478 : */
479 : void setNeutralSpeciesActivity();
480 :
481 : /// Database filename
482 : const FileName _filename;
483 : /// JSON data
484 : nlohmann::json _root;
485 : /// List of basis (primary) species names read from database
486 : std::vector<std::string> _bs_names;
487 : /// List of secondary equilibrium species to read from database
488 : std::vector<std::string> _es_names;
489 : /// List of secondary mineral species to read from database
490 : std::vector<std::string> _ms_names;
491 : /// Temperature points in database
492 : std::vector<Real> _temperature_points;
493 : /// Pressure points in database
494 : std::vector<Real> _pressure_points;
495 : /// Elements and their molecular weight read from the database
496 : std::map<std::string, GeochemistryElements> _elements;
497 : /// Basis species data read from the database
498 : std::map<std::string, GeochemistryBasisSpecies> _basis_species;
499 : /// Secondary equilibrium species and free electron data read from the database
500 : std::map<std::string, GeochemistryEquilibriumSpecies> _equilibrium_species;
501 : /// Mineral species data read from the database
502 : std::map<std::string, GeochemistryMineralSpecies> _mineral_species;
503 : /// Gas species data read from the database
504 : std::map<std::string, GeochemistryGasSpecies> _gas_species;
505 : /// Redox species (couples) data read from the database
506 : std::map<std::string, GeochemistryRedoxSpecies> _redox_species;
507 : /// Oxide species data read from the database
508 : std::map<std::string, GeochemistryOxideSpecies> _oxide_species;
509 : /// Surface sorbing species data read from the database
510 : std::map<std::string, GeochemistrySurfaceSpecies> _surface_species;
511 : /// Debye-Huckel activity coefficients
512 : GeochemistryDebyeHuckel _debye_huckel;
513 : /// Neutral species activity coefficients
514 : std::map<std::string, GeochemistryNeutralSpeciesActivity> _neutral_species_activity;
515 : // Helper for converting json node to Real from string
516 : static Real getReal(const nlohmann::json & node);
517 :
518 : private:
519 : /**
520 : * Generates a formatted vector of strings representing all reactions
521 : * @param names list of reaction species
522 : * @param basis species list of basis species for each reaction species (this vector must be of
523 : * same size as names)
524 : * @return formatted reaction equations
525 : */
526 : std::vector<std::string>
527 : printReactions(const std::vector<std::string> & names,
528 : const std::vector<std::map<std::string, Real>> & basis_species) const;
529 : };
|
