https://mooseframework.inl.gov
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 
15 const GeochemicalDatabaseReader db_solver("database/moose_testdb.json", true, true, false);
16 const GeochemicalDatabaseReader db_full("../database/moose_geochemdb.json", true, true, false);
17 const GeochemicalDatabaseReader
18  db_ferric("../test/database/ferric_hydroxide_sorption.json", true, true);
19 // Following model only has OH- as an equilibrium species
20 const PertinentGeochemicalSystem
21  model_simplest(db_solver, {"H2O", "H+"}, {}, {}, {}, {}, {}, "O2(aq)", "e-");
22 GeochemistrySpeciesSwapper swapper2(2, 1E-6);
23 GeochemistrySpeciesSwapper swapper_kin(4, 1E-6);
24 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm2 = {
25  GeochemicalSystem::ConstraintUserMeaningEnum::KG_SOLVENT_WATER,
26  GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION};
27 const std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu2 = {
28  GeochemistryUnitConverter::GeochemistryUnit::KG,
29  GeochemistryUnitConverter::GeochemistryUnit::MOLES};
30 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm4 = {
31  GeochemicalSystem::ConstraintUserMeaningEnum::KG_SOLVENT_WATER,
32  GeochemicalSystem::ConstraintUserMeaningEnum::ACTIVITY,
33  GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
34  GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION_WITH_KINETIC};
35 const std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu4 = {
36  GeochemistryUnitConverter::GeochemistryUnit::KG,
37  GeochemistryUnitConverter::GeochemistryUnit::DIMENSIONLESS,
38  GeochemistryUnitConverter::GeochemistryUnit::MOLES,
39  GeochemistryUnitConverter::GeochemistryUnit::MOLES};
40 GeochemistryIonicStrength is_solver(3.0, 3.0, false, false);
41 GeochemistryActivityCoefficientsDebyeHuckel ac_solver(is_solver, db_solver);
42 
44 TEST(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 
240 TEST(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 
279 TEST(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 
346 TEST(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 
649 TEST(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 
962 TEST(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 
1171 TEST(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 
1513 TEST(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 
1692 TEST(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 
1964 TEST(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 
2127 TEST(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 
2218 TEST(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 
2337 TEST(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 
2389 TEST(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 
2562 TEST(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 
2797 TEST(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 }
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:242
GeochemistryConstants::MOLES_PER_KG_WATER
constexpr Real MOLES_PER_KG_WATER
Definition: GeochemistryConstants.h:16
GeochemistryConstants::CELSIUS_TO_KELVIN
constexpr Real CELSIUS_TO_KELVIN
Definition: GeochemistryConstants.h:22
is_solver
GeochemistryIonicStrength is_solver(3.0, 3.0, false, false)
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:199
GeochemistryConstants::DIELECTRIC_CONSTANT_WATER
constexpr Real DIELECTRIC_CONSTANT_WATER
Definition: GeochemistryConstants.h:24
cm4
const std::vector< GeochemicalSystem::ConstraintUserMeaningEnum > cm4
Definition: GeochemicalSolverTest.C:30
mgd
const ModelGeochemicalDatabase mgd
Definition: GeochemistryActivityCoefficientsDebyeHuckelTest.C:28
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:330
GeochemicalSolver::setMaxInitialResidual
void setMaxInitialResidual(Real max_initial_residual)
Set value for max_initial_residual.
Definition: GeochemicalSolver.C:64
swapper2
GeochemistrySpeciesSwapper swapper2(2, 1E-6)
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:300
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:220
GeochemistrySpeciesSwapper
Class to swap basis species with equilibrium species.
Definition: GeochemistrySpeciesSwapper.h:19
GeochemistryConstants::PERMITTIVITY_FREE_SPACE
constexpr Real PERMITTIVITY_FREE_SPACE
Definition: GeochemistryConstants.h:23
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:251
GeochemistryActivityCoefficientsDebyeHuckel
Computes activity coefficients for non-minerals and non-gases (since these species do not have activi...
Definition: GeochemistryActivityCoefficientsDebyeHuckel.h:58
PertinentGeochemicalSystem::modelGeochemicalDatabase
const ModelGeochemicalDatabase & modelGeochemicalDatabase() const
Return a reference to the ModelGeochemicalDatabase structure.
Definition: PertinentGeochemicalSystem.C:850
db_full
const GeochemicalDatabaseReader db_full("../database/moose_geochemdb.json", true, true, false)
PertinentGeochemicalSystem::addKineticRate
void addKineticRate(const KineticRateUserDescription &description)
Adds a rate description for kinetic_species.
Definition: PertinentGeochemicalSystem.C:856
db_ferric
const GeochemicalDatabaseReader db_ferric("../test/database/ferric_hydroxide_sorption.json", true, true)
PertinentGeochemicalSystem
Constructs and stores a minimal amount of information that is pertinent to the user-defined geochemic...
Definition: PertinentGeochemicalSystem.h:382
GeochemistryIonicStrength
Calculators to compute ionic strength and stoichiometric ionic strength.
Definition: GeochemistryIonicStrength.h:17
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
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-")
swapper_kin
GeochemistrySpeciesSwapper swapper_kin(4, 1E-6)
TEST
TEST(GeochemicalSolverTest, exception)
Check exception.
Definition: GeochemicalSolverTest.C:44
Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
cu2
const std::vector< GeochemistryUnitConverter::GeochemistryUnit > cu2
Definition: GeochemicalSolverTest.C:27
ac_solver
GeochemistryActivityCoefficientsDebyeHuckel ac_solver(is_solver, db_solver)
GeochemicalSystem
This class holds information about bulk composition, molalities, activities, activity coefficients...
Definition: GeochemicalSystem.h:39
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:173
KineticRateUserDescription
Holds a user-specified description of a kinetic rate.
Definition: GeochemistryKineticRateCalculator.h:75
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:167
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:164
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:205
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:303
ModelGeochemicalDatabase
Data structure to hold all relevant information from the database file.
Definition: PertinentGeochemicalSystem.h:99
EM::j
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
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:196
model_simplest
const PertinentGeochemicalSystem model_simplest(db_solver, {"H2O", "H+"}, {}, {}, {}, {}, {}, "O2(aq)", "e-")
db_solver
const GeochemicalDatabaseReader db_solver("database/moose_testdb.json", true, true, false)
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
cm2
const std::vector< GeochemicalSystem::ConstraintUserMeaningEnum > cm2
Definition: GeochemicalSolverTest.C:24
GeochemicalSolver
This class contains methods to solve the algebraic system in GeochemicalSystem.
Definition: GeochemicalSolver.h:18
std::pow
MooseUnits pow(const MooseUnits &, int)
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:193
NS::k
static const std::string k
Definition: NS.h:130
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:161
GeochemicalDatabaseReader
Class for reading geochemical reactions from a MOOSE geochemical database.
Definition: GeochemicalDatabaseReader.h:213
cu4
const std::vector< GeochemistryUnitConverter::GeochemistryUnit > cu4
Definition: GeochemicalSolverTest.C:35
Generated on Thu Jul 17 2025 01:43:52 for https://mooseframework.inl.gov by   doxygen 1.8.14