Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "PertinentGeochemicalSystem.h"
11 :
12 1178 : PertinentGeochemicalSystem::PertinentGeochemicalSystem(
13 : const GeochemicalDatabaseReader & db,
14 : const std::vector<std::string> & basis_species,
15 : const std::vector<std::string> & minerals,
16 : const std::vector<std::string> & gases,
17 : const std::vector<std::string> & kinetic_minerals,
18 : const std::vector<std::string> & kinetic_redox,
19 : const std::vector<std::string> & kinetic_surface_species,
20 : const std::string & redox_ox,
21 1178 : const std::string & redox_e)
22 1178 : : _db(db),
23 1178 : _basis_index(),
24 1178 : _basis_info(),
25 1178 : _mineral_index(),
26 1178 : _mineral_info(),
27 1178 : _gas_index(),
28 1178 : _gas_info(),
29 1178 : _kinetic_mineral_index(),
30 1178 : _kinetic_mineral_info(),
31 1178 : _kinetic_redox_index(),
32 1178 : _kinetic_redox_info(),
33 1178 : _kinetic_surface_index(),
34 1178 : _kinetic_surface_info(),
35 1178 : _secondary_index(),
36 1178 : _secondary_info(),
37 1178 : _redox_ox(redox_ox),
38 1178 : _redox_e(redox_e),
39 1178 : _model(_db)
40 : {
41 : // Use the constructor info to build the "index" and "info" structures
42 1178 : buildBasis(basis_species);
43 1175 : buildMinerals(minerals);
44 1173 : buildGases(gases);
45 1171 : buildKineticMinerals(kinetic_minerals);
46 1168 : buildKineticRedox(kinetic_redox);
47 1165 : buildKineticSurface(kinetic_surface_species);
48 :
49 : // Pull out all secondary equilibrium species: these are all relevant "redox couples" and
50 : // "secondary species" and "surface species"
51 1163 : buildSecondarySpecies();
52 :
53 : // Build minerals in the case that minerals = {"*"}
54 1163 : buildAllMinerals(minerals);
55 :
56 : // Check that everything can be expressed in terms of the basis, possibly via redox and secondary
57 : // species
58 1163 : checkMinerals(_mineral_info);
59 1160 : checkGases();
60 1159 : checkMinerals(_kinetic_mineral_info);
61 1157 : checkKineticRedox();
62 1156 : checkKineticSurfaceSpecies();
63 :
64 : // Populate _model
65 1155 : createModel();
66 1378 : }
67 :
68 : unsigned
69 196 : PertinentGeochemicalSystem::getIndexOfOriginalBasisSpecies(const std::string & name) const
70 : {
71 : try
72 : {
73 195 : return _basis_index.at(name);
74 : }
75 1 : catch (const std::out_of_range &)
76 : {
77 1 : mooseError("species ", name, " is not in the original basis");
78 1 : }
79 0 : catch (...)
80 : {
81 0 : throw;
82 0 : }
83 : }
84 :
85 : std::vector<std::string>
86 318 : PertinentGeochemicalSystem::originalBasisNames() const
87 : {
88 318 : std::vector<std::string> names(_basis_info.size());
89 2091 : for (const auto & name_ind : _basis_index)
90 1773 : names[name_ind.second] = name_ind.first;
91 318 : return names;
92 0 : }
93 :
94 : void
95 1178 : PertinentGeochemicalSystem::buildBasis(const std::vector<std::string> & basis_species)
96 : {
97 1178 : unsigned ind = 0;
98 8148 : for (const auto & name : basis_species)
99 : {
100 8151 : if (ind == 0 and name != "H2O")
101 1 : mooseError("First member of basis species list must be H2O");
102 : if (_basis_index.count(name) == 1)
103 1 : mooseError(name, " exists more than once in the basis species list");
104 : _basis_index.emplace(name, ind);
105 6971 : ind += 1;
106 6971 : if (_db.isBasisSpecies(name))
107 18828 : _basis_info.push_back(_db.getBasisSpecies({name})[name]);
108 695 : else if (_db.isRedoxSpecies(name))
109 : {
110 2082 : const GeochemistryRedoxSpecies rs = _db.getRedoxSpecies({name})[name];
111 : GeochemistryBasisSpecies bs;
112 : bs.name = rs.name;
113 694 : bs.radius = rs.radius;
114 694 : bs.charge = rs.charge;
115 694 : bs.molecular_weight = rs.molecular_weight;
116 694 : _basis_info.push_back(bs);
117 694 : }
118 : else
119 1 : mooseError(name + " does not exist in the basis species or redox species in " +
120 1 : _db.filename());
121 : }
122 8145 : }
123 :
124 : void
125 1175 : PertinentGeochemicalSystem::buildMinerals(const std::vector<std::string> & minerals)
126 : {
127 1175 : unsigned ind = 0;
128 1591 : if (minerals.size() == 1 && minerals[0] == "*")
129 1 : return; // buildAllMinerals is called later
130 2817 : for (const auto & name : minerals)
131 : {
132 : if (_mineral_index.count(name) == 1)
133 1 : mooseError(name + " exists more than once in the minerals list");
134 : _mineral_index.emplace(name, ind);
135 1644 : ind += 1;
136 4932 : _mineral_info.push_back(_db.getMineralSpecies({name})[name]);
137 : }
138 1644 : }
139 :
140 : void
141 1173 : PertinentGeochemicalSystem::buildGases(const std::vector<std::string> & gases)
142 : {
143 1173 : unsigned ind = 0;
144 1514 : for (const auto & name : gases)
145 : {
146 : if (_gas_index.count(name) == 1)
147 1 : mooseError(name + " exists more than once in the gases list");
148 : _gas_index.emplace(name, ind);
149 342 : ind += 1;
150 1026 : const GeochemistryGasSpecies gas = _db.getGasSpecies({name})[name];
151 341 : _gas_info.push_back(gas);
152 341 : }
153 1513 : }
154 :
155 : void
156 1171 : PertinentGeochemicalSystem::buildKineticMinerals(const std::vector<std::string> & kinetic_minerals)
157 : {
158 1171 : unsigned ind = 0;
159 1423 : for (const auto & name : kinetic_minerals)
160 : {
161 : if (_kinetic_mineral_index.count(name) == 1)
162 1 : mooseError(name + " exists more than once in the kinetic_minerals list");
163 : if (_mineral_index.count(name) == 1)
164 1 : mooseError(name + " exists in both the minerals and kinetic_minerals lists");
165 : _kinetic_mineral_index.emplace(name, ind);
166 253 : ind += 1;
167 759 : _kinetic_mineral_info.push_back(_db.getMineralSpecies({name})[name]);
168 : }
169 1421 : }
170 :
171 : void
172 1168 : PertinentGeochemicalSystem::buildKineticRedox(const std::vector<std::string> & kinetic_redox)
173 : {
174 1168 : unsigned ind = 0;
175 1218 : for (const auto & name : kinetic_redox)
176 : {
177 : if (_kinetic_redox_index.count(name) == 1)
178 1 : mooseError(name + " exists more than once in the kinetic_redox list");
179 : if (_basis_index.count(name) == 1)
180 1 : mooseError(name + " exists in both the basis_species and kinetic_redox lists");
181 : _kinetic_redox_index.emplace(name, ind);
182 51 : ind += 1;
183 153 : _kinetic_redox_info.push_back(_db.getRedoxSpecies({name})[name]);
184 : }
185 1216 : }
186 :
187 : void
188 1165 : PertinentGeochemicalSystem::buildKineticSurface(
189 : const std::vector<std::string> & kinetic_surface_species)
190 : {
191 1165 : unsigned ind = 0;
192 1188 : for (const auto & name : kinetic_surface_species)
193 : {
194 : if (_kinetic_surface_index.count(name) == 1)
195 1 : mooseError(name + " exists more than once in the kinetic_surface_species list");
196 : _kinetic_surface_index.emplace(name, ind);
197 24 : ind += 1;
198 72 : _kinetic_surface_info.push_back(_db.getSurfaceSpecies({name})[name]);
199 : }
200 1187 : }
201 :
202 : void
203 1163 : PertinentGeochemicalSystem::buildSecondarySpecies()
204 : {
205 : // run through all redox couples, including them if:
206 : // - they are not part of the kinetic_redox list
207 : // - they are not part of the basis_species list
208 : // - their reaction involves only basis_species or secondary species already encountered
209 1163 : unsigned ind = 0;
210 30362 : for (const auto & name : _db.redoxCoupleNames())
211 : {
212 : if (_kinetic_redox_index.count(name) == 0 && _basis_index.count(name) == 0)
213 : {
214 85401 : const GeochemistryRedoxSpecies rs = _db.getRedoxSpecies({name})[name];
215 : // check all reaction species are in the basis
216 : bool all_species_in_basis_or_sec = true;
217 67341 : for (const auto & element : rs.basis_species)
218 66041 : if (_basis_index.count(element.first) == 0 && _secondary_index.count(element.first) == 0)
219 : {
220 : all_species_in_basis_or_sec = false;
221 : break;
222 : }
223 28467 : if (all_species_in_basis_or_sec)
224 : {
225 : _secondary_index.emplace(name, ind);
226 1300 : ind += 1;
227 : GeochemistryEquilibriumSpecies ss;
228 : ss.name = rs.name;
229 : ss.basis_species = rs.basis_species;
230 1300 : ss.equilibrium_const = rs.equilibrium_const;
231 1300 : ss.radius = rs.radius;
232 1300 : ss.charge = rs.charge;
233 1300 : ss.molecular_weight = rs.molecular_weight;
234 1300 : _secondary_info.push_back(ss);
235 1300 : }
236 28467 : }
237 1163 : }
238 :
239 : // run through all secondary species, including them if:
240 : // - their reaction involves only basis_species, or secondary species already encountered
241 : // - the name is not _redox_e (which is usually "e-" the free electron)
242 300486 : for (const auto & name : _db.secondarySpeciesNames())
243 : {
244 299323 : if (name == _redox_e)
245 989 : continue;
246 895002 : const GeochemistryEquilibriumSpecies ss = _db.getEquilibriumSpecies({name})[name];
247 : // check all reaction species are in the basis
248 : bool all_species_in_basis_or_sec = true;
249 516298 : for (const auto & element : ss.basis_species)
250 499334 : if (_basis_index.count(element.first) == 0 && _secondary_index.count(element.first) == 0)
251 : {
252 : all_species_in_basis_or_sec = false;
253 : break;
254 : }
255 298334 : if (all_species_in_basis_or_sec)
256 : {
257 : _secondary_index.emplace(name, ind);
258 16964 : ind += 1;
259 16964 : _secondary_info.push_back(ss);
260 : }
261 299497 : }
262 :
263 : // run through all surface species, including them if:
264 : // - their reaction involves only basis_species or secondary species encountered so far
265 : // - they are not in the kinetic_surface_species list
266 12083 : for (const auto & name : _db.surfaceSpeciesNames())
267 : {
268 : if (_kinetic_surface_index.count(name) == 0)
269 : {
270 32703 : const GeochemistrySurfaceSpecies ss = _db.getSurfaceSpecies({name})[name];
271 : // check all reaction species are in the basis
272 : bool all_species_in_basis_or_sec = true;
273 19125 : for (const auto & element : ss.basis_species)
274 18687 : if (_basis_index.count(element.first) == 0 && _secondary_index.count(element.first) == 0)
275 : {
276 : all_species_in_basis_or_sec = false;
277 : break;
278 : }
279 10901 : if (all_species_in_basis_or_sec)
280 : {
281 : GeochemistryEquilibriumSpecies to_sec;
282 : to_sec.name = ss.name;
283 : to_sec.basis_species = ss.basis_species;
284 438 : to_sec.radius = -1.5; // flag to activity calculators that activity coefficient = 1
285 438 : to_sec.charge = ss.charge;
286 438 : to_sec.molecular_weight = ss.molecular_weight;
287 438 : const Real T0 = _db.getTemperatures()[0];
288 3942 : for (const auto & temp : _db.getTemperatures())
289 3504 : to_sec.equilibrium_const.push_back(ss.log10K + ss.dlog10KdT * (temp - T0));
290 :
291 : _secondary_index.emplace(name, ind);
292 438 : ind += 1;
293 438 : _secondary_info.push_back(to_sec);
294 438 : }
295 10901 : }
296 1163 : }
297 338865 : }
298 :
299 : void
300 1163 : PertinentGeochemicalSystem::buildAllMinerals(const std::vector<std::string> & minerals)
301 : {
302 1575 : if (!(minerals.size() == 1 && minerals[0] == "*"))
303 1162 : return; // buildMinerals has done its job of building _mineral_info and _mineral_index
304 1 : unsigned ind = 0;
305 7 : for (const auto & name_ms : _db.getMineralSpecies(_db.mineralSpeciesNames()))
306 : {
307 6 : if (_kinetic_mineral_index.count(name_ms.first) == 1)
308 1 : continue;
309 : bool known_basis_only = true;
310 15 : for (const auto & basis_stoi : name_ms.second.basis_species)
311 : {
312 12 : if (_basis_index.count(basis_stoi.first) == 0 &&
313 : _secondary_index.count(basis_stoi.first) == 0)
314 : {
315 : known_basis_only = false;
316 : break;
317 : }
318 : }
319 5 : if (known_basis_only)
320 : {
321 : _mineral_index.emplace(name_ms.first, ind);
322 3 : ind += 1;
323 3 : _mineral_info.push_back(name_ms.second);
324 : }
325 : }
326 : }
327 :
328 : bool
329 335 : PertinentGeochemicalSystem::checkRedoxe()
330 : {
331 : bool found = false;
332 89496 : for (const auto & name : _db.secondarySpeciesNames())
333 89162 : if (name == _redox_e)
334 : {
335 : found = true;
336 1005 : const GeochemistryEquilibriumSpecies ss = _db.getEquilibriumSpecies({name})[name];
337 1339 : for (const auto & element : ss.basis_species)
338 1005 : if (_basis_index.count(element.first) == 0 && _secondary_index.count(element.first) == 0)
339 : return false;
340 669 : }
341 334 : return found;
342 335 : }
343 :
344 : void
345 334 : PertinentGeochemicalSystem::buildRedoxeInfo(std::vector<Real> & redox_e_stoichiometry,
346 : std::vector<Real> & redox_e_log10K)
347 : {
348 334 : const unsigned num_basis = _basis_info.size();
349 334 : const unsigned numT = _db.getTemperatures().size();
350 334 : redox_e_stoichiometry.assign(num_basis, 0.0);
351 334 : redox_e_log10K.assign(numT, 0.0);
352 89490 : for (const auto & name : _db.secondarySpeciesNames())
353 89156 : if (name == _redox_e)
354 : {
355 1002 : const GeochemistryEquilibriumSpecies ss = _db.getEquilibriumSpecies({name})[name];
356 3006 : for (unsigned i = 0; i < numT; ++i)
357 2672 : redox_e_log10K[i] = ss.equilibrium_const[i];
358 1336 : for (const auto & react : ss.basis_species)
359 : {
360 1002 : const Real stoi_coeff = react.second;
361 1002 : if (_model.basis_species_index.count(react.first) == 1)
362 : {
363 1000 : const unsigned col = _model.basis_species_index[react.first];
364 1000 : redox_e_stoichiometry[col] += react.second;
365 : }
366 : else if (_secondary_index.count(react.first) == 1)
367 : {
368 : // reaction species is not a basis component, but a secondary component.
369 : // So express stoichiometry in terms of the secondary component's reaction
370 2 : const unsigned sec_row = _model.eqm_species_index[react.first];
371 18 : for (unsigned i = 0; i < numT; ++i)
372 16 : redox_e_log10K[i] += stoi_coeff * _model.eqm_log10K(sec_row, i);
373 12 : for (unsigned col = 0; col < num_basis; ++col)
374 10 : redox_e_stoichiometry[col] += stoi_coeff * _model.eqm_stoichiometry(sec_row, col);
375 : }
376 : }
377 668 : }
378 668 : }
379 :
380 : void
381 2322 : PertinentGeochemicalSystem::checkMinerals(
382 : const std::vector<GeochemistryMineralSpecies> & mineral_info) const
383 : {
384 4203 : for (const auto & mineral : mineral_info)
385 : {
386 8312 : for (const auto & element : mineral.basis_species)
387 6429 : if (_basis_index.count(element.first) == 0 && _secondary_index.count(element.first) == 0)
388 9 : mooseError("The reaction for " + mineral.name + " depends on " + element.first +
389 : " which is not reducable to a set of basis species");
390 2101 : for (const auto & element : mineral.sorption_sites)
391 220 : if (_basis_index.count(element.first) == 0)
392 6 : mooseError("The sorbing sites for " + mineral.name + " include " + element.first +
393 : " which is not in the basis_species list");
394 : }
395 2317 : }
396 :
397 : void
398 1160 : PertinentGeochemicalSystem::checkGases() const
399 : {
400 1498 : for (const auto & gas : _gas_info)
401 1037 : for (const auto & element : gas.basis_species)
402 699 : if (_basis_index.count(element.first) == 0 && _secondary_index.count(element.first) == 0)
403 3 : mooseError("The reaction for " + gas.name + " depends on " + element.first +
404 : " which is not reducable to a set of basis species");
405 1159 : }
406 :
407 : void
408 1157 : PertinentGeochemicalSystem::checkKineticRedox() const
409 : {
410 1202 : for (const auto & kr : _kinetic_redox_info)
411 234 : for (const auto & element : kr.basis_species)
412 189 : if (_basis_index.count(element.first) == 0 && _secondary_index.count(element.first) == 0)
413 3 : mooseError("The reaction for " + kr.name + " depends on " + element.first +
414 : " which is not reducable to a set of basis species");
415 1156 : }
416 :
417 : void
418 1156 : PertinentGeochemicalSystem::checkKineticSurfaceSpecies() const
419 : {
420 1174 : for (const auto & kr : _kinetic_surface_info)
421 63 : for (const auto & element : kr.basis_species)
422 45 : if (_basis_index.count(element.first) == 0 && _secondary_index.count(element.first) == 0)
423 3 : mooseError("The reaction for " + kr.name + " depends on " + element.first +
424 : " which is not reducable to a set of basis species");
425 1155 : }
426 :
427 : void
428 1155 : PertinentGeochemicalSystem::createModel()
429 : {
430 1155 : const unsigned num_rows = _secondary_info.size() + _mineral_info.size() + _gas_info.size();
431 1155 : const unsigned num_cols = _basis_info.size();
432 1155 : const unsigned num_temperatures = _db.getTemperatures().size();
433 : unsigned ind = 0;
434 :
435 : // create basis_species_index map
436 : _model.basis_species_index = _basis_index;
437 :
438 : // extract the names
439 1155 : _model.basis_species_name = std::vector<std::string>(num_cols);
440 8053 : for (const auto & species : _basis_info)
441 6898 : _model.basis_species_name[_model.basis_species_index[species.name]] = species.name;
442 :
443 : // initially no basis species are minerals
444 1155 : _model.basis_species_mineral = std::vector<bool>(num_cols, false);
445 :
446 : // initially no basis species are gases
447 1155 : _model.basis_species_gas = std::vector<bool>(num_cols, false);
448 :
449 : // initially all basis species are involved in transport (this gets modified for surface species
450 : // below)
451 1155 : _model.basis_species_transported = std::vector<bool>(num_cols, true);
452 :
453 : // record the charge
454 2310 : _model.basis_species_charge = std::vector<Real>(num_cols, 0.0);
455 8053 : for (const auto & species : _basis_info)
456 6898 : _model.basis_species_charge[_model.basis_species_index[species.name]] = species.charge;
457 :
458 : // record the ionic radius
459 2310 : _model.basis_species_radius = std::vector<Real>(num_cols, 0.0);
460 8053 : for (const auto & species : _basis_info)
461 6898 : _model.basis_species_radius[_model.basis_species_index[species.name]] = species.radius;
462 :
463 : // record the molecular weight
464 2310 : _model.basis_species_molecular_weight = std::vector<Real>(num_cols, 0.0);
465 8053 : for (const auto & species : _basis_info)
466 6898 : _model.basis_species_molecular_weight[_model.basis_species_index[species.name]] =
467 6898 : species.molecular_weight;
468 :
469 : // record the molecular weight (zero for all species except minerals)
470 1155 : _model.basis_species_molecular_volume = std::vector<Real>(num_cols, 0.0);
471 :
472 : // the "eqm_species" stuff is rather long-winded because of the different data structures used
473 : // to hold secondary, mineral and gas info. There is a bit of an overlap, however, so let's
474 : // create a new data structure that contains that overlapping info
475 :
476 1155 : std::vector<GeochemistryEquilibriumSpecies> overlap(_secondary_info);
477 2789 : for (const auto & species : _mineral_info)
478 : {
479 : GeochemistryEquilibriumSpecies es;
480 1634 : es.name = species.name;
481 1634 : es.molecular_weight = species.molecular_weight;
482 1634 : es.equilibrium_const = species.equilibrium_const;
483 : es.basis_species = species.basis_species;
484 1634 : overlap.push_back(es);
485 1634 : }
486 1493 : for (const auto & species : _gas_info)
487 : {
488 : GeochemistryEquilibriumSpecies es;
489 338 : es.name = species.name;
490 338 : es.molecular_weight = species.molecular_weight;
491 338 : es.equilibrium_const = species.equilibrium_const;
492 : es.basis_species = species.basis_species;
493 338 : overlap.push_back(es);
494 338 : }
495 :
496 : // create the eqm_species_index map
497 : ind = 0;
498 21809 : for (const auto & species : overlap)
499 20654 : _model.eqm_species_index[species.name] = ind++;
500 :
501 : // extract the names
502 1155 : _model.eqm_species_name = std::vector<std::string>(num_rows);
503 21809 : for (const auto & species : _model.eqm_species_index)
504 20654 : _model.eqm_species_name[species.second] = species.first;
505 :
506 : // create the eqm_species_mineral vector
507 2310 : _model.eqm_species_mineral = std::vector<bool>(num_rows, false);
508 2789 : for (const auto & species : _mineral_info)
509 1634 : _model.eqm_species_mineral[_model.eqm_species_index[species.name]] = true;
510 1493 : for (const auto & species : _gas_info)
511 338 : _model.eqm_species_mineral[_model.eqm_species_index[species.name]] = false;
512 :
513 : // create the eqm_species_gas vector
514 2310 : _model.eqm_species_gas = std::vector<bool>(num_rows, false);
515 2789 : for (const auto & species : _mineral_info)
516 1634 : _model.eqm_species_gas[_model.eqm_species_index[species.name]] = false;
517 1493 : for (const auto & species : _gas_info)
518 338 : _model.eqm_species_gas[_model.eqm_species_index[species.name]] = true;
519 :
520 : // create the eqm_species_transported vector (true for non-minerals) - gets modified below for
521 : // surface species
522 2310 : _model.eqm_species_transported = std::vector<bool>(num_rows, true);
523 2789 : for (const auto & species : _mineral_info)
524 1634 : _model.eqm_species_transported[_model.eqm_species_index.at(species.name)] = false;
525 :
526 : // record the charge
527 : _model.eqm_species_charge =
528 2310 : std::vector<Real>(num_rows, 0.0); // charge of gases and minerals is zero
529 19837 : for (const auto & species : _secondary_info)
530 18682 : _model.eqm_species_charge[_model.eqm_species_index[species.name]] = species.charge;
531 :
532 : // record the radius
533 : _model.eqm_species_radius =
534 2310 : std::vector<Real>(num_rows, 0.0); // ionic radius of gases and minerals is zero
535 19837 : for (const auto & species : _secondary_info)
536 18682 : _model.eqm_species_radius[_model.eqm_species_index[species.name]] = species.radius;
537 :
538 : // record the molecular weight
539 1155 : _model.eqm_species_molecular_weight = std::vector<Real>(num_rows, 0.0);
540 21809 : for (const auto & species : overlap)
541 20654 : _model.eqm_species_molecular_weight[_model.eqm_species_index[species.name]] =
542 20654 : species.molecular_weight;
543 :
544 : // record the molecular volume (zero for all species except minerals)
545 2310 : _model.eqm_species_molecular_volume = std::vector<Real>(num_rows, 0.0);
546 2789 : for (const auto & species : _mineral_info)
547 1634 : _model.eqm_species_molecular_volume[_model.eqm_species_index[species.name]] =
548 1634 : species.molecular_volume;
549 :
550 : // record surface-complexation info
551 2789 : for (const auto & species : _mineral_info)
552 1634 : if (species.surface_area != 0.0)
553 : {
554 : SurfaceComplexationInfo sci;
555 60 : sci.surface_area = species.surface_area;
556 : sci.sorption_sites = species.sorption_sites;
557 60 : _model.surface_complexation_info[species.name] = sci;
558 : }
559 1401 : for (const auto & species : _kinetic_mineral_info)
560 246 : if (species.surface_area != 0.0)
561 : {
562 : SurfaceComplexationInfo sci;
563 48 : sci.surface_area = species.surface_area;
564 : sci.sorption_sites = species.sorption_sites;
565 48 : _model.surface_complexation_info[species.name] = sci;
566 : }
567 :
568 : // record gas fugacity info
569 1493 : for (const auto & species : _gas_info)
570 338 : _model.gas_chi[species.name] = species.chi;
571 :
572 : // create the stoichiometry matrix
573 1155 : _model.eqm_stoichiometry.resize(num_rows, num_cols);
574 1155 : _model.eqm_log10K.resize(num_rows, num_temperatures);
575 :
576 : // populate the stoichiometry
577 21809 : for (const auto & species : overlap)
578 : {
579 20654 : const unsigned row = _model.eqm_species_index[species.name];
580 185886 : for (unsigned i = 0; i < num_temperatures; ++i)
581 165232 : _model.eqm_log10K(row, i) = species.equilibrium_const[i];
582 75315 : for (const auto & react : species.basis_species)
583 : {
584 54661 : const Real stoi_coeff = react.second;
585 54661 : if (_model.basis_species_index.count(react.first) == 1)
586 : {
587 49956 : const unsigned col = _model.basis_species_index[react.first];
588 49956 : _model.eqm_stoichiometry(row, col) += react.second;
589 : }
590 : else if (_secondary_index.count(react.first) == 1)
591 : {
592 : // reaction species is not a basis component, but a secondary component.
593 : // So express stoichiometry in terms of the secondary component's reaction
594 4705 : const unsigned sec_row = _model.eqm_species_index[react.first];
595 42345 : for (unsigned i = 0; i < num_temperatures; ++i)
596 37640 : _model.eqm_log10K(row, i) += stoi_coeff * _model.eqm_log10K(sec_row, i);
597 65525 : for (unsigned col = 0; col < num_cols; ++col)
598 60820 : _model.eqm_stoichiometry(row, col) += stoi_coeff * _model.eqm_stoichiometry(sec_row, col);
599 : }
600 : else
601 0 : mooseError("Species " + species.name + " includes " + react.first +
602 : ", which cannot be expressed in terms of the basis. Previous checks must be "
603 : "erroneous!");
604 : }
605 : }
606 :
607 : // Build the redox information, if any. Here we express any O2(aq) in the redox equations in
608 : // terms of redox_e (which is usually e-)
609 1155 : _model.redox_lhs = _redox_e;
610 1155 : std::vector<Real> redox_e_stoichiometry(num_cols, 0.0);
611 1157 : std::vector<Real> redox_e_log10K(num_temperatures, 0.0);
612 : std::vector<Real> redox_stoi;
613 : std::vector<Real> redox_log10K;
614 1490 : if ((_model.basis_species_index.count(_redox_ox) == 1) && checkRedoxe())
615 : {
616 : // construct the stoichiometry and log10K for _redox_e and put it
617 334 : buildRedoxeInfo(redox_e_stoichiometry, redox_e_log10K);
618 334 : redox_stoi.insert(redox_stoi.end(), redox_e_stoichiometry.begin(), redox_e_stoichiometry.end());
619 334 : redox_log10K.insert(redox_log10K.end(), redox_e_log10K.begin(), redox_e_log10K.end());
620 : // the electron reaction is
621 : // e- = nuw_i * basis_i + beta * O2(aq), where we've pulled out the O2(aq) because it's
622 : // special
623 334 : const unsigned o2_index = _model.basis_species_index.at(_redox_ox);
624 334 : const Real beta = redox_e_stoichiometry[o2_index];
625 334 : if (beta != 0.0)
626 : {
627 3422 : for (const auto & bs : _model.basis_species_index)
628 3088 : if (_db.isRedoxSpecies(bs.first))
629 : {
630 : // this basis species is a redox couple in disequilibrium
631 1550 : const GeochemistryRedoxSpecies rs = _db.getRedoxSpecies({bs.first})[bs.first];
632 : // check that its reaction involves only basis species, and record the stoichiometry in
633 : // the current basis
634 516 : std::vector<Real> stoi(num_cols, 0.0);
635 : bool only_involves_basis_species = true;
636 1858 : for (const auto & name_stoi : rs.basis_species)
637 : {
638 1563 : if (_model.basis_species_index.count(name_stoi.first) == 1)
639 1342 : stoi[_model.basis_species_index.at(name_stoi.first)] = name_stoi.second;
640 : else
641 : {
642 : only_involves_basis_species = false;
643 : break;
644 : }
645 : }
646 516 : if (!only_involves_basis_species)
647 221 : continue;
648 : // Reaction is now
649 : // redox = nu_i * basis_i + alpha * O2(aq), where we've pulled the O2(aq) out because
650 : // it's special. Now pull the redox couple to the RHS of the reaction, so we have 0 =
651 : // -redox + nu_i * basis_i + alpha * O2(aq)
652 295 : stoi[bs.second] = -1.0;
653 : // check that the stoichiometry involves O2(aq)
654 295 : const Real alpha = stoi[o2_index];
655 295 : if (alpha == 0.0)
656 1 : continue;
657 : // multiply equation -beta/alpha so it reads
658 : // 0 = -beta/alpha * (-redox + nu_i * basis_i) - beta * O2(aq)
659 3280 : for (unsigned basis_i = 0; basis_i < num_cols; ++basis_i)
660 2986 : stoi[basis_i] *= -beta / alpha;
661 : // add the equation to e- = nuw_i * basis_i + beta * O2(aq)
662 3280 : for (unsigned basis_i = 0; basis_i < num_cols; ++basis_i)
663 2986 : stoi[basis_i] += redox_e_stoichiometry[basis_i];
664 : // now the reation is e- = nuw_i * basis_i - beta/alpha * (-redox + nu_i * basis_i)
665 294 : redox_stoi.insert(redox_stoi.end(), stoi.begin(), stoi.end());
666 :
667 : // record the equilibrium constants
668 2646 : for (unsigned temp = 0; temp < num_temperatures; ++temp)
669 2352 : redox_log10K.push_back((-beta / alpha) * rs.equilibrium_const[temp] +
670 : redox_e_log10K[temp]);
671 516 : }
672 : }
673 : }
674 : // record the above in the model.redox_stoichiometry and model.redox_log10K DenseMatrices
675 1155 : const unsigned num_redox = redox_stoi.size() / num_cols;
676 1155 : _model.redox_stoichiometry.resize(num_redox, num_cols);
677 1783 : for (unsigned red = 0; red < num_redox; ++red)
678 6702 : for (unsigned basis_i = 0; basis_i < num_cols; ++basis_i)
679 6074 : _model.redox_stoichiometry(red, basis_i) = redox_stoi[red * num_cols + basis_i];
680 1155 : _model.redox_log10K.resize(num_redox, num_temperatures);
681 1783 : for (unsigned red = 0; red < num_redox; ++red)
682 5652 : for (unsigned temp = 0; temp < num_temperatures; ++temp)
683 5024 : _model.redox_log10K(red, temp) = redox_log10K[red * num_temperatures + temp];
684 :
685 : // To build the kin_species_index, kin_species_name, etc, let's build an "overlap", similar to
686 : // above
687 1155 : std::vector<GeochemistryMineralSpecies> overlap_kin(_kinetic_mineral_info);
688 1200 : for (const auto & species : _kinetic_redox_info)
689 : {
690 : GeochemistryMineralSpecies ms;
691 45 : ms.name = species.name;
692 45 : ms.molecular_volume = 0.0;
693 : ms.basis_species = species.basis_species;
694 45 : ms.molecular_weight = species.molecular_weight;
695 45 : ms.equilibrium_const = species.equilibrium_const;
696 45 : overlap_kin.push_back(ms);
697 45 : }
698 1173 : for (const auto & species : _kinetic_surface_info)
699 : {
700 : GeochemistryMineralSpecies ms;
701 18 : ms.name = species.name;
702 18 : ms.molecular_volume = 0.0;
703 : ms.basis_species = species.basis_species;
704 18 : ms.molecular_weight = species.molecular_weight;
705 18 : const Real T0 = _db.getTemperatures()[0];
706 162 : for (const auto & temp : _db.getTemperatures())
707 144 : ms.equilibrium_const.push_back(species.log10K + species.dlog10KdT * (temp - T0));
708 18 : overlap_kin.push_back(ms);
709 18 : }
710 1155 : const unsigned num_kin = overlap_kin.size();
711 :
712 : // create the kin_species_index map
713 : ind = 0;
714 1464 : for (const auto & species : overlap_kin)
715 309 : _model.kin_species_index[species.name] = ind++;
716 :
717 : // extract the names
718 1155 : _model.kin_species_name = std::vector<std::string>(num_kin);
719 1464 : for (const auto & species : _model.kin_species_index)
720 309 : _model.kin_species_name[species.second] = species.first;
721 :
722 : // build the kin_species_mineral info
723 2310 : _model.kin_species_mineral = std::vector<bool>(num_kin, true);
724 1200 : for (const auto & species : _kinetic_redox_info)
725 45 : _model.kin_species_mineral[_model.kin_species_index[species.name]] = false;
726 1173 : for (const auto & species : _kinetic_surface_info)
727 18 : _model.kin_species_mineral[_model.kin_species_index[species.name]] = false;
728 :
729 : // create the kin_species_transported vector (false for minerals and surface species)
730 2310 : _model.kin_species_transported = std::vector<bool>(num_kin, true);
731 1401 : for (const auto & species : _kinetic_mineral_info)
732 246 : _model.kin_species_transported[_model.kin_species_index.at(species.name)] = false;
733 1173 : for (const auto & species : _kinetic_surface_info)
734 18 : _model.kin_species_transported[_model.kin_species_index.at(species.name)] = false;
735 :
736 : // build the kin_species_charge info
737 2310 : _model.kin_species_charge = std::vector<Real>(num_kin, 0.0);
738 1200 : for (const auto & species : _kinetic_redox_info)
739 45 : _model.kin_species_charge[_model.kin_species_index[species.name]] = species.charge;
740 1173 : for (const auto & species : _kinetic_surface_info)
741 18 : _model.kin_species_charge[_model.kin_species_index[species.name]] = species.charge;
742 :
743 : // extract the molecular weight
744 1155 : _model.kin_species_molecular_weight = std::vector<Real>(num_kin, 0.0);
745 1464 : for (const auto & species : overlap_kin)
746 309 : _model.kin_species_molecular_weight[_model.kin_species_index[species.name]] =
747 309 : species.molecular_weight;
748 :
749 : // extract the molecular volume
750 1155 : _model.kin_species_molecular_volume = std::vector<Real>(num_kin, 0.0);
751 1464 : for (const auto & species : overlap_kin)
752 309 : _model.kin_species_molecular_volume[_model.kin_species_index[species.name]] =
753 309 : species.molecular_volume;
754 :
755 : // extract the stoichiometry
756 1155 : _model.kin_stoichiometry.resize(num_kin, num_cols);
757 1155 : _model.kin_log10K.resize(num_kin, num_temperatures);
758 :
759 : // populate the stoichiometry
760 1464 : for (const auto & species : overlap_kin)
761 : {
762 309 : const unsigned row = _model.kin_species_index[species.name];
763 2781 : for (unsigned i = 0; i < num_temperatures; ++i)
764 2472 : _model.kin_log10K(row, i) = species.equilibrium_const[i];
765 1267 : for (const auto & react : species.basis_species)
766 : {
767 958 : const Real stoi_coeff = react.second;
768 958 : if (_model.basis_species_index.count(react.first) == 1)
769 : {
770 922 : const unsigned col = _model.basis_species_index[react.first];
771 922 : _model.kin_stoichiometry(row, col) += react.second;
772 : }
773 : else if (_model.eqm_species_index.count(react.first) == 1)
774 : {
775 : // reaction species is not a basis component, but a secondary component.
776 : // So express stoichiometry in terms of the secondary component's reaction
777 36 : const unsigned sec_row = _model.eqm_species_index[react.first];
778 324 : for (unsigned i = 0; i < num_temperatures; ++i)
779 288 : _model.kin_log10K(row, i) += stoi_coeff * _model.eqm_log10K(sec_row, i);
780 344 : for (unsigned col = 0; col < num_cols; ++col)
781 308 : _model.kin_stoichiometry(row, col) += stoi_coeff * _model.eqm_stoichiometry(sec_row, col);
782 : }
783 : else
784 0 : mooseError("Kinetic species " + species.name + " includes " + react.first +
785 : ", which cannot be expressed in terms of the basis. Previous checks must be "
786 : "erroneous!");
787 : }
788 : }
789 :
790 : // check that there are no repeated sorbing sites in the SurfaceComplexationInfo
791 : std::vector<std::string> all_sorbing_sites;
792 1262 : for (const auto & name_info : _model.surface_complexation_info)
793 322 : for (const auto & name_frac : name_info.second.sorption_sites)
794 215 : if (std::find(all_sorbing_sites.begin(), all_sorbing_sites.end(), name_frac.first) !=
795 : all_sorbing_sites.end())
796 1 : mooseError(
797 : "The sorbing site ", name_frac.first, " appears in more than one sorbing mineral");
798 : else
799 214 : all_sorbing_sites.push_back(name_frac.first);
800 :
801 : // build the information related to surface sorption, and modify the species_transported vectors
802 1154 : _model.surface_sorption_related.assign(num_rows, false);
803 1154 : _model.surface_sorption_number.assign(num_rows, 99);
804 1154 : for (const auto & name_info :
805 1259 : _model.surface_complexation_info) // all minerals involved in surface complexation
806 : {
807 106 : for (const auto & name_frac :
808 318 : name_info.second.sorption_sites) // all sorption sites on the given mineral
809 : {
810 212 : const unsigned basis_index_of_sorption_site = _model.basis_species_index.at(name_frac.first);
811 : _model.basis_species_transported[basis_index_of_sorption_site] = false;
812 : }
813 : bool mineral_involved_in_eqm = false;
814 106 : for (const auto & name_frac :
815 318 : name_info.second.sorption_sites) // all sorption sites on the given mineral
816 : {
817 212 : const unsigned basis_index_of_sorption_site = _model.basis_species_index.at(name_frac.first);
818 1253 : for (unsigned j = 0; j < num_rows; ++j) // all equilibrium species
819 1225 : if (_model.eqm_stoichiometry(j, basis_index_of_sorption_site) != 0.0)
820 : {
821 : mineral_involved_in_eqm = true;
822 : break;
823 : }
824 : }
825 106 : if (!mineral_involved_in_eqm)
826 4 : continue;
827 102 : const unsigned num_surface_sorption = _model.surface_sorption_name.size();
828 102 : _model.surface_sorption_name.push_back(name_info.first);
829 102 : _model.surface_sorption_area.push_back(name_info.second.surface_area);
830 102 : for (const auto & name_frac :
831 305 : name_info.second.sorption_sites) // all sorption sites on the given mineral
832 : {
833 204 : const unsigned basis_index_of_sorption_site = _model.basis_species_index.at(name_frac.first);
834 1920 : for (unsigned j = 0; j < num_rows; ++j) // all equilibrium species
835 1717 : if (_model.eqm_stoichiometry(j, basis_index_of_sorption_site) != 0.0)
836 : {
837 413 : if (_model.surface_sorption_related[j])
838 1 : mooseError("It is an error for any equilibrium species (such as ",
839 : _model.eqm_species_name[j],
840 : ") to have a reaction involving more than one sorbing site");
841 : _model.surface_sorption_related[j] = true;
842 412 : _model.surface_sorption_number[j] = num_surface_sorption;
843 : _model.eqm_species_transported[j] = false;
844 : }
845 : }
846 : }
847 2828 : }
848 :
849 : const ModelGeochemicalDatabase &
850 112153 : PertinentGeochemicalSystem::modelGeochemicalDatabase() const
851 : {
852 112153 : return _model;
853 : }
854 :
855 : void
856 111 : PertinentGeochemicalSystem::addKineticRate(const KineticRateUserDescription & description)
857 : {
858 111 : const std::string kinetic_species = description.kinetic_species_name;
859 : if (_model.kin_species_index.count(kinetic_species) == 0)
860 1 : mooseError("Cannot prescribe a kinetic rate to species ",
861 : kinetic_species,
862 : " since it is not a kinetic species");
863 110 : const unsigned kinetic_species_index = _model.kin_species_index.at(kinetic_species);
864 :
865 : // build the promoting index list
866 110 : const unsigned num_pro = description.promoting_species.size();
867 110 : const unsigned num_basis = _model.basis_species_name.size();
868 110 : const unsigned num_eqm = _model.eqm_species_name.size();
869 113 : std::vector<Real> promoting_ind(num_basis + num_eqm, 0.0);
870 112 : std::vector<Real> promoting_m_ind(num_basis + num_eqm, 0.0);
871 112 : std::vector<Real> promoting_k(num_basis + num_eqm, 0.0);
872 247 : for (unsigned i = 0; i < num_pro; ++i)
873 : {
874 : unsigned index = 0;
875 138 : const std::string promoting_species = description.promoting_species[i];
876 : if (_model.basis_species_index.count(promoting_species) == 1)
877 95 : index = _model.basis_species_index.at(promoting_species);
878 : else if (_model.eqm_species_index.count(promoting_species) == 1)
879 42 : index = num_basis + _model.eqm_species_index.at(promoting_species);
880 : else
881 1 : mooseError(
882 : "Promoting species ", promoting_species, " must be a basis or a secondary species");
883 137 : promoting_ind[index] = description.promoting_indices[i];
884 137 : promoting_m_ind[index] = description.promoting_monod_indices[i];
885 137 : promoting_k[index] = description.promoting_half_saturation[i];
886 : }
887 : unsigned progeny_num = 0;
888 109 : if (_model.basis_species_index.count(description.progeny) == 1)
889 105 : progeny_num = _model.basis_species_index.at(description.progeny);
890 : else if (_model.eqm_species_index.count(description.progeny) == 1)
891 3 : progeny_num = num_basis + _model.eqm_species_index.at(description.progeny);
892 : else
893 1 : mooseError("Progeny ", description.progeny, " must be a basis or a secondary species");
894 :
895 : // append the result to kin_rate
896 218 : _model.kin_rate.push_back(KineticRateDefinition(kinetic_species_index,
897 : promoting_ind,
898 : promoting_m_ind,
899 : promoting_k,
900 : progeny_num,
901 : description));
902 108 : }
|
