https://mooseframework.inl.gov
Functions | Variables
GeochemicalSolverTest.C File Reference

Go to the source code of this file.

Functions

 TEST (GeochemicalSolverTest, exception)
 Check exception. More...
 
 TEST (GeochemicalSolverTest, setgetMaxInitialResidual)
 Check set/getMaxInitialResidual. More...
 
 TEST (GeochemicalSolverTest, solve1)
 Solve super-simple case. More...
 
 TEST (GeochemicalSolverTest, solve2)
 Solve realistic case with minerals and gases, but no precipitation, no redox, no sorption. More...
 
 TEST (GeochemicalSolverTest, solve3)
 Solve realistic case with minerals that precipitate, no redox, no sorption. More...
 
 TEST (GeochemicalSolverTest, solve3_restore)
 Solve realistic case with minerals that precipitate, no redox, no sorption, and then "restore" to check it works correctly, and then check copy-assignment of GeochemicalSystem works properly. More...
 
 TEST (GeochemicalSolverTest, solve4)
 Solve realistic case with redox disequilibrium, no minerals, no sorption. More...
 
 TEST (GeochemicalSolverTest, solve4_restore)
 Solve realistic case with redox disequilibrium, no minerals, no sorption, then "restore" the solution and check. More...
 
 TEST (GeochemicalSolverTest, solve5)
 Solve realistic case with sorption and minerals, no redox. More...
 
 TEST (GeochemicalSolverTest, solve5_restore)
 Repeat the realistic case with sorption and minerals, no redox, this time with a "restore", and then use the copy assignment operator of GeochemicalSystem to create a new geochemical system. More...
 
 TEST (GeochemicalSolverTest, solve_addH)
 Test progressively adding H+ to a system. More...
 
 TEST (GeochemicalSolverTest, maxSwapsException)
 Test that the max swaps allowed works OK. More...
 
 TEST (GeochemicalSolverTest, setRampMaxIonicStrength)
 Check setRampMaxIonicStrength. More...
 
 TEST (GeochemicalSolverTest, solve_kinetic1)
 Solve case that involves kinetic species with zero rates (so kinetic species should have no impact except to modify the bulk composition) More...
 
 TEST (GeochemicalSolverTest, solve_kinetic2)
 Solve case that involves kinetic species with constant rates. More...
 
 TEST (GeochemicalSolverTest, solve_kinetic3)
 Solve case that involves kinetic species with promoting indices and implicit solve. More...
 

Variables

const GeochemicalDatabaseReader db_solver ("database/moose_testdb.json", true, true, false)
 
const GeochemicalDatabaseReader db_full ("../database/moose_geochemdb.json", true, true, false)
 
const GeochemicalDatabaseReader db_ferric ("../test/database/ferric_hydroxide_sorption.json", true, true)
 
const PertinentGeochemicalSystem model_simplest (db_solver, {"H2O", "H+"}, {}, {}, {}, {}, {}, "O2(aq)", "e-")
 
GeochemistrySpeciesSwapper swapper2 (2, 1E-6)
 
GeochemistrySpeciesSwapper swapper_kin (4, 1E-6)
 
const std::vector< GeochemicalSystem::ConstraintUserMeaningEnum > cm2
 
const std::vector< GeochemistryUnitConverter::GeochemistryUnit > cu2
 
const std::vector< GeochemicalSystem::ConstraintUserMeaningEnum > cm4
 
const std::vector< GeochemistryUnitConverter::GeochemistryUnit > cu4
 
GeochemistryIonicStrength is_solver (3.0, 3.0, false, false)
 
GeochemistryActivityCoefficientsDebyeHuckel ac_solver (is_solver, db_solver)
 

Function Documentation

◆ TEST() [1/16]

TEST ( GeochemicalSolverTest  ,
exception   
)

Check exception.

Definition at line 44 of file GeochemicalSolverTest.C.

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 }
is_solver
GeochemistryIonicStrength is_solver(3.0, 3.0, false, false)
mgd
const ModelGeochemicalDatabase mgd
Definition: GeochemistryActivityCoefficientsDebyeHuckelTest.C:28
GeochemicalSolver::setMaxInitialResidual
void setMaxInitialResidual(Real max_initial_residual)
Set value for max_initial_residual.
Definition: GeochemicalSolver.C:64
swapper2
GeochemistrySpeciesSwapper swapper2(2, 1E-6)
PertinentGeochemicalSystem::modelGeochemicalDatabase
const ModelGeochemicalDatabase & modelGeochemicalDatabase() const
Return a reference to the ModelGeochemicalDatabase structure.
Definition: PertinentGeochemicalSystem.C:850
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::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
model_simplest
const PertinentGeochemicalSystem model_simplest(db_solver, {"H2O", "H+"}, {}, {}, {}, {}, {}, "O2(aq)", "e-")
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
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

