https://mooseframework.inl.gov
Loading...
Searching...
No Matches
GeochemicalSolverTest.C
Go to the documentation of this file.
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 "gtest/gtest.h"
11
12#include "GeochemicalSolver.h"
13#include "GeochemistryKineticRateCalculator.h"
14
15const GeochemicalDatabaseReader db_solver("database/moose_testdb.json", true, true, false);
16const GeochemicalDatabaseReader db_full("../database/moose_geochemdb.json", true, true, false);
17const GeochemicalDatabaseReader
18 db_ferric("../test/database/ferric_hydroxide_sorption.json", true, true);
19// Following model only has OH- as an equilibrium species
20const PertinentGeochemicalSystem
21 model_simplest(db_solver, {"H2O", "H+"}, {}, {}, {}, {}, {}, "O2(aq)", "e-");
22GeochemistrySpeciesSwapper swapper2(2, 1E-6);
23GeochemistrySpeciesSwapper swapper_kin(4, 1E-6);
27const std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu2 = {
28 GeochemistryUnitConverter::GeochemistryUnit::KG,
29 GeochemistryUnitConverter::GeochemistryUnit::MOLES};
40GeochemistryIonicStrength is_solver(3.0, 3.0, false, false);
41GeochemistryActivityCoefficientsDebyeHuckel ac_solver(is_solver, db_solver);
42
44TEST(GeochemicalSolverTest, exception)
45{
46 ModelGeochemicalDatabase mgd = model_simplest.modelGeochemicalDatabase();
47 GeochemicalSystem egs(mgd,
48 ac_solver,
49 is_solver,
50 swapper2,
51 {},
52 {},
53 "H+",
54 {"H2O", "H+"},
55 {1.75, 3.0},
56 cu2,
57 cm2,
58 25,
59 0,
60 1E-20,
61 {},
62 {},
63 {});
64 try
65 {
66 const GeochemicalSolver solver(mgd.basis_species_name.size(),
67 mgd.kin_species_name.size(),
68 is_solver,
69 1.0,
70 0.1,
71 1,
72 1E100,
73 0.1,
74 1,
75 {},
76 3.0,
77 10,
78 true);
79 FAIL() << "Missing expected exception.";
80 }
81 catch (const std::exception & e)
82 {
83 std::string msg(e.what());
84 ASSERT_TRUE(msg.find("GeochemicalSolver: ramp_max_ionic_strength must be less than max_iter") !=
85 std::string::npos)
86 << "Failed with unexpected error message: " << msg;
87 }
88
89 try
90 {
91 const GeochemicalSolver solver(mgd.basis_species_name.size(),
92 mgd.kin_species_name.size(),
93 is_solver,
94 1.0,
95 0.1,
96 100,
97 0.0,
98 0.1,
99 1,
100 {},
101 3.0,
102 10,
103 true);
104 FAIL() << "Missing expected exception.";
105 }
106 catch (const std::exception & e)
107 {
108 std::string msg(e.what());
109 ASSERT_TRUE(msg.find("GeochemicalSolver: max_initial_residual must be positive") !=
110 std::string::npos)
111 << "Failed with unexpected error message: " << msg;
112 }
113
114 try
115 {
116 const GeochemicalSolver solver(mgd.basis_species_name.size(),
117 mgd.kin_species_name.size(),
118 is_solver,
119 1.0,
120 0.1,
121 100,
122 1E100,
123 0.1,
124 1,
125 {},
126 -1.0,
127 10,
128 true);
129 FAIL() << "Missing expected exception.";
130 }
131 catch (const std::exception & e)
132 {
133 std::string msg(e.what());
134 ASSERT_TRUE(msg.find("GeochemicalSolver: max_ionic_strength must not be negative") !=
135 std::string::npos)
136 << "Failed with unexpected error message: " << msg;
137 }
138
139 try
140 {
141 const GeochemicalSolver solver(mgd.basis_species_name.size(),
142 mgd.kin_species_name.size(),
143 is_solver,
144 1.0,
145 -0.1,
146 100,
147 1E100,
148 0.1,
149 1,
150 {},
151 1.0,
152 10,
153 true);
154 FAIL() << "Missing expected exception.";
155 }
156 catch (const std::exception & e)
157 {
158 std::string msg(e.what());
159 ASSERT_TRUE(msg.find("GeochemicalSolver: rel_tol must not be negative") != std::string::npos)
160 << "Failed with unexpected error message: " << msg;
161 }
162
163 try
164 {
165 const GeochemicalSolver solver(mgd.basis_species_name.size(),
166 mgd.kin_species_name.size(),
167 is_solver,
168 -1.0,
169 0.1,
170 100,
171 1E100,
172 0.1,
173 1,
174 {},
175 1.0,
176 10,
177 true);
178 FAIL() << "Missing expected exception.";
179 }
180 catch (const std::exception & e)
181 {
182 std::string msg(e.what());
183 ASSERT_TRUE(msg.find("GeochemicalSolver: abs_tol must not be negative") != std::string::npos)
184 << "Failed with unexpected error message: " << msg;
185 }
186
187 try
188 {
189 const GeochemicalSolver solver(mgd.basis_species_name.size(),
190 mgd.kin_species_name.size(),
191 is_solver,
192 0.0,
193 0.0,
194 100,
195 1E100,
196 0.1,
197 1,
198 {},
199 1.0,
200 10,
201 true);
202 FAIL() << "Missing expected exception.";
203 }
204 catch (const std::exception & e)
205 {
206 std::string msg(e.what());
207 ASSERT_TRUE(msg.find("GeochemicalSolver: either rel_tol or abs_tol must be positive") !=
208 std::string::npos)
209 << "Failed with unexpected error message: " << msg;
210 }
211
212 try
213 {
214 GeochemicalSolver solver(mgd.basis_species_name.size(),
215 mgd.kin_species_name.size(),
216 is_solver,
217 1.0,
218 0.1,
219 100,
220 1E100,
221 0.1,
222 1,
223 {},
224 1.0,
225 10,
226 true);
227 solver.setMaxInitialResidual(0.0);
228 FAIL() << "Missing expected exception.";
229 }
230 catch (const std::exception & e)
231 {
232 std::string msg(e.what());
233 ASSERT_TRUE(msg.find("GeochemicalSolver: max_initial_residual must be positive") !=
234 std::string::npos)
235 << "Failed with unexpected error message: " << msg;
236 }
237}
238
240TEST(GeochemicalSolverTest, setgetMaxInitialResidual)
241{
242 ModelGeochemicalDatabase mgd = model_simplest.modelGeochemicalDatabase();
243 GeochemicalSystem egs(mgd,
244 ac_solver,
245 is_solver,
246 swapper2,
247 {},
248 {},
249 "H+",
250 {"H2O", "H+"},
251 {1.75, 3.0},
252 cu2,
253 cm2,
254 25,
255 0,
256 1E-20,
257 {},
258 {},
259 {});
260 GeochemicalSolver solver(mgd.basis_species_name.size(),
261 mgd.kin_species_name.size(),
262 is_solver,
263 1.0,
264 0.1,
265 100,
266 99.0,
267 0.1,
268 1,
269 {},
270 1.0,
271 10,
272 true);
273 EXPECT_EQ(solver.getMaxInitialResidual(), 99.0);
274 solver.setMaxInitialResidual(123.0);
275 EXPECT_EQ(solver.getMaxInitialResidual(), 123.0);
276}
277
279TEST(GeochemicalSolverTest, solve1)
280{
281 ModelGeochemicalDatabase mgd = model_simplest.modelGeochemicalDatabase();
282 GeochemicalSystem egs(mgd,
283 ac_solver,
284 is_solver,
285 swapper2,
286 {},
287 {},
288 "H+",
289 {"H2O", "H+"},
290 {1.75, 1},
291 cu2,
292 cm2,
293 25,
294 0,
295 1E-20,
296 {},
297 {},
298 {});
299 GeochemicalSolver solver(mgd.basis_species_name.size(),
300 mgd.kin_species_name.size(),
301 is_solver,
302 1.0E-15,
303 1.0E-200,
304 100,
305 10.0,
306 0.1,
307 1,
308 {},
309 3.0,
310 0,
311 false);
312
313 std::stringstream ss;
314 unsigned tot_iter;
315 Real abs_residual;
316 DenseVector<Real> mole_additions(egs.getNumInBasis() + egs.getNumKinetic());
317 DenseMatrix<Real> dmole_additions(egs.getNumInBasis() + egs.getNumKinetic(),
318 egs.getNumInBasis() + egs.getNumKinetic());
319 solver.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
320
321 // check Newton has converged
322 EXPECT_LE(tot_iter, (unsigned)100);
323 EXPECT_LE(abs_residual, 1.0E-15);
324
325 // check that equilibrium molality is set correctly
326 EXPECT_NEAR(egs.getEquilibriumMolality(0) /
327 (egs.getBasisActivity(0) / egs.getBasisActivity(1) /
328 std::pow(10.0, egs.getLog10K(0)) / egs.getEquilibriumActivityCoefficient(0)),
329 1.0,
330 1.0E-9);
331
332 // check that basis molality is correct
333 EXPECT_EQ(egs.getSolventWaterMass(), 1.75);
334 EXPECT_NEAR(egs.getSolventMassAndFreeMolalityAndMineralMoles()[1] / egs.getEquilibriumMolality(0),
335 1.0,
336 1.0E-9);
337
338 // check that basis bulk composition is correct
339 EXPECT_NEAR(egs.getBulkMolesOld()[0],
340 1.75 * (GeochemistryConstants::MOLES_PER_KG_WATER + egs.getEquilibriumMolality(0)),
341 1.0E-9);
342 EXPECT_NEAR(egs.getBulkMolesOld()[1], 0.0, 1.0E-9);
343}
344
346TEST(GeochemicalSolverTest, solve2)
347{
348 // build the model
349 const PertinentGeochemicalSystem model(
350 db_full,
351 {"H2O", "H+", "Cl-", "Na+", "SO4--", "Mg++", "Ca++", "K+", "HCO3-", "SiO2(aq)", "O2(aq)"},
352 {"Antigorite",
353 "Tremolite",
354 "Talc",
355 "Chrysotile",
356 "Sepiolite",
357 "Anthophyllite",
358 "Dolomite",
359 "Dolomite-ord",
360 "Huntite",
361 "Dolomite-dis",
362 "Magnesite",
363 "Calcite",
364 "Aragonite",
365 "Quartz"},
366 {"O2(g)", "CO2(g)"},
367 {},
368 {},
369 {},
370 "O2(aq)",
371 "e-");
372 ModelGeochemicalDatabase mgd = model.modelGeochemicalDatabase();
373
374 // build the equilibrium system
375 GeochemistrySpeciesSwapper swapper(11, 1E-6);
376 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
377 GeochemicalSystem::ConstraintUserMeaningEnum::KG_SOLVENT_WATER,
378 GeochemicalSystem::ConstraintUserMeaningEnum::FUGACITY,
379 GeochemicalSystem::ConstraintUserMeaningEnum::FUGACITY,
380 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
381 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
382 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
383 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
384 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
385 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
386 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
387 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION};
388 const std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu = {
389 GeochemistryUnitConverter::GeochemistryUnit::KG,
390 GeochemistryUnitConverter::GeochemistryUnit::DIMENSIONLESS,
391 GeochemistryUnitConverter::GeochemistryUnit::DIMENSIONLESS,
392 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
393 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
394 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
395 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
396 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
397 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
398 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
399 GeochemistryUnitConverter::GeochemistryUnit::MOLES};
400 GeochemicalSystem egs(
401 mgd,
402 ac_solver,
403 is_solver,
404 swapper,
405 {"H+", "O2(aq)"},
406 {"CO2(g)", "O2(g)"},
407 "Cl-",
408 {"H2O", "CO2(g)", "O2(g)", "Cl-", "Na+", "SO4--", "Mg++", "Ca++", "K+", "HCO3-", "SiO2(aq)"},
409 {1.0,
410 0.0003162278,
411 0.2,
412 0.5656,
413 0.4850,
414 0.02924,
415 0.05501,
416 0.01063,
417 0.010576055,
418 0.002412,
419 0.00010349},
420 cu,
421 cm,
422 25,
423 0,
424 1E-20,
425 {},
426 {},
427 {});
428
429 // build solver
430 GeochemicalSolver solver(mgd.basis_species_name.size(),
431 mgd.kin_species_name.size(),
432 is_solver,
433 1.0E-15,
434 1.0E-200,
435 100,
436 1E3,
437 0.1,
438 1,
439 {"Antigorite",
440 "Tremolite",
441 "Talc",
442 "Chrysotile",
443 "Sepiolite",
444 "Anthophyllite",
445 "Dolomite",
446 "Dolomite-ord",
447 "Huntite",
448 "Dolomite-dis",
449 "Magnesite",
450 "Calcite",
451 "Aragonite",
452 "Quartz"},
453 3.0,
454 10,
455 false);
456
457 // solve
458 std::stringstream ss;
459 unsigned tot_iter;
460 Real abs_residual;
461 DenseVector<Real> mole_additions(egs.getNumInBasis() + egs.getNumKinetic());
462 DenseMatrix<Real> dmole_additions(egs.getNumInBasis() + egs.getNumKinetic(),
463 egs.getNumInBasis() + egs.getNumKinetic());
464 solver.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
465
466 // check converged
467 EXPECT_LE(tot_iter, (unsigned)100);
468 EXPECT_LE(abs_residual, 1.0E-15);
469
470 // check number in basis, number in redox disequilibrium and number of surface potentials
471 EXPECT_EQ(egs.getNumInBasis(), (unsigned)11);
472 EXPECT_EQ(egs.getNumRedox(), (unsigned)1);
473 EXPECT_EQ(egs.getNumSurfacePotentials(), (unsigned)0);
474 EXPECT_EQ(egs.getNumInEquilibrium(), mgd.eqm_species_name.size());
475 EXPECT_EQ(egs.getNumKinetic(), (unsigned)0);
476
477 // check that the constraints are satisfied
478 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
479 if (mgd.basis_species_name[i] == "H2O")
480 {
481 EXPECT_FALSE(mgd.basis_species_gas[i]);
482 EXPECT_FALSE(mgd.basis_species_mineral[i]);
483 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
484 EXPECT_NEAR(egs.getSolventWaterMass(), 1.0, 1.0E-15);
485 EXPECT_NEAR(egs.getSolventMassAndFreeMolalityAndMineralMoles()[0], 1.0, 1.0E-15);
486 }
487 else if (mgd.basis_species_name[i] == "CO2(g)")
488 {
489 EXPECT_TRUE(mgd.basis_species_gas[i]);
490 EXPECT_FALSE(mgd.basis_species_mineral[i]);
491 EXPECT_TRUE(egs.getBasisActivityKnown()[i]);
492 EXPECT_NEAR(egs.getBasisActivity(i), 0.0003162278, 1.0E-15);
493 }
494 else if (mgd.basis_species_name[i] == "O2(g)")
495 {
496 EXPECT_TRUE(mgd.basis_species_gas[i]);
497 EXPECT_FALSE(mgd.basis_species_mineral[i]);
498 EXPECT_TRUE(egs.getBasisActivityKnown()[i]);
499 EXPECT_NEAR(egs.getBasisActivity(i), 0.2, 1.0E-15);
500 }
501 else if (mgd.basis_species_name[i] == "Cl-")
502 {
503 EXPECT_FALSE(mgd.basis_species_gas[i]);
504 EXPECT_FALSE(mgd.basis_species_mineral[i]);
505 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
506 // do not know the bulk composition as it is dictated by charge neutrality
507 EXPECT_EQ(egs.getChargeBalanceBasisIndex(), i);
508 }
509 else if (mgd.basis_species_name[i] == "Na+")
510 {
511 EXPECT_FALSE(mgd.basis_species_gas[i]);
512 EXPECT_FALSE(mgd.basis_species_mineral[i]);
513 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
514 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.4850, 1.0E-15);
515 }
516 else if (mgd.basis_species_name[i] == "SO4--")
517 {
518 EXPECT_FALSE(mgd.basis_species_gas[i]);
519 EXPECT_FALSE(mgd.basis_species_mineral[i]);
520 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
521 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.02924, 1.0E-15);
522 }
523 else if (mgd.basis_species_name[i] == "Mg++")
524 {
525 EXPECT_FALSE(mgd.basis_species_gas[i]);
526 EXPECT_FALSE(mgd.basis_species_mineral[i]);
527 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
528 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.05501, 1.0E-15);
529 }
530 else if (mgd.basis_species_name[i] == "Ca++")
531 {
532 EXPECT_FALSE(mgd.basis_species_gas[i]);
533 EXPECT_FALSE(mgd.basis_species_mineral[i]);
534 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
535 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.01063, 1.0E-15);
536 }
537 else if (mgd.basis_species_name[i] == "K+")
538 {
539 EXPECT_FALSE(mgd.basis_species_gas[i]);
540 EXPECT_FALSE(mgd.basis_species_mineral[i]);
541 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
542 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.010576055, 1.0E-15);
543 }
544 else if (mgd.basis_species_name[i] == "HCO3-")
545 {
546 EXPECT_FALSE(mgd.basis_species_gas[i]);
547 EXPECT_FALSE(mgd.basis_species_mineral[i]);
548 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
549 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.002412, 1.0E-15);
550 }
551 else if (mgd.basis_species_name[i] == "SiO2(aq)")
552 {
553 EXPECT_FALSE(mgd.basis_species_gas[i]);
554 EXPECT_FALSE(mgd.basis_species_mineral[i]);
555 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
556 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.00010349, 1.0E-15);
557 }
558 else
559 FAIL() << "Incorrect basis species";
560
561 // check total charge
562 EXPECT_NEAR(egs.getTotalChargeOld(), 0.0, 1E-15);
563 // check total charge by summing up basis charges
564 Real tot_charge = 0.0;
565 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
566 tot_charge += egs.getBulkMolesOld()[i] * mgd.basis_species_charge[i];
567 EXPECT_NEAR(tot_charge, 0.0, 1E-15);
568
569 // surface charges
570 for (unsigned sp = 0; sp < egs.getNumSurfacePotentials(); ++sp)
571 {
572 Real charge = 0.0;
573 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
574 if (mgd.surface_sorption_related[j] && mgd.surface_sorption_number[j] == sp)
575 charge += mgd.eqm_species_charge[j] * egs.getEquilibriumMolality(j);
576 EXPECT_NEAR(charge,
577 egs.getSorbingSurfaceArea()[sp] * egs.getSurfaceCharge(sp) /
578 GeochemistryConstants::FARADAY,
579 1E-15);
580 }
581
582 // check basis activity = activity_coefficient * molality
583 for (unsigned i = 1; i < egs.getNumInBasis(); ++i) // don't loop over water
584 if (mgd.basis_species_gas[i] || mgd.basis_species_mineral[i] || egs.getBasisActivityKnown()[i])
585 continue;
586 else
587 EXPECT_NEAR(egs.getBasisActivity(i),
588 egs.getBasisActivityCoefficient(i) *
589 egs.getSolventMassAndFreeMolalityAndMineralMoles()[i],
590 1E-15);
591
592 // check residuals are zero
593 for (unsigned a = 0; a < egs.getNumInAlgebraicSystem(); ++a)
594 EXPECT_LE(std::abs(egs.getResidualComponent(a, mole_additions)), 1E-15);
595 // check residuals are zero by summing molalities, bulk compositions, etc
596 Real nw = egs.getSolventWaterMass();
597 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
598 {
599 Real res = -egs.getBulkMolesOld()[i];
600 if (i == 0)
601 res += nw * GeochemistryConstants::MOLES_PER_KG_WATER;
602 else if (mgd.basis_species_mineral[i])
603 res += egs.getSolventMassAndFreeMolalityAndMineralMoles()[i];
604 else if (mgd.basis_species_gas[i])
605 res += 0.0;
606 else
607 res += nw * egs.getSolventMassAndFreeMolalityAndMineralMoles()[i];
608 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
609 res += nw * mgd.eqm_stoichiometry(j, i) * egs.getEquilibriumMolality(j);
610 EXPECT_LE(std::abs(res), 1E-13);
611 }
612 // surface potentials
613 const Real prefactor = std::sqrt(GeochemistryConstants::GAS_CONSTANT *
614 (25.0 + GeochemistryConstants::CELSIUS_TO_KELVIN) *
615 GeochemistryConstants::PERMITTIVITY_FREE_SPACE *
616 GeochemistryConstants::DIELECTRIC_CONSTANT_WATER *
617 GeochemistryConstants::DENSITY_WATER * egs.getIonicStrength()) /
618 nw;
619 for (unsigned sp = 0; sp < egs.getNumSurfacePotentials(); ++sp)
620 {
621 EXPECT_NEAR(prefactor * std::sinh(egs.getSurfacePotential(sp) * GeochemistryConstants::FARADAY /
622 2.0 / GeochemistryConstants::GAS_CONSTANT /
623 (25.0 + GeochemistryConstants::CELSIUS_TO_KELVIN)),
624 egs.getSurfaceCharge(sp),
625 1.0E-13);
626 }
627
628 // check equilibrium mass balance
629 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
630 {
631 if (mgd.eqm_species_mineral[j] || mgd.eqm_species_gas[j])
632 continue;
633 Real log10ap = 0.0;
634 for (unsigned i = 0; i < mgd.basis_species_name.size(); ++i)
635 log10ap += mgd.eqm_stoichiometry(j, i) * std::log10(egs.getBasisActivity(i));
636 log10ap -= egs.getLog10K(j);
637 if (mgd.surface_sorption_related[j])
638 log10ap -= std::log10(egs.getSurfacePotential(mgd.surface_sorption_number[j]));
639 else
640 log10ap -= std::log10(egs.getEquilibriumActivityCoefficient(j));
641 if (log10ap < -300.0)
642 EXPECT_LE(std::abs(egs.getEquilibriumMolality(j)), 1E-25);
643 else
644 EXPECT_LE(std::abs(egs.getEquilibriumMolality(j) - std::pow(10.0, log10ap)), 1E-15);
645 }
646}
647
649TEST(GeochemicalSolverTest, solve3)
650{
651 // build the model
652 const PertinentGeochemicalSystem model(
653 db_full,
654 {"H2O", "H+", "Cl-", "Na+", "SO4--", "Mg++", "Ca++", "K+", "HCO3-", "SiO2(aq)", "O2(aq)"},
655 {"Antigorite",
656 "Tremolite",
657 "Talc",
658 "Chrysotile",
659 "Sepiolite",
660 "Anthophyllite",
661 "Dolomite",
662 "Dolomite-ord",
663 "Huntite",
664 "Dolomite-dis",
665 "Magnesite",
666 "Calcite",
667 "Aragonite",
668 "Quartz"},
669 {},
670 {},
671 {},
672 {},
673 "O2(aq)",
674 "e-");
675 ModelGeochemicalDatabase mgd = model.modelGeochemicalDatabase();
676
677 // build the equilibrium system
678 GeochemistrySpeciesSwapper swapper(11, 1E-6);
679 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
680 GeochemicalSystem::ConstraintUserMeaningEnum::KG_SOLVENT_WATER,
681 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
682 GeochemicalSystem::ConstraintUserMeaningEnum::FREE_CONCENTRATION,
683 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
684 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
685 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
686 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
687 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
688 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
689 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
690 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION};
691 const std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu = {
692 GeochemistryUnitConverter::GeochemistryUnit::KG,
693 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
694 GeochemistryUnitConverter::GeochemistryUnit::MOLAL,
695 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
696 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
697 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
698 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
699 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
700 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
701 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
702 GeochemistryUnitConverter::GeochemistryUnit::MOLES};
703 GeochemicalSystem egs(
704 mgd,
705 ac_solver,
706 is_solver,
707 swapper,
708 {"H+"},
709 {"MgCO3"},
710 "Cl-",
711 {"H2O", "MgCO3", "O2(aq)", "Cl-", "Na+", "SO4--", "Mg++", "Ca++", "K+", "HCO3-", "SiO2(aq)"},
712 {1.0,
713 0.196E-3,
714 0.2151E-3,
715 0.5656,
716 0.4850,
717 0.02924,
718 0.05501,
719 0.01063,
720 0.010576055,
721 0.002412,
722 0.00010349},
723 cu,
724 cm,
725 25,
726 0,
727 1E-20,
728 {},
729 {},
730 {});
731
732 // build solver (4 swaps are needed to solve this system)
733 GeochemicalSolver solver(mgd.basis_species_name.size(),
734 mgd.kin_species_name.size(),
735 is_solver,
736 1.0E-15,
737 1.0E-200,
738 100,
739 1E3,
740 0.1,
741 4,
742 {"Dolomite-ord", "Dolomite-dis"},
743 3.0,
744 10,
745 false);
746
747 // solve
748 std::stringstream ss;
749 unsigned tot_iter;
750 Real abs_residual;
751 DenseVector<Real> mole_additions(egs.getNumInBasis() + egs.getNumKinetic());
752 DenseMatrix<Real> dmole_additions(egs.getNumInBasis() + egs.getNumKinetic(),
753 egs.getNumInBasis() + egs.getNumKinetic());
754 solver.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
755
756 // check converged
757 EXPECT_LE(tot_iter, (unsigned)100);
758 EXPECT_LE(abs_residual, 1.0E-15);
759
760 // check number in basis, number in redox disequilibrium and number of surface potentials
761 EXPECT_EQ(egs.getNumInBasis(), (unsigned)11);
762 EXPECT_EQ(egs.getNumRedox(), (unsigned)1);
763 EXPECT_EQ(egs.getNumSurfacePotentials(), (unsigned)0);
764 EXPECT_EQ(egs.getNumInEquilibrium(), mgd.eqm_species_name.size());
765
766 // check that the constraints are satisfied
767 Real bulk_dolo = 0.0;
768 Real bulk_ca = 0.0;
769 Real bulk_mg = 0.0;
770 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
771 if (mgd.basis_species_name[i] == "H2O")
772 {
773 EXPECT_FALSE(mgd.basis_species_gas[i]);
774 EXPECT_FALSE(mgd.basis_species_mineral[i]);
775 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
776 EXPECT_NEAR(egs.getSolventWaterMass(), 1.0, 1.0E-15);
777 EXPECT_NEAR(egs.getSolventMassAndFreeMolalityAndMineralMoles()[0], 1.0, 1.0E-15);
778 }
779 else if (mgd.basis_species_name[i] == "O2(aq)")
780 {
781 EXPECT_FALSE(mgd.basis_species_gas[i]);
782 EXPECT_FALSE(mgd.basis_species_mineral[i]);
783 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
784 EXPECT_NEAR(egs.getSolventMassAndFreeMolalityAndMineralMoles()[i], 0.2151E-3, 1.0E-15);
785 }
786 else if (mgd.basis_species_name[i] == "Cl-")
787 {
788 EXPECT_FALSE(mgd.basis_species_gas[i]);
789 EXPECT_FALSE(mgd.basis_species_mineral[i]);
790 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
791 // do not know the bulk composition as it is dictated by charge neutrality
792 EXPECT_EQ(egs.getChargeBalanceBasisIndex(), i);
793 }
794 else if (mgd.basis_species_name[i] == "Na+")
795 {
796 EXPECT_FALSE(mgd.basis_species_gas[i]);
797 EXPECT_FALSE(mgd.basis_species_mineral[i]);
798 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
799 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.4850, 1.0E-15);
800 }
801 else if (mgd.basis_species_name[i] == "SO4--")
802 {
803 EXPECT_FALSE(mgd.basis_species_gas[i]);
804 EXPECT_FALSE(mgd.basis_species_mineral[i]);
805 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
806 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.02924, 1.0E-15);
807 }
808 else if (mgd.basis_species_name[i] == "Mg++")
809 {
810 EXPECT_FALSE(mgd.basis_species_gas[i]);
811 EXPECT_FALSE(mgd.basis_species_mineral[i]);
812 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
813 bulk_mg = egs.getBulkMolesOld()[i];
814 }
815 else if (mgd.basis_species_name[i] == "Ca++")
816 {
817 EXPECT_FALSE(mgd.basis_species_gas[i]);
818 EXPECT_FALSE(mgd.basis_species_mineral[i]);
819 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
820 bulk_ca = egs.getBulkMolesOld()[i];
821 }
822 else if (mgd.basis_species_name[i] == "K+")
823 {
824 EXPECT_FALSE(mgd.basis_species_gas[i]);
825 EXPECT_FALSE(mgd.basis_species_mineral[i]);
826 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
827 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.010576055, 1.0E-15);
828 }
829 else if (mgd.basis_species_name[i] == "Quartz")
830 {
831 EXPECT_FALSE(mgd.basis_species_gas[i]);
832 EXPECT_TRUE(mgd.basis_species_mineral[i]);
833 EXPECT_TRUE(egs.getBasisActivityKnown()[i]);
834 EXPECT_GE(egs.getSolventMassAndFreeMolalityAndMineralMoles()[i], 0.0);
835 // Quartz = SiO2(aq)
836 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.00010349, 1.0);
837 }
838 else if (mgd.basis_species_name[i] == "Dolomite")
839 {
840 EXPECT_FALSE(mgd.basis_species_gas[i]);
841 EXPECT_TRUE(mgd.basis_species_mineral[i]);
842 EXPECT_TRUE(egs.getBasisActivityKnown()[i]);
843 EXPECT_GE(egs.getSolventMassAndFreeMolalityAndMineralMoles()[i], 0.0);
844 bulk_dolo = egs.getBulkMolesOld()[i];
845 }
846 else if (mgd.basis_species_name[i] == "HCO3-")
847 {
848 EXPECT_FALSE(mgd.basis_species_gas[i]);
849 EXPECT_FALSE(mgd.basis_species_mineral[i]);
850 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
851 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.002412, 1.0E-13);
852 }
853 else
854 FAIL() << "Incorrect basis species";
855 // dolomite = 2MgCO3 - Mg + Ca
856 EXPECT_NEAR(bulk_dolo + bulk_ca, 0.01063, 1E-13);
857 EXPECT_NEAR(-bulk_dolo + bulk_mg, 0.05501, 1E-13);
858 EXPECT_NEAR(2.0 * bulk_dolo, 0.196E-3, 1E-13);
859
860 // check total charge
861 EXPECT_NEAR(egs.getTotalChargeOld(), 0.0, 1E-15);
862 // check total charge by summing up basis charges
863 Real tot_charge = 0.0;
864 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
865 tot_charge += egs.getBulkMolesOld()[i] * mgd.basis_species_charge[i];
866 EXPECT_NEAR(tot_charge, 0.0, 1E-15);
867
868 // surface charges
869 for (unsigned sp = 0; sp < egs.getNumSurfacePotentials(); ++sp)
870 {
871 Real charge = 0.0;
872 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
873 if (mgd.surface_sorption_related[j] && mgd.surface_sorption_number[j] == sp)
874 charge += mgd.eqm_species_charge[j] * egs.getEquilibriumMolality(j);
875 EXPECT_NEAR(charge,
876 egs.getSorbingSurfaceArea()[sp] * egs.getSurfaceCharge(sp) /
877 GeochemistryConstants::FARADAY,
878 1E-15);
879 }
880
881 // check basis activity = activity_coefficient * molality
882 for (unsigned i = 1; i < egs.getNumInBasis(); ++i) // don't loop over water
883 if (mgd.basis_species_gas[i] || mgd.basis_species_mineral[i] || egs.getBasisActivityKnown()[i])
884 continue;
885 else
886 EXPECT_NEAR(egs.getBasisActivity(i),
887 egs.getBasisActivityCoefficient(i) *
888 egs.getSolventMassAndFreeMolalityAndMineralMoles()[i],
889 1E-15);
890
891 // check residuals are zero
892 for (unsigned a = 0; a < egs.getNumInAlgebraicSystem(); ++a)
893 EXPECT_LE(std::abs(egs.getResidualComponent(a, mole_additions)), 1E-15);
894 // check residuals are zero by summing molalities, bulk compositions, etc
895 Real nw = egs.getSolventWaterMass();
896 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
897 {
898 Real res = -egs.getBulkMolesOld()[i];
899 if (i == 0)
900 res += nw * GeochemistryConstants::MOLES_PER_KG_WATER;
901 else if (mgd.basis_species_mineral[i])
902 res += egs.getSolventMassAndFreeMolalityAndMineralMoles()[i];
903 else if (mgd.basis_species_gas[i])
904 res += 0.0;
905 else
906 res += nw * egs.getSolventMassAndFreeMolalityAndMineralMoles()[i];
907 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
908 res += nw * mgd.eqm_stoichiometry(j, i) * egs.getEquilibriumMolality(j);
909 EXPECT_LE(std::abs(res), 1E-13);
910 }
911 // surface potentials
912 const Real prefactor = std::sqrt(GeochemistryConstants::GAS_CONSTANT *
913 (25.0 + GeochemistryConstants::CELSIUS_TO_KELVIN) *
914 GeochemistryConstants::PERMITTIVITY_FREE_SPACE *
915 GeochemistryConstants::DIELECTRIC_CONSTANT_WATER *
916 GeochemistryConstants::DENSITY_WATER * egs.getIonicStrength()) /
917 nw;
918 for (unsigned sp = 0; sp < egs.getNumSurfacePotentials(); ++sp)
919 {
920 EXPECT_NEAR(prefactor * std::sinh(egs.getSurfacePotential(sp) * GeochemistryConstants::FARADAY /
921 2.0 / GeochemistryConstants::GAS_CONSTANT /
922 (25.0 + GeochemistryConstants::CELSIUS_TO_KELVIN)),
923 egs.getSurfaceCharge(sp),
924 1.0E-13);
925 }
926
927 // check equilibrium mass balance
928 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
929 {
930 if (mgd.eqm_species_mineral[j] || mgd.eqm_species_gas[j])
931 continue;
932 Real log10ap = 0.0;
933 for (unsigned i = 0; i < mgd.basis_species_name.size(); ++i)
934 log10ap += mgd.eqm_stoichiometry(j, i) * std::log10(egs.getBasisActivity(i));
935 log10ap -= egs.getLog10K(j);
936 if (mgd.surface_sorption_related[j])
937 log10ap -= std::log10(egs.getSurfacePotential(mgd.surface_sorption_number[j]));
938 else
939 log10ap -= std::log10(egs.getEquilibriumActivityCoefficient(j));
940 if (log10ap < -300.0)
941 EXPECT_LE(std::abs(egs.getEquilibriumMolality(j)), 1E-25);
942 else
943 EXPECT_LE(std::abs(egs.getEquilibriumMolality(j) - std::pow(10.0, log10ap)), 1E-15);
944 }
945
946 // check equilibrium species have negative saturation indices, except for prevent_precipitation
947 // minerals
948 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
949 {
950 if (!mgd.eqm_species_mineral[j])
951 continue;
952 if (mgd.eqm_species_name[j] == "Dolomite-ord" || mgd.eqm_species_name[j] == "Dolomite-dis")
953 continue;
954 Real log10ap = 0.0;
955 for (unsigned i = 0; i < mgd.basis_species_name.size(); ++i)
956 log10ap += mgd.eqm_stoichiometry(j, i) * std::log10(egs.getBasisActivity(i));
957 EXPECT_LE(log10ap, egs.getLog10K(j));
958 }
959}
960
962TEST(GeochemicalSolverTest, solve3_restore)
963{
964 // build the model
965 const PertinentGeochemicalSystem model(
966 db_full,
967 {"H2O", "H+", "Cl-", "Na+", "SO4--", "Mg++", "Ca++", "K+", "HCO3-", "SiO2(aq)", "O2(aq)"},
968 {"Antigorite",
969 "Tremolite",
970 "Talc",
971 "Chrysotile",
972 "Sepiolite",
973 "Anthophyllite",
974 "Dolomite",
975 "Dolomite-ord",
976 "Huntite",
977 "Dolomite-dis",
978 "Magnesite",
979 "Calcite",
980 "Aragonite",
981 "Quartz"},
982 {},
983 {},
984 {},
985 {},
986 "O2(aq)",
987 "e-");
988 ModelGeochemicalDatabase mgd = model.modelGeochemicalDatabase();
989
990 // build the equilibrium system
991 GeochemistrySpeciesSwapper swapper(11, 1E-6);
992 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
993 GeochemicalSystem::ConstraintUserMeaningEnum::KG_SOLVENT_WATER,
994 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
995 GeochemicalSystem::ConstraintUserMeaningEnum::FREE_CONCENTRATION,
996 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
997 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
998 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
999 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1000 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1001 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1002 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1003 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION};
1004 const std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu = {
1005 GeochemistryUnitConverter::GeochemistryUnit::KG,
1006 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1007 GeochemistryUnitConverter::GeochemistryUnit::MOLAL,
1008 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1009 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1010 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1011 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1012 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1013 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1014 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1015 GeochemistryUnitConverter::GeochemistryUnit::MOLES};
1016 GeochemicalSystem egs(
1017 mgd,
1018 ac_solver,
1019 is_solver,
1020 swapper,
1021 {"H+"},
1022 {"MgCO3"},
1023 "Cl-",
1024 {"H2O", "MgCO3", "O2(aq)", "Cl-", "Na+", "SO4--", "Mg++", "Ca++", "K+", "HCO3-", "SiO2(aq)"},
1025 {1.0,
1026 0.196E-3,
1027 0.2151E-3,
1028 0.5656,
1029 0.4850,
1030 0.02924,
1031 0.05501,
1032 0.01063,
1033 0.010576055,
1034 0.002412,
1035 0.00010349},
1036 cu,
1037 cm,
1038 25,
1039 0,
1040 1E-20,
1041 {},
1042 {},
1043 {});
1044
1045 // build solver
1046 GeochemicalSolver solver(mgd.basis_species_name.size(),
1047 mgd.kin_species_name.size(),
1048 is_solver,
1049 1.0E-15,
1050 1.0E-200,
1051 100,
1052 1E3,
1053 0.1,
1054 5,
1055 {"Dolomite-ord", "Dolomite-dis"},
1056 3.0,
1057 10,
1058 false);
1059
1060 // solve
1061 std::stringstream ss;
1062 unsigned tot_iter;
1063 Real abs_residual;
1064 DenseVector<Real> mole_additions(egs.getNumInBasis() + egs.getNumKinetic());
1065 DenseMatrix<Real> dmole_additions(egs.getNumInBasis() + egs.getNumKinetic(),
1066 egs.getNumInBasis() + egs.getNumKinetic());
1067 solver.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
1068 const Real old_residual = abs_residual;
1069
1070 // check converged
1071 EXPECT_LE(tot_iter, (unsigned)100);
1072 EXPECT_LE(abs_residual, 1.0E-15);
1073
1074 // retrieve the molalities, and set them into egs: this should result in no change if the
1075 // "restore" is working correctly
1076 std::vector<std::string> names = mgd.basis_species_name;
1077 names.insert(names.end(), mgd.eqm_species_name.begin(), mgd.eqm_species_name.end());
1078 std::vector<Real> molal = egs.getSolventMassAndFreeMolalityAndMineralMoles();
1079 for (unsigned j = 0; j < egs.getNumInEquilibrium(); ++j)
1080 molal.push_back(egs.getEquilibriumMolality(j));
1081 // form constraints_from_molalities, which are all false
1082 const std::vector<bool> com_false(egs.getNumInBasis(), false);
1083 egs.setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(
1084 names, molal, com_false);
1085
1086 // build solver with no ramping of the maximum ionic strength
1087 GeochemicalSolver solver0(mgd.basis_species_name.size(),
1088 mgd.kin_species_name.size(),
1089 is_solver,
1090 1.0E-15,
1091 1.0E-200,
1092 100,
1093 1E3,
1094 0.1,
1095 0,
1096 {"Dolomite-ord", "Dolomite-dis"},
1097 3.0,
1098 0,
1099 false);
1100
1101 // solve this: no swaps should be necessary
1102 solver0.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
1103
1104 // check that the soler thinks this is truly a solution and the residual has not changed
1105 EXPECT_EQ(tot_iter, (unsigned)0);
1106 EXPECT_EQ(abs_residual, old_residual);
1107
1108 // Now use constraints_from_molalities = true
1109 const std::vector<bool> com_true(egs.getNumInBasis(), true);
1110 egs.setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(
1111 names, molal, com_true);
1112
1113 solver0.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
1114
1115 // check that the soler thinks this is truly a solution and the residual has not increased
1116 EXPECT_EQ(tot_iter, (unsigned)0);
1117 EXPECT_LE(abs_residual, old_residual);
1118
1119 // now check the copy-assignment of GeochemicalSystem
1120 const PertinentGeochemicalSystem model_dest(
1121 db_full,
1122 {"H2O", "H+", "Cl-", "Na+", "SO4--", "Mg++", "Ca++", "K+", "HCO3-", "SiO2(aq)", "O2(aq)"},
1123 {"Antigorite",
1124 "Tremolite",
1125 "Talc",
1126 "Chrysotile",
1127 "Sepiolite",
1128 "Anthophyllite",
1129 "Dolomite",
1130 "Dolomite-ord",
1131 "Huntite",
1132 "Dolomite-dis",
1133 "Magnesite",
1134 "Calcite",
1135 "Aragonite",
1136 "Quartz"},
1137 {},
1138 {},
1139 {},
1140 {},
1141 "O2(aq)",
1142 "e-");
1143 ModelGeochemicalDatabase mgd_dest = model_dest.modelGeochemicalDatabase();
1144 GeochemicalSystem egs_dest(
1145 mgd_dest,
1146 ac_solver,
1147 is_solver,
1148 swapper,
1149 {"H+"},
1150 {"MgCO3"},
1151 "Cl-",
1152 {"H2O", "MgCO3", "O2(aq)", "Cl-", "Na+", "SO4--", "Mg++", "Ca++", "K+", "HCO3-", "SiO2(aq)"},
1153 {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0},
1154 cu,
1155 cm,
1156 25,
1157 0,
1158 1E-20,
1159 {},
1160 {},
1161 {});
1162
1163 // change constraint meanings in egs to provide a more thorough check of the copy assignment
1164 egs.closeSystem();
1165 // copy assignment
1166 egs_dest = egs;
1167 EXPECT_EQ(egs_dest, egs);
1168}
1169
1171TEST(GeochemicalSolverTest, solve4)
1172{
1173 // build the model
1174 const PertinentGeochemicalSystem model(db_full,
1175 {"H2O",
1176 "H+",
1177 "Cl-",
1178 "O2(aq)",
1179 "HCO3-",
1180 "Ca++",
1181 "Mg++",
1182 "Na+",
1183 "K+",
1184 "Fe++",
1185 "Fe+++",
1186 "Mn++",
1187 "Zn++",
1188 "SO4--"},
1189 {},
1190 {},
1191 {},
1192 {},
1193 {},
1194 "O2(aq)",
1195 "e-");
1196 ModelGeochemicalDatabase mgd = model.modelGeochemicalDatabase();
1197
1198 // build the equilibrium system
1199 GeochemistrySpeciesSwapper swapper(14, 1E-6);
1200 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
1201 GeochemicalSystem::ConstraintUserMeaningEnum::KG_SOLVENT_WATER,
1202 GeochemicalSystem::ConstraintUserMeaningEnum::ACTIVITY,
1203 GeochemicalSystem::ConstraintUserMeaningEnum::FREE_CONCENTRATION,
1204 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1205 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1206 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1207 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1208 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1209 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1210 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1211 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1212 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1213 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1214 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION};
1215 const std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu = {
1216 GeochemistryUnitConverter::GeochemistryUnit::KG,
1217 GeochemistryUnitConverter::GeochemistryUnit::DIMENSIONLESS,
1218 GeochemistryUnitConverter::GeochemistryUnit::MOLAL,
1219 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1220 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1221 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1222 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1223 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1224 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1225 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1226 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1227 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1228 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1229 GeochemistryUnitConverter::GeochemistryUnit::MOLES};
1230 GeochemicalSystem egs(mgd,
1231 ac_solver,
1232 is_solver,
1233 swapper,
1234 {},
1235 {},
1236 "Cl-",
1237 {"H2O",
1238 "H+",
1239 "O2(aq)",
1240 "Cl-",
1241 "HCO3-",
1242 "Ca++",
1243 "Mg++",
1244 "Na+",
1245 "K+",
1246 "Fe++",
1247 "Fe+++",
1248 "Mn++",
1249 "Zn++",
1250 "SO4--"},
1251 {1.0,
1252 0.8913E-6,
1253 0.13438E-3,
1254 3.041E-5,
1255 0.0295E-3,
1256 0.005938E-3,
1257 0.01448E-3,
1258 0.0018704E-3,
1259 0.005115E-3,
1260 0.012534E-3,
1261 0.0005372E-3,
1262 0.005042E-3,
1263 0.001897E-3,
1264 0.01562E-3},
1265 cu,
1266 cm,
1267 25,
1268 0,
1269 1E-20,
1270 {},
1271 {},
1272 {});
1273
1274 // build solver
1275 GeochemicalSolver solver(mgd.basis_species_name.size(),
1276 mgd.kin_species_name.size(),
1277 is_solver,
1278 1.0E-15,
1279 1.0E-200,
1280 100,
1281 1E-2,
1282 0.1,
1283 1,
1284 {},
1285 3.0,
1286 10,
1287 false);
1288
1289 // solve
1290 std::stringstream ss;
1291 unsigned tot_iter;
1292 Real abs_residual;
1293 DenseVector<Real> mole_additions(egs.getNumInBasis() + egs.getNumKinetic());
1294 DenseMatrix<Real> dmole_additions(egs.getNumInBasis() + egs.getNumKinetic(),
1295 egs.getNumInBasis() + egs.getNumKinetic());
1296 solver.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
1297
1298 // check converged
1299 EXPECT_LE(tot_iter, (unsigned)100);
1300 EXPECT_LE(abs_residual, 1.0E-15);
1301
1302 // check number in basis, number in redox disequilibrium and number of surface potentials
1303 EXPECT_EQ(egs.getNumInBasis(), (unsigned)14);
1304 EXPECT_EQ(egs.getNumRedox(), (unsigned)2);
1305 EXPECT_EQ(egs.getNumSurfacePotentials(), (unsigned)0);
1306 EXPECT_EQ(egs.getNumInEquilibrium(), mgd.eqm_species_name.size());
1307
1308 // check that the constraints are satisfied
1309 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
1310 if (mgd.basis_species_name[i] == "H2O")
1311 {
1312 EXPECT_FALSE(mgd.basis_species_gas[i]);
1313 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1314 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1315 EXPECT_NEAR(egs.getSolventWaterMass(), 1.0, 1.0E-15);
1316 EXPECT_NEAR(egs.getSolventMassAndFreeMolalityAndMineralMoles()[0], 1.0, 1.0E-15);
1317 }
1318 else if (mgd.basis_species_name[i] == "H+")
1319 {
1320 EXPECT_FALSE(mgd.basis_species_gas[i]);
1321 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1322 EXPECT_TRUE(egs.getBasisActivityKnown()[i]);
1323 EXPECT_NEAR(egs.getBasisActivity(i), 0.8913E-6, 1.0E-15);
1324 }
1325 else if (mgd.basis_species_name[i] == "O2(aq)")
1326 {
1327 EXPECT_FALSE(mgd.basis_species_gas[i]);
1328 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1329 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1330 EXPECT_NEAR(egs.getSolventMassAndFreeMolalityAndMineralMoles()[i], 0.13438E-3, 1.0E-15);
1331 }
1332 else if (mgd.basis_species_name[i] == "Cl-")
1333 {
1334 EXPECT_FALSE(mgd.basis_species_gas[i]);
1335 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1336 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1337 // do not know the bulk composition as it is dictated by charge neutrality
1338 EXPECT_EQ(egs.getChargeBalanceBasisIndex(), i);
1339 }
1340 else if (mgd.basis_species_name[i] == "HCO3-")
1341 {
1342 EXPECT_FALSE(mgd.basis_species_gas[i]);
1343 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1344 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1345 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.0295E-3, 1.0E-15);
1346 }
1347 else if (mgd.basis_species_name[i] == "Ca++")
1348 {
1349 EXPECT_FALSE(mgd.basis_species_gas[i]);
1350 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1351 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1352 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.005938E-3, 1.0E-15);
1353 }
1354 else if (mgd.basis_species_name[i] == "Mg++")
1355 {
1356 EXPECT_FALSE(mgd.basis_species_gas[i]);
1357 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1358 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1359 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.01448E-3, 1.0E-15);
1360 }
1361 else if (mgd.basis_species_name[i] == "Na+")
1362 {
1363 EXPECT_FALSE(mgd.basis_species_gas[i]);
1364 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1365 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1366 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.0018704E-3, 1.0E-15);
1367 }
1368 else if (mgd.basis_species_name[i] == "K+")
1369 {
1370 EXPECT_FALSE(mgd.basis_species_gas[i]);
1371 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1372 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1373 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.005115E-3, 1.0E-15);
1374 }
1375 else if (mgd.basis_species_name[i] == "Fe++")
1376 {
1377 EXPECT_FALSE(mgd.basis_species_gas[i]);
1378 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1379 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1380 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.012534E-3, 1.0E-15);
1381 }
1382 else if (mgd.basis_species_name[i] == "Fe+++")
1383 {
1384 EXPECT_FALSE(mgd.basis_species_gas[i]);
1385 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1386 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1387 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.0005372E-3, 1.0E-15);
1388 }
1389 else if (mgd.basis_species_name[i] == "Mn++")
1390 {
1391 EXPECT_FALSE(mgd.basis_species_gas[i]);
1392 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1393 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1394 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.005042E-3, 1.0E-15);
1395 }
1396 else if (mgd.basis_species_name[i] == "Zn++")
1397 {
1398 EXPECT_FALSE(mgd.basis_species_gas[i]);
1399 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1400 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1401 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.001897E-3, 1.0E-15);
1402 }
1403 else if (mgd.basis_species_name[i] == "SO4--")
1404 {
1405 EXPECT_FALSE(mgd.basis_species_gas[i]);
1406 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1407 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1408 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.01562E-3, 1.0E-15);
1409 }
1410 else
1411 FAIL() << "Incorrect basis species";
1412
1413 // check total charge
1414 EXPECT_NEAR(egs.getTotalChargeOld(), 0.0, 1E-15);
1415 // check total charge by summing up basis charges
1416 Real tot_charge = 0.0;
1417 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
1418 tot_charge += egs.getBulkMolesOld()[i] * mgd.basis_species_charge[i];
1419 EXPECT_NEAR(tot_charge, 0.0, 1E-15);
1420
1421 // surface charges
1422 for (unsigned sp = 0; sp < egs.getNumSurfacePotentials(); ++sp)
1423 {
1424 Real charge = 0.0;
1425 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
1426 if (mgd.surface_sorption_related[j] && mgd.surface_sorption_number[j] == sp)
1427 charge += mgd.eqm_species_charge[j] * egs.getEquilibriumMolality(j);
1428 EXPECT_NEAR(charge,
1429 egs.getSorbingSurfaceArea()[sp] * egs.getSurfaceCharge(sp) /
1430 GeochemistryConstants::FARADAY,
1431 1E-15);
1432 }
1433
1434 // check basis activity = activity_coefficient * molality
1435 for (unsigned i = 1; i < egs.getNumInBasis(); ++i) // don't loop over water
1436 if (mgd.basis_species_gas[i] || mgd.basis_species_mineral[i] || egs.getBasisActivityKnown()[i])
1437 continue;
1438 else
1439 EXPECT_NEAR(egs.getBasisActivity(i),
1440 egs.getBasisActivityCoefficient(i) *
1441 egs.getSolventMassAndFreeMolalityAndMineralMoles()[i],
1442 1E-15);
1443
1444 // check residuals are zero
1445 for (unsigned a = 0; a < egs.getNumInAlgebraicSystem(); ++a)
1446 EXPECT_LE(std::abs(egs.getResidualComponent(a, mole_additions)), 1E-15);
1447 // check residuals are zero by summing molalities, bulk compositions, etc
1448 Real nw = egs.getSolventWaterMass();
1449 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
1450 {
1451 Real res = -egs.getBulkMolesOld()[i];
1452 if (i == 0)
1453 res += nw * GeochemistryConstants::MOLES_PER_KG_WATER;
1454 else if (mgd.basis_species_mineral[i])
1455 res += egs.getSolventMassAndFreeMolalityAndMineralMoles()[i];
1456 else if (mgd.basis_species_gas[i])
1457 res += 0.0;
1458 else
1459 res += nw * egs.getSolventMassAndFreeMolalityAndMineralMoles()[i];
1460 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
1461 res += nw * mgd.eqm_stoichiometry(j, i) * egs.getEquilibriumMolality(j);
1462 EXPECT_LE(std::abs(res), 1E-14);
1463 }
1464 // surface potentials
1465 const Real prefactor = std::sqrt(GeochemistryConstants::GAS_CONSTANT *
1466 (25.0 + GeochemistryConstants::CELSIUS_TO_KELVIN) *
1467 GeochemistryConstants::PERMITTIVITY_FREE_SPACE *
1468 GeochemistryConstants::DIELECTRIC_CONSTANT_WATER *
1469 GeochemistryConstants::DENSITY_WATER * egs.getIonicStrength()) /
1470 nw;
1471 for (unsigned sp = 0; sp < egs.getNumSurfacePotentials(); ++sp)
1472 {
1473 EXPECT_NEAR(prefactor * std::sinh(egs.getSurfacePotential(sp) * GeochemistryConstants::FARADAY /
1474 2.0 / GeochemistryConstants::GAS_CONSTANT /
1475 (25.0 + GeochemistryConstants::CELSIUS_TO_KELVIN)),
1476 egs.getSurfaceCharge(sp),
1477 1.0E-13);
1478 }
1479
1480 // check equilibrium mass balance
1481 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
1482 {
1483 if (mgd.eqm_species_mineral[j] || mgd.eqm_species_gas[j])
1484 continue;
1485 Real log10ap = 0.0;
1486 for (unsigned i = 0; i < mgd.basis_species_name.size(); ++i)
1487 log10ap += mgd.eqm_stoichiometry(j, i) * std::log10(egs.getBasisActivity(i));
1488 log10ap -= egs.getLog10K(j);
1489 if (mgd.surface_sorption_related[j])
1490 log10ap -= std::log10(egs.getSurfacePotential(mgd.surface_sorption_number[j]));
1491 else
1492 log10ap -= std::log10(egs.getEquilibriumActivityCoefficient(j));
1493 if (log10ap < -300.0)
1494 EXPECT_LE(std::abs(egs.getEquilibriumMolality(j)), 1E-25);
1495 else
1496 EXPECT_LE(std::abs(egs.getEquilibriumMolality(j) - std::pow(10.0, log10ap)), 1E-15);
1497 }
1498
1499 // check equilibrium species have negative saturation indices, except for prevent_precipitation
1500 // minerals
1501 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
1502 {
1503 if (!mgd.eqm_species_mineral[j])
1504 continue;
1505 Real log10ap = 0.0;
1506 for (unsigned i = 0; i < mgd.basis_species_name.size(); ++i)
1507 log10ap += mgd.eqm_stoichiometry(j, i) * std::log10(egs.getBasisActivity(i));
1508 EXPECT_LE(log10ap, egs.getLog10K(j));
1509 }
1510}
1511
1513TEST(GeochemicalSolverTest, solve4_restore)
1514{
1515 // build the model
1516 const PertinentGeochemicalSystem model(db_full,
1517 {"H2O",
1518 "H+",
1519 "Cl-",
1520 "O2(aq)",
1521 "HCO3-",
1522 "Ca++",
1523 "Mg++",
1524 "Na+",
1525 "K+",
1526 "Fe++",
1527 "Fe+++",
1528 "Mn++",
1529 "Zn++",
1530 "SO4--"},
1531 {},
1532 {},
1533 {},
1534 {},
1535 {},
1536 "O2(aq)",
1537 "e-");
1538 ModelGeochemicalDatabase mgd = model.modelGeochemicalDatabase();
1539
1540 // build the equilibrium system
1541 GeochemistrySpeciesSwapper swapper(14, 1E-6);
1542 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
1543 GeochemicalSystem::ConstraintUserMeaningEnum::KG_SOLVENT_WATER,
1544 GeochemicalSystem::ConstraintUserMeaningEnum::ACTIVITY,
1545 GeochemicalSystem::ConstraintUserMeaningEnum::FREE_CONCENTRATION,
1546 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1547 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1548 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1549 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1550 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1551 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1552 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1553 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1554 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1555 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1556 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION};
1557 const std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu = {
1558 GeochemistryUnitConverter::GeochemistryUnit::KG,
1559 GeochemistryUnitConverter::GeochemistryUnit::DIMENSIONLESS,
1560 GeochemistryUnitConverter::GeochemistryUnit::MOLAL,
1561 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1562 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1563 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1564 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1565 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1566 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1567 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1568 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1569 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1570 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1571 GeochemistryUnitConverter::GeochemistryUnit::MOLES};
1572 GeochemicalSystem egs(mgd,
1573 ac_solver,
1574 is_solver,
1575 swapper,
1576 {},
1577 {},
1578 "Cl-",
1579 {"H2O",
1580 "H+",
1581 "O2(aq)",
1582 "Cl-",
1583 "HCO3-",
1584 "Ca++",
1585 "Mg++",
1586 "Na+",
1587 "K+",
1588 "Fe++",
1589 "Fe+++",
1590 "Mn++",
1591 "Zn++",
1592 "SO4--"},
1593 {1.0,
1594 0.8913E-6,
1595 0.13438E-3,
1596 3.041E-5,
1597 0.0295E-3,
1598 0.005938E-3,
1599 0.01448E-3,
1600 0.0018704E-3,
1601 0.005115E-3,
1602 0.012534E-3,
1603 0.0005372E-3,
1604 0.005042E-3,
1605 0.001897E-3,
1606 0.01562E-3},
1607 cu,
1608 cm,
1609 25,
1610 0,
1611 1E-20,
1612 {},
1613 {},
1614 {});
1615
1616 // build solver
1617 GeochemicalSolver solver(mgd.basis_species_name.size(),
1618 mgd.kin_species_name.size(),
1619 is_solver,
1620 1.0E-15,
1621 1.0E-200,
1622 100,
1623 1E-2,
1624 0.1,
1625 1,
1626 {},
1627 3.0,
1628 10,
1629 false);
1630
1631 // solve
1632 std::stringstream ss;
1633 unsigned tot_iter;
1634 Real abs_residual;
1635 DenseVector<Real> mole_additions(egs.getNumInBasis() + egs.getNumKinetic());
1636 DenseMatrix<Real> dmole_additions(egs.getNumInBasis() + egs.getNumKinetic(),
1637 egs.getNumInBasis() + egs.getNumKinetic());
1638 solver.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
1639 const Real old_residual = abs_residual;
1640
1641 // check converged
1642 EXPECT_LE(tot_iter, (unsigned)100);
1643 EXPECT_LE(abs_residual, 1.0E-15);
1644
1645 // retrieve the molalities, and set them into egs: this should result in no change if the
1646 // "restore" is working correctly
1647 std::vector<std::string> names = mgd.basis_species_name;
1648 names.insert(names.end(), mgd.eqm_species_name.begin(), mgd.eqm_species_name.end());
1649 std::vector<Real> molal = egs.getSolventMassAndFreeMolalityAndMineralMoles();
1650 for (unsigned j = 0; j < egs.getNumInEquilibrium(); ++j)
1651 molal.push_back(egs.getEquilibriumMolality(j));
1652 // form constraints_from_molalities, which are all false
1653 const std::vector<bool> com_false(egs.getNumInBasis(), false);
1654 egs.setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(
1655 names, molal, com_false);
1656
1657 // build solver with no ramping of the maximum ionic strength
1658 GeochemicalSolver solver0(mgd.basis_species_name.size(),
1659 mgd.kin_species_name.size(),
1660 is_solver,
1661 1.0E-15,
1662 1.0E-200,
1663 100,
1664 1E-2,
1665 0.1,
1666 1,
1667 {},
1668 3.0,
1669 0,
1670 false);
1671
1672 solver0.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
1673
1674 // check that the soler thinks this is truly a solution and the residual has not changed (up to
1675 // precision-loss)
1676 EXPECT_EQ(tot_iter, (unsigned)0);
1677 EXPECT_LE(abs_residual, old_residual);
1678
1679 // Now use constraints_from_molalities = true
1680 const std::vector<bool> com_true(egs.getNumInBasis(), true);
1681 egs.setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(
1682 names, molal, com_true);
1683
1684 solver0.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
1685
1686 // check that the soler thinks this is truly a solution and the residual has not increased
1687 EXPECT_EQ(tot_iter, (unsigned)0);
1688 EXPECT_LE(abs_residual, old_residual);
1689}
1690
1692TEST(GeochemicalSolverTest, solve5)
1693{
1694 // build the model
1695 const PertinentGeochemicalSystem model(
1696 db_ferric,
1697 {"H2O", "H+", "Na+", "Cl-", "Hg++", "Pb++", "SO4--", "Fe+++", ">(s)FeOH", ">(w)FeOH"},
1698 {"Fe(OH)3(ppd)"},
1699 {},
1700 {},
1701 {},
1702 {},
1703 "O2(aq)",
1704 "e-");
1705 ModelGeochemicalDatabase mgd = model.modelGeochemicalDatabase();
1706
1707 // build the equilibrium system
1708 GeochemistrySpeciesSwapper swapper(10, 1E-6);
1709 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
1710 GeochemicalSystem::ConstraintUserMeaningEnum::KG_SOLVENT_WATER,
1711 GeochemicalSystem::ConstraintUserMeaningEnum::ACTIVITY,
1712 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1713 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1714 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1715 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1716 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1717 GeochemicalSystem::ConstraintUserMeaningEnum::FREE_MINERAL,
1718 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1719 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION};
1720 const std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu = {
1721 GeochemistryUnitConverter::GeochemistryUnit::KG,
1722 GeochemistryUnitConverter::GeochemistryUnit::DIMENSIONLESS,
1723 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1724 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1725 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1726 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1727 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1728 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1729 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1730 GeochemistryUnitConverter::GeochemistryUnit::MOLES};
1731 GeochemicalSystem egs(
1732 mgd,
1733 ac_solver,
1734 is_solver,
1735 swapper,
1736 {"Fe+++"},
1737 {"Fe(OH)3(ppd)"},
1738 "Cl-",
1739 {"H2O", "H+", "Na+", "Cl-", "Hg++", "Pb++", "SO4--", "Fe(OH)3(ppd)", ">(s)FeOH", ">(w)FeOH"},
1740 {1.0, 1.0E-4, 10E-3, 10E-3, 0.1E-3, 0.1E-3, 0.2E-3, 9.3573E-3, 4.6786E-5, 1.87145E-3},
1741 cu,
1742 cm,
1743 25,
1744 0,
1745 1E-20,
1746 {},
1747 {},
1748 {});
1749
1750 // build solver
1751 GeochemicalSolver solver(mgd.basis_species_name.size(),
1752 mgd.kin_species_name.size(),
1753 is_solver,
1754 1.0E-15,
1755 1.0E-200,
1756 100,
1757 1.0,
1758 0.1,
1759 1,
1760 {},
1761 3.0,
1762 10,
1763 false);
1764
1765 // solve
1766 std::stringstream ss;
1767 unsigned tot_iter;
1768 Real abs_residual;
1769 DenseVector<Real> mole_additions(egs.getNumInBasis() + egs.getNumKinetic());
1770 DenseMatrix<Real> dmole_additions(egs.getNumInBasis() + egs.getNumKinetic(),
1771 egs.getNumInBasis() + egs.getNumKinetic());
1772 solver.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
1773
1774 // check converged
1775 EXPECT_LE(tot_iter, (unsigned)100);
1776 EXPECT_LE(abs_residual, 1.0E-15);
1777
1778 // check number in basis, number in redox disequilibrium and number of surface potentials
1779 EXPECT_EQ(egs.getNumInBasis(), (unsigned)10);
1780 EXPECT_EQ(egs.getNumRedox(), (unsigned)0);
1781 EXPECT_EQ(egs.getNumSurfacePotentials(), (unsigned)1);
1782 EXPECT_EQ(egs.getNumInEquilibrium(), mgd.eqm_species_name.size());
1783
1784 // check that the constraints are satisfied
1785 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
1786 if (mgd.basis_species_name[i] == "H2O")
1787 {
1788 EXPECT_FALSE(mgd.basis_species_gas[i]);
1789 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1790 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1791 EXPECT_NEAR(egs.getSolventWaterMass(), 1.0, 1.0E-15);
1792 EXPECT_NEAR(egs.getSolventMassAndFreeMolalityAndMineralMoles()[0], 1.0, 1.0E-15);
1793 }
1794 else if (mgd.basis_species_name[i] == "H+")
1795 {
1796 EXPECT_FALSE(mgd.basis_species_gas[i]);
1797 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1798 EXPECT_TRUE(egs.getBasisActivityKnown()[i]);
1799 EXPECT_NEAR(egs.getBasisActivity(i), 1.0E-4, 1.0E-15);
1800 }
1801 else if (mgd.basis_species_name[i] == "Na+")
1802 {
1803 EXPECT_FALSE(mgd.basis_species_gas[i]);
1804 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1805 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1806 EXPECT_NEAR(egs.getBulkMolesOld()[i], 10E-3, 1.0E-15);
1807 }
1808 else if (mgd.basis_species_name[i] == "Cl-")
1809 {
1810 EXPECT_FALSE(mgd.basis_species_gas[i]);
1811 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1812 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1813 // do not know the bulk composition as it is dictated by charge neutrality
1814 EXPECT_EQ(egs.getChargeBalanceBasisIndex(), i);
1815 }
1816 else if (mgd.basis_species_name[i] == "Hg++")
1817 {
1818 EXPECT_FALSE(mgd.basis_species_gas[i]);
1819 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1820 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1821 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.1E-3, 1.0E-15);
1822 }
1823 else if (mgd.basis_species_name[i] == "Pb++")
1824 {
1825 EXPECT_FALSE(mgd.basis_species_gas[i]);
1826 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1827 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1828 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.1E-3, 1.0E-15);
1829 }
1830 else if (mgd.basis_species_name[i] == "SO4--")
1831 {
1832 EXPECT_FALSE(mgd.basis_species_gas[i]);
1833 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1834 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1835 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.2E-3, 1.0E-15);
1836 }
1837 else if (mgd.basis_species_name[i] == "Fe(OH)3(ppd)")
1838 {
1839 EXPECT_FALSE(mgd.basis_species_gas[i]);
1840 EXPECT_TRUE(mgd.basis_species_mineral[i]);
1841 EXPECT_TRUE(egs.getBasisActivityKnown()[i]);
1842 EXPECT_NEAR(egs.getSolventMassAndFreeMolalityAndMineralMoles()[i], 9.3573E-3, 1.0E-15);
1843 }
1844 else if (mgd.basis_species_name[i] == ">(s)FeOH")
1845 {
1846 EXPECT_FALSE(mgd.basis_species_gas[i]);
1847 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1848 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1849 EXPECT_NEAR(egs.getBulkMolesOld()[i], 4.6786E-5, 1.0E-15);
1850 }
1851 else if (mgd.basis_species_name[i] == ">(w)FeOH")
1852 {
1853 EXPECT_FALSE(mgd.basis_species_gas[i]);
1854 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1855 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1856 EXPECT_NEAR(egs.getBulkMolesOld()[i], 1.87145E-3, 1.0E-15);
1857 }
1858 else
1859 FAIL() << "Incorrect basis species";
1860
1861 // check total charge
1862 EXPECT_NEAR(egs.getTotalChargeOld(), 0.0, 1E-15);
1863 // check total charge by summing up basis charges
1864 Real tot_charge = 0.0;
1865 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
1866 tot_charge += egs.getBulkMolesOld()[i] * mgd.basis_species_charge[i];
1867 EXPECT_NEAR(tot_charge, 0.0, 1E-15);
1868
1869 // surface charges
1870 for (unsigned sp = 0; sp < egs.getNumSurfacePotentials(); ++sp)
1871 {
1872 Real charge = 0.0;
1873 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
1874 if (mgd.surface_sorption_related[j] && mgd.surface_sorption_number[j] == sp)
1875 charge += mgd.eqm_species_charge[j] * egs.getEquilibriumMolality(j);
1876 EXPECT_NEAR(charge,
1877 egs.getSorbingSurfaceArea()[sp] * egs.getSurfaceCharge(sp) /
1878 GeochemistryConstants::FARADAY,
1879 1E-15);
1880 }
1881
1882 // check basis activity = activity_coefficient * molality
1883 for (unsigned i = 1; i < egs.getNumInBasis(); ++i) // don't loop over water
1884 if (mgd.basis_species_gas[i] || mgd.basis_species_mineral[i] || egs.getBasisActivityKnown()[i])
1885 continue;
1886 else
1887 EXPECT_NEAR(egs.getBasisActivity(i),
1888 egs.getBasisActivityCoefficient(i) *
1889 egs.getSolventMassAndFreeMolalityAndMineralMoles()[i],
1890 1E-15);
1891
1892 // check residuals are zero
1893 for (unsigned a = 0; a < egs.getNumInAlgebraicSystem(); ++a)
1894 EXPECT_LE(std::abs(egs.getResidualComponent(a, mole_additions)), 1E-15);
1895 // check residuals are zero by summing molalities, bulk compositions, etc
1896 Real nw = egs.getSolventWaterMass();
1897 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
1898 {
1899 Real res = -egs.getBulkMolesOld()[i];
1900 if (i == 0)
1901 res += nw * GeochemistryConstants::MOLES_PER_KG_WATER;
1902 else if (mgd.basis_species_mineral[i])
1903 res += egs.getSolventMassAndFreeMolalityAndMineralMoles()[i];
1904 else if (mgd.basis_species_gas[i])
1905 res += 0.0;
1906 else
1907 res += nw * egs.getSolventMassAndFreeMolalityAndMineralMoles()[i];
1908 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
1909 res += nw * mgd.eqm_stoichiometry(j, i) * egs.getEquilibriumMolality(j);
1910 EXPECT_LE(std::abs(res), 1E-14);
1911 }
1912 // surface potentials
1913 const Real prefactor = std::sqrt(8.0 * GeochemistryConstants::GAS_CONSTANT *
1914 (25.0 + GeochemistryConstants::CELSIUS_TO_KELVIN) *
1915 GeochemistryConstants::PERMITTIVITY_FREE_SPACE *
1916 GeochemistryConstants::DIELECTRIC_CONSTANT_WATER *
1917 GeochemistryConstants::DENSITY_WATER * egs.getIonicStrength()) /
1918 nw;
1919 for (unsigned sp = 0; sp < egs.getNumSurfacePotentials(); ++sp)
1920 {
1921 EXPECT_NEAR(prefactor * std::sinh(egs.getSurfacePotential(sp) * GeochemistryConstants::FARADAY /
1922 2.0 / GeochemistryConstants::GAS_CONSTANT /
1923 (25.0 + GeochemistryConstants::CELSIUS_TO_KELVIN)),
1924 egs.getSurfaceCharge(sp),
1925 1.0E-13);
1926 }
1927
1928 // check equilibrium mass balance
1929 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
1930 {
1931 if (mgd.eqm_species_mineral[j] || mgd.eqm_species_gas[j])
1932 continue;
1933 Real log10ap = 0.0;
1934 for (unsigned i = 0; i < mgd.basis_species_name.size(); ++i)
1935 log10ap += mgd.eqm_stoichiometry(j, i) * std::log10(egs.getBasisActivity(i));
1936 log10ap -= egs.getLog10K(j);
1937 if (mgd.surface_sorption_related[j])
1938 log10ap -= std::log10(std::exp(mgd.eqm_species_charge[j] * GeochemistryConstants::FARADAY *
1939 egs.getSurfacePotential(mgd.surface_sorption_number[j]) /
1940 GeochemistryConstants::GAS_CONSTANT /
1941 (25.0 + GeochemistryConstants::CELSIUS_TO_KELVIN)));
1942 else
1943 log10ap -= std::log10(egs.getEquilibriumActivityCoefficient(j));
1944 if (log10ap < -300.0)
1945 EXPECT_LE(std::abs(egs.getEquilibriumMolality(j)), 1E-25);
1946 else
1947 EXPECT_LE(std::abs(egs.getEquilibriumMolality(j) - std::pow(10.0, log10ap)), 1E-15);
1948 }
1949
1950 // check equilibrium species have negative saturation indices, except for prevent_precipitation
1951 // minerals
1952 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
1953 {
1954 if (!mgd.eqm_species_mineral[j])
1955 continue;
1956 Real log10ap = 0.0;
1957 for (unsigned i = 0; i < mgd.basis_species_name.size(); ++i)
1958 log10ap += mgd.eqm_stoichiometry(j, i) * std::log10(egs.getBasisActivity(i));
1959 EXPECT_LE(log10ap, egs.getLog10K(j));
1960 }
1961}
1962
1964TEST(GeochemicalSolverTest, solve5_restore)
1965{
1966 // build the model
1967 const PertinentGeochemicalSystem model(
1968 db_ferric,
1969 {"H2O", "H+", "Na+", "Cl-", "Hg++", "Pb++", "SO4--", "Fe+++", ">(s)FeOH", ">(w)FeOH"},
1970 {"Fe(OH)3(ppd)"},
1971 {},
1972 {},
1973 {},
1974 {},
1975 "O2(aq)",
1976 "e-");
1977 ModelGeochemicalDatabase mgd = model.modelGeochemicalDatabase();
1978
1979 // build the equilibrium system
1980 GeochemistrySpeciesSwapper swapper(10, 1E-6);
1981 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
1982 GeochemicalSystem::ConstraintUserMeaningEnum::KG_SOLVENT_WATER,
1983 GeochemicalSystem::ConstraintUserMeaningEnum::ACTIVITY,
1984 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1985 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1986 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1987 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1988 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1989 GeochemicalSystem::ConstraintUserMeaningEnum::FREE_MINERAL,
1990 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1991 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION};
1992 const std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu = {
1993 GeochemistryUnitConverter::GeochemistryUnit::KG,
1994 GeochemistryUnitConverter::GeochemistryUnit::DIMENSIONLESS,
1995 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1996 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1997 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1998 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1999 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
2000 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
2001 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
2002 GeochemistryUnitConverter::GeochemistryUnit::MOLES};
2003 GeochemicalSystem egs(
2004 mgd,
2005 ac_solver,
2006 is_solver,
2007 swapper,
2008 {"Fe+++"},
2009 {"Fe(OH)3(ppd)"},
2010 "Cl-",
2011 {"H2O", "H+", "Na+", "Cl-", "Hg++", "Pb++", "SO4--", "Fe(OH)3(ppd)", ">(s)FeOH", ">(w)FeOH"},
2012 {1.0, 1.0E-4, 10E-3, 10E-3, 0.1E-3, 0.1E-3, 0.2E-3, 9.3573E-3, 4.6786E-5, 1.87145E-3},
2013 cu,
2014 cm,
2015 25,
2016 0,
2017 1E-20,
2018 {},
2019 {},
2020 {});
2021
2022 // build solver
2023 GeochemicalSolver solver(mgd.basis_species_name.size(),
2024 mgd.kin_species_name.size(),
2025 is_solver,
2026 1.0E-15,
2027 1.0E-200,
2028 100,
2029 1.0,
2030 0.1,
2031 1,
2032 {},
2033 3.0,
2034 0,
2035 false);
2036
2037 // solve
2038 std::stringstream ss;
2039 unsigned tot_iter;
2040 Real abs_residual;
2041 DenseVector<Real> mole_additions(egs.getNumInBasis() + egs.getNumKinetic());
2042 DenseMatrix<Real> dmole_additions(egs.getNumInBasis() + egs.getNumKinetic(),
2043 egs.getNumInBasis() + egs.getNumKinetic());
2044 solver.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
2045 const Real old_residual = abs_residual;
2046
2047 // check converged
2048 EXPECT_LE(tot_iter, (unsigned)100);
2049 EXPECT_LE(abs_residual, 1.0E-15);
2050
2051 // now setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles to the solution
2052 // obtained This should actually do nothing: we are just checking the "restore" doesn't muck
2053 // something up! Of course, activity coefficients will slightly change due to making the
2054 // configuration slightly more consistent (activities and molalities match) so there will be very
2055 // minor changes
2056 std::vector<std::string> names = mgd.basis_species_name;
2057 names.insert(names.end(), mgd.eqm_species_name.begin(), mgd.eqm_species_name.end());
2058 names.push_back("Fe(OH)3(ppd)_surface_potential_expr");
2059 std::vector<Real> molal = egs.getSolventMassAndFreeMolalityAndMineralMoles();
2060 for (unsigned j = 0; j < egs.getNumInEquilibrium(); ++j)
2061 molal.push_back(egs.getEquilibriumMolality(j));
2062 molal.push_back(std::exp(egs.getSurfacePotential(0) * GeochemistryConstants::FARADAY /
2063 (-2.0 * GeochemistryConstants::GAS_CONSTANT *
2064 (25.0 + GeochemistryConstants::CELSIUS_TO_KELVIN))));
2065 // form constraints_from_molalities, which are all false
2066 const std::vector<bool> com_false(egs.getNumInBasis(), false);
2067 egs.setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(
2068 names, molal, com_false);
2069
2070 solver.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
2071 EXPECT_EQ(tot_iter, (unsigned)0);
2072 EXPECT_LE(abs_residual, 2.0 * old_residual);
2073
2074 // Now use constraints_from_molalities = true
2075 const std::vector<bool> com_true(egs.getNumInBasis(), true);
2076 egs.setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(
2077 names, molal, com_true);
2078
2079 solver.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
2080
2081 EXPECT_EQ(tot_iter, (unsigned)0);
2082 EXPECT_LE(abs_residual, 2.0 * old_residual);
2083
2084 // now check the copy-assignment operator of GeochemicalSystem
2085
2086 // build the destination model
2087 const PertinentGeochemicalSystem dest_model(
2088 db_ferric,
2089 {"H2O", "H+", "Na+", "Cl-", "Hg++", "Pb++", "SO4--", "Fe+++", ">(s)FeOH", ">(w)FeOH"},
2090 {"Fe(OH)3(ppd)"},
2091 {},
2092 {},
2093 {},
2094 {},
2095 "O2(aq)",
2096 "e-");
2097 ModelGeochemicalDatabase dest_mgd = dest_model.modelGeochemicalDatabase();
2098
2099 // build the destination equilibrium system
2100 GeochemicalSystem dest_egs(
2101 dest_mgd,
2102 ac_solver,
2103 is_solver,
2104 swapper,
2105 {"Fe+++"},
2106 {"Fe(OH)3(ppd)"},
2107 "Cl-",
2108 {"H2O", "H+", "Na+", "Cl-", "Hg++", "Pb++", "SO4--", "Fe(OH)3(ppd)", ">(s)FeOH", ">(w)FeOH"},
2109 {2.0, 2.0E-4, 20E-3, 20E-3, 0.2E-3, 0.2E-3, 0.4E-3, 19.3573E-3, 14.6786E-5, 2.87145E-3},
2110 cu,
2111 cm,
2112 250,
2113 0,
2114 1E-20,
2115 {},
2116 {},
2117 {});
2118
2119 // change constraint meanings in egs to provide a more thorough check of the copy assignment
2120 egs.closeSystem();
2121 // copy assignment
2122 dest_egs = egs;
2123 EXPECT_EQ(dest_egs, egs);
2124}
2125
2127TEST(GeochemicalSolverTest, solve_addH)
2128{
2129 ModelGeochemicalDatabase mgd = model_simplest.modelGeochemicalDatabase();
2130 GeochemicalSystem egs(mgd,
2131 ac_solver,
2132 is_solver,
2133 swapper2,
2134 {},
2135 {},
2136 "H+",
2137 {"H2O", "H+"},
2138 {1.75, 1},
2139 cu2,
2140 cm2,
2141 25,
2142 0,
2143 1E-20,
2144 {},
2145 {},
2146 {});
2147 GeochemicalSolver solver(mgd.basis_species_name.size(),
2148 mgd.kin_species_name.size(),
2149 is_solver,
2150 1.0E-12,
2151 1.0E-200,
2152 100,
2153 10.0,
2154 0.1,
2155 1,
2156 {},
2157 3.0,
2158 0,
2159 false);
2160
2161 std::stringstream ss;
2162 unsigned tot_iter;
2163 Real abs_residual;
2164 DenseVector<Real> mole_additions(egs.getNumInBasis() + egs.getNumKinetic());
2165 DenseMatrix<Real> dmole_additions(egs.getNumInBasis() + egs.getNumKinetic(),
2166 egs.getNumInBasis() + egs.getNumKinetic());
2167 solver.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
2168
2169 // check Newton has converged
2170 EXPECT_LE(tot_iter, (unsigned)100);
2171 EXPECT_LE(abs_residual, 1.0E-12);
2172
2173 // check Bulk is as expected
2174 const std::vector<Real> bulk = egs.getBulkMolesOld();
2175 EXPECT_NEAR(bulk[0],
2176 1.75 * (GeochemistryConstants::MOLES_PER_KG_WATER + egs.getEquilibriumMolality(0)),
2177 1.0E-9);
2178 EXPECT_NEAR(bulk[1], 0.0, 1.0E-9);
2179
2180 // close the system by setting to fixed number of moles
2181 egs.changeConstraintToBulk(0);
2182
2183 // check Bulk is as expected
2184 std::vector<Real> bulk_new = egs.getBulkMolesOld();
2185 for (unsigned i = 0; i < 2; ++i)
2186 EXPECT_EQ(bulk_new[i], bulk[i]);
2187
2188 // add 1 mol of water
2189 mole_additions(0) = 1.0;
2190 solver.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
2191
2192 // check Newton has converged
2193 EXPECT_LE(tot_iter, (unsigned)100);
2194 EXPECT_LE(abs_residual, 1.0E-12);
2195
2196 // check Bulk is as expected
2197 bulk_new = egs.getBulkMolesOld();
2198 EXPECT_EQ(bulk_new[0], bulk[0] + 1.0);
2199 EXPECT_EQ(bulk_new[1], bulk[1]);
2200
2201 // add 1 mol of H+ : this shouldn't make any difference because charge-balance instantly removes
2202 // it!
2203 mole_additions(0) = 0.0;
2204 mole_additions(1) = 1.0;
2205 solver.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
2206
2207 // check Newton has converged
2208 EXPECT_LE(tot_iter, (unsigned)100);
2209 EXPECT_LE(abs_residual, 1.0E-12);
2210
2211 // check Bulk is as expected
2212 bulk_new = egs.getBulkMolesOld();
2213 EXPECT_EQ(bulk_new[0], bulk[0] + 1.0);
2214 EXPECT_EQ(bulk_new[1], bulk[1]);
2215}
2216
2218TEST(GeochemicalSolverTest, maxSwapsException)
2219{
2220 // build the model
2221 const PertinentGeochemicalSystem model(
2222 db_full,
2223 {"H2O", "H+", "Cl-", "Na+", "SO4--", "Mg++", "Ca++", "K+", "HCO3-", "SiO2(aq)", "O2(aq)"},
2224 {"Antigorite",
2225 "Tremolite",
2226 "Talc",
2227 "Chrysotile",
2228 "Sepiolite",
2229 "Anthophyllite",
2230 "Dolomite",
2231 "Dolomite-ord",
2232 "Huntite",
2233 "Dolomite-dis",
2234 "Magnesite",
2235 "Calcite",
2236 "Aragonite",
2237 "Quartz"},
2238 {},
2239 {},
2240 {},
2241 {},
2242 "O2(aq)",
2243 "e-");
2244 ModelGeochemicalDatabase mgd = model.modelGeochemicalDatabase();
2245
2246 // build the equilibrium system
2247 GeochemistrySpeciesSwapper swapper(11, 1E-6);
2248 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
2249 GeochemicalSystem::ConstraintUserMeaningEnum::KG_SOLVENT_WATER,
2250 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
2251 GeochemicalSystem::ConstraintUserMeaningEnum::FREE_CONCENTRATION,
2252 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
2253 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
2254 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
2255 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
2256 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
2257 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
2258 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
2259 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION};
2260 const std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu = {
2261 GeochemistryUnitConverter::GeochemistryUnit::KG,
2262 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
2263 GeochemistryUnitConverter::GeochemistryUnit::MOLAL,
2264 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
2265 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
2266 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
2267 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
2268 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
2269 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
2270 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
2271 GeochemistryUnitConverter::GeochemistryUnit::MOLES};
2272 GeochemicalSystem egs(
2273 mgd,
2274 ac_solver,
2275 is_solver,
2276 swapper,
2277 {"H+"},
2278 {"MgCO3"},
2279 "Cl-",
2280 {"H2O", "MgCO3", "O2(aq)", "Cl-", "Na+", "SO4--", "Mg++", "Ca++", "K+", "HCO3-", "SiO2(aq)"},
2281 {1.0,
2282 0.196E-3,
2283 0.2151E-3,
2284 0.5656,
2285 0.4850,
2286 0.02924,
2287 0.05501,
2288 0.01063,
2289 0.010576055,
2290 0.002412,
2291 0.00010349},
2292 cu,
2293 cm,
2294 25,
2295 0,
2296 1E-20,
2297 {},
2298 {},
2299 {});
2300
2301 // build solver (4 swaps are needed to solve this system)
2302 GeochemicalSolver solver(mgd.basis_species_name.size(),
2303 mgd.kin_species_name.size(),
2304 is_solver,
2305 1.0E-15,
2306 1.0E-200,
2307 100,
2308 1E3,
2309 0.1,
2310 3,
2311 {"Dolomite-ord", "Dolomite-dis"},
2312 3.0,
2313 10,
2314 false);
2315
2316 // solve
2317 std::stringstream ss;
2318 unsigned tot_iter;
2319 Real abs_residual;
2320 try
2321 {
2322 DenseVector<Real> mole_additions(egs.getNumInBasis() + egs.getNumKinetic());
2323 DenseMatrix<Real> dmole_additions(egs.getNumInBasis() + egs.getNumKinetic(),
2324 egs.getNumInBasis() + egs.getNumKinetic());
2325 solver.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
2326 FAIL() << "Missing expected exception.";
2327 }
2328 catch (const std::exception & e)
2329 {
2330 std::string msg(e.what());
2331 ASSERT_TRUE(msg.find("Maximum number of swaps performed during solve") != std::string::npos)
2332 << "Failed with unexpected error message: " << msg;
2333 }
2334}
2335
2337TEST(GeochemicalSolverTest, setRampMaxIonicStrength)
2338{
2339 ModelGeochemicalDatabase mgd = model_simplest.modelGeochemicalDatabase();
2340 GeochemicalSystem egs(mgd,
2341 ac_solver,
2342 is_solver,
2343 swapper2,
2344 {},
2345 {},
2346 "H+",
2347 {"H2O", "H+"},
2348 {1.75, 3.0},
2349 cu2,
2350 cm2,
2351 25,
2352 0,
2353 1E-20,
2354 {},
2355 {},
2356 {});
2357 GeochemicalSolver solver(mgd.basis_species_name.size(),
2358 mgd.kin_species_name.size(),
2359 is_solver,
2360 1.0,
2361 0.1,
2362 100,
2363 1E100,
2364 0.1,
2365 1,
2366 {},
2367 3.0,
2368 10,
2369 true);
2370
2371 ASSERT_EQ(solver.getRampMaxIonicStrength(), (unsigned)10);
2372 try
2373 {
2374 solver.setRampMaxIonicStrength(101);
2375 FAIL() << "Missing expected exception.";
2376 }
2377 catch (const std::exception & e)
2378 {
2379 std::string msg(e.what());
2380 ASSERT_TRUE(msg.find("GeochemicalSolver: ramp_max_ionic_strength must be less than max_iter") !=
2381 std::string::npos)
2382 << "Failed with unexpected error message: " << msg;
2383 }
2384 solver.setRampMaxIonicStrength(21);
2385 ASSERT_EQ(solver.getRampMaxIonicStrength(), (unsigned)21);
2386}
2387
2389TEST(GeochemicalSolverTest, solve_kinetic1)
2390{
2391 const PertinentGeochemicalSystem model(db_solver,
2392 {"H2O", "H+", "Fe+++", "HCO3-"},
2393 {},
2394 {},
2395 {"Something", "Fe(OH)3(ppd)fake"},
2396 {},
2397 {},
2398 "O2(aq)",
2399 "e-");
2400 ModelGeochemicalDatabase mgd = model.modelGeochemicalDatabase();
2401 std::vector<GeochemistryUnitConverter::GeochemistryUnit> ku;
2402 ku.push_back(GeochemistryUnitConverter::GeochemistryUnit::MOLES);
2403 ku.push_back(GeochemistryUnitConverter::GeochemistryUnit::MOLES);
2404 GeochemicalSystem egs(mgd,
2405 ac_solver,
2406 is_solver,
2407 swapper_kin,
2408 {},
2409 {},
2410 "HCO3-",
2411 {"H2O", "H+", "Fe+++", "HCO3-"},
2412 {1.75, 1E-5, 1E-5, 4E-5},
2413 cu4,
2414 cm4,
2415 25,
2416 0,
2417 1E-20,
2418 {"Something", "Fe(OH)3(ppd)fake"},
2419 {1.0E-6, 2.0E-6},
2420 ku);
2421 GeochemicalSolver solver(mgd.basis_species_name.size(),
2422 mgd.kin_species_name.size(),
2423 is_solver,
2424 1.0E-15,
2425 1.0E-200,
2426 100,
2427 1E-5,
2428 1E-16,
2429 1,
2430 {},
2431 3.0,
2432 0,
2433 true);
2434
2435 std::stringstream ss;
2436 unsigned tot_iter;
2437 Real abs_residual;
2438 DenseVector<Real> mole_additions(egs.getNumInBasis() + egs.getNumKinetic());
2439 DenseMatrix<Real> dmole_additions(egs.getNumInBasis() + egs.getNumKinetic(),
2440 egs.getNumInBasis() + egs.getNumKinetic());
2441 solver.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
2442
2443 // check Newton has converged
2444 EXPECT_LE(tot_iter, (unsigned)100);
2445 EXPECT_LE(abs_residual, 1.0E-15);
2446
2447 // check numbers are correct
2448 EXPECT_EQ(egs.getNumInBasis(), (unsigned)4);
2449 EXPECT_EQ(egs.getNumRedox(), (unsigned)0);
2450 EXPECT_EQ(egs.getNumSurfacePotentials(), (unsigned)0);
2451 EXPECT_EQ(egs.getNumInEquilibrium(), mgd.eqm_species_name.size());
2452 EXPECT_EQ(egs.getNumKinetic(), (unsigned)2);
2453
2454 // check that kinetic moles have not changed (kinetic rates are zero)
2455 for (unsigned k = 0; k < egs.getNumKinetic(); ++k)
2456 if (mgd.kin_species_name[k] == "Something")
2457 EXPECT_EQ(egs.getKineticMoles(k), 1.0E-6);
2458 else if (mgd.kin_species_name[k] == "Fe(OH)3(ppd)fake")
2459 EXPECT_EQ(egs.getKineticMoles(k), 2.0E-6);
2460 else
2461 FAIL() << "Incorrect kinetic species";
2462
2463 // check that the constraints are satisfied
2464 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
2465 if (mgd.basis_species_name[i] == "H2O")
2466 {
2467 EXPECT_FALSE(mgd.basis_species_gas[i]);
2468 EXPECT_FALSE(mgd.basis_species_mineral[i]);
2469 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
2470 EXPECT_NEAR(egs.getSolventWaterMass(), 1.75, 1.0E-15);
2471 EXPECT_NEAR(egs.getSolventMassAndFreeMolalityAndMineralMoles()[0], 1.75, 1.0E-15);
2472 }
2473 else if (mgd.basis_species_name[i] == "H+")
2474 {
2475 EXPECT_FALSE(mgd.basis_species_gas[i]);
2476 EXPECT_FALSE(mgd.basis_species_mineral[i]);
2477 EXPECT_TRUE(egs.getBasisActivityKnown()[i]);
2478 EXPECT_NEAR(egs.getBasisActivity(i), 1E-5, 1.0E-15);
2479 }
2480 else if (mgd.basis_species_name[i] == "Fe+++")
2481 {
2482 EXPECT_FALSE(mgd.basis_species_gas[i]);
2483 EXPECT_FALSE(mgd.basis_species_mineral[i]);
2484 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
2485 EXPECT_NEAR(egs.getBulkMolesOld()[i], 1E-5 + 1E-6 + 2 * 2E-6, 1.0E-15);
2486 }
2487 else if (mgd.basis_species_name[i] == "HCO3-")
2488 {
2489 EXPECT_FALSE(mgd.basis_species_gas[i]);
2490 EXPECT_FALSE(mgd.basis_species_mineral[i]);
2491 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
2492 // do not know the bulk composition as it is dictated by charge neutrality
2493 EXPECT_EQ(egs.getChargeBalanceBasisIndex(), i);
2494 }
2495 else
2496 FAIL() << "Incorrect basis species";
2497
2498 // check total charge
2499 EXPECT_NEAR(egs.getTotalChargeOld(), 0.0, 1E-15);
2500 // check total charge by summing up basis charges
2501 Real tot_charge = 0.0;
2502 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
2503 tot_charge += egs.getBulkMolesOld()[i] * mgd.basis_species_charge[i];
2504 EXPECT_NEAR(tot_charge, 0.0, 1E-15);
2505
2506 // check basis activity = activity_coefficient * molality
2507 for (unsigned i = 1; i < egs.getNumInBasis(); ++i) // don't loop over water
2508 if (mgd.basis_species_gas[i] || mgd.basis_species_mineral[i] || egs.getBasisActivityKnown()[i])
2509 continue;
2510 else
2511 EXPECT_NEAR(egs.getBasisActivity(i),
2512 egs.getBasisActivityCoefficient(i) *
2513 egs.getSolventMassAndFreeMolalityAndMineralMoles()[i],
2514 1E-15);
2515
2516 // check residuals are zero
2517 for (unsigned a = 0; a < egs.getNumInAlgebraicSystem(); ++a)
2518 EXPECT_LE(std::abs(egs.getResidualComponent(a, mole_additions)), 1E-15);
2519 // check residuals are zero by summing molalities, bulk compositions, etc
2520 Real nw = egs.getSolventWaterMass();
2521 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
2522 {
2523 Real res = -egs.getBulkMolesOld()[i];
2524 if (i == 0)
2525 res += nw * GeochemistryConstants::MOLES_PER_KG_WATER;
2526 else if (mgd.basis_species_mineral[i])
2527 res += egs.getSolventMassAndFreeMolalityAndMineralMoles()[i];
2528 else if (mgd.basis_species_gas[i])
2529 res += 0.0;
2530 else
2531 res += nw * egs.getSolventMassAndFreeMolalityAndMineralMoles()[i];
2532 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
2533 res += nw * mgd.eqm_stoichiometry(j, i) * egs.getEquilibriumMolality(j);
2534 for (unsigned k = 0; k < mgd.kin_species_name.size(); ++k)
2535 res += mgd.kin_stoichiometry(k, i) *
2536 egs.getKineticMoles(
2537 k); // this should be the only important kinetic contribution to this test!
2538 EXPECT_LE(std::abs(res), 1E-13);
2539 }
2540
2541 // check equilibrium mass balance
2542 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
2543 {
2544 if (mgd.eqm_species_mineral[j] || mgd.eqm_species_gas[j])
2545 continue;
2546 Real log10ap = 0.0;
2547 for (unsigned i = 0; i < mgd.basis_species_name.size(); ++i)
2548 log10ap += mgd.eqm_stoichiometry(j, i) * std::log10(egs.getBasisActivity(i));
2549 log10ap -= egs.getLog10K(j);
2550 if (mgd.surface_sorption_related[j])
2551 log10ap -= std::log10(egs.getSurfacePotential(mgd.surface_sorption_number[j]));
2552 else
2553 log10ap -= std::log10(egs.getEquilibriumActivityCoefficient(j));
2554 if (log10ap < -300.0)
2555 EXPECT_LE(std::abs(egs.getEquilibriumMolality(j)), 1E-25);
2556 else
2557 EXPECT_LE(std::abs(egs.getEquilibriumMolality(j) - std::pow(10.0, log10ap)), 1E-15);
2558 }
2559}
2560
2562TEST(GeochemicalSolverTest, solve_kinetic2)
2563{
2564 PertinentGeochemicalSystem model(db_solver,
2565 {"H2O", "H+", "Fe+++", "HCO3-"},
2566 {},
2567 {},
2568 {"Something", "Fe(OH)3(ppd)fake"},
2569 {},
2570 {},
2571 "O2(aq)",
2572 "e-");
2573 // Define constant rates
2574 // Something = -3H+ + Fe+++ + 2H2O + 1.5HCO3-, so Q=1E15*1E-5*1E-7=1E3 (approx) >> K, so
2575 // rate_Something produces a negative rate, so mole_additions > 0, so Something will increase in
2576 // mole number
2577 KineticRateUserDescription rate_Something("Something",
2578 1.0E-7,
2579 1.0,
2580 false,
2581 0.0,
2582 0.0,
2583 0.0,
2584 {},
2585 {},
2586 {},
2587 {},
2588 1.0,
2589 0.0,
2590 0.0,
2591 0.0,
2592 DirectionChoiceEnum::BOTH,
2593 "H2O",
2594 0.0,
2595 -1.0,
2596 0.0);
2597 // Fe(OH)3(ppd)fake = -3H+ + 2Fe+++ + 3H2O, so Q=1E15*1E-10=1E5 (approx) < K, so rate_Fe produces
2598 // a positive rate, so mole_additions < 0, so Something will decrease in mole number
2599 KineticRateUserDescription rate_Fe("Fe(OH)3(ppd)fake",
2600 1.0E-7,
2601 1.0,
2602 false,
2603 0.0,
2604 0.0,
2605 0.0,
2606 {},
2607 {},
2608 {},
2609 {},
2610 1.0,
2611 0.0,
2612 0.0,
2613 0.0,
2614 DirectionChoiceEnum::BOTH,
2615 "H2O",
2616 0.0,
2617 -1.0,
2618 0.0);
2619 model.addKineticRate(rate_Something);
2620 model.addKineticRate(rate_Fe);
2621
2622 ModelGeochemicalDatabase mgd = model.modelGeochemicalDatabase();
2623 std::vector<GeochemistryUnitConverter::GeochemistryUnit> ku;
2624 ku.push_back(GeochemistryUnitConverter::GeochemistryUnit::MOLES);
2625 ku.push_back(GeochemistryUnitConverter::GeochemistryUnit::MOLES);
2626 GeochemicalSystem egs(mgd,
2627 ac_solver,
2628 is_solver,
2629 swapper_kin,
2630 {},
2631 {},
2632 "HCO3-",
2633 {"H2O", "H+", "Fe+++", "HCO3-"},
2634 {1.75, 1E-5, 1E-5, 4E-5},
2635 cu4,
2636 cm4,
2637 25,
2638 0,
2639 1E-20,
2640 {"Something", "Fe(OH)3(ppd)fake"},
2641 {1.0E-6, 2.0E-6},
2642 ku);
2643 GeochemicalSolver solver(mgd.basis_species_name.size(),
2644 mgd.kin_species_name.size(),
2645 is_solver,
2646 1.0E-15,
2647 1.0E-200,
2648 100,
2649 1E-5,
2650 1E-16,
2651 1,
2652 {},
2653 3.0,
2654 0,
2655 true);
2656
2657 const unsigned index_something = mgd.kin_species_index.at("Something");
2658 const unsigned index_fe = mgd.kin_species_index.at("Fe(OH)3(ppd)fake");
2659
2660 std::stringstream ss;
2661 unsigned tot_iter;
2662 Real abs_residual;
2663 DenseVector<Real> mole_additions(egs.getNumInBasis() + egs.getNumKinetic());
2664 mole_additions(egs.getNumInBasis() + index_something) = 1.5E-6; // external input to "Something"
2665 DenseMatrix<Real> dmole_additions(egs.getNumInBasis() + egs.getNumKinetic(),
2666 egs.getNumInBasis() + egs.getNumKinetic());
2667 // solve with a time-step size of 10.0
2668 solver.solveSystem(egs, ss, tot_iter, abs_residual, 10.0, mole_additions, dmole_additions);
2669
2670 // check Newton has converged
2671 EXPECT_LE(tot_iter, (unsigned)100);
2672 EXPECT_LE(abs_residual, 1.0E-15);
2673
2674 // check numbers are correct
2675 EXPECT_EQ(egs.getNumInBasis(), (unsigned)4);
2676 EXPECT_EQ(egs.getNumRedox(), (unsigned)0);
2677 EXPECT_EQ(egs.getNumSurfacePotentials(), (unsigned)0);
2678 EXPECT_EQ(egs.getNumInEquilibrium(), mgd.eqm_species_name.size());
2679 EXPECT_EQ(egs.getNumKinetic(), (unsigned)2);
2680
2681 // check activity products and log10K are ordered as calculated above
2682 EXPECT_GT(egs.log10KineticActivityProduct(index_something),
2683 egs.getKineticLog10K(index_something));
2684 EXPECT_LT(egs.log10KineticActivityProduct(index_fe), egs.getKineticLog10K(index_fe));
2685
2686 // check that mole additions is as calculated above
2687 EXPECT_EQ(mole_additions(egs.getNumInBasis() + index_something), 1.5E-6 + 1.0E-6);
2688 EXPECT_EQ(mole_additions(egs.getNumInBasis() + index_fe), -1.0E-6);
2689
2690 // check that kinetic moles have changed according to the constant rates
2691 for (unsigned k = 0; k < egs.getNumKinetic(); ++k)
2692 if (mgd.kin_species_name[k] == "Something")
2693 EXPECT_NEAR(egs.getKineticMoles(k), 2.0E-6 + 1.5E-6, 1.0E-12);
2694 else if (mgd.kin_species_name[k] == "Fe(OH)3(ppd)fake")
2695 EXPECT_EQ(egs.getKineticMoles(k), 1.0E-6);
2696 else
2697 FAIL() << "Incorrect kinetic species";
2698
2699 // check that the constraints are satisfied
2700 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
2701 if (mgd.basis_species_name[i] == "H2O")
2702 {
2703 EXPECT_FALSE(mgd.basis_species_gas[i]);
2704 EXPECT_FALSE(mgd.basis_species_mineral[i]);
2705 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
2706 EXPECT_NEAR(egs.getSolventWaterMass(), 1.75, 1.0E-15);
2707 EXPECT_NEAR(egs.getSolventMassAndFreeMolalityAndMineralMoles()[0], 1.75, 1.0E-15);
2708 }
2709 else if (mgd.basis_species_name[i] == "H+")
2710 {
2711 EXPECT_FALSE(mgd.basis_species_gas[i]);
2712 EXPECT_FALSE(mgd.basis_species_mineral[i]);
2713 EXPECT_TRUE(egs.getBasisActivityKnown()[i]);
2714 EXPECT_NEAR(egs.getBasisActivity(i), 1E-5, 1.0E-15);
2715 }
2716 else if (mgd.basis_species_name[i] == "Fe+++")
2717 {
2718 EXPECT_FALSE(mgd.basis_species_gas[i]);
2719 EXPECT_FALSE(mgd.basis_species_mineral[i]);
2720 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
2721 EXPECT_NEAR(egs.getBulkMolesOld()[i], 1E-5 + 1E-6 + 2 * 2E-6, 1.0E-15);
2722 }
2723 else if (mgd.basis_species_name[i] == "HCO3-")
2724 {
2725 EXPECT_FALSE(mgd.basis_species_gas[i]);
2726 EXPECT_FALSE(mgd.basis_species_mineral[i]);
2727 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
2728 // do not know the bulk composition as it is dictated by charge neutrality
2729 EXPECT_EQ(egs.getChargeBalanceBasisIndex(), i);
2730 }
2731 else
2732 FAIL() << "Incorrect basis species";
2733
2734 // check total charge
2735 EXPECT_NEAR(egs.getTotalChargeOld(), 0.0, 1E-15);
2736 // check total charge by summing up basis charges
2737 Real tot_charge = 0.0;
2738 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
2739 tot_charge += egs.getBulkMolesOld()[i] * mgd.basis_species_charge[i];
2740 EXPECT_NEAR(tot_charge, 0.0, 1E-15);
2741
2742 // check basis activity = activity_coefficient * molality
2743 for (unsigned i = 1; i < egs.getNumInBasis(); ++i) // don't loop over water
2744 if (mgd.basis_species_gas[i] || mgd.basis_species_mineral[i] || egs.getBasisActivityKnown()[i])
2745 continue;
2746 else
2747 EXPECT_NEAR(egs.getBasisActivity(i),
2748 egs.getBasisActivityCoefficient(i) *
2749 egs.getSolventMassAndFreeMolalityAndMineralMoles()[i],
2750 1E-15);
2751
2752 // check residuals are zero
2753 mole_additions.zero(); // to remove the kinetic rate contributions
2754 for (unsigned a = 0; a < egs.getNumInAlgebraicSystem(); ++a)
2755 EXPECT_LE(std::abs(egs.getResidualComponent(a, mole_additions)), 1E-15);
2756 // check residuals are zero by summing molalities, bulk compositions, etc
2757 Real nw = egs.getSolventWaterMass();
2758 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
2759 {
2760 Real res = -egs.getBulkMolesOld()[i];
2761 if (i == 0)
2762 res += nw * GeochemistryConstants::MOLES_PER_KG_WATER;
2763 else if (mgd.basis_species_mineral[i])
2764 res += egs.getSolventMassAndFreeMolalityAndMineralMoles()[i];
2765 else if (mgd.basis_species_gas[i])
2766 res += 0.0;
2767 else
2768 res += nw * egs.getSolventMassAndFreeMolalityAndMineralMoles()[i];
2769 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
2770 res += nw * mgd.eqm_stoichiometry(j, i) * egs.getEquilibriumMolality(j);
2771 for (unsigned k = 0; k < mgd.kin_species_name.size(); ++k)
2772 res += mgd.kin_stoichiometry(k, i) * egs.getKineticMoles(k);
2773 EXPECT_LE(std::abs(res), 1E-13);
2774 }
2775
2776 // check equilibrium mass balance
2777 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
2778 {
2779 if (mgd.eqm_species_mineral[j] || mgd.eqm_species_gas[j])
2780 continue;
2781 Real log10ap = 0.0;
2782 for (unsigned i = 0; i < mgd.basis_species_name.size(); ++i)
2783 log10ap += mgd.eqm_stoichiometry(j, i) * std::log10(egs.getBasisActivity(i));
2784 log10ap -= egs.getLog10K(j);
2785 if (mgd.surface_sorption_related[j])
2786 log10ap -= std::log10(egs.getSurfacePotential(mgd.surface_sorption_number[j]));
2787 else
2788 log10ap -= std::log10(egs.getEquilibriumActivityCoefficient(j));
2789 if (log10ap < -300.0)
2790 EXPECT_LE(std::abs(egs.getEquilibriumMolality(j)), 1E-25);
2791 else
2792 EXPECT_LE(std::abs(egs.getEquilibriumMolality(j) - std::pow(10.0, log10ap)), 1E-15);
2793 }
2794}
2795
2797TEST(GeochemicalSolverTest, solve_kinetic3)
2798{
2799 PertinentGeochemicalSystem model(
2800 db_solver, {"H2O", "H+", "Fe+++", "HCO3-"}, {}, {}, {"Something"}, {}, {}, "O2(aq)", "e-");
2801 // Define rate. The following produces
2802 // rate = 1.0E-2 * moles_something * 88.8537 * a_{H+} * exp(1E5 / R * (1/303.15 - 1/T))
2803 // = 4.56796E-06 * moles_something
2804 // Something = -3H+ + Fe+++ + 2H2O + 1.5HCO3-, so Q=1E15*1E-5*1E-7=1E3 (approx) >> K, so
2805 // rate_Something produces a negative rate, so mole_additions > 0, so Something will increase in
2806 // mole number
2807 KineticRateUserDescription rate_Something("Something",
2808 1.0E-2,
2809 1.0,
2810 true,
2811 0.0,
2812 0.0,
2813 0.0,
2814 {"H+"},
2815 {1.0},
2816 {0.0},
2817 {0.0},
2818 1.0,
2819 0.0,
2820 1E5,
2821 1.0 / 303.15,
2822 DirectionChoiceEnum::BOTH,
2823 "H2O",
2824 0.0,
2825 -1.0,
2826 0.0);
2827 model.addKineticRate(rate_Something);
2828
2829 ModelGeochemicalDatabase mgd = model.modelGeochemicalDatabase();
2830 std::vector<GeochemistryUnitConverter::GeochemistryUnit> ku;
2831 ku.push_back(GeochemistryUnitConverter::GeochemistryUnit::MOLES);
2832 GeochemicalSystem egs(mgd,
2833 ac_solver,
2834 is_solver,
2835 swapper_kin,
2836 {},
2837 {},
2838 "HCO3-",
2839 {"H2O", "H+", "Fe+++", "HCO3-"},
2840 {1.75, 1E-5, 1E-5, 4E-5},
2841 cu4,
2842 cm4,
2843 25,
2844 0,
2845 1E-20,
2846 {"Something"},
2847 {1.0E-6},
2848 ku);
2849 GeochemicalSolver solver(mgd.basis_species_name.size(),
2850 mgd.kin_species_name.size(),
2851 is_solver,
2852 1.0E-15,
2853 1.0E-200,
2854 100,
2855 1E-5,
2856 1E-16,
2857 1,
2858 {},
2859 3.0,
2860 0,
2861 true);
2862
2863 std::stringstream ss;
2864 unsigned tot_iter;
2865 Real abs_residual;
2866 DenseVector<Real> mole_additions(egs.getNumInBasis() + egs.getNumKinetic());
2867 DenseMatrix<Real> dmole_additions(egs.getNumInBasis() + egs.getNumKinetic(),
2868 egs.getNumInBasis() + egs.getNumKinetic());
2869 // solve with a time-step size of 1E4
2870 solver.solveSystem(egs, ss, tot_iter, abs_residual, 1.0E4, mole_additions, dmole_additions);
2871
2872 // check Newton has converged
2873 EXPECT_LE(tot_iter, (unsigned)100);
2874 EXPECT_LE(abs_residual, 1.0E-15);
2875
2876 // check numbers are correct
2877 EXPECT_EQ(egs.getNumInBasis(), (unsigned)4);
2878 EXPECT_EQ(egs.getNumRedox(), (unsigned)0);
2879 EXPECT_EQ(egs.getNumSurfacePotentials(), (unsigned)0);
2880 EXPECT_EQ(egs.getNumInEquilibrium(), mgd.eqm_species_name.size());
2881 EXPECT_EQ(egs.getNumKinetic(), (unsigned)1);
2882
2883 // check activity products and log10K are ordered as calculated above
2884 const unsigned index_something = 0;
2885 EXPECT_GT(egs.log10KineticActivityProduct(index_something),
2886 egs.getKineticLog10K(index_something));
2887
2888 // check that mole additions is as calculated above
2889 EXPECT_NEAR(mole_additions(egs.getNumInBasis() + index_something) /
2890 (1.0E4 * 4.56796204026978E-06 * egs.getKineticMoles(0)),
2891 1.0,
2892 1.0E-8);
2893
2894 // check that kinetic moles have changed according to the constant rates
2895 for (unsigned k = 0; k < egs.getNumKinetic(); ++k)
2896 if (mgd.kin_species_name[k] == "Something")
2897 EXPECT_NEAR(egs.getKineticMoles(k),
2898 1.0E-6 + mole_additions(egs.getNumInBasis() + index_something),
2899 1.0E-12);
2900 else
2901 FAIL() << "Incorrect kinetic species";
2902
2903 // check that the constraints are satisfied
2904 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
2905 if (mgd.basis_species_name[i] == "H2O")
2906 {
2907 EXPECT_FALSE(mgd.basis_species_gas[i]);
2908 EXPECT_FALSE(mgd.basis_species_mineral[i]);
2909 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
2910 EXPECT_NEAR(egs.getSolventWaterMass(), 1.75, 1.0E-15);
2911 EXPECT_NEAR(egs.getSolventMassAndFreeMolalityAndMineralMoles()[0], 1.75, 1.0E-15);
2912 }
2913 else if (mgd.basis_species_name[i] == "H+")
2914 {
2915 EXPECT_FALSE(mgd.basis_species_gas[i]);
2916 EXPECT_FALSE(mgd.basis_species_mineral[i]);
2917 EXPECT_TRUE(egs.getBasisActivityKnown()[i]);
2918 EXPECT_NEAR(egs.getBasisActivity(i), 1E-5, 1.0E-15);
2919 }
2920 else if (mgd.basis_species_name[i] == "Fe+++")
2921 {
2922 EXPECT_FALSE(mgd.basis_species_gas[i]);
2923 EXPECT_FALSE(mgd.basis_species_mineral[i]);
2924 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
2925 EXPECT_NEAR(egs.getBulkMolesOld()[i], 1E-5 + 1E-6, 1E-15); // 1E-6 from Something
2926 }
2927 else if (mgd.basis_species_name[i] == "HCO3-")
2928 {
2929 EXPECT_FALSE(mgd.basis_species_gas[i]);
2930 EXPECT_FALSE(mgd.basis_species_mineral[i]);
2931 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
2932 // do not know the bulk composition as it is dictated by charge neutrality
2933 EXPECT_EQ(egs.getChargeBalanceBasisIndex(), i);
2934 }
2935 else
2936 FAIL() << "Incorrect basis species";
2937
2938 // check total charge
2939 EXPECT_NEAR(egs.getTotalChargeOld(), 0.0, 1E-15);
2940 // check total charge by summing up basis charges
2941 Real tot_charge = 0.0;
2942 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
2943 tot_charge += egs.getBulkMolesOld()[i] * mgd.basis_species_charge[i];
2944 EXPECT_NEAR(tot_charge, 0.0, 1E-15);
2945
2946 // check basis activity = activity_coefficient * molality
2947 for (unsigned i = 1; i < egs.getNumInBasis(); ++i) // don't loop over water
2948 if (mgd.basis_species_gas[i] || mgd.basis_species_mineral[i] || egs.getBasisActivityKnown()[i])
2949 continue;
2950 else
2951 EXPECT_NEAR(egs.getBasisActivity(i),
2952 egs.getBasisActivityCoefficient(i) *
2953 egs.getSolventMassAndFreeMolalityAndMineralMoles()[i],
2954 1E-15);
2955
2956 // check residuals are zero
2957 mole_additions.zero(); // to remove the kinetic rate contributions
2958 for (unsigned a = 0; a < egs.getNumInAlgebraicSystem(); ++a)
2959 EXPECT_LE(std::abs(egs.getResidualComponent(a, mole_additions)), 1E-15);
2960 // check residuals are zero by summing molalities, bulk compositions, etc
2961 Real nw = egs.getSolventWaterMass();
2962 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
2963 {
2964 Real res = -egs.getBulkMolesOld()[i];
2965 if (i == 0)
2966 res += nw * GeochemistryConstants::MOLES_PER_KG_WATER;
2967 else if (mgd.basis_species_mineral[i])
2968 res += egs.getSolventMassAndFreeMolalityAndMineralMoles()[i];
2969 else if (mgd.basis_species_gas[i])
2970 res += 0.0;
2971 else
2972 res += nw * egs.getSolventMassAndFreeMolalityAndMineralMoles()[i];
2973 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
2974 res += nw * mgd.eqm_stoichiometry(j, i) * egs.getEquilibriumMolality(j);
2975 for (unsigned k = 0; k < mgd.kin_species_name.size(); ++k)
2976 res += mgd.kin_stoichiometry(k, i) * egs.getKineticMoles(k);
2977 EXPECT_LE(std::abs(res), 1E-13);
2978 }
2979
2980 // check equilibrium mass balance
2981 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
2982 {
2983 if (mgd.eqm_species_mineral[j] || mgd.eqm_species_gas[j])
2984 continue;
2985 Real log10ap = 0.0;
2986 for (unsigned i = 0; i < mgd.basis_species_name.size(); ++i)
2987 log10ap += mgd.eqm_stoichiometry(j, i) * std::log10(egs.getBasisActivity(i));
2988 log10ap -= egs.getLog10K(j);
2989 if (mgd.surface_sorption_related[j])
2990 log10ap -= std::log10(egs.getSurfacePotential(mgd.surface_sorption_number[j]));
2991 else
2992 log10ap -= std::log10(egs.getEquilibriumActivityCoefficient(j));
2993 if (log10ap < -300.0)
2994 EXPECT_LE(std::abs(egs.getEquilibriumMolality(j)), 1E-25);
2995 else
2996 EXPECT_LE(std::abs(egs.getEquilibriumMolality(j) - std::pow(10.0, log10ap)), 1E-15);
2997 }
2998}
db_full
const GeochemicalDatabaseReader db_full("../database/moose_geochemdb.json", true, true, false)
model_simplest
const PertinentGeochemicalSystem model_simplest(db_solver, {"H2O", "H+"}, {}, {}, {}, {}, {}, "O2(aq)", "e-")
TEST
TEST(GeochemicalSolverTest, exception)
Check exception.
Definition GeochemicalSolverTest.C:44
swapper2
GeochemistrySpeciesSwapper swapper2(2, 1E-6)
is_solver
GeochemistryIonicStrength is_solver(3.0, 3.0, false, false)
cm2
const std::vector< GeochemicalSystem::ConstraintUserMeaningEnum > cm2
Definition GeochemicalSolverTest.C:24
cu4
const std::vector< GeochemistryUnitConverter::GeochemistryUnit > cu4
Definition GeochemicalSolverTest.C:35
db_ferric
const GeochemicalDatabaseReader db_ferric("../test/database/ferric_hydroxide_sorption.json", true, true)
cm4
const std::vector< GeochemicalSystem::ConstraintUserMeaningEnum > cm4
Definition GeochemicalSolverTest.C:30
ac_solver
GeochemistryActivityCoefficientsDebyeHuckel ac_solver(is_solver, db_solver)
cu2
const std::vector< GeochemistryUnitConverter::GeochemistryUnit > cu2
Definition GeochemicalSolverTest.C:27
db_solver
const GeochemicalDatabaseReader db_solver("database/moose_testdb.json", true, true, false)
swapper_kin
GeochemistrySpeciesSwapper swapper_kin(4, 1E-6)
model
const PertinentGeochemicalSystem model(database, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++", ">(s)FeOH", "radius_neg1", "radius_neg1.5"}, {"Calcite"}, {}, {"Calcite_asdf"}, {"CH4(aq)"}, {">(s)FeOCa+"}, "O2(aq)", "e-")
mgd
const ModelGeochemicalDatabase mgd
Definition GeochemistryActivityCoefficientsDebyeHuckelTest.C:28
GeochemicalDatabaseReader
Class for reading geochemical reactions from a MOOSE geochemical database.
Definition GeochemicalDatabaseReader.h:216
GeochemicalSolver
This class contains methods to solve the algebraic system in GeochemicalSystem.
Definition GeochemicalSolver.h:19
GeochemicalSolver::solveSystem
void solveSystem(GeochemicalSystem &egs, std::stringstream &ss, unsigned &tot_iter, Real &abs_residual, Real dt, DenseVector< Real > &mole_additions, DenseMatrix< Real > &dmole_additions)
Solve the system.
Definition GeochemicalSolver.C:316
GeochemicalSolver::setMaxInitialResidual
void setMaxInitialResidual(Real max_initial_residual)
Set value for max_initial_residual.
Definition GeochemicalSolver.C:64
GeochemicalSystem
This class holds information about bulk composition, molalities, activities, activity coefficients,...
Definition GeochemicalSystem.h:40
GeochemicalSystem::closeSystem
void closeSystem()
Changes a KG_SOLVENT_WATER constraint to MOLES_BULK_WATER (if there is such a constraint) and all FRE...
Definition GeochemicalSystem.C:2100
GeochemistryActivityCoefficientsDebyeHuckel
Computes activity coefficients for non-minerals and non-gases (since these species do not have activi...
Definition GeochemistryActivityCoefficientsDebyeHuckel.h:59
GeochemistryIonicStrength
Calculators to compute ionic strength and stoichiometric ionic strength.
Definition GeochemistryIonicStrength.h:18
GeochemistrySpeciesSwapper
Class to swap basis species with equilibrium species.
Definition GeochemistrySpeciesSwapper.h:20
PertinentGeochemicalSystem
Constructs and stores a minimal amount of information that is pertinent to the user-defined geochemic...
Definition PertinentGeochemicalSystem.h:387
PertinentGeochemicalSystem::modelGeochemicalDatabase
const ModelGeochemicalDatabase & modelGeochemicalDatabase() const
Return a reference to the ModelGeochemicalDatabase structure.
Definition PertinentGeochemicalSystem.C:850
PertinentGeochemicalSystem::addKineticRate
void addKineticRate(const KineticRateUserDescription &description)
Adds a rate description for kinetic_species.
Definition PertinentGeochemicalSystem.C:856
GeochemistryConstants::DIELECTRIC_CONSTANT_WATER
constexpr Real DIELECTRIC_CONSTANT_WATER
Definition GeochemistryConstants.h:24
GeochemistryConstants::PERMITTIVITY_FREE_SPACE
constexpr Real PERMITTIVITY_FREE_SPACE
Definition GeochemistryConstants.h:23
KineticRateUserDescription
Holds a user-specified description of a kinetic rate.
Definition GeochemistryKineticRateCalculator.h:76
ModelGeochemicalDatabase
Data structure to hold all relevant information from the database file.
Definition PertinentGeochemicalSystem.h:102
ModelGeochemicalDatabase::eqm_species_mineral
std::vector< bool > eqm_species_mineral
eqm_species_mineral[i] = true iff the i^th equilibrium species is a mineral
Definition PertinentGeochemicalSystem.h:200
ModelGeochemicalDatabase::eqm_stoichiometry
DenseMatrix< Real > eqm_stoichiometry
eqm_stoichiometry(i, j) = stoichiometric coefficient for equilibrium species "i" in terms of the basi...
Definition PertinentGeochemicalSystem.h:224
ModelGeochemicalDatabase::kin_stoichiometry
DenseMatrix< Real > kin_stoichiometry
kin_stoichiometry(i, j) = stoichiometric coefficient for kinetic species "i" in terms of the basis sp...
Definition PertinentGeochemicalSystem.h:334
ModelGeochemicalDatabase::basis_species_charge
std::vector< Real > basis_species_charge
all quantities have a charge (mineral charge = 0, gas charge = 0, oxide charge = 0)
Definition PertinentGeochemicalSystem.h:177
ModelGeochemicalDatabase::surface_sorption_number
std::vector< unsigned > surface_sorption_number
surface_sorption_number[j] = the index of the surface potential that should be used to modify the equ...
Definition PertinentGeochemicalSystem.h:255
ModelGeochemicalDatabase::kin_species_name
std::vector< std::string > kin_species_name
kin_species_name[j] = name of the j^th kinetic species
Definition PertinentGeochemicalSystem.h:307
ModelGeochemicalDatabase::basis_species_mineral
std::vector< bool > basis_species_mineral
basis_species_mineral[j] = true iff the j^th basis species is a mineral
Definition PertinentGeochemicalSystem.h:168
ModelGeochemicalDatabase::surface_sorption_related
std::vector< bool > surface_sorption_related
surface_sorption_related[j] = true iff the j^th equilibrium species is involved in surface sorption
Definition PertinentGeochemicalSystem.h:246
ModelGeochemicalDatabase::basis_species_name
std::vector< std::string > basis_species_name
basis_species_name[j] = name of the j^th basis species
Definition PertinentGeochemicalSystem.h:165
ModelGeochemicalDatabase::basis_species_gas
std::vector< bool > basis_species_gas
basis_species_gas[j] = true iff the j^th basis species is a gas
Definition PertinentGeochemicalSystem.h:171
ModelGeochemicalDatabase::kin_species_index
std::unordered_map< std::string, unsigned > kin_species_index
kin_species_index[name] = index of the kinetic species, within all ModelGeochemicalDatabase internal ...
Definition PertinentGeochemicalSystem.h:304
ModelGeochemicalDatabase::eqm_species_name
std::vector< std::string > eqm_species_name
eqm_species_name[i] = name of the i^th eqm species
Definition PertinentGeochemicalSystem.h:197
ModelGeochemicalDatabase::eqm_species_charge
std::vector< Real > eqm_species_charge
all quantities have a charge (mineral charge = 0, gas charge = 0, oxide charge = 0)
Definition PertinentGeochemicalSystem.h:209
ModelGeochemicalDatabase::eqm_species_gas
std::vector< bool > eqm_species_gas
eqm_species_gas[i] = true iff the i^th equilibrium species is a gas
Definition PertinentGeochemicalSystem.h:203
Generated on Fri Aug 21 2026 13:41:03 for https://mooseframework.inl.gov by doxygen 1.9.8