◆ TEST() [2/16]

TEST ( GeochemicalSolverTest  ,
setgetMaxInitialResidual   
)

Check set/getMaxInitialResidual.

Definition at line 240 of file GeochemicalSolverTest.C.

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 }
is_solver
GeochemistryIonicStrength is_solver(3.0, 3.0, false, false)
mgd
const ModelGeochemicalDatabase mgd
Definition: GeochemistryActivityCoefficientsDebyeHuckelTest.C:28
swapper2
GeochemistrySpeciesSwapper swapper2(2, 1E-6)
PertinentGeochemicalSystem::modelGeochemicalDatabase
const ModelGeochemicalDatabase & modelGeochemicalDatabase() const
Return a reference to the ModelGeochemicalDatabase structure.
Definition: PertinentGeochemicalSystem.C:850
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::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
model_simplest
const PertinentGeochemicalSystem model_simplest(db_solver, {"H2O", "H+"}, {}, {}, {}, {}, {}, "O2(aq)", "e-")
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
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

◆ TEST() [3/16]

TEST ( GeochemicalSolverTest  ,
solve1   
)

Solve super-simple case.

Definition at line 279 of file GeochemicalSolverTest.C.

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 }
GeochemistryConstants::MOLES_PER_KG_WATER
constexpr Real MOLES_PER_KG_WATER
Definition: GeochemistryConstants.h:16
is_solver
GeochemistryIonicStrength is_solver(3.0, 3.0, false, false)
mgd
const ModelGeochemicalDatabase mgd
Definition: GeochemistryActivityCoefficientsDebyeHuckelTest.C:28
swapper2
GeochemistrySpeciesSwapper swapper2(2, 1E-6)
PertinentGeochemicalSystem::modelGeochemicalDatabase
const ModelGeochemicalDatabase & modelGeochemicalDatabase() const
Return a reference to the ModelGeochemicalDatabase structure.
Definition: PertinentGeochemicalSystem.C:850
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::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
model_simplest
const PertinentGeochemicalSystem model_simplest(db_solver, {"H2O", "H+"}, {}, {}, {}, {}, {}, "O2(aq)", "e-")
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::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

◆ TEST() [4/16]

TEST ( GeochemicalSolverTest  ,
solve2   
)

Solve realistic case with minerals and gases, but no precipitation, no redox, no sorption.

Definition at line 346 of file GeochemicalSolverTest.C.

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 }
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
mgd
const ModelGeochemicalDatabase mgd
Definition: GeochemistryActivityCoefficientsDebyeHuckelTest.C:28
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
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
Constructs and stores a minimal amount of information that is pertinent to the user-defined geochemic...
Definition: PertinentGeochemicalSystem.h:382
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-")
Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
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
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
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
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

◆ TEST() [5/16]

TEST ( GeochemicalSolverTest  ,
solve3   
)

Solve realistic case with minerals that precipitate, no redox, no sorption.

Definition at line 649 of file GeochemicalSolverTest.C.

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 }
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
mgd
const ModelGeochemicalDatabase mgd
Definition: GeochemistryActivityCoefficientsDebyeHuckelTest.C:28
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
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
Constructs and stores a minimal amount of information that is pertinent to the user-defined geochemic...
Definition: PertinentGeochemicalSystem.h:382
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-")
Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
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
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
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
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

◆ TEST() [6/16]

TEST ( GeochemicalSolverTest  ,
solve3_restore   
)

Solve realistic case with minerals that precipitate, no redox, no sorption, and then "restore" to check it works correctly, and then check copy-assignment of GeochemicalSystem works properly.

Definition at line 962 of file GeochemicalSolverTest.C.

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 }
is_solver
GeochemistryIonicStrength is_solver(3.0, 3.0, false, false)
mgd
const ModelGeochemicalDatabase mgd
Definition: GeochemistryActivityCoefficientsDebyeHuckelTest.C:28
GeochemistrySpeciesSwapper
Class to swap basis species with equilibrium species.
Definition: GeochemistrySpeciesSwapper.h:19
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
Constructs and stores a minimal amount of information that is pertinent to the user-defined geochemic...
Definition: PertinentGeochemicalSystem.h:382
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-")
Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
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::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")
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
GeochemicalSolver
This class contains methods to solve the algebraic system in GeochemicalSystem.
Definition: GeochemicalSolver.h:18
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
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

◆ TEST() [7/16]

TEST ( GeochemicalSolverTest  ,
solve4   
)

Solve realistic case with redox disequilibrium, no minerals, no sorption.

Definition at line 1171 of file GeochemicalSolverTest.C.

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 }
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
mgd
const ModelGeochemicalDatabase mgd
Definition: GeochemistryActivityCoefficientsDebyeHuckelTest.C:28
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
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
Constructs and stores a minimal amount of information that is pertinent to the user-defined geochemic...
Definition: PertinentGeochemicalSystem.h:382
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-")
Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
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
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
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
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

◆ TEST() [8/16]

TEST ( GeochemicalSolverTest  ,
solve4_restore   
)

Solve realistic case with redox disequilibrium, no minerals, no sorption, then "restore" the solution and check.

Definition at line 1513 of file GeochemicalSolverTest.C.

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 }
is_solver
GeochemistryIonicStrength is_solver(3.0, 3.0, false, false)
mgd
const ModelGeochemicalDatabase mgd
Definition: GeochemistryActivityCoefficientsDebyeHuckelTest.C:28
GeochemistrySpeciesSwapper
Class to swap basis species with equilibrium species.
Definition: GeochemistrySpeciesSwapper.h:19
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
Constructs and stores a minimal amount of information that is pertinent to the user-defined geochemic...
Definition: PertinentGeochemicalSystem.h:382
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-")
Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
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::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")
GeochemicalSolver
This class contains methods to solve the algebraic system in GeochemicalSystem.
Definition: GeochemicalSolver.h:18
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
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

◆ TEST() [9/16]

TEST ( GeochemicalSolverTest  ,
solve5   
)

Solve realistic case with sorption and minerals, no redox.

Definition at line 1692 of file GeochemicalSolverTest.C.

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 }
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
mgd
const ModelGeochemicalDatabase mgd
Definition: GeochemistryActivityCoefficientsDebyeHuckelTest.C:28
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
PertinentGeochemicalSystem::modelGeochemicalDatabase
const ModelGeochemicalDatabase & modelGeochemicalDatabase() const
Return a reference to the ModelGeochemicalDatabase structure.
Definition: PertinentGeochemicalSystem.C:850
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
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-")
Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
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
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
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
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

◆ TEST() [10/16]

TEST ( GeochemicalSolverTest  ,
solve5_restore   
)

Repeat the realistic case with sorption and minerals, no redox, this time with a "restore", and then use the copy assignment operator of GeochemicalSystem to create a new geochemical system.

Definition at line 1964 of file GeochemicalSolverTest.C.

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 }
GeochemistryConstants::CELSIUS_TO_KELVIN
constexpr Real CELSIUS_TO_KELVIN
Definition: GeochemistryConstants.h:22
is_solver
GeochemistryIonicStrength is_solver(3.0, 3.0, false, false)
mgd
const ModelGeochemicalDatabase mgd
Definition: GeochemistryActivityCoefficientsDebyeHuckelTest.C:28
GeochemistrySpeciesSwapper
Class to swap basis species with equilibrium species.
Definition: GeochemistrySpeciesSwapper.h:19
PertinentGeochemicalSystem::modelGeochemicalDatabase
const ModelGeochemicalDatabase & modelGeochemicalDatabase() const
Return a reference to the ModelGeochemicalDatabase structure.
Definition: PertinentGeochemicalSystem.C:850
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
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-")
Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
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::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")
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
GeochemicalSolver
This class contains methods to solve the algebraic system in GeochemicalSystem.
Definition: GeochemicalSolver.h:18
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
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

◆ TEST() [11/16]

TEST ( GeochemicalSolverTest  ,
solve_addH   
)

Test progressively adding H+ to a system.

Definition at line 2127 of file GeochemicalSolverTest.C.

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 }
GeochemistryConstants::MOLES_PER_KG_WATER
constexpr Real MOLES_PER_KG_WATER
Definition: GeochemistryConstants.h:16
is_solver
GeochemistryIonicStrength is_solver(3.0, 3.0, false, false)
mgd
const ModelGeochemicalDatabase mgd
Definition: GeochemistryActivityCoefficientsDebyeHuckelTest.C:28
swapper2
GeochemistrySpeciesSwapper swapper2(2, 1E-6)
PertinentGeochemicalSystem::modelGeochemicalDatabase
const ModelGeochemicalDatabase & modelGeochemicalDatabase() const
Return a reference to the ModelGeochemicalDatabase structure.
Definition: PertinentGeochemicalSystem.C:850
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::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
model_simplest
const PertinentGeochemicalSystem model_simplest(db_solver, {"H2O", "H+"}, {}, {}, {}, {}, {}, "O2(aq)", "e-")
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
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

◆ TEST() [12/16]

TEST ( GeochemicalSolverTest  ,
maxSwapsException   
)

Test that the max swaps allowed works OK.

Definition at line 2218 of file GeochemicalSolverTest.C.

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 }
is_solver
GeochemistryIonicStrength is_solver(3.0, 3.0, false, false)
mgd
const ModelGeochemicalDatabase mgd
Definition: GeochemistryActivityCoefficientsDebyeHuckelTest.C:28
GeochemistrySpeciesSwapper
Class to swap basis species with equilibrium species.
Definition: GeochemistrySpeciesSwapper.h:19
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
Constructs and stores a minimal amount of information that is pertinent to the user-defined geochemic...
Definition: PertinentGeochemicalSystem.h:382
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-")
Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
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::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
GeochemicalSolver
This class contains methods to solve the algebraic system in GeochemicalSystem.
Definition: GeochemicalSolver.h:18
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

◆ TEST() [13/16]

TEST ( GeochemicalSolverTest  ,
setRampMaxIonicStrength   
)

Check setRampMaxIonicStrength.

Definition at line 2337 of file GeochemicalSolverTest.C.

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 }
is_solver
GeochemistryIonicStrength is_solver(3.0, 3.0, false, false)
mgd
const ModelGeochemicalDatabase mgd
Definition: GeochemistryActivityCoefficientsDebyeHuckelTest.C:28
swapper2
GeochemistrySpeciesSwapper swapper2(2, 1E-6)
PertinentGeochemicalSystem::modelGeochemicalDatabase
const ModelGeochemicalDatabase & modelGeochemicalDatabase() const
Return a reference to the ModelGeochemicalDatabase structure.
Definition: PertinentGeochemicalSystem.C:850
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::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
model_simplest
const PertinentGeochemicalSystem model_simplest(db_solver, {"H2O", "H+"}, {}, {}, {}, {}, {}, "O2(aq)", "e-")
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
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

◆ TEST() [14/16]

TEST ( GeochemicalSolverTest  ,
solve_kinetic1   
)

Solve case that involves kinetic species with zero rates (so kinetic species should have no impact except to modify the bulk composition)

Definition at line 2389 of file GeochemicalSolverTest.C.

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 }
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
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
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
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
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
PertinentGeochemicalSystem::modelGeochemicalDatabase
const ModelGeochemicalDatabase & modelGeochemicalDatabase() const
Return a reference to the ModelGeochemicalDatabase structure.
Definition: PertinentGeochemicalSystem.C:850
PertinentGeochemicalSystem
Constructs and stores a minimal amount of information that is pertinent to the user-defined geochemic...
Definition: PertinentGeochemicalSystem.h:382
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)
Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
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
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::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
db_solver
const GeochemicalDatabaseReader db_solver("database/moose_testdb.json", true, true, false)
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
cu4
const std::vector< GeochemistryUnitConverter::GeochemistryUnit > cu4
Definition: GeochemicalSolverTest.C:35

◆ TEST() [15/16]

TEST ( GeochemicalSolverTest  ,
solve_kinetic2   
)

Solve case that involves kinetic species with constant rates.

Definition at line 2562 of file GeochemicalSolverTest.C.

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 }
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
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
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
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
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
PertinentGeochemicalSystem::modelGeochemicalDatabase
const ModelGeochemicalDatabase & modelGeochemicalDatabase() const
Return a reference to the ModelGeochemicalDatabase structure.
Definition: PertinentGeochemicalSystem.C:850
PertinentGeochemicalSystem::addKineticRate
void addKineticRate(const KineticRateUserDescription &description)
Adds a rate description for kinetic_species.
Definition: PertinentGeochemicalSystem.C:856
PertinentGeochemicalSystem
Constructs and stores a minimal amount of information that is pertinent to the user-defined geochemic...
Definition: PertinentGeochemicalSystem.h:382
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)
Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
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::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
db_solver
const GeochemicalDatabaseReader db_solver("database/moose_testdb.json", true, true, false)
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
cu4
const std::vector< GeochemistryUnitConverter::GeochemistryUnit > cu4
Definition: GeochemicalSolverTest.C:35

◆ TEST() [16/16]

TEST ( GeochemicalSolverTest  ,
solve_kinetic3   
)

Solve case that involves kinetic species with promoting indices and implicit solve.

Definition at line 2797 of file GeochemicalSolverTest.C.

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
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
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
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
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
PertinentGeochemicalSystem::modelGeochemicalDatabase
const ModelGeochemicalDatabase & modelGeochemicalDatabase() const
Return a reference to the ModelGeochemicalDatabase structure.
Definition: PertinentGeochemicalSystem.C:850
PertinentGeochemicalSystem::addKineticRate
void addKineticRate(const KineticRateUserDescription &description)
Adds a rate description for kinetic_species.
Definition: PertinentGeochemicalSystem.C:856
PertinentGeochemicalSystem
Constructs and stores a minimal amount of information that is pertinent to the user-defined geochemic...
Definition: PertinentGeochemicalSystem.h:382
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)
Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
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::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
db_solver
const GeochemicalDatabaseReader db_solver("database/moose_testdb.json", true, true, false)
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
cu4
const std::vector< GeochemistryUnitConverter::GeochemistryUnit > cu4
Definition: GeochemicalSolverTest.C:35

Variable Documentation

◆ ac_solver

GeochemistryActivityCoefficientsDebyeHuckel ac_solver(is_solver, db_solver)

Referenced by TEST().

◆ cm2

const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm2
Initial value:
= {
GeochemicalSystem::ConstraintUserMeaningEnum::KG_SOLVENT_WATER,
GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION}

Definition at line 24 of file GeochemicalSolverTest.C.

Referenced by TEST(), and TEST_F().

◆ cm4

const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm4
Initial value:
= {
GeochemicalSystem::ConstraintUserMeaningEnum::KG_SOLVENT_WATER,
GeochemicalSystem::ConstraintUserMeaningEnum::ACTIVITY,
GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION_WITH_KINETIC}

Definition at line 30 of file GeochemicalSolverTest.C.

Referenced by TEST().

◆ cu2

const std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu2
Initial value:
= {
GeochemistryUnitConverter::GeochemistryUnit::KG,
GeochemistryUnitConverter::GeochemistryUnit::MOLES}

Definition at line 27 of file GeochemicalSolverTest.C.

Referenced by TEST(), and TEST_F().

◆ cu4

const std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu4
Initial value:
= {
GeochemistryUnitConverter::GeochemistryUnit::KG,
GeochemistryUnitConverter::GeochemistryUnit::DIMENSIONLESS,
GeochemistryUnitConverter::GeochemistryUnit::MOLES,
GeochemistryUnitConverter::GeochemistryUnit::MOLES}

Definition at line 35 of file GeochemicalSolverTest.C.

Referenced by TEST().

◆ db_ferric

const GeochemicalDatabaseReader db_ferric("../test/database/ferric_hydroxide_sorption.json", true, true)

Referenced by TEST().

◆ db_full

const GeochemicalDatabaseReader db_full("../database/moose_geochemdb.json", true, true, false)

Referenced by TEST().

◆ db_solver

const GeochemicalDatabaseReader db_solver("database/moose_testdb.json", true, true, false)

Referenced by TEST().

◆ is_solver

GeochemistryIonicStrength is_solver(3.0, 3.0, false, false)

Referenced by TEST().

◆ model_simplest

const PertinentGeochemicalSystem model_simplest(db_solver, {"H2O", "H+"}, {}, {}, {}, {}, {}, "O2(aq)", "e-")

Referenced by TEST().

◆ swapper2

GeochemistrySpeciesSwapper swapper2(2, 1E-6)

Referenced by TEST().

◆ swapper_kin

GeochemistrySpeciesSwapper swapper_kin(4, 1E-6)

Referenced by TEST().

Generated on Wed Nov 5 2025 20:23:13 for https://mooseframework.inl.gov by   doxygen 1.8.14