https://mooseframework.inl.gov
GeochemicalSystemTest.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://mooseframework.inl.gov
3 //*
4 //* All rights reserved, see COPYRIGHT for full restrictions
5 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 //*
7 //* Licensed under LGPL 2.1, please see LICENSE for details
8 //* https://www.gnu.org/licenses/lgpl-2.1.html
9 
10 #include "GeochemicalSystemTest.h"
11 
13 TEST_F(GeochemicalSystemTest, constructWithMultiMooseEnum)
14 {
15  MultiMooseEnum constraint_user_meaning(
16  "kg_solvent_water bulk_composition bulk_composition_with_kinetic free_concentration "
17  "free_mineral activity log10activity fugacity log10fugacity");
18  constraint_user_meaning.setAdditionalValue(
19  "activity bulk_composition bulk_composition free_concentration");
20  MultiMooseEnum constraint_unit("dimensionless moles molal kg g mg ug kg_per_kg_solvent "
21  "g_per_kg_solvent mg_per_kg_solvent ug_per_kg_solvent cm3");
22  constraint_unit.setAdditionalValue("dimensionless moles moles molal");
23  ModelGeochemicalDatabase mgd_calcite2 = _model_calcite.modelGeochemicalDatabase();
24  const GeochemicalSystem egs(mgd_calcite2,
25  _ac3,
26  _is3,
27  _swapper4,
28  {"Ca++"},
29  {"Calcite"},
30  "H+",
31  {"H2O", "Calcite", "H+", "HCO3-"},
32  {1.75, 3.0, 2.0, 1.0},
33  constraint_unit,
34  constraint_user_meaning,
35  25,
36  0,
37  1E-20,
38  {},
39  {},
40  MultiMooseEnum(""));
41 }
42 
44 TEST_F(GeochemicalSystemTest, numInBasis) { EXPECT_EQ(_egs_calcite.getNumInBasis(), (unsigned)4); }
45 
47 TEST_F(GeochemicalSystemTest, numInEquilibrium)
48 {
49  EXPECT_EQ(_egs_calcite.getNumInEquilibrium(), (unsigned)6);
50 }
51 
53 TEST_F(GeochemicalSystemTest, swapSizeException)
54 {
55  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
57  {"H2O", "H+", ">(s)FeOH", ">(w)FeOH", "Fe++", "HCO3-", "O2(aq)"},
58  {},
59  {},
60  {},
61  {},
62  {},
63  "O2(aq)",
64  "e-");
66 
67  try
68  {
70  _ac3,
71  _is3,
72  _swapper7,
73  {"H+"},
74  {},
75  "H+",
76  {},
77  {},
78  _cu_dummy,
79  _cm_dummy,
80  25,
81  0,
82  1E-20,
83  {},
84  {},
85  {});
86  FAIL() << "Missing expected exception.";
87  }
88  catch (const std::exception & e)
89  {
90  std::string msg(e.what());
91  ASSERT_TRUE(msg.find("swap_out_of_basis must have same length as swap_into_basis") !=
92  std::string::npos)
93  << "Failed with unexpected error message: " << msg;
94  }
95 }
96 
98 TEST_F(GeochemicalSystemTest, swapChargeBalanceException)
99 {
100  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
102  {"H2O", "H+", ">(s)FeOH", ">(w)FeOH", "Fe++", "HCO3-", "O2(aq)"},
103  {},
104  {},
105  {},
106  {},
107  {},
108  "O2(aq)",
109  "e-");
111 
112  try
113  {
115  _ac3,
116  _is3,
117  _swapper7,
118  {"H+"},
119  {"FeOH"},
120  "H+",
121  {},
122  {},
123  _cu_dummy,
124  _cm_dummy,
125  25,
126  0,
127  1E-20,
128  {},
129  {},
130  {});
131  FAIL() << "Missing expected exception.";
132  }
133  catch (const std::exception & e)
134  {
135  std::string msg(e.what());
136  ASSERT_TRUE(msg.find("Cannot swap out H+ because it is the charge-balance species") !=
137  std::string::npos)
138  << "Failed with unexpected error message: " << msg;
139  }
140 }
143 {
144  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
146  {"H2O", "H+", ">(s)FeOH", ">(w)FeOH", "Fe++", "HCO3-", "O2(aq)"},
147  {},
148  {},
149  {},
150  {},
151  {},
152  "O2(aq)",
153  "e-");
155 
156  try
157  {
159  _ac3,
160  _is3,
161  _swapper7,
162  {"CO2(aq)"},
163  {"CO2(aq)"},
164  "H+",
165  {},
166  {},
167  _cu_dummy,
168  _cm_dummy,
169  25,
170  0,
171  1E-20,
172  {},
173  {},
174  {});
175  FAIL() << "Missing expected exception.";
176  }
177  catch (const std::exception & e)
178  {
179  std::string msg(e.what());
180  ASSERT_TRUE(msg.find("CO2(aq) is not in the basis, so cannot be removed "
181  "from the basis") != std::string::npos)
182  << "Failed with unexpected error message: " << msg;
183  }
184 }
185 
187 TEST_F(GeochemicalSystemTest, constraintSizeExceptions)
188 {
189  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
191  database, {"H2O", "H+", "HCO3-"}, {}, {}, {}, {}, {}, "O2(aq)", "e-");
193 
194  std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
198  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
202 
203  std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm_poor;
206  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu_poor;
209 
210  try
211  {
213  _ac3,
214  _is3,
215  _swapper3,
216  {},
217  {},
218  "H+",
219  {"H2O", "H+", "HCO3-"},
220  {1.0, 2.0},
221  cu,
222  cm,
223  25,
224  0,
225  1E-20,
226  {},
227  {},
228  {});
229  FAIL() << "Missing expected exception.";
230  }
231  catch (const std::exception & e)
232  {
233  std::string msg(e.what());
234  ASSERT_TRUE(msg.find("Constrained species names must have same length as constraint values") !=
235  std::string::npos)
236  << "Failed with unexpected error message: " << msg;
237  }
238 
239  try
240  {
242  _ac3,
243  _is3,
244  _swapper3,
245  {},
246  {},
247  "H+",
248  {"H2O", "H+", "HCO3-"},
249  {1.0, 2.0, 3.0},
250  cu,
251  cm_poor,
252  25,
253  0,
254  1E-20,
255  {},
256  {},
257  {});
258  FAIL() << "Missing expected exception.";
259  }
260  catch (const std::exception & e)
261  {
262  std::string msg(e.what());
263  ASSERT_TRUE(
264  msg.find("Constrained species names must have same length as constraint meanings") !=
265  std::string::npos)
266  << "Failed with unexpected error message: " << msg;
267  }
268 
269  try
270  {
272  _ac3,
273  _is3,
274  _swapper3,
275  {},
276  {},
277  "H+",
278  {"H2O", "H+", "HCO3-"},
279  {1.0, 2.0, 3.0},
280  cu_poor,
281  cm,
282  25,
283  0,
284  1E-20,
285  {},
286  {},
287  {});
288  FAIL() << "Missing expected exception.";
289  }
290  catch (const std::exception & e)
291  {
292  std::string msg(e.what());
293  ASSERT_TRUE(msg.find("Constrained species names must have same length as constraint units") !=
294  std::string::npos)
295  << "Failed with unexpected error message: " << msg;
296  }
297 
298  try
299  {
301  _ac3,
302  _is3,
303  _swapper3,
304  {},
305  {},
306  "H+",
307  {"H2O", "H+"},
308  {1.0, 2.0},
309  cu_poor,
310  cm_poor,
311  25,
312  0,
313  1E-20,
314  {},
315  {},
316  {});
317  FAIL() << "Missing expected exception.";
318  }
319  catch (const std::exception & e)
320  {
321  std::string msg(e.what());
322  ASSERT_TRUE(
323  msg.find("Constrained species names must have same length as the number of species in the "
324  "basis (each component must be provided with a single constraint") !=
325  std::string::npos)
326  << "Failed with unexpected error message: " << msg;
327  }
328 }
330 TEST_F(GeochemicalSystemTest, constraintNameExceptions)
331 {
332  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
334  database, {"H2O", "H+", "HCO3-"}, {}, {}, {}, {}, {}, "O2(aq)", "e-");
336 
337  std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
341  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
345 
346  try
347  {
349  _ac3,
350  _is3,
351  _swapper3,
352  {},
353  {},
354  "H+",
355  {"H2O", "H+", "H20"},
356  {1.0, 2.0, 3.0},
357  cu,
358  cm,
359  25,
360  0,
361  1E-20,
362  {},
363  {},
364  {});
365  FAIL() << "Missing expected exception.";
366  }
367  catch (const std::exception & e)
368  {
369  std::string msg(e.what());
370  ASSERT_TRUE(msg.find("The basis species HCO3- must appear in the constrained species list") !=
371  std::string::npos)
372  << "Failed with unexpected error message: " << msg;
373  }
374 
375  try
376  {
378  _ac3,
379  _is3,
380  _swapper3,
381  {"HCO3-"},
382  {"CO2(aq)"},
383  "H+",
384  {"H2O", "H+", "HCO3-"},
385  {1.0, 2.0, 3.0},
386  cu,
387  cm,
388  25,
389  0,
390  1E-20,
391  {},
392  {},
393  {});
394  FAIL() << "Missing expected exception.";
395  }
396  catch (const std::exception & e)
397  {
398  std::string msg(e.what());
399  ASSERT_TRUE(msg.find("The basis species CO2(aq) must appear in the constrained species list") !=
400  std::string::npos)
401  << "Failed with unexpected error message: " << msg;
402  }
403 }
404 
406 TEST_F(GeochemicalSystemTest, positiveException1)
407 {
408  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
410  database, {"H2O", "H+", "HCO3-"}, {}, {}, {}, {}, {}, "O2(aq)", "e-");
412 
413  std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
417  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
421 
422  try
423  {
425  _ac3,
426  _is3,
427  _swapper3,
428  {"HCO3-"},
429  {"CO2(aq)"},
430  "H+",
431  {"H2O", "H+", "CO2(aq)"},
432  {-1.0, 2.0, 3.0},
433  cu,
434  cm,
435  25,
436  0,
437  1E-20,
438  {},
439  {},
440  {});
441  FAIL() << "Missing expected exception.";
442  }
443  catch (const std::exception & e)
444  {
445  std::string msg(e.what());
446  ASSERT_TRUE(msg.find("Specified mass of solvent water must be positive: you entered -1") !=
447  std::string::npos)
448  << "Failed with unexpected error message: " << msg;
449  }
450 }
452 TEST_F(GeochemicalSystemTest, positiveException2)
453 {
454  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
456  database, {"H2O", "H+", "HCO3-"}, {}, {}, {}, {}, {}, "O2(aq)", "e-");
458 
459  std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
463  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
467 
468  try
469  {
471  _ac3,
472  _is3,
473  _swapper3,
474  {"HCO3-"},
475  {"CO2(aq)"},
476  "H+",
477  {"H2O", "H+", "CO2(aq)"},
478  {-1.0, 2.0, 3.0},
479  cu,
480  cm,
481  25,
482  0,
483  1E-20,
484  {},
485  {},
486  {});
487  FAIL() << "Missing expected exception.";
488  }
489  catch (const std::exception & e)
490  {
491  std::string msg(e.what());
492  ASSERT_TRUE(msg.find("Specified activity values must be positive: you entered -1") !=
493  std::string::npos)
494  << "Failed with unexpected error message: " << msg;
495  }
496 }
497 
499 TEST_F(GeochemicalSystemTest, positiveException3)
500 {
501  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
503  database, {"H2O", "H+", "HCO3-", "O2(aq)"}, {}, {"O2(g)"}, {}, {}, {}, "O2(aq)", "e-");
505 
506  std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
511  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
516 
517  try
518  {
520  _ac3,
521  _is3,
522  _swapper4,
523  {"O2(aq)"},
524  {"O2(g)"},
525  "H+",
526  {"H2O", "H+", "O2(g)", "HCO3-"},
527  {1.0, 2.0, -3.0, 4.0},
528  cu,
529  cm,
530  25,
531  0,
532  1E-20,
533  {},
534  {},
535  {});
536  FAIL() << "Missing expected exception.";
537  }
538  catch (const std::exception & e)
539  {
540  std::string msg(e.what());
541  ASSERT_TRUE(msg.find("Specified fugacity values must be positive: you entered -3") !=
542  std::string::npos)
543  << "Failed with unexpected error message: " << msg;
544  }
545 }
546 
548 TEST_F(GeochemicalSystemTest, positiveException4)
549 {
550  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
552  database, {"H2O", "H+", "HCO3-"}, {}, {}, {}, {}, {}, "O2(aq)", "e-");
554 
555  std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
559  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
563 
564  try
565  {
567  _ac3,
568  _is3,
569  _swapper3,
570  {"HCO3-"},
571  {"CO2(aq)"},
572  "H+",
573  {"H2O", "H+", "CO2(aq)"},
574  {1.0, 2.0, -2.0},
575  cu,
576  cm,
577  25,
578  0,
579  1E-20,
580  {},
581  {},
582  {});
583  FAIL() << "Missing expected exception.";
584  }
585  catch (const std::exception & e)
586  {
587  std::string msg(e.what());
588  ASSERT_TRUE(msg.find("Specified free concentration values must be positive: you entered -2") !=
589  std::string::npos)
590  << "Failed with unexpected error message: " << msg;
591  }
592 }
593 
595 TEST_F(GeochemicalSystemTest, positiveException5)
596 {
597  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
599  database, {"H2O", "H+", "HCO3-", "Ca++"}, {"Calcite"}, {}, {}, {}, {}, "O2(aq)", "e-");
601 
602  std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
607  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
612 
613  try
614  {
616  _ac3,
617  _is3,
618  _swapper4,
619  {"Ca++"},
620  {"Calcite"},
621  "H+",
622  {"H2O", "Calcite", "H+", "HCO3-"},
623  {1.0, -3.0, 2.0, 1.0},
624  cu,
625  cm,
626  25,
627  0,
628  1E-20,
629  {},
630  {},
631  {});
632  FAIL() << "Missing expected exception.";
633  }
634  catch (const std::exception & e)
635  {
636  std::string msg(e.what());
637  ASSERT_TRUE(msg.find("Specified free mineral values must be positive: you entered -3") !=
638  std::string::npos)
639  << "Failed with unexpected error message: " << msg;
640  }
641 }
642 
645 {
646  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
648  database, {"H2O", "H+", "HCO3-"}, {}, {}, {}, {}, {}, "O2(aq)", "e-");
650 
651  std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
655  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
659 
660  try
661  {
663  _ac3,
664  _is3,
665  _swapper3,
666  {"HCO3-"},
667  {"CO2(aq)"},
668  "H+",
669  {"H2O", "H+", "CO2(aq)"},
670  {1.0, 2.0, 3.0},
671  cu,
672  cm,
673  25,
674  0,
675  1E-20,
676  {},
677  {},
678  {});
679  FAIL() << "Missing expected exception.";
680  }
681  catch (const std::exception & e)
682  {
683  std::string msg(e.what());
684  ASSERT_TRUE(msg.find("H2O must be provided with either a mass of solvent water, a bulk "
685  "composition in moles or mass, or an activity") != std::string::npos)
686  << "Failed with unexpected error message: " << msg;
687  }
688 }
689 
691 TEST_F(GeochemicalSystemTest, fugacityException)
692 {
693  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
695  database, {"H2O", "H+", "HCO3-", "O2(aq)"}, {}, {"O2(g)"}, {}, {}, {}, "O2(aq)", "e-");
697 
698  std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
703  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
708 
709  try
710  {
712  _ac3,
713  _is3,
714  _swapper4,
715  {"O2(aq)"},
716  {"O2(g)"},
717  "H+",
718  {"H2O", "H+", "O2(g)", "HCO3-"},
719  {1.0, 2.0, 3.0, 4.0},
720  cu,
721  cm,
722  25,
723  0,
724  1E-20,
725  {},
726  {},
727  {});
728  FAIL() << "Missing expected exception.";
729  }
730  catch (const std::exception & e)
731  {
732  std::string msg(e.what());
733  ASSERT_TRUE(msg.find("The gas O2(g) must be provided with a fugacity") != std::string::npos)
734  << "Failed with unexpected error message: " << msg;
735  }
736 }
737 
739 TEST_F(GeochemicalSystemTest, mineralMeaningExecption)
740 {
741  std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
746  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
751 
752  try
753  {
754  GeochemicalSystem egs(_mgd_calcite,
755  _ac3,
756  _is3,
757  _swapper4,
758  {},
759  {},
760  "H+",
761  {"H2O", "Calcite", "H+", "HCO3-"},
762  {1.0, 3.0, 2.0, 1.0},
763  cu,
764  cm,
765  25,
766  0,
767  1E-20,
768  {},
769  {},
770  {});
771  FAIL() << "Missing expected exception.";
772  }
773  catch (const std::exception & e)
774  {
775  std::string msg(e.what());
776  ASSERT_TRUE(
777  msg.find("The mineral Calcite must be provided with either: a free number of moles, a free "
778  "mass or a free volume; or a bulk composition of moles or mass") !=
779  std::string::npos)
780  << "Failed with unexpected error message: " << msg;
781  }
782 
783  std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm2;
788  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu2;
793 
794  try
795  {
796  GeochemicalSystem egs(_mgd_calcite,
797  _ac3,
798  _is3,
799  _swapper4,
800  {},
801  {},
802  "H+",
803  {"H2O", "Calcite", "H+", "HCO3-"},
804  {1.0, 3.0, 2.0, 1.0},
805  cu2,
806  cm2,
807  25,
808  0,
809  1E-20,
810  {},
811  {},
812  {});
813  FAIL() << "Missing expected exception.";
814  }
815  catch (const std::exception & e)
816  {
817  std::string msg(e.what());
818  ASSERT_TRUE(
819  msg.find("The mineral Calcite must be provided with either: a free number of "
820  "moles, a free mass or a free volume; or a bulk composition of moles or mass") !=
821  std::string::npos)
822  << "Failed with unexpected error message: " << msg;
823  }
824 }
825 
827 TEST_F(GeochemicalSystemTest, aqueousSpeciesMeaningExecption)
828 {
829  std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
834  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
839 
840  try
841  {
842  GeochemicalSystem egs(_mgd_calcite,
843  _ac3,
844  _is3,
845  _swapper4,
846  {},
847  {},
848  "HCO3-",
849  {"H2O", "Calcite", "H+", "HCO3-"},
850  {1.0, 3.0, 2.0, 1.0},
851  cu,
852  cm,
853  25,
854  0,
855  1E-20,
856  {},
857  {},
858  {});
859  FAIL() << "Missing expected exception.";
860  }
861  catch (const std::exception & e)
862  {
863  std::string msg(e.what());
864  ASSERT_TRUE(msg.find("The basis species H+ must be provided with a free concentration, bulk "
865  "composition or an activity") != std::string::npos)
866  << "Failed with unexpected error message: " << msg;
867  }
868 
869  std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm2;
874  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu2;
879 
880  try
881  {
882  GeochemicalSystem egs(_mgd_calcite,
883  _ac3,
884  _is3,
885  _swapper4,
886  {},
887  {},
888  "H+",
889  {"H2O", "Calcite", "H+", "HCO3-"},
890  {1.0, 3.0, 2.0, 1.0},
891  cu2,
892  cm2,
893  25,
894  0,
895  1E-20,
896  {},
897  {},
898  {});
899  FAIL() << "Missing expected exception.";
900  }
901  catch (const std::exception & e)
902  {
903  std::string msg(e.what());
904  ASSERT_TRUE(msg.find("The basis species HCO3- must be provided with a free concentration, bulk "
905  "composition or an activity") != std::string::npos)
906  << "Failed with unexpected error message: " << msg;
907  }
908 }
909 
911 TEST_F(GeochemicalSystemTest, waterUnitsException)
912 {
913  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
915  database, {"H2O", "H+", "HCO3-"}, {}, {}, {}, {}, {}, "O2(aq)", "e-");
917 
918  std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
922  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
926 
927  try
928  {
930  _ac3,
931  _is3,
932  _swapper3,
933  {},
934  {},
935  "H+",
936  {"H2O", "H+", "HCO3-"},
937  {1.0, 2.0, 3.0},
938  cu,
939  cm,
940  25,
941  0,
942  1E-20,
943  {},
944  {},
945  {});
946  FAIL() << "Missing expected exception.";
947  }
948  catch (const std::exception & e)
949  {
950  std::string msg(e.what());
951  ASSERT_TRUE(msg.find("Units for kg_solvent_water must be kg") != std::string::npos)
952  << "Failed with unexpected error message: " << msg;
953  }
954 }
955 
957 TEST_F(GeochemicalSystemTest, bulkUnitsException)
958 {
959  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
961  database, {"H2O", "H+", "HCO3-"}, {}, {}, {}, {}, {}, "O2(aq)", "e-");
963 
964  std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
968  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
972 
973  try
974  {
976  _ac3,
977  _is3,
978  _swapper3,
979  {},
980  {},
981  "H+",
982  {"H2O", "H+", "HCO3-"},
983  {1.0, 2.0, 3.0},
984  cu,
985  cm,
986  25,
987  0,
988  1E-20,
989  {},
990  {},
991  {});
992  FAIL() << "Missing expected exception.";
993  }
994  catch (const std::exception & e)
995  {
996  std::string msg(e.what());
997  ASSERT_TRUE(msg.find("Species H+: units for bulk composition must be moles or mass") !=
998  std::string::npos)
999  << "Failed with unexpected error message: " << msg;
1000  }
1001 }
1002 
1004 TEST_F(GeochemicalSystemTest, freeConcUnitsException)
1005 {
1006  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
1008  database, {"H2O", "H+", "HCO3-"}, {}, {}, {}, {}, {}, "O2(aq)", "e-");
1010 
1011  std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
1015  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
1019 
1020  try
1021  {
1022  GeochemicalSystem egs(mgd,
1023  _ac3,
1024  _is3,
1025  _swapper3,
1026  {},
1027  {},
1028  "H+",
1029  {"H2O", "H+", "HCO3-"},
1030  {1.0, 2.0, 3.0},
1031  cu,
1032  cm,
1033  25,
1034  0,
1035  1E-20,
1036  {},
1037  {},
1038  {});
1039  FAIL() << "Missing expected exception.";
1040  }
1041  catch (const std::exception & e)
1042  {
1043  std::string msg(e.what());
1044  ASSERT_TRUE(msg.find("Species HCO3-: units for free concentration quantities must be molal or "
1045  "mass_per_kg_solvent") != std::string::npos)
1046  << "Failed with unexpected error message: " << msg;
1047  }
1048 }
1049 
1051 TEST_F(GeochemicalSystemTest, freeMineralUnitsException)
1052 {
1053  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
1055  database, {"H2O", "H+", "HCO3-", "Ca++"}, {"Calcite"}, {}, {}, {}, {}, "O2(aq)", "e-");
1057 
1058  std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
1063  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
1068 
1069  try
1070  {
1071  GeochemicalSystem egs(mgd,
1072  _ac3,
1073  _is3,
1074  _swapper4,
1075  {"Ca++"},
1076  {"Calcite"},
1077  "H+",
1078  {"H2O", "Calcite", "H+", "HCO3-"},
1079  {1.0, 3.0, 2.0, 1.0},
1080  cu,
1081  cm,
1082  25,
1083  0,
1084  1E-20,
1085  {},
1086  {},
1087  {});
1088  FAIL() << "Missing expected exception.";
1089  }
1090  catch (const std::exception & e)
1091  {
1092  std::string msg(e.what());
1093  ASSERT_TRUE(
1094  msg.find(
1095  "Species Calcite: units for free mineral quantities must be moles, mass or volume") !=
1096  std::string::npos)
1097  << "Failed with unexpected error message: " << msg;
1098  }
1099 }
1100 
1102 TEST_F(GeochemicalSystemTest, activityUnitsException)
1103 {
1104  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
1106  database, {"H2O", "H+", "HCO3-"}, {}, {}, {}, {}, {}, "O2(aq)", "e-");
1108 
1109  std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
1113  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
1117 
1118  try
1119  {
1120  GeochemicalSystem egs(mgd,
1121  _ac3,
1122  _is3,
1123  _swapper3,
1124  {},
1125  {},
1126  "H+",
1127  {"H2O", "H+", "HCO3-"},
1128  {1.0, 2.0, 3.0},
1129  cu,
1130  cm,
1131  25,
1132  0,
1133  1E-20,
1134  {},
1135  {},
1136  {});
1137  FAIL() << "Missing expected exception.";
1138  }
1139  catch (const std::exception & e)
1140  {
1141  std::string msg(e.what());
1142  ASSERT_TRUE(
1143  msg.find("Species HCO3-: dimensionless units must be used when specifying activity") !=
1144  std::string::npos)
1145  << "Failed with unexpected error message: " << msg;
1146  }
1147 
1149  try
1150  {
1151  GeochemicalSystem egs(mgd,
1152  _ac3,
1153  _is3,
1154  _swapper3,
1155  {},
1156  {},
1157  "H+",
1158  {"H2O", "H+", "HCO3-"},
1159  {1.0, 2.0, 3.0},
1160  cu,
1161  cm,
1162  25,
1163  0,
1164  1E-20,
1165  {},
1166  {},
1167  {});
1168  FAIL() << "Missing expected exception.";
1169  }
1170  catch (const std::exception & e)
1171  {
1172  std::string msg(e.what());
1173  ASSERT_TRUE(
1174  msg.find("Species HCO3-: dimensionless units must be used when specifying log10activity") !=
1175  std::string::npos)
1176  << "Failed with unexpected error message: " << msg;
1177  }
1178 }
1179 
1181 TEST_F(GeochemicalSystemTest, fugacityUnitsException)
1182 {
1183  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
1185  database, {"H2O", "H+", "HCO3-", "O2(aq)"}, {}, {"O2(g)"}, {}, {}, {}, "O2(aq)", "e-");
1187 
1188  std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
1193  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
1198 
1199  try
1200  {
1201  GeochemicalSystem egs(mgd,
1202  _ac3,
1203  _is3,
1204  _swapper4,
1205  {"O2(aq)"},
1206  {"O2(g)"},
1207  "H+",
1208  {"H2O", "H+", "O2(g)", "HCO3-"},
1209  {1.0, 2.0, 3.0, 4.0},
1210  cu,
1211  cm,
1212  25,
1213  0,
1214  1E-20,
1215  {},
1216  {},
1217  {});
1218  FAIL() << "Missing expected exception.";
1219  }
1220  catch (const std::exception & e)
1221  {
1222  std::string msg(e.what());
1223  ASSERT_TRUE(
1224  msg.find("Species O2(g): dimensionless units must be used when specifying fugacity") !=
1225  std::string::npos)
1226  << "Failed with unexpected error message: " << msg;
1227  }
1228 
1230  try
1231  {
1232  // no need to swap, since the above constructor of egs will have performed the swaps in mgd
1233  GeochemicalSystem egs(mgd,
1234  _ac3,
1235  _is3,
1236  _swapper4,
1237  {},
1238  {},
1239  "H+",
1240  {"H2O", "H+", "O2(g)", "HCO3-"},
1241  {1.0, 2.0, 3.0, 4.0},
1242  cu,
1243  cm,
1244  25,
1245  0,
1246  1E-20,
1247  {},
1248  {},
1249  {});
1250  FAIL() << "Missing expected exception.";
1251  }
1252  catch (const std::exception & e)
1253  {
1254  std::string msg(e.what());
1255  ASSERT_TRUE(
1256  msg.find("Species O2(g): dimensionless units must be used when specifying log10fugacity") !=
1257  std::string::npos)
1258  << "Failed with unexpected error message: " << msg;
1259  }
1260 }
1261 
1263 TEST_F(GeochemicalSystemTest, unitsConversion1)
1264 {
1265  std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
1270  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
1275 
1276  const GeochemicalSystem egs(_mgd_calcite,
1277  _ac3,
1278  _is3,
1279  _swapper4,
1280  {},
1281  {},
1282  "H+",
1283  {"H2O", "H+", "HCO3-", "Calcite"},
1284  {100.0, 3.0, 1.0, 2.0},
1285  cu,
1286  cm,
1287  25,
1288  0,
1289  1E-20,
1290  {},
1291  {},
1292  {});
1293  const std::vector<GeochemicalSystem::ConstraintMeaningEnum> cim = {
1298  for (unsigned i = 0; i < 4; ++i)
1299  EXPECT_EQ(egs.getConstraintMeaning()[i], cim[i]);
1300  EXPECT_NEAR(egs.getBulkMolesOld()[0], 100.0 / 18.0152, 1.0E-6); // H2O
1301  EXPECT_NEAR(egs.getBulkMolesOld()[1], 3.0, 1.0E-6); // H+
1302  EXPECT_NEAR(
1303  egs.getSolventMassAndFreeMolalityAndMineralMoles()[2], 1.0 / 61.0171, 1.0E-6); // HCO3-
1304  EXPECT_NEAR(
1305  egs.getSolventMassAndFreeMolalityAndMineralMoles()[3], 2.0 / 36.9340, 1.0E-6); // Calcite
1306 }
1307 
1309 TEST_F(GeochemicalSystemTest, unitsConversion2)
1310 {
1311  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
1313  database, {"H2O", "H+", "HCO3-", "O2(aq)"}, {}, {"O2(g)"}, {}, {}, {}, "O2(aq)", "e-");
1315 
1316  std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
1321  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
1326 
1327  const GeochemicalSystem egs(mgd,
1328  _ac3,
1329  _is3,
1330  _swapper4,
1331  {"O2(aq)"},
1332  {"O2(g)"},
1333  "HCO3-",
1334  {"H2O", "H+", "HCO3-", "O2(g)"},
1335  {1.0, 2.0, 3.0, 4.0},
1336  cu,
1337  cm,
1338  25,
1339  0,
1340  1E-20,
1341  {},
1342  {},
1343  {});
1344  const std::vector<GeochemicalSystem::ConstraintMeaningEnum> cim = {
1349  for (unsigned i = 0; i < 4; ++i)
1350  EXPECT_EQ(egs.getConstraintMeaning()[i], cim[i]);
1351  EXPECT_NEAR(egs.getSolventMassAndFreeMolalityAndMineralMoles()[0], 1.0, 1.0E-6); // H2O
1352  EXPECT_NEAR(egs.getBasisActivity(1), 1E2, 1.0E-6); // H+
1353  EXPECT_NEAR(egs.getBulkMolesOld()[2], 2.0E-6 / 61.0171, 1.0E-6); // HCO3-
1354  EXPECT_NEAR(egs.getBasisActivity(3), 1E4, 1.0E-6); // O2(g)
1355 }
1356 
1358 TEST_F(GeochemicalSystemTest, chargeBalanceMeaningExecption)
1359 {
1360  std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
1365  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
1370 
1371  try
1372  {
1373  GeochemicalSystem egs(_mgd_calcite,
1374  _ac3,
1375  _is3,
1376  _swapper4,
1377  {},
1378  {},
1379  "H+",
1380  {"H2O", "Calcite", "H+", "HCO3-"},
1381  {1.0, 3.0, 2.0, 1.0},
1382  cu,
1383  cm,
1384  25,
1385  0,
1386  1E-20,
1387  {},
1388  {},
1389  {});
1390  FAIL() << "Missing expected exception.";
1391  }
1392  catch (const std::exception & e)
1393  {
1394  std::string msg(e.what());
1395  ASSERT_TRUE(
1396  msg.find(
1397  "For code consistency, the species H+ must be provided with a bulk composition "
1398  "because it is the charge-balance species. The value provided should be a reasonable "
1399  "estimate of the mole number, but will be overridden as the solve progresses") !=
1400  std::string::npos)
1401  << "Failed with unexpected error message: " << msg;
1402  }
1403 
1404  std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm2;
1409  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu2;
1414 
1415  try
1416  {
1417  GeochemicalSystem egs(_mgd_calcite,
1418  _ac3,
1419  _is3,
1420  _swapper4,
1421  {},
1422  {},
1423  "H+",
1424  {"H2O", "Calcite", "H+", "HCO3-"},
1425  {1.0, 3.0, 2.0, 1.0},
1426  cu2,
1427  cm2,
1428  25,
1429  0,
1430  1E-20,
1431  {},
1432  {},
1433  {});
1434  FAIL() << "Missing expected exception.";
1435  }
1436  catch (const std::exception & e)
1437  {
1438  std::string msg(e.what());
1439  ASSERT_TRUE(
1440  msg.find(
1441  "For code consistency, the species H+ must be provided with a bulk composition "
1442  "because it is the charge-balance species. The value provided should be a reasonable "
1443  "estimate of the mole number, but will be overridden as the solve progresses") !=
1444  std::string::npos)
1445  << "Failed with unexpected error message: " << msg;
1446  }
1447 }
1448 
1450 TEST_F(GeochemicalSystemTest, chargeBalanceInBasisExecption)
1451 {
1452  try
1453  {
1454  GeochemicalSystem egs(_mgd_calcite,
1455  _ac3,
1456  _is3,
1457  _swapper4,
1458  {},
1459  {},
1460  "OH-",
1461  {"H2O", "Calcite", "H+", "HCO3-"},
1462  {1.0, 3.0, 2.0, 1.0},
1463  _cu_calcite,
1464  _cm_calcite,
1465  25,
1466  0,
1467  1E-20,
1468  {},
1469  {},
1470  {});
1471  FAIL() << "Missing expected exception.";
1472  }
1473  catch (const std::exception & e)
1474  {
1475  std::string msg(e.what());
1476  ASSERT_TRUE(
1477  msg.find("Cannot enforce charge balance using OH- because it is not in the basis") !=
1478  std::string::npos)
1479  << "Failed with unexpected error message: " << msg;
1480  }
1481 }
1482 
1484 TEST_F(GeochemicalSystemTest, chargeBalanceHasChargeExecption)
1485 {
1486  try
1487  {
1488  GeochemicalSystem egs(_mgd_calcite,
1489  _ac3,
1490  _is3,
1491  _swapper4,
1492  {},
1493  {},
1494  "Calcite",
1495  {"H2O", "Calcite", "H+", "HCO3-"},
1496  {1.0, 3.0, 2.0, 1.0},
1497  _cu_calcite,
1498  _cm_calcite,
1499  25,
1500  0,
1501  1E-20,
1502  {},
1503  {},
1504  {});
1505  FAIL() << "Missing expected exception.";
1506  }
1507  catch (const std::exception & e)
1508  {
1509  std::string msg(e.what());
1510  ASSERT_TRUE(
1511  msg.find("Cannot enforce charge balance using Calcite because it has zero charge") !=
1512  std::string::npos)
1513  << "Failed with unexpected error message: " << msg;
1514  }
1515 }
1516 
1518 TEST_F(GeochemicalSystemTest, chargeBalanceIndex)
1519 {
1520  EXPECT_EQ(_egs_calcite.getChargeBalanceBasisIndex(), (unsigned)1);
1521 }
1522 
1525 {
1526  try
1527  {
1528  EXPECT_EQ(_egs_calcite.getLog10K(123), 1);
1529  FAIL() << "Missing expected exception.";
1530  }
1531  catch (const std::exception & e)
1532  {
1533  std::string msg(e.what());
1534  ASSERT_TRUE(msg.find("Cannot retrieve log10K for equilibrium species 123 since there are only "
1535  "6 equilibrium species") != std::string::npos)
1536  << "Failed with unexpected error message: " << msg;
1537  }
1538  EXPECT_EQ(_egs_calcite.getLog10K(0), -6.3660); // CO2(aq) at 25degC
1539 }
1540 
1542 TEST_F(GeochemicalSystemTest, getInAlgebraicSystem)
1543 {
1544  const std::vector<bool> as_gold = {false, true, false, false};
1545  EXPECT_EQ(_egs_calcite.getInAlgebraicSystem().size(), (std::size_t)4);
1546  for (unsigned i = 0; i < as_gold.size(); ++i)
1547  EXPECT_EQ(_egs_calcite.getInAlgebraicSystem()[i], as_gold[i]);
1548 }
1549 
1551 TEST_F(GeochemicalSystemTest, NumInAlgebraicSystem)
1552 {
1553  EXPECT_EQ(_egs_calcite.getNumInAlgebraicSystem(), (unsigned)1);
1554  EXPECT_EQ(_egs_calcite.getNumBasisInAlgebraicSystem(), (unsigned)1);
1555  EXPECT_EQ(_egs_calcite.getNumSurfacePotentials(), (unsigned)0);
1556  EXPECT_EQ(_egs_kinetic_calcite.getNumInAlgebraicSystem(), (unsigned)3);
1557  EXPECT_EQ(_egs_kinetic_calcite.getNumBasisInAlgebraicSystem(), (unsigned)2);
1558  EXPECT_EQ(_egs_kinetic_calcite.getNumSurfacePotentials(), (unsigned)0);
1559 }
1560 
1562 TEST_F(GeochemicalSystemTest, getBasisIndexOfAlgebraicSystem)
1563 {
1564  EXPECT_EQ(_egs_calcite.getBasisIndexOfAlgebraicSystem()[0], (unsigned)1);
1565  EXPECT_EQ(_egs_kinetic_calcite.getBasisIndexOfAlgebraicSystem()[0], (unsigned)1);
1566  EXPECT_EQ(_egs_kinetic_calcite.getBasisIndexOfAlgebraicSystem()[1], (unsigned)3);
1567 }
1568 
1570 TEST_F(GeochemicalSystemTest, getAlgebraicIndexOfBasisSystem)
1571 {
1572  EXPECT_EQ(_egs_calcite.getAlgebraicIndexOfBasisSystem()[1], (unsigned)0);
1573  EXPECT_EQ(_egs_kinetic_calcite.getAlgebraicIndexOfBasisSystem()[1], (unsigned)0);
1574  EXPECT_EQ(_egs_kinetic_calcite.getAlgebraicIndexOfBasisSystem()[3], (unsigned)1);
1575 }
1576 
1578 TEST_F(GeochemicalSystemTest, getSolventWaterMass)
1579 {
1580  EXPECT_EQ(_egs_calcite.getSolventWaterMass(), 1.0);
1581 
1582  const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
1587  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
1592  const GeochemicalSystem egs(_mgd_calcite,
1593  _ac3,
1594  _is3,
1595  _swapper4,
1596  {},
1597  {},
1598  "H+",
1599  {"H2O", "Calcite", "H+", "HCO3-"},
1600  {2.5, 3.0, 2.0, 1.0},
1601  cu,
1602  cm,
1603  25,
1604  0,
1605  1E-20,
1606  {},
1607  {},
1608  {});
1609  EXPECT_EQ(egs.getSolventWaterMass(), 2.5);
1610 }
1611 
1614 {
1615  EXPECT_EQ(_egs_calcite.getBulkMolesOld()[3], 3.0); // Calcite
1616  EXPECT_EQ(_egs_calcite.getBulkMolesOld()[1], 2.0); // H+
1617 
1618  const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
1623  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
1628  const GeochemicalSystem egs(_mgd_calcite,
1629  _ac3,
1630  _is3,
1631  _swapper4,
1632  {},
1633  {},
1634  "H+",
1635  {"H2O", "Calcite", "H+", "HCO3-"},
1636  {-0.5, 3.5, 2.5, 1.0},
1637  cu,
1638  cm,
1639  25,
1640  0,
1641  1E-20,
1642  {},
1643  {},
1644  {});
1645  EXPECT_EQ(egs.getBulkMolesOld()[0], -0.5); // H2O
1646  EXPECT_EQ(egs.getBulkMolesOld()[1], 1.0); // H+ : note that charge neutrality has forced this
1647  EXPECT_EQ(egs.getBulkMolesOld()[2], 1.0); // HCO3-
1648  EXPECT_EQ(egs.getBulkMolesOld()[3], 3.5); // Calcite
1649 }
1650 
1652 TEST_F(GeochemicalSystemTest, getSolventMassAndFreeMolalityAndMineralMoles)
1653 {
1654  EXPECT_EQ(_egs_calcite.getSolventMassAndFreeMolalityAndMineralMoles()[2], 1.0); // HCO3-
1655  EXPECT_EQ(_egs_calcite.getAlgebraicVariableValues()[0],
1656  _egs_calcite.getSolventMassAndFreeMolalityAndMineralMoles()[1]); // H+
1657  EXPECT_EQ(_egs_calcite.getAlgebraicBasisValues()[0],
1658  _egs_calcite.getSolventMassAndFreeMolalityAndMineralMoles()[1]); // H+
1659  EXPECT_EQ(_egs_calcite.getAlgebraicVariableDenseValues()(0),
1660  _egs_calcite.getSolventMassAndFreeMolalityAndMineralMoles()[1]); // H+
1661 
1662  EXPECT_EQ(_egs_kinetic_calcite.getSolventMassAndFreeMolalityAndMineralMoles()[2], 1.0); // HCO3-
1663  EXPECT_EQ(_egs_kinetic_calcite.getAlgebraicVariableValues()[0],
1664  _egs_kinetic_calcite.getSolventMassAndFreeMolalityAndMineralMoles()[1]); // H+
1665  EXPECT_EQ(_egs_kinetic_calcite.getAlgebraicBasisValues()[0],
1666  _egs_kinetic_calcite.getSolventMassAndFreeMolalityAndMineralMoles()[1]); // H+
1667  EXPECT_EQ(_egs_kinetic_calcite.getAlgebraicVariableDenseValues()(0),
1668  _egs_kinetic_calcite.getSolventMassAndFreeMolalityAndMineralMoles()[1]); // H+
1669  EXPECT_EQ(_egs_kinetic_calcite.getAlgebraicVariableValues()[1],
1670  _egs_kinetic_calcite.getSolventMassAndFreeMolalityAndMineralMoles()[3]); // Ca++
1671  EXPECT_EQ(_egs_kinetic_calcite.getAlgebraicBasisValues()[1],
1672  _egs_kinetic_calcite.getSolventMassAndFreeMolalityAndMineralMoles()[3]); // Ca++
1673  EXPECT_EQ(_egs_kinetic_calcite.getAlgebraicVariableDenseValues()(1),
1674  _egs_kinetic_calcite.getSolventMassAndFreeMolalityAndMineralMoles()[3]); // Ca++
1675  EXPECT_EQ(_egs_kinetic_calcite.getAlgebraicVariableValues()[2], 1.1); // Calcite
1676  EXPECT_EQ(_egs_kinetic_calcite.getAlgebraicVariableDenseValues()(2), 1.1); // Calcite
1677 
1678  const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
1683  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
1688  const GeochemicalSystem egs(_mgd_calcite,
1689  _ac3,
1690  _is3,
1691  _swapper4,
1692  {},
1693  {},
1694  "H+",
1695  {"H2O", "Calcite", "H+", "HCO3-"},
1696  {0.5, 3.5, 2.5, 1.0},
1697  cu,
1698  cm,
1699  25,
1700  0,
1701  1E-20,
1702  {},
1703  {},
1704  {});
1705  EXPECT_EQ(egs.getSolventMassAndFreeMolalityAndMineralMoles()[0], 0.5); // H2O
1706  EXPECT_EQ(egs.getSolventMassAndFreeMolalityAndMineralMoles()[2], 1.0); // HCO3-
1707  EXPECT_EQ(egs.getSolventMassAndFreeMolalityAndMineralMoles()[3], 3.5); // Calcite
1708  EXPECT_EQ(_egs_calcite.getAlgebraicVariableValues()[0],
1709  _egs_calcite.getSolventMassAndFreeMolalityAndMineralMoles()[1]); // H+
1710  EXPECT_EQ(_egs_calcite.getAlgebraicBasisValues()[0],
1711  _egs_calcite.getSolventMassAndFreeMolalityAndMineralMoles()[1]); // H+
1712  EXPECT_EQ(_egs_calcite.getAlgebraicVariableDenseValues()(0),
1713  _egs_calcite.getSolventMassAndFreeMolalityAndMineralMoles()[1]); // H+
1714 }
1715 
1717 TEST_F(GeochemicalSystemTest, getBasisActivityKnown)
1718 {
1719  EXPECT_EQ(_egs_calcite.getBasisActivityKnown()[0], true); // H2O
1720  EXPECT_EQ(_egs_calcite.getBasisActivityKnown()[1], false); // H+
1721  EXPECT_EQ(_egs_calcite.getBasisActivityKnown()[2], false); // HCO3-
1722  EXPECT_EQ(_egs_calcite.getBasisActivityKnown()[3], true); // Calcite
1723 
1724  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
1726  database, {"H2O", "H+", "HCO3-", "O2(aq)"}, {}, {"O2(g)"}, {}, {}, {}, "O2(aq)", "e-");
1728 
1729  std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
1734  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
1739 
1740  GeochemicalSystem egs(mgd,
1741  _ac3,
1742  _is3,
1743  _swapper4,
1744  {"O2(aq)"},
1745  {"O2(g)"},
1746  "H+",
1747  {"H2O", "H+", "O2(g)", "HCO3-"},
1748  {1.0, 2.0, 3.0, 4.0},
1749  cu,
1750  cm,
1751  25,
1752  0,
1753  1E-20,
1754  {},
1755  {},
1756  {});
1757 
1758  EXPECT_EQ(egs.getBasisActivityKnown()[0], false); // H2O
1759  EXPECT_EQ(egs.getBasisActivityKnown()[1], false); // H+
1760  EXPECT_EQ(egs.getBasisActivityKnown()[2], false); // HCO3-
1761  EXPECT_EQ(egs.getBasisActivityKnown()[3], true); // O2(g)
1762 }
1763 
1765 TEST_F(GeochemicalSystemTest, getBasisActivity)
1766 {
1767  EXPECT_EQ(_egs_calcite.getBasisActivity(0), 1.75); // H2O
1768  EXPECT_EQ(_egs_calcite.getBasisActivity(3), 1.0); // Calcite
1769  EXPECT_EQ(_egs_calcite.getBasisActivity()[0], 1.75); // H2O
1770  EXPECT_EQ(_egs_calcite.getBasisActivity()[3], 1.0); // Calcite
1771  try
1772  {
1773  EXPECT_EQ(_egs_calcite.getBasisActivity(4), 1);
1774  FAIL() << "Missing expected exception.";
1775  }
1776  catch (const std::exception & e)
1777  {
1778  std::string msg(e.what());
1779  ASSERT_TRUE(msg.find("Cannot retrieve basis activity for species 4 since there are only 4 "
1780  "basis species") != std::string::npos)
1781  << "Failed with unexpected error message: " << msg;
1782  }
1783 
1784  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
1786  database, {"H2O", "H+", "HCO3-", "O2(aq)"}, {}, {"O2(g)"}, {}, {}, {}, "O2(aq)", "e-");
1788 
1789  std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
1794  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
1799 
1800  GeochemicalSystem egs(mgd,
1801  _ac0,
1802  _is0,
1803  _swapper4,
1804  {"O2(aq)"},
1805  {"O2(g)"},
1806  "H+",
1807  {"H2O", "H+", "O2(g)", "HCO3-"},
1808  {1.0, 1.5, 3.0, 2.5},
1809  cu,
1810  cm,
1811  25,
1812  0,
1813  1E-20,
1814  {},
1815  {},
1816  {}); // maximum ionic_strength = 0, so all activity_coefficients = 1
1817 
1818  EXPECT_EQ(egs.getBasisActivity(1), egs.getSolventMassAndFreeMolalityAndMineralMoles()[1]); // H+
1819  EXPECT_EQ(egs.getBasisActivity(2), 2.5); // HCO3-
1820  EXPECT_EQ(egs.getBasisActivity(3), 3.0); // O2(g)
1821  EXPECT_EQ(egs.getBasisActivity()[2], 2.5); // HCO3-
1822  EXPECT_EQ(egs.getBasisActivity()[3], 3.0); // O2(g)
1823 }
1824 
1826 TEST_F(GeochemicalSystemTest, getBasisActivityCoeff)
1827 {
1828  try
1829  {
1830  EXPECT_EQ(_egs_calcite.getBasisActivityCoefficient(4), 1);
1831  FAIL() << "Missing expected exception.";
1832  }
1833  catch (const std::exception & e)
1834  {
1835  std::string msg(e.what());
1836  ASSERT_TRUE(
1837  msg.find("Cannot retrieve basis activity coefficient for species 4 since there are only 4 "
1838  "basis species") != std::string::npos)
1839  << "Failed with unexpected error message: " << msg;
1840  }
1841 
1842  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
1844  database, {"H2O", "H+", "HCO3-", "O2(aq)"}, {}, {"O2(g)"}, {}, {}, {}, "O2(aq)", "e-");
1846 
1847  std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
1852  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
1857 
1858  GeochemicalSystem egs(mgd,
1859  _ac0,
1860  _is0,
1861  _swapper4,
1862  {"O2(aq)"},
1863  {"O2(g)"},
1864  "H+",
1865  {"H2O", "H+", "O2(g)", "HCO3-"},
1866  {1.0, 1.5, 3.0, 2.5},
1867  cu,
1868  cm,
1869  25,
1870  0,
1871  1E-20,
1872  {},
1873  {},
1874  {}); // maximum ionic_strength = 0, so all activity_coefficients = 1
1875 
1876  EXPECT_EQ(egs.getBasisActivityCoefficient(1), 1.0); // H+
1877  EXPECT_EQ(egs.getBasisActivityCoefficient(2), 1.0); // HCO3-
1878  EXPECT_EQ(egs.getBasisActivityCoefficient()[1], 1.0); // H+
1879  EXPECT_EQ(egs.getBasisActivityCoefficient()[2], 1.0); // HCO3-
1880 }
1881 
1883 TEST_F(GeochemicalSystemTest, getEqmActivityCoefficientException)
1884 {
1885  try
1886  {
1887  EXPECT_EQ(_egs_calcite.getEquilibriumActivityCoefficient(44), 1);
1888  FAIL() << "Missing expected exception.";
1889  }
1890  catch (const std::exception & e)
1891  {
1892  std::string msg(e.what());
1893  ASSERT_TRUE(msg.find("Cannot retrieve activity coefficient for equilibrium species 44 since "
1894  "there are only 6 "
1895  "equilibrium species") != std::string::npos)
1896  << "Failed with unexpected error message: " << msg;
1897  }
1898 }
1899 
1902 {
1903  try
1904  {
1905  EXPECT_EQ(_egs_calcite.getEquilibriumMolality(44), 1);
1906  FAIL() << "Missing expected exception.";
1907  }
1908  catch (const std::exception & e)
1909  {
1910  std::string msg(e.what());
1911  ASSERT_TRUE(
1912  msg.find("Cannot retrieve molality for equilibrium species 44 since there are only 6 "
1913  "equilibrium species") != std::string::npos)
1914  << "Failed with unexpected error message: " << msg;
1915  }
1916 
1917  // CO2(aq) = -H2O + H+ + HCO3-
1918  EXPECT_NEAR(_egs_calcite.getEquilibriumMolality(0) /
1919  (_egs_calcite.getBasisActivity(1) * _egs_calcite.getBasisActivity(2) /
1920  _egs_calcite.getBasisActivity(0) /
1921  _egs_calcite.getEquilibriumActivityCoefficient(0) /
1922  std::pow(10.0, _egs_calcite.getLog10K(0))),
1923  1.0,
1924  1.0E-8);
1925  EXPECT_NEAR(_egs_calcite.getEquilibriumMolality()[0] /
1926  (_egs_calcite.getBasisActivity(1) * _egs_calcite.getBasisActivity(2) /
1927  _egs_calcite.getBasisActivity(0) /
1928  _egs_calcite.getEquilibriumActivityCoefficient()[0] /
1929  std::pow(10.0, _egs_calcite.getLog10K(0))),
1930  1.0,
1931  1.0E-8);
1932  // CO3-- = HCO3- - H+
1933  EXPECT_NEAR(_egs_calcite.getEquilibriumMolality(1) /
1934  (_egs_calcite.getBasisActivity(2) / _egs_calcite.getBasisActivity(1) /
1935  _egs_calcite.getEquilibriumActivityCoefficient(1) /
1936  std::pow(10.0, _egs_calcite.getLog10K(1))),
1937  1.0,
1938  1.0E-8);
1939  EXPECT_NEAR(_egs_calcite.getEquilibriumMolality()[1] /
1940  (_egs_calcite.getBasisActivity(2) / _egs_calcite.getBasisActivity(1) /
1941  _egs_calcite.getEquilibriumActivityCoefficient()[1] /
1942  std::pow(10.0, _egs_calcite.getLog10K(1))),
1943  1.0,
1944  1.0E-8);
1945  // CaCO3 = Calcite
1946  EXPECT_NEAR(_egs_calcite.getEquilibriumMolality(2) /
1947  (1.0 / _egs_calcite.getEquilibriumActivityCoefficient(2) /
1948  std::pow(10.0, _egs_calcite.getLog10K(2))),
1949  1.0,
1950  1.0E-8);
1951  EXPECT_NEAR(_egs_calcite.getEquilibriumMolality()[2] /
1952  (1.0 / _egs_calcite.getEquilibriumActivityCoefficient()[2] /
1953  std::pow(10.0, _egs_calcite.getLog10K(2))),
1954  1.0,
1955  1.0E-8);
1956  // CaOH+ = Calcite - HCO3- + H20
1957  EXPECT_NEAR(_egs_calcite.getEquilibriumMolality(3) /
1958  (_egs_calcite.getBasisActivity(0) / _egs_calcite.getBasisActivity(2) /
1959  _egs_calcite.getEquilibriumActivityCoefficient()[3] /
1960  std::pow(10.0, _egs_calcite.getLog10K(3))),
1961  1.0,
1962  1.0E-8);
1963  EXPECT_NEAR(_egs_calcite.getEquilibriumMolality()[3] /
1964  (_egs_calcite.getBasisActivity(0) / _egs_calcite.getBasisActivity(2) /
1965  _egs_calcite.getEquilibriumActivityCoefficient(3) /
1966  std::pow(10.0, _egs_calcite.getLog10K(3))),
1967  1.0,
1968  1.0E-8);
1969  // OH- = H20 - H+
1970  EXPECT_NEAR(_egs_calcite.getEquilibriumMolality(4) /
1971  (_egs_calcite.getBasisActivity(0) / _egs_calcite.getBasisActivity(1) /
1972  _egs_calcite.getEquilibriumActivityCoefficient(4) /
1973  std::pow(10.0, _egs_calcite.getLog10K(4))),
1974  1.0,
1975  1.0E-8);
1976  EXPECT_NEAR(_egs_calcite.getEquilibriumMolality()[4] /
1977  (_egs_calcite.getBasisActivity(0) / _egs_calcite.getBasisActivity(1) /
1978  _egs_calcite.getEquilibriumActivityCoefficient()[4] /
1979  std::pow(10.0, _egs_calcite.getLog10K(4))),
1980  1.0,
1981  1.0E-8);
1982  // Ca++ = Calcite - HCO3- + H+
1983  EXPECT_NEAR(_egs_calcite.getEquilibriumMolality(5) /
1984  (_egs_calcite.getBasisActivity(1) / _egs_calcite.getBasisActivity(2) /
1985  _egs_calcite.getEquilibriumActivityCoefficient(5) /
1986  std::pow(10.0, _egs_calcite.getLog10K(5))),
1987  1.0,
1988  1.0E-8);
1989  EXPECT_NEAR(_egs_calcite.getEquilibriumMolality()[5] /
1990  (_egs_calcite.getBasisActivity(1) / _egs_calcite.getBasisActivity(2) /
1991  _egs_calcite.getEquilibriumActivityCoefficient()[5] /
1992  std::pow(10.0, _egs_calcite.getLog10K(5))),
1993  1.0,
1994  1.0E-8);
1995 }
1996 
1998 TEST_F(GeochemicalSystemTest, computeAndGetEquilibriumActivity)
1999 {
2000  GeochemicalSystem nonconst = _egs_redox;
2002  const std::vector<Real> & eqm_m = nonconst.getEquilibriumMolality();
2003  const std::vector<Real> & eqm_ga = nonconst.getEquilibriumActivityCoefficient();
2004  const std::vector<Real> & eqm_act = nonconst.computeAndGetEquilibriumActivity();
2005  for (unsigned j = 0; j < nonconst.getNumInEquilibrium(); ++j)
2006  {
2007  EXPECT_EQ(eqm_act[j], nonconst.getEquilibriumActivity(j));
2008  if (mgd.eqm_species_mineral[j])
2009  EXPECT_EQ(eqm_act[j], 1.0);
2010  else if (mgd.eqm_species_name[j] == "O2(g)")
2011  EXPECT_NEAR(eqm_act[j],
2012  nonconst.getBasisActivity(mgd.basis_species_index.at("O2(aq)")) /
2013  std::pow(10.0, nonconst.getLog10K(j)),
2014  1E-3);
2015  else if (mgd.eqm_species_name[j] == "CH4(g)fake")
2016  EXPECT_NEAR(
2017  eqm_act[j],
2018  std::pow(nonconst.getBasisActivity(mgd.basis_species_index.at("CH4(aq)")), 2.0) *
2019  std::pow(nonconst.getBasisActivity(mgd.basis_species_index.at("Fe+++")), -2.0) *
2020  std::pow(nonconst.getBasisActivity(mgd.basis_species_index.at("HCO3-")), 1.5) /
2021  std::pow(10.0, nonconst.getLog10K(j)),
2022  1E-3);
2023  else
2024  EXPECT_EQ(eqm_act[j], eqm_m[j] * eqm_ga[j]);
2025  }
2026 }
2027 
2029 TEST_F(GeochemicalSystemTest, getEquilibriumActivityExeception)
2030 {
2031  try
2032  {
2033  EXPECT_EQ(_egs_calcite.getEquilibriumActivity(6), 1);
2034  FAIL() << "Missing expected exception.";
2035  }
2036  catch (const std::exception & e)
2037  {
2038  std::string msg(e.what());
2039  ASSERT_TRUE(
2040  msg.find("Cannot retrieve activity for equilibrium species 6 since there are only 6 "
2041  "equilibrium species") != std::string::npos)
2042  << "Failed with unexpected error message: " << msg;
2043  }
2044 }
2045 
2047 TEST_F(GeochemicalSystemTest, getTotalChargeOld)
2048 {
2049  const GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
2051  database, {"H2O", "StoiCheckBasis", "Fe+++"}, {}, {}, {}, {}, {}, "O2(aq)", "e-");
2053  const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
2057  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
2061  const GeochemicalSystem egs(mgd,
2062  _ac3,
2063  _is3,
2064  _swapper3,
2065  {},
2066  {},
2067  "StoiCheckBasis",
2068  {"H2O", "Fe+++", "StoiCheckBasis"},
2069  {1.75, 5.0, 1.0},
2070  cu,
2071  cm,
2072  25,
2073  0,
2074  1E-20,
2075  {},
2076  {},
2077  {});
2078  EXPECT_EQ(egs.getBulkMolesOld()[2], 5.0); // Fe+++
2079  EXPECT_NEAR(egs.getBulkMolesOld()[1],
2080  -6.0,
2081  1.0E-12); // StoiCheckBasis (charge=2.5) computed to ensure charge neutrality
2082  EXPECT_NEAR(egs.getTotalChargeOld(), 0.0, 1.0E-12);
2083 }
2084 
2087 {
2088  const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
2093  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
2098  const GeochemicalSystem egs(_mgd_calcite,
2099  _ac3,
2100  _is3,
2101  _swapper4,
2102  {},
2103  {},
2104  "HCO3-",
2105  {"H2O", "Calcite", "HCO3-", "H+"},
2106  {175.0, 3.0, 2.0, 1.0},
2107  cu,
2108  cm,
2109  25,
2110  0,
2111  1E-20,
2112  {},
2113  {},
2114  {});
2115  // Here, water and HCO3- are the components in the algebraic system
2116  EXPECT_EQ(egs.getNumInAlgebraicSystem(), (unsigned)2);
2117  EXPECT_EQ(egs.getNumBasisInAlgebraicSystem(), (unsigned)2);
2118  EXPECT_EQ(egs.getNumSurfacePotentials(), (unsigned)0);
2119 
2120  const DenseVector<Real> mole_additions(4);
2121 
2122  try
2123  {
2124  EXPECT_EQ(egs.getResidualComponent(3, mole_additions), 2);
2125  FAIL() << "Missing expected exception.";
2126  }
2127  catch (const std::exception & e)
2128  {
2129  std::string msg(e.what());
2130  ASSERT_TRUE(
2131  msg.find(
2132  "Cannot retrieve residual for algebraic index 3 because there are only 2 "
2133  "molalities in the algebraic system and 0 surface potentials and 0 kinetic species") !=
2134  std::string::npos)
2135  << "Failed with unexpected error message: " << msg;
2136  }
2137 
2138  try
2139  {
2140  const DenseVector<Real> bad(5);
2141  EXPECT_EQ(egs.getResidualComponent(0, bad), 2);
2142  FAIL() << "Missing expected exception.";
2143  }
2144  catch (const std::exception & e)
2145  {
2146  std::string msg(e.what());
2147  ASSERT_TRUE(msg.find("The increment in mole numbers (mole_additions) needs to be of size 4 + 0 "
2148  "but it is of size 5") != std::string::npos)
2149  << "Failed with unexpected error message: " << msg;
2150  }
2151 
2152  const Real nw = egs.getSolventMassAndFreeMolalityAndMineralMoles()[0];
2153  const Real mHCO3 = egs.getSolventMassAndFreeMolalityAndMineralMoles()[2];
2154  // H2O has algebraic index = 0 and basis index = 0
2155  Real res = -175.0 + nw * GeochemistryConstants::MOLES_PER_KG_WATER;
2156  for (unsigned j = 0; j < _mgd_calcite.eqm_species_name.size(); ++j)
2157  res += nw * _mgd_calcite.eqm_stoichiometry(j, 0) * egs.getEquilibriumMolality(j);
2158  EXPECT_NEAR(egs.getResidualComponent(0, mole_additions) / res, 1.0, 1.0E-12);
2159  // HCO3- has algebraic index = 1 and basis index = 2
2160  res = -egs.getBulkMolesOld()[1] +
2161  nw * mHCO3; // first term is because of charge balance (H+ has basis index 1)
2162  for (unsigned j = 0; j < _mgd_calcite.eqm_species_name.size(); ++j)
2163  res += nw * _mgd_calcite.eqm_stoichiometry(j, 2) * egs.getEquilibriumMolality(j);
2164  EXPECT_NEAR(egs.getResidualComponent(1, mole_additions) / res, 1.0, 1.0E-12);
2165 
2166  // The above is a good test of the indexing, but due to crazy log10K it is difficult to
2167  // demonstrate correct residual, so instead:
2169  _db_calcite, {"H2O", "H+", "StoiCheckBasis"}, {}, {}, {}, {}, {}, "O2(aq)", "e-");
2170  ModelGeochemicalDatabase mgd2 = model2.modelGeochemicalDatabase();
2171  const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm2 = {
2175  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu2;
2179  const GeochemicalSystem egs2(mgd2,
2180  _ac3,
2181  _is3,
2182  _swapper3,
2183  {},
2184  {},
2185  "StoiCheckBasis",
2186  {"H2O", "H+", "StoiCheckBasis"},
2187  {175.0, 8.0, 2.0},
2188  cu2,
2189  cm2,
2190  25,
2191  0,
2192  1E-20,
2193  {},
2194  {},
2195  {});
2196  DenseVector<Real> mole_add3(3);
2197  // Here, water and StoiCheckBasis are the components in the algebraic system
2198  EXPECT_EQ(egs2.getNumInAlgebraicSystem(), (unsigned)2);
2199  EXPECT_EQ(egs2.getNumBasisInAlgebraicSystem(), (unsigned)2);
2200  EXPECT_EQ(egs2.getNumSurfacePotentials(), (unsigned)0);
2201  // The equilibrium species are StoiCheckRedox (m = 9.82934e+06) and OH- (m = 4.78251e-15)
2202  const Real nw2 = egs2.getSolventMassAndFreeMolalityAndMineralMoles()[0];
2203  const Real mscb = egs2.getSolventMassAndFreeMolalityAndMineralMoles()[2];
2204  const Real mscr = egs2.getEquilibriumMolality(0);
2205  // H2O
2206  const Real resh2o = -175.0 + nw2 * (GeochemistryConstants::MOLES_PER_KG_WATER + (-0.5) * mscr);
2207  EXPECT_NEAR(egs2.getResidualComponent(0, mole_add3) / resh2o, 1.0, 1.0E-12);
2208  // StoiCheckBasis (algebraic index = 1, basis index = 2)
2209  const Real resscb = (1.0 / 2.5) * egs2.getBulkMolesOld()[1] +
2210  nw2 * (mscb + 1.5 * mscr); // first term is from charge balance
2211  EXPECT_NEAR(egs2.getResidualComponent(1, mole_add3) / resscb, 1.0, 1.0E-12);
2212 
2213  for (unsigned i = 0; i < 3; ++i)
2214  mole_add3(i) =
2215  (i + 1) *
2216  1E6; // although moles of StoiCheckBasis are added, they are removed by charge balance
2217  EXPECT_NEAR((egs2.getResidualComponent(0, mole_add3) + 1E6) / resh2o, 1.0, 1.0E-12);
2218  EXPECT_NEAR(egs2.getResidualComponent(1, mole_add3) / resscb, 1.0, 1.0E-12);
2219 }
2220 
2222 TEST_F(GeochemicalSystemTest, enforceChargeBalance)
2223 {
2224  GeochemicalSystem nonconst = _egs_calcite;
2225  ASSERT_TRUE(std::abs(nonconst.getTotalChargeOld()) > 1.0); // initially there is nonzero charge
2226  nonconst.enforceChargeBalance();
2227  EXPECT_NEAR(nonconst.getTotalChargeOld(), 0.0, 1E-12);
2228 }
2229 
2232 {
2233  GeochemicalSystem egs(_mgd_calcite,
2234  _ac3,
2235  _is3,
2236  _swapper4,
2237  {},
2238  {},
2239  "H+",
2240  {"H2O", "Calcite", "H+", "HCO3-"},
2241  {1.75, 3.0, 2.0, 1.0},
2242  _cu_calcite,
2243  _cm_calcite,
2244  25,
2245  0,
2246  1E-20,
2247  {},
2248  {},
2249  {});
2250 
2251  DenseVector<Real> var_incorrect_size(123);
2252  try
2253  {
2254  egs.setAlgebraicVariables(var_incorrect_size);
2255  FAIL() << "Missing expected exception.";
2256  }
2257  catch (const std::exception & e)
2258  {
2259  std::string msg(e.what());
2260  ASSERT_TRUE(msg.find("Incorrect size in setAlgebraicVariables") != std::string::npos)
2261  << "Failed with unexpected error message: " << msg;
2262  }
2263 
2264  DenseVector<Real> var_neg(egs.getNumInAlgebraicSystem());
2265  var_neg(0) = -1.25;
2266  try
2267  {
2268  egs.setAlgebraicVariables(var_neg);
2269  FAIL() << "Missing expected exception.";
2270  }
2271  catch (const std::exception & e)
2272  {
2273  std::string msg(e.what());
2274  ASSERT_TRUE(msg.find("Cannot set algebraic variables to non-positive values such as -1.25") !=
2275  std::string::npos)
2276  << "Failed with unexpected error message: " << msg;
2277  }
2278 }
2279 
2286 {
2287  const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
2292  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
2297  GeochemicalSystem egs(
2298  _mgd_calcite,
2299  _ac8,
2300  _is8,
2301  _swapper4,
2302  {},
2303  {},
2304  "HCO3-",
2305  {"H2O", "Calcite", "HCO3-", "H+"},
2306  {175.0, 3.0, 3.2E-2, 1E-8}, // the molality of H+ is carefully chosen so the secondary species
2307  // do not contribute much to the ionic strengths
2308  cu,
2309  cm,
2310  25,
2311  4,
2312  1E-20,
2313  {},
2314  {},
2315  {});
2316  DenseVector<Real> mole_add(4);
2317  DenseMatrix<Real> dmole_add(4, 4);
2318 
2319  // Here, water and HCO3- are the components in the algebraic system
2320  EXPECT_EQ(egs.getNumInAlgebraicSystem(), (unsigned)2);
2321  EXPECT_EQ(egs.getNumBasisInAlgebraicSystem(), (unsigned)2);
2322  EXPECT_EQ(egs.getNumSurfacePotentials(), (unsigned)0);
2323  DenseVector<Real> res_orig(2);
2324  res_orig(0) = egs.getResidualComponent(0, mole_add);
2325  res_orig(1) = egs.getResidualComponent(1, mole_add);
2326  const std::vector<Real> var_orig =
2327  egs.getAlgebraicVariableValues(); // var[0] ~ 175/55 = 3.2 (kg of solvent water). var[1]
2328  // ~ 3.2E-2/3.2 = 0.01 (molality of HCO3-)
2329  DenseMatrix<Real> jac(1, 1); // deliberately size incorrectly to check resizing done correctly
2330  egs.computeJacobian(res_orig, jac, mole_add, dmole_add);
2331 
2332  const std::vector<Real> mol_orig = egs.getAlgebraicBasisValues();
2333  EXPECT_EQ(mol_orig.size(), var_orig.size());
2334  for (unsigned a = 0; a < mol_orig.size(); ++a)
2335  EXPECT_EQ(mol_orig[a], var_orig[a]);
2336 
2337  const Real eps = 1E-9;
2338 
2339  for (unsigned var_num = 0; var_num < 2; ++var_num)
2340  {
2341  std::vector<Real> var_new = var_orig;
2342  var_new[var_num] += eps;
2343  egs.setAlgebraicVariables(var_new);
2344  for (unsigned res_comp = 0; res_comp < 2; ++res_comp)
2345  {
2346  const Real expected =
2347  (egs.getResidualComponent(res_comp, mole_add) - res_orig(res_comp)) / eps;
2348  EXPECT_NEAR(jac(res_comp, var_num) / expected, 1.0, 1E-5);
2349  }
2350  }
2351 
2352  // now check the derivatives including mole_additions
2353  for (unsigned i = 0; i < 4; ++i)
2354  mole_add(i) = i + 1; // add 1 mole of solvent water and of 3 of HCO3-. The HCO3- addition gets
2355  // removed by charge-balance
2356  egs.setAlgebraicVariables(var_orig);
2357  res_orig(0) = egs.getResidualComponent(0, mole_add);
2358  res_orig(1) = egs.getResidualComponent(1, mole_add);
2359  for (unsigned i = 0; i < 4; ++i)
2360  for (unsigned j = 0; j < 4; ++j)
2361  dmole_add(i, j) = std::sin(i + 1.0) * (j + 1); // just some randomish derivatives
2362  egs.computeJacobian(res_orig, jac, mole_add, dmole_add);
2363  for (unsigned var_num = 0; var_num < 2; ++var_num)
2364  {
2365  std::vector<Real> var_new = var_orig;
2366  var_new[var_num] += eps;
2367  egs.setAlgebraicVariables(var_new);
2368  for (unsigned i = 0; i < 4; ++i)
2369  {
2370  if (var_num == 0)
2371  mole_add(i) = i + 1 + dmole_add(i, 0) * eps;
2372  else
2373  mole_add(i) = i + 1 + dmole_add(i, 2) * eps;
2374  }
2375  for (unsigned res_comp = 0; res_comp < 2; ++res_comp)
2376  {
2377  const Real expected =
2378  (egs.getResidualComponent(res_comp, mole_add) - res_orig(res_comp)) / eps;
2379  EXPECT_NEAR(jac(res_comp, var_num) / expected, 1.0, 1E-5);
2380  }
2381  }
2382 }
2383 
2390 {
2391  const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
2396  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
2401  GeochemicalSystem egs(_mgd_calcite,
2402  _ac8,
2403  _is8,
2404  _swapper4,
2405  {},
2406  {},
2407  "HCO3-",
2408  {"H2O", "Calcite", "HCO3-", "H+"},
2409  {1.0, 3.0, 3.2E-4, 1E-4},
2410  cu,
2411  cm,
2412  25,
2413  2,
2414  1E-20,
2415  {},
2416  {},
2417  {});
2418  DenseVector<Real> mole_add(4);
2419  DenseMatrix<Real> dmole_add(4, 4);
2420 
2421  // Here, H+ and HCO3- are the components in the algebraic system
2422  EXPECT_EQ(egs.getNumInAlgebraicSystem(), (unsigned)2);
2423  EXPECT_EQ(egs.getNumBasisInAlgebraicSystem(), (unsigned)2);
2424  EXPECT_EQ(egs.getNumSurfacePotentials(), (unsigned)0);
2425  DenseVector<Real> res_orig(2);
2426  res_orig(0) = egs.getResidualComponent(0, mole_add);
2427  res_orig(1) = egs.getResidualComponent(1, mole_add);
2428  const std::vector<Real> var_orig = egs.getAlgebraicVariableValues(); // both are 9E-5
2429  DenseMatrix<Real> jac(1, 1); // deliberately size incorrectly to check resizing done correctly
2430  egs.computeJacobian(res_orig, jac, mole_add, dmole_add);
2431 
2432  const Real eps = 1E-9;
2433  for (unsigned var_num = 0; var_num < 2; ++var_num)
2434  {
2435  std::vector<Real> var_new = var_orig;
2436  var_new[var_num] += eps;
2437  egs.setAlgebraicVariables(var_new);
2438  for (unsigned res_comp = 0; res_comp < 2; ++res_comp)
2439  {
2440  const Real expected =
2441  (egs.getResidualComponent(res_comp, mole_add) - res_orig(res_comp)) / eps;
2442  if (var_num == 0)
2443  EXPECT_NEAR(jac(res_comp, var_num) / expected,
2444  1.0,
2445  1.0E-10); // everything happens to be linear in H+
2446  else
2447  EXPECT_NEAR(jac(res_comp, var_num) / expected, 1.0, 1.0E-4); // nonlinear in HCO3-
2448  }
2449  }
2450 
2451  // now check the derivatives including mole_additions
2452  for (unsigned i = 0; i < 4; ++i)
2453  mole_add(i) = (i + 1) * 1E4;
2454  egs.setAlgebraicVariables(var_orig);
2455  res_orig(0) = egs.getResidualComponent(0, mole_add);
2456  res_orig(1) = egs.getResidualComponent(1, mole_add);
2457  for (unsigned i = 0; i < 4; ++i)
2458  for (unsigned j = 0; j < 4; ++j)
2459  dmole_add(i, j) = std::sin(i + 1.0) * (j + 1) * 1E5; // just some randomish derivatives
2460  egs.computeJacobian(res_orig, jac, mole_add, dmole_add);
2461  for (unsigned var_num = 0; var_num < 2; ++var_num)
2462  {
2463  std::vector<Real> var_new = var_orig;
2464  var_new[var_num] += eps;
2465  egs.setAlgebraicVariables(var_new);
2466  for (unsigned i = 0; i < 4; ++i)
2467  {
2468  if (var_num == 0)
2469  mole_add(i) = (i + 1) * 1E4 + dmole_add(i, 1) * eps; // change H+
2470  else
2471  mole_add(i) = (i + 1) * 1E4 + dmole_add(i, 2) * eps; // change HCO3-
2472  }
2473  for (unsigned res_comp = 0; res_comp < 2; ++res_comp)
2474  {
2475  const Real expected =
2476  (egs.getResidualComponent(res_comp, mole_add) - res_orig(res_comp)) / eps;
2477  if (var_num == 0)
2478  EXPECT_NEAR(jac(res_comp, var_num) / expected,
2479  1.0,
2480  1.0E-7); // everything happens to be linear in H+
2481  else
2482  EXPECT_NEAR(jac(res_comp, var_num) / expected, 1.0, 1.0E-4); // nonlinear in HCO3-
2483  }
2484  }
2485 }
2486 
2493 {
2494  const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
2499  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
2504  GeochemicalSystem egs(_mgd_calcite,
2505  _ac8,
2506  _is8,
2507  _swapper4,
2508  {},
2509  {},
2510  "H+",
2511  {"H2O", "Calcite", "H+", "HCO3-"},
2512  {1.0, 3.0, 3.2E-4, 1E-4},
2513  cu,
2514  cm,
2515  25,
2516  0,
2517  1E-20,
2518  {},
2519  {},
2520  {});
2521  DenseVector<Real> mole_add(4);
2522  DenseMatrix<Real> dmole_add(4, 4);
2523 
2524  // Here H+ is the only component in the algebraic system
2525  EXPECT_EQ(egs.getNumInAlgebraicSystem(), (unsigned)1);
2526  EXPECT_EQ(egs.getNumBasisInAlgebraicSystem(), (unsigned)1);
2527  EXPECT_EQ(egs.getNumSurfacePotentials(), (unsigned)0);
2528  DenseVector<Real> res_orig(1);
2529  res_orig(0) = egs.getResidualComponent(0, mole_add);
2530  const std::vector<Real> var_orig = egs.getAlgebraicVariableValues(); // around 3E-4
2531  DenseMatrix<Real> jac(1, 1);
2532  egs.computeJacobian(res_orig, jac, mole_add, dmole_add);
2533 
2534  const Real eps = 1E-11;
2535  std::vector<Real> var_new = var_orig;
2536  var_new[0] += eps;
2537  egs.setAlgebraicVariables(var_new);
2538  Real expected = (egs.getResidualComponent(0, mole_add) - res_orig(0)) / eps;
2539  EXPECT_NEAR(jac(0, 0) / expected, 1.0, 1.0E-7);
2540 
2541  // now check the derivatives including mole_additions
2542  // because of charge balance there is actually no dependence on mole_add
2543  for (unsigned i = 0; i < 4; ++i)
2544  mole_add(i) = i + 1;
2545  egs.setAlgebraicVariables(var_orig);
2546  res_orig(0) = egs.getResidualComponent(0, mole_add);
2547  for (unsigned i = 0; i < 4; ++i)
2548  for (unsigned j = 0; j < 4; ++j)
2549  dmole_add(i, j) = std::sin(i + 1.0) * (j + 1) * 1E10; // huge randomish derivatives
2550  egs.computeJacobian(res_orig, jac, mole_add, dmole_add);
2551  egs.setAlgebraicVariables(var_new);
2552  for (unsigned i = 0; i < 4; ++i)
2553  mole_add(i) = (i + 1) + dmole_add(i, 1) * eps; // change H+
2554  expected = (egs.getResidualComponent(0, mole_add) - res_orig(0)) / eps;
2555  EXPECT_NEAR(jac(0, 0) / expected, 1.0, 1.0E-7);
2556 }
2557 
2564 {
2566  _db_calcite, {"H2O", "H+", "Ca++", "HCO3-"}, {"Calcite"}, {}, {}, {}, {}, "O2(aq)", "e-");
2568  const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
2573  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
2578  GeochemicalSystem egs(mgd,
2579  _ac8,
2580  _is8,
2581  _swapper4,
2582  {},
2583  {},
2584  "H+",
2585  {"H2O", "H+", "Ca++", "HCO3-"},
2586  {175.0, 1E1, 4E1, 1E-4},
2587  cu,
2588  cm,
2589  25,
2590  0,
2591  1E-20,
2592  {},
2593  {},
2594  {});
2595  DenseVector<Real> mole_add(4);
2596  DenseMatrix<Real> dmole_add(4, 4);
2597 
2598  EXPECT_EQ(egs.getNumInAlgebraicSystem(), (unsigned)3);
2599  EXPECT_EQ(egs.getNumBasisInAlgebraicSystem(), (unsigned)3);
2600  EXPECT_EQ(egs.getNumSurfacePotentials(), (unsigned)0);
2601  DenseVector<Real> res_orig(3);
2602  res_orig(0) = egs.getResidualComponent(0, mole_add);
2603  res_orig(1) = egs.getResidualComponent(1, mole_add);
2604  res_orig(2) = egs.getResidualComponent(2, mole_add);
2605  const std::vector<Real> var_orig = egs.getAlgebraicVariableValues(); // around {3.1, 2.9, 11.4}
2606  DenseMatrix<Real> jac(1, 1);
2607  egs.computeJacobian(res_orig, jac, mole_add, dmole_add);
2608 
2609  const std::vector<Real> mol_orig = egs.getAlgebraicBasisValues();
2610  EXPECT_EQ(mol_orig.size(), var_orig.size());
2611  for (unsigned a = 0; a < mol_orig.size(); ++a)
2612  EXPECT_EQ(mol_orig[a], var_orig[a]);
2613 
2614  const Real eps = 1E-6;
2615  for (unsigned var_num = 0; var_num < 3; ++var_num)
2616  {
2617  std::vector<Real> var_new = var_orig;
2618  var_new[var_num] += eps;
2619  egs.setAlgebraicVariables(var_new);
2620  for (unsigned res_comp = 0; res_comp < 3; ++res_comp)
2621  {
2622  const Real expected =
2623  (egs.getResidualComponent(res_comp, mole_add) - res_orig(res_comp)) / eps;
2624  if (std::abs(expected) < 1.0E-6)
2625  EXPECT_NEAR(jac(res_comp, var_num), 0.0, 1.0E-6);
2626  else
2627  EXPECT_NEAR(expected / jac(res_comp, var_num), 1.0, 1.0E-6);
2628  }
2629  }
2630 
2631  // now check the derivatives including mole_additions
2632  for (unsigned i = 0; i < 4; ++i)
2633  mole_add(i) = i + 1;
2634  egs.setAlgebraicVariables(var_orig);
2635  res_orig(0) = egs.getResidualComponent(0, mole_add);
2636  res_orig(1) = egs.getResidualComponent(1, mole_add);
2637  res_orig(2) = egs.getResidualComponent(2, mole_add);
2638  for (unsigned i = 0; i < 4; ++i)
2639  for (unsigned j = 0; j < 4; ++j)
2640  dmole_add(i, j) = 10 * std::sin(i + 1.0) * (j + 1); // just some randomish derivatives
2641  egs.computeJacobian(res_orig, jac, mole_add, dmole_add);
2642  for (unsigned var_num = 0; var_num < 3; ++var_num)
2643  {
2644  std::vector<Real> var_new = var_orig;
2645  var_new[var_num] += eps;
2646  egs.setAlgebraicVariables(var_new);
2647  for (unsigned i = 0; i < 4; ++i)
2648  mole_add(i) = i + 1 + dmole_add(i, var_num) * eps;
2649  for (unsigned res_comp = 0; res_comp < 3; ++res_comp)
2650  {
2651  const Real expected =
2652  (egs.getResidualComponent(res_comp, mole_add) - res_orig(res_comp)) / eps;
2653  EXPECT_NEAR(jac(res_comp, var_num) / expected, 1.0, 1E-7);
2654  }
2655  }
2656 }
2657 
2659 TEST_F(GeochemicalSystemTest, saturationIndices)
2660 {
2662  _db_calcite, {"H2O", "H+", "HCO3-", "Ca++"}, {"Calcite"}, {}, {}, {}, {}, "O2(aq)", "e-");
2664  const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm_fixed_activity = {
2669  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu_fixed_activity;
2670  cu_fixed_activity.push_back(GeochemistryUnitConverter::GeochemistryUnit::KG);
2672  cu_fixed_activity.push_back(GeochemistryUnitConverter::GeochemistryUnit::MOLES);
2674  const GeochemicalSystem egs(mgd,
2675  _ac3,
2676  _is3,
2677  _swapper4,
2678  {},
2679  {},
2680  "H+",
2681  {"H2O", "Ca++", "H+", "HCO3-"},
2682  {1.11, 3.0, 2.0, 1.5},
2683  cu_fixed_activity,
2684  cm_fixed_activity,
2685  25,
2686  0,
2687  1E-20,
2688  {},
2689  {},
2690  {});
2691  std::vector<Real> si = egs.getSaturationIndices();
2692  for (const auto & name_index : mgd.eqm_species_index)
2693  if (name_index.first == "Calcite")
2694  EXPECT_NEAR(
2695  si[name_index.second],
2696  std::log10(egs.getBasisActivity(3) * egs.getBasisActivity(2) / egs.getBasisActivity(1)) -
2697  1.7130,
2698  1E-9);
2699  else
2700  EXPECT_EQ(si[name_index.second], 0.0);
2701 }
2702 
2705 {
2706  GeochemicalSystem nonconst = _egs_calcite;
2707  try
2708  {
2709  nonconst.performSwap(0, 0);
2710  FAIL() << "Missing expected exception.";
2711  }
2712  catch (const std::exception & e)
2713  {
2714  std::string msg(e.what());
2715  ASSERT_TRUE(
2716  msg.find("GeochemicalSystem: attempting to swap out water and replace it by CO2(aq)."
2717  " This could be because the algorithm would like to "
2718  "swap out the charge-balance species, in which case you should choose a "
2719  "different charge-balance species") != std::string::npos)
2720  << "Failed with unexpected error message: " << msg;
2721  }
2722 
2723  try
2724  {
2725  nonconst.performSwap(1, 0);
2726  FAIL() << "Missing expected exception.";
2727  }
2728  catch (const std::exception & e)
2729  {
2730  std::string msg(e.what());
2731  ASSERT_TRUE(msg.find("GeochemicalSystem: attempting to swap the charge-balance "
2732  "species out of the basis") != std::string::npos)
2733  << "Failed with unexpected error message: " << msg;
2734  }
2735 
2736  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
2738  database, {"H2O", "H+", "HCO3-", "O2(aq)"}, {}, {"O2(g)"}, {}, {}, {}, "O2(aq)", "e-");
2740 
2741  std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
2746  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
2751 
2752  GeochemicalSystem egs(mgd,
2753  _ac3,
2754  _is3,
2755  _swapper4,
2756  {},
2757  {},
2758  "H+",
2759  {"H2O", "H+", "O2(aq)", "HCO3-"},
2760  {1.0, 1.5, 3.0, 2.5},
2761  cu,
2762  cm,
2763  25,
2764  0,
2765  1E-20,
2766  {},
2767  {},
2768  {});
2769 
2770  try
2771  {
2772  egs.performSwap(3, 5);
2773  FAIL() << "Missing expected exception.";
2774  }
2775  catch (const std::exception & e)
2776  {
2777  std::string msg(e.what());
2778  ASSERT_TRUE(msg.find("GeochemicalSystem: attempting to swap a gas into the basis") !=
2779  std::string::npos)
2780  << "Failed with unexpected error message: " << msg;
2781  }
2782 
2783  std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm2;
2788  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu2;
2793  GeochemicalSystem egs2(mgd,
2794  _ac3,
2795  _is3,
2796  _swapper4,
2797  {"O2(aq)"},
2798  {"O2(g)"},
2799  "H+",
2800  {"H2O", "H+", "O2(g)", "HCO3-"},
2801  {1.0, 1.5, 3.0, 2.5},
2802  cu2,
2803  cm2,
2804  25,
2805  0,
2806  1E-20,
2807  {},
2808  {},
2809  {});
2810 
2811  try
2812  {
2813  egs2.performSwap(3, 0);
2814  FAIL() << "Missing expected exception.";
2815  }
2816  catch (const std::exception & e)
2817  {
2818  std::string msg(e.what());
2819  ASSERT_TRUE(msg.find("GeochemicalSystem: attempting to swap a gas out of the basis") !=
2820  std::string::npos)
2821  << "Failed with unexpected error message: " << msg;
2822  }
2823 }
2824 
2827 {
2828  const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
2833  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
2838  ModelGeochemicalDatabase mgd = _model_calcite.modelGeochemicalDatabase();
2839  GeochemicalSystem egs(mgd,
2840  _ac3,
2841  _is3,
2842  _swapper4,
2843  {"Ca++"},
2844  {"Calcite"},
2845  "H+",
2846  {"H2O", "Calcite", "H+", "HCO3-"},
2847  {1.75, 3.0, 2.0, 1.0},
2848  cu,
2849  cm,
2850  25,
2851  0,
2852  1E-20,
2853  {},
2854  {},
2855  {});
2856 
2857  EXPECT_EQ(egs.getNumInAlgebraicSystem(), (unsigned)1);
2858  EXPECT_EQ(egs.getNumBasisInAlgebraicSystem(), (unsigned)1);
2859  EXPECT_EQ(egs.getNumSurfacePotentials(), (unsigned)0);
2860  EXPECT_EQ(egs.getInAlgebraicSystem()[0], false);
2861  EXPECT_EQ(egs.getInAlgebraicSystem()[1], true);
2862  EXPECT_EQ(egs.getInAlgebraicSystem()[2], false);
2863  EXPECT_EQ(egs.getInAlgebraicSystem()[3], false);
2864  EXPECT_EQ(egs.getChargeBalanceBasisIndex(), (unsigned)1);
2865  const std::vector<Real> bm = egs.getBulkMolesOld();
2866 
2867  egs.performSwap(3, 5); // swap out Calcite and replace by Ca++
2868 
2869  EXPECT_EQ(egs.getNumInAlgebraicSystem(), (unsigned)2); // now Ca++ has a bulk moles attached to it
2870  EXPECT_EQ(egs.getNumBasisInAlgebraicSystem(),
2871  (unsigned)2); // now Ca++ has a bulk moles attached to it
2872  EXPECT_EQ(egs.getNumSurfacePotentials(), (unsigned)0);
2873  EXPECT_EQ(egs.getInAlgebraicSystem()[0], false);
2874  EXPECT_EQ(egs.getInAlgebraicSystem()[1], true);
2875  EXPECT_EQ(egs.getInAlgebraicSystem()[2], false);
2876  EXPECT_EQ(egs.getInAlgebraicSystem()[3], true);
2877  // Ca++ = Calcite - HCO3- + H+. So, same number of bulk moles of Ca++ as Calcite, but number of
2878  // bulk moles of H+ = (original bulk moles of H+) - bulk moles of Calcite
2879  EXPECT_EQ(egs.getChargeBalanceBasisIndex(), (unsigned)1);
2880  EXPECT_NEAR(egs.getBulkMolesOld()[1], bm[1] - bm[3], 1.0E-7);
2881  EXPECT_NEAR(egs.getBulkMolesOld()[3], bm[3], 1.0E-7);
2882 }
2883 
2886 {
2887  const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
2892  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
2897  GeochemistryIonicStrength is(0.0078125, 0.0078125, false, false);
2899  ModelGeochemicalDatabase mgd = _model_calcite.modelGeochemicalDatabase();
2900  const GeochemicalSystem egs_small_IS(mgd,
2901  ac,
2902  is,
2903  _swapper4,
2904  {},
2905  {},
2906  "H+",
2907  {"H2O", "Ca++", "H+", "HCO3-"},
2908  {1.75, 3.0E-1, 2.0E-1, 1.0E-1},
2909  cu,
2910  cm,
2911  25,
2912  0,
2913  1E-20,
2914  {},
2915  {},
2916  {});
2917  EXPECT_EQ(egs_small_IS.getIonicStrength(), 0.0078125);
2918  EXPECT_EQ(egs_small_IS.getStoichiometricIonicStrength(), 0.0078125);
2919 
2920  GeochemistryIonicStrength is_false(3.0, 3.0, true, false);
2921  GeochemistryActivityCoefficientsDebyeHuckel ac_false(is_false, _db_calcite);
2922  const GeochemicalSystem egs(
2923  mgd,
2924  ac_false,
2925  is_false,
2926  _swapper4,
2927  {},
2928  {},
2929  "H+",
2930  {"H2O", "Ca++", "H+", "HCO3-"},
2931  {1.75, 1E-10, 1E-10, 1.0}, // up to 1E-10, the only contributor to IS is HCO3-
2932  cu,
2933  cm,
2934  25,
2935  0,
2936  1E-20,
2937  {},
2938  {},
2939  {});
2940  EXPECT_NEAR(egs.getIonicStrength(), 0.5, 1.0E-8);
2941  EXPECT_NEAR(egs.getStoichiometricIonicStrength(), 0.5, 1.0E-8);
2942 }
2943 
2945 TEST_F(GeochemicalSystemTest, alterAndRevertChargeBalance)
2946 {
2947  const GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
2949  database, {"H2O", "StoiCheckBasis", "HCO3-"}, {}, {}, {}, {}, {}, "O2(aq)", "e-");
2951  const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
2955  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
2959  GeochemicalSystem egs(mgd,
2960  _ac3,
2961  _is3,
2962  _swapper3,
2963  {},
2964  {},
2965  "StoiCheckBasis",
2966  {"H2O", "StoiCheckBasis", "HCO3-"},
2967  {1.75, 5.0, 1.0},
2968  cu,
2969  cm,
2970  25,
2971  0,
2972  1E-20,
2973  {},
2974  {},
2975  {});
2976  // in this case, charge-neutrality can be enforced immediately
2977  EXPECT_NEAR(egs.getTotalChargeOld(), 0.0, 1.0E-12);
2978  EXPECT_NEAR(egs.getBulkMolesOld()[1],
2979  1.0 / 2.5,
2980  1E-9); // this has been changed from 5.0 by the charge-balancing
2981  EXPECT_EQ(egs.getBulkMolesOld()[2], 1.0);
2982  EXPECT_EQ(egs.getChargeBalanceBasisIndex(), (unsigned)1);
2983 
2984  // StoiCheckBasis has molality 0.2, while HCO3- has molality 0.5, so can change the charge-balance
2985  // species via:
2986  egs.alterChargeBalanceSpecies(0.3); // now the charge-balance species should be HCO3-
2987  EXPECT_EQ(egs.getChargeBalanceBasisIndex(), (unsigned)2);
2988  EXPECT_EQ(egs.getBulkMolesOld()[1], 5.0);
2989  EXPECT_EQ(egs.getBulkMolesOld()[2],
2990  5 * 2.5); // this has been changed from 1.0 by the charge-balancing
2991 
2992  // revert to the original
2993  EXPECT_TRUE(egs.revertToOriginalChargeBalanceSpecies());
2994  EXPECT_EQ(egs.getChargeBalanceBasisIndex(), (unsigned)1);
2995  EXPECT_NEAR(egs.getBulkMolesOld()[1], 1.0 / 2.5, 1E-9);
2996  EXPECT_EQ(egs.getBulkMolesOld()[2], 1.0);
2997 }
2998 
3000 TEST_F(GeochemicalSystemTest, getNumRedox) { EXPECT_EQ(_egs_redox.getNumRedox(), (unsigned)3); }
3001 
3003 TEST_F(GeochemicalSystemTest, getOriginalRedoxLHS)
3004 {
3005  EXPECT_EQ(_egs_redox.getOriginalRedoxLHS(), "e-");
3006 }
3007 
3010 {
3011  // not sure which order the redox has been ordered in. The reactions are:
3012  // e- = (1/4/7.5)(O-phth)-- + (1/2 + 5/4/7.5)H2O + (-1 - 6/4/7.5)H+ - 8/4/7.5HCO3-
3013  // e- = (1/8)CH4(aq) + (1/2 - 1/8)H2O - (1+1/8)H+ - (1/8)HCO3-
3014  const bool ophth_is_slot_one = (_mgd_redox.redox_stoichiometry(1, 4) > 1.0E-6);
3015  const unsigned ophth_slot = (ophth_is_slot_one ? 1 : 2);
3016  const unsigned ch4_slot = (ophth_is_slot_one ? 2 : 1);
3017  Real boa = 1.0 / 4.0 / 7.5;
3018  EXPECT_NEAR(
3019  _egs_redox.getRedoxLog10K(ophth_slot), -boa * 542.8292 + 20.7757 - 0.25 * (-2.8990), 1E-8);
3020  boa = 1.0 / 8.0;
3021  EXPECT_NEAR(
3022  _egs_redox.getRedoxLog10K(ch4_slot), -boa * 144.1080 + 20.7757 - 0.25 * (-2.8990), 1E-8);
3023  try
3024  {
3025  _egs_redox.getRedoxLog10K(3);
3026  FAIL() << "Missing expected exception.";
3027  }
3028  catch (const std::exception & e)
3029  {
3030  std::string msg(e.what());
3031  ASSERT_TRUE(
3032  msg.find(
3033  "Cannot retrieve log10K for redox species 3 since there are only 3 redox species") !=
3034  std::string::npos)
3035  << "Failed with unexpected error message: " << msg;
3036  }
3037 }
3038 
3040 TEST_F(GeochemicalSystemTest, log10RedoxActivityProduct)
3041 {
3042  // not sure which order the redox has been ordered in. The reactions are:
3043  // e- = (1/4/7.5)(O-phth)-- + (1/2 + 5/4/7.5)H2O + (-1 - 6/4/7.5)H+ - 8/4/7.5HCO3-
3044  // e- = (1/8)CH4(aq) + (1/2 - 1/8)H2O - (1+1/8)H+ - (1/8)HCO3-
3045  const bool ophth_is_slot_one = (_mgd_redox.redox_stoichiometry(1, 4) > 1.0E-6);
3046  const unsigned ophth_slot = (ophth_is_slot_one ? 1 : 2);
3047  const unsigned ch4_slot = (ophth_is_slot_one ? 2 : 1);
3048 
3049  Real boa = 1.0 / 4.0 / 7.5;
3050  const Real log10ap_o =
3051  boa * std::log10(5.0) + (-1.0 - 6.0 * boa) * std::log10(2.0) - 8 * boa * std::log10(3.0);
3052  EXPECT_NEAR(_egs_redox.log10RedoxActivityProduct(ophth_slot), log10ap_o, 1E-8);
3053 
3054  boa = 1.0 / 8.0;
3055  const Real log10ap_c =
3056  boa * std::log10(6.0) - (1.0 + boa) * std::log10(2.0) - boa * std::log10(3.0);
3057  EXPECT_NEAR(_egs_redox.log10RedoxActivityProduct(ch4_slot), log10ap_c, 1E-8);
3058  try
3059  {
3060  _egs_redox.log10RedoxActivityProduct(3);
3061  FAIL() << "Missing expected exception.";
3062  }
3063  catch (const std::exception & e)
3064  {
3065  std::string msg(e.what());
3066  ASSERT_TRUE(msg.find("Cannot retrieve activity product for redox species 3 since there are "
3067  "only 3 redox species") != std::string::npos)
3068  << "Failed with unexpected error message: " << msg;
3069  }
3070 }
3071 
3073 TEST_F(GeochemicalSystemTest, getSurfacePotentialException)
3074 {
3075  try
3076  {
3077  _egs_calcite.getSurfacePotential(0);
3078  FAIL() << "Missing expected exception.";
3079  }
3080  catch (const std::exception & e)
3081  {
3082  std::string msg(e.what());
3083  ASSERT_TRUE(msg.find("Cannot retrieve the surface potential for surface 0 since there are only "
3084  "0 surfaces involved in surface complexation") != std::string::npos)
3085  << "Failed with unexpected error message: " << msg;
3086  }
3087 }
3088 
3090 TEST_F(GeochemicalSystemTest, getSurfaceChargeException)
3091 {
3092  try
3093  {
3094  _egs_calcite.getSurfaceCharge(0);
3095  FAIL() << "Missing expected exception.";
3096  }
3097  catch (const std::exception & e)
3098  {
3099  std::string msg(e.what());
3100  ASSERT_TRUE(msg.find("Cannot retrieve the surface charge for surface 0 since there are only 0 "
3101  "surfaces involved in surface complexation") != std::string::npos)
3102  << "Failed with unexpected error message: " << msg;
3103  }
3104 }
3105 
3107 TEST_F(GeochemicalSystemTest, sorbingSurfaceArea)
3108 {
3109  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
3111  {"H2O", "H+", ">(s)FeOH", ">(w)FeOH", "Fe+++"},
3112  {"Fe(OH)3(ppd)"},
3113  {},
3114  {},
3115  {},
3116  {},
3117  "O2(aq)",
3118  "e-");
3120  const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
3126  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
3132  GeochemicalSystem egs(mgd,
3133  _ac3,
3134  _is3,
3135  _swapper5,
3136  {"Fe+++"},
3137  {"Fe(OH)3(ppd)"},
3138  "H+",
3139  {"H2O", "H+", ">(s)FeOH", ">(w)FeOH", "Fe(OH)3(ppd)"},
3140  {1.75, 1.0, 2.0, 1.0, 1.5},
3141  cu,
3142  cm,
3143  25.0,
3144  0,
3145  1E-20,
3146  {},
3147  {},
3148  {});
3149  EXPECT_EQ(egs.getSorbingSurfaceArea().size(), (std::size_t)1);
3150  EXPECT_EQ(egs.getSorbingSurfaceArea()[0], 1.5 * 106.8689 * 600.0);
3151 }
3152 
3155 {
3156  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
3158  {"H2O", "H+", ">(s)FeOH", ">(w)FeOH", "Fe+++", "HCO3-"},
3159  {"Fe(OH)3(ppd)"},
3160  {},
3161  {},
3162  {},
3163  {},
3164  "O2(aq)",
3165  "e-");
3167  const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
3174  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
3181  const Real temp = 45.0;
3182  GeochemicalSystem egs(mgd,
3183  _ac3,
3184  _is3,
3185  _swapper6,
3186  {},
3187  {},
3188  "H+",
3189  {"H2O", "H+", ">(s)FeOH", ">(w)FeOH", "Fe+++", "HCO3-"},
3190  {1.75, 1.0, 2.0, 3.0, 1.0, 1.0},
3191  cu,
3192  cm,
3193  temp,
3194  0,
3195  1E-20,
3196  {},
3197  {},
3198  {});
3199  const DenseVector<Real> mole_add(6);
3200  const DenseMatrix<Real> dmole_add(6, 6);
3201 
3202  // test basic things like sizes and that the algebraic variables are correctly set
3203  EXPECT_EQ(egs.getNumInAlgebraicSystem(), (unsigned)3);
3204  EXPECT_EQ(egs.getNumBasisInAlgebraicSystem(), (unsigned)2);
3205  EXPECT_EQ(egs.getNumSurfacePotentials(), (unsigned)1);
3206  EXPECT_EQ(egs.getBasisIndexOfAlgebraicSystem()[0],
3207  (unsigned)1); // H+ is slot=0 in algebraic system and slot=1 in basis
3208  EXPECT_EQ(egs.getBasisIndexOfAlgebraicSystem()[1],
3209  (unsigned)5); // HCO3-
3210  EXPECT_EQ(egs.getAlgebraicIndexOfBasisSystem()[1], (unsigned)0);
3211  EXPECT_EQ(egs.getAlgebraicIndexOfBasisSystem()[5], (unsigned)1);
3212  DenseVector<Real> alg_vars = egs.getAlgebraicVariableDenseValues();
3213  EXPECT_EQ(alg_vars.size(), (std::size_t)3);
3214  std::vector<Real> mols = egs.getAlgebraicBasisValues();
3215  EXPECT_EQ(mols.size(), (std::size_t)2);
3216  EXPECT_EQ(mols[0], alg_vars(0));
3217  EXPECT_EQ(mols[1], alg_vars(1));
3218  EXPECT_EQ(egs.getAlgebraicVariableValues().size(), (std::size_t)3);
3219 
3220  // try setting a _surface_pot_expr to a non-positive quantity
3221  try
3222  {
3223  alg_vars(2) = -1.0;
3224  egs.setAlgebraicVariables(alg_vars);
3225  FAIL() << "Missing expected exception.";
3226  }
3227  catch (const std::exception & e)
3228  {
3229  std::string msg(e.what());
3230  ASSERT_TRUE(msg.find("Cannot set algebraic variables to non-positive values such as -1") !=
3231  std::string::npos)
3232  << "Failed with unexpected error message: " << msg;
3233  }
3234 
3235  // set the _surface_pot_expr to 2.0
3236  alg_vars(2) = 2.0;
3237  egs.setAlgebraicVariables(alg_vars);
3238 
3239  // check that the setting worked
3240  const std::vector<Real> av = egs.getAlgebraicVariableValues();
3241  EXPECT_EQ(av.size(), (std::size_t)3);
3242  for (unsigned i = 0; i < 3; ++i)
3243  EXPECT_EQ(av[i], alg_vars(i));
3244 
3245  // check the surface potential is correctly computed
3246  const Real surf_pot_gold = -std::log(alg_vars(2)) * 2.0 * 8.314472 * 318.15 / 96485.3415;
3247  EXPECT_NEAR(egs.getSurfacePotential(0), surf_pot_gold, 1.0E-9);
3248  try
3249  {
3250  egs.getSurfacePotential(1);
3251  FAIL() << "Missing expected exception.";
3252  }
3253  catch (const std::exception & e)
3254  {
3255  std::string msg(e.what());
3256  ASSERT_TRUE(msg.find("Cannot retrieve the surface potential for surface 1 since there are only "
3257  "1 surfaces involved in surface complexation") != std::string::npos)
3258  << "Failed with unexpected error message: " << msg;
3259  }
3260 
3261  // check the specific surface charge is correctly computed
3262  const Real pref = 0.5 / 96485.3415 *
3263  std::sqrt(8.0 * 8.314472 * 318.15 * 8.8541878128E-12 * 78.5 * 1000.0 *
3264  egs.getIonicStrength()) *
3265  (-alg_vars(2) + 1.0 / alg_vars(2));
3266  const Real surf_charge_gold = pref * 96485.3415;
3267  EXPECT_NEAR(egs.getSurfaceCharge(0), surf_charge_gold, 1.0E-9);
3268  try
3269  {
3270  egs.getSurfaceCharge(1);
3271  FAIL() << "Missing expected exception.";
3272  }
3273  catch (const std::exception & e)
3274  {
3275  std::string msg(e.what());
3276  ASSERT_TRUE(msg.find("Cannot retrieve the surface charge for surface 1 since there are only 1 "
3277  "surfaces involved in surface complexation") != std::string::npos)
3278  << "Failed with unexpected error message: " << msg;
3279  }
3280 
3281  // check that equilibrium molality for the surface species is correctly computed
3282  const Real act_h = egs.getBasisActivity(1);
3283  const Real act_bicarb = egs.getBasisActivity(2);
3284  const unsigned ind_surf = model.modelGeochemicalDatabase().eqm_species_index.at(">(s)FeO-");
3285  const Real mol_surf = egs.getEquilibriumMolality(ind_surf);
3286  EXPECT_NEAR(mol_surf,
3287  act_bicarb / act_h / std::pow(10.0, egs.getLog10K(ind_surf)) *
3288  std::pow(alg_vars(2), -2.0),
3289  1.0E-9);
3290 
3291  // check that residuals are correctly calculated (non-surface stuff is checked in other tests)
3292  EXPECT_NEAR(
3293  egs.getResidualComponent(2, mole_add), -pref * 600.0 + 1.75 * (-1.0) * mol_surf, 1.0E-9);
3294  try
3295  {
3296  egs.getResidualComponent(3, mole_add);
3297  FAIL() << "Missing expected exception.";
3298  }
3299  catch (const std::exception & e)
3300  {
3301  std::string msg(e.what());
3302  ASSERT_TRUE(msg.find("Cannot retrieve residual for algebraic index 3 because there are only 2 "
3303  "molalities in the algebraic system and 1 surface potentials") !=
3304  std::string::npos)
3305  << "Failed with unexpected error message: " << msg;
3306  }
3307 }
3308 
3311 {
3312  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
3314  {"H2O", "H+", ">(s)FeOH", ">(w)FeOH", "Fe+++", "HCO3-"},
3315  {"Fe(OH)3(ppd)"},
3316  {},
3317  {},
3318  {},
3319  {},
3320  "O2(aq)",
3321  "e-");
3323  const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
3330  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
3337  const Real temp = 45.0;
3338  GeochemistryIonicStrength is2(1E-2, 1E-2, false, false);
3340  GeochemicalSystem egs(mgd,
3341  ac2,
3342  is2,
3343  _swapper6,
3344  {},
3345  {},
3346  "H+",
3347  {"H2O", "H+", ">(s)FeOH", ">(w)FeOH", "Fe+++", "HCO3-"},
3348  {1.75, 1.0, 2.0, 3.0, 0.5, 1.0},
3349  cu,
3350  cm,
3351  temp,
3352  0,
3353  1E-20,
3354  {},
3355  {},
3356  {});
3357  const DenseVector<Real> mole_add(6);
3358  const DenseMatrix<Real> dmole_add(6, 6);
3359 
3360  const unsigned num_alg = egs.getNumInAlgebraicSystem();
3361  DenseVector<Real> res_orig(num_alg);
3362  for (unsigned i = 0; i < num_alg; ++i)
3363  res_orig(i) = egs.getResidualComponent(i, mole_add);
3364  const std::vector<Real> var_orig = egs.getAlgebraicVariableValues();
3365  DenseMatrix<Real> jac(1, 1);
3366  egs.computeJacobian(res_orig, jac, mole_add, dmole_add);
3367 
3368  const Real eps = 1E-3;
3369  for (unsigned var_num = 0; var_num < num_alg; ++var_num)
3370  {
3371  std::vector<Real> var_new = var_orig;
3372  var_new[var_num] *= (1.0 + eps);
3373  egs.setAlgebraicVariables(var_new);
3374  for (unsigned res_comp = 0; res_comp < num_alg; ++res_comp)
3375  {
3376  const Real expected = (egs.getResidualComponent(res_comp, mole_add) - res_orig(res_comp)) /
3377  (var_new[var_num] - var_orig[var_num]);
3378  if (std::abs(expected) < 1.0E-7)
3379  EXPECT_NEAR(jac(res_comp, var_num), 0.0, 1.0E-7);
3380  else
3381  EXPECT_NEAR(expected / jac(res_comp, var_num), 1.0, 1.0E-2);
3382  }
3383  }
3384 }
3385 
3388 {
3389  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
3391  {"H2O", "H+", ">(s)FeOH", ">(w)FeOH", "Fe+++", "HCO3-"},
3392  {"Fe(OH)3(ppd)"},
3393  {},
3394  {"Something", "Fe(OH)3(ppd)fake"},
3395  {},
3396  {},
3397  "O2(aq)",
3398  "e-");
3400  const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
3407  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
3414  std::vector<GeochemistryUnitConverter::GeochemistryUnit> ku;
3417 
3418  const Real temp = 45.0;
3419  GeochemistryIonicStrength is2(1E-2, 1E-2, false, false);
3421  GeochemicalSystem egs(mgd,
3422  ac2,
3423  is2,
3424  _swapper6,
3425  {},
3426  {},
3427  "H+",
3428  {"H2O", "H+", ">(s)FeOH", ">(w)FeOH", "Fe+++", "HCO3-"},
3429  {1.75, 1.0, 2.0, 3.0, 0.5, 1.0},
3430  cu,
3431  cm,
3432  temp,
3433  0,
3434  1E-20,
3435  {"Something", "Fe(OH)3(ppd)fake"},
3436  {4.4, 5.5},
3437  ku);
3438  const DenseVector<Real> mole_add(8);
3439  const DenseMatrix<Real> dmole_add(8, 8);
3440 
3441  const unsigned num_alg = egs.getNumInAlgebraicSystem();
3442  DenseVector<Real> res_orig(num_alg);
3443  for (unsigned i = 0; i < num_alg; ++i)
3444  res_orig(i) = egs.getResidualComponent(i, mole_add);
3445  const std::vector<Real> var_orig = egs.getAlgebraicVariableValues();
3446  DenseMatrix<Real> jac(1, 1);
3447  egs.computeJacobian(res_orig, jac, mole_add, dmole_add);
3448 
3449  const Real eps = 1E-3;
3450  for (unsigned var_num = 0; var_num < num_alg; ++var_num)
3451  {
3452  std::vector<Real> var_new = var_orig;
3453  var_new[var_num] *= (1.0 + eps);
3454  egs.setAlgebraicVariables(var_new);
3455  for (unsigned res_comp = 0; res_comp < num_alg; ++res_comp)
3456  {
3457  const Real expected = (egs.getResidualComponent(res_comp, mole_add) - res_orig(res_comp)) /
3458  (var_new[var_num] - var_orig[var_num]);
3459  if (std::abs(expected) < 1.0E-7)
3460  EXPECT_NEAR(jac(res_comp, var_num), 0.0, 1.0E-7);
3461  else
3462  EXPECT_NEAR(expected / jac(res_comp, var_num), 1.0, 1.0E-2);
3463  }
3464  }
3465 }
3466 
3469 {
3470  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
3472  {"H2O", "H+", ">(s)FeOH", ">(w)FeOH", "Fe+++", "HCO3-"},
3473  {"Fe(OH)3(ppd)"},
3474  {},
3475  {"Something", "Fe(OH)3(ppd)fake"},
3476  {},
3477  {},
3478  "O2(aq)",
3479  "e-");
3481  const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
3488  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
3495  std::vector<GeochemistryUnitConverter::GeochemistryUnit> ku;
3498  const Real temp = 45.0;
3499  GeochemistryIonicStrength is2(1E-2, 1E-2, false, false);
3501  GeochemicalSystem egs(mgd,
3502  ac2,
3503  is2,
3504  _swapper6,
3505  {},
3506  {},
3507  "H+",
3508  {"H2O", "H+", ">(s)FeOH", ">(w)FeOH", "Fe+++", "HCO3-"},
3509  {1.75, 1.0, 2.0, 3.0, 0.5, 1.0},
3510  cu,
3511  cm,
3512  temp,
3513  0,
3514  1E-20,
3515  {"Something", "Fe(OH)3(ppd)fake"},
3516  {4.4, 5.5},
3517  ku);
3518  DenseVector<Real> mole_add(8);
3519  for (unsigned i = 0; i < 8; ++i)
3520  mole_add(i) = 1.0 + i;
3521  DenseMatrix<Real> dmole_add(8, 8);
3522  dmole_add(0, 0) = 1E15;
3523  for (unsigned i = 0; i < 8; ++i)
3524  for (unsigned j = 0; j < 8; ++j)
3525  dmole_add(i, j) = std::sin(i + 1.0) * (j + 1); // just some randomish derivatives
3526 
3527  const unsigned num_alg = egs.getNumInAlgebraicSystem();
3528  DenseVector<Real> res_orig(num_alg);
3529  for (unsigned i = 0; i < num_alg; ++i)
3530  res_orig(i) = egs.getResidualComponent(i, mole_add);
3531  const std::vector<Real> var_orig = egs.getAlgebraicVariableValues();
3532  DenseMatrix<Real> jac(1, 1);
3533  egs.computeJacobian(res_orig, jac, mole_add, dmole_add);
3534 
3535  const Real eps = 1E-3;
3536  for (unsigned var_num = 0; var_num < num_alg; ++var_num)
3537  {
3538  std::vector<Real> var_new = var_orig;
3539  var_new[var_num] *= (1.0 + eps);
3540  const Real del = var_new[var_num] - var_orig[var_num];
3541  egs.setAlgebraicVariables(var_new);
3542  for (unsigned i = 0; i < 8; ++i)
3543  {
3544  if (var_num < 3) // H2O, H+, >(s)FeOH
3545  mole_add(i) = 1.0 + i + dmole_add(i, var_num) * del;
3546  else if (var_num == 3) // HCO3-
3547  mole_add(i) = 1.0 + i + dmole_add(i, 5) * del;
3548  else if (var_num == 4) // surface potential
3549  mole_add(i) = 1.0 + i;
3550  else if (var_num == 5 || var_num == 6) // kinetic species
3551  mole_add(i) = 1.0 + i + dmole_add(i, var_num + 1) * del;
3552  }
3553 
3554  for (unsigned res_comp = 0; res_comp < num_alg; ++res_comp)
3555  {
3556  const Real expected =
3557  (egs.getResidualComponent(res_comp, mole_add) - res_orig(res_comp)) / del;
3558  if (std::abs(expected) < 1.0E-7)
3559  EXPECT_NEAR(jac(res_comp, var_num), 0.0, 1.0E-7);
3560  else
3561  EXPECT_NEAR(expected / jac(res_comp, var_num), 1.0, 1.0E-2);
3562  }
3563  }
3564 }
3565 
3567 TEST_F(GeochemicalSystemTest, setGetTemperature)
3568 {
3569  GeochemicalSystem nonconst = _egs_calcite;
3570  EXPECT_EQ(nonconst.getTemperature(), 25.0);
3571  nonconst.setTemperature(40.0);
3572  EXPECT_EQ(nonconst.getTemperature(), 40.0);
3573 }
3574 
3576 TEST_F(GeochemicalSystemTest, setMolalitiesExcept1)
3577 {
3578  GeochemicalSystem nonconst = _egs_calcite;
3579 
3580  try
3581  {
3583  {"H2O", "H+", "HCO3-", "Ca++", "Calcite", "CO2(aq)", "CO3--", "CaCO3", "CaOH+", "OH-"},
3584  {1.1, 2.2, 3.3, 4.4, 5.5, 6.6, 7.7, 8.8, 9.9},
3585  {false, true, true, true});
3586  FAIL() << "Missing expected exception.";
3587  }
3588  catch (const std::exception & e)
3589  {
3590  std::string msg(e.what());
3591  ASSERT_TRUE(msg.find("When setting all molalities, names and values must have same size") !=
3592  std::string::npos)
3593  << "Failed with unexpected error message: " << msg;
3594  }
3595 
3596  try
3597  {
3599  {"H2O", "H+", "HCO3-", "Ca++", "Calcite", "CO2(aq)", "CO3--", "CaCO3", "CaOH+"},
3600  {1.1, 2.2, 3.3, 4.4, 5.5, 6.6, 7.7, 8.8, 9.9},
3601  {false, true, true, true});
3602  FAIL() << "Missing expected exception.";
3603  }
3604  catch (const std::exception & e)
3605  {
3606  std::string msg(e.what());
3607  ASSERT_TRUE(
3608  msg.find(
3609  "When setting all molalities, values must be provided for every species and surface "
3610  "potentials") != std::string::npos)
3611  << "Failed with unexpected error message: " << msg;
3612  }
3613 
3614  try
3615  {
3617  {"H2O", "H+", "HCO3-", "Ca++", "Calcite", "CO2(aq)", "CO3--", "CaCO3", "CaOH+", "OH-"},
3618  {1.1, 2.2, 3.3, 4.4, 5.5, 6.6, 7.7, 8.8, 9.9, 10.1},
3619  {false, true, true});
3620  FAIL() << "Missing expected exception.";
3621  }
3622  catch (const std::exception & e)
3623  {
3624  std::string msg(e.what());
3625  ASSERT_TRUE(
3626  msg.find("constraints_from_molalities has size 3 while number of basis species is 4") !=
3627  std::string::npos)
3628  << "Failed with unexpected error message: " << msg;
3629  }
3630 
3631  try
3632  {
3634  {"H2O", "H+", "HCO3-", "Ca++", "Calcite", "CO2(aq)", "CO3--", "CaCO3", "CaOH+", "OH-"},
3635  {1.1, 2.2, 3.3, 4.4, -1.0, 6.6, 7.7, 8.8, 9.9, 10.1},
3636  {false, true, true, true});
3637  FAIL() << "Missing expected exception.";
3638  }
3639  catch (const std::exception & e)
3640  {
3641  std::string msg(e.what());
3642  ASSERT_TRUE(msg.find("Molality for mineral Calcite cannot be -1: it must be non-negative") !=
3643  std::string::npos)
3644  << "Failed with unexpected error message: " << msg;
3645  }
3646 
3647  try
3648  {
3650  {"H2O", "H+", "HCO3-", "Ca++", "Calcite", "CO2(aq)", "CO3--", "CaCO3", "CaOH+", "OH-"},
3651  {1.1, 2.2, 3.3, -1.0, 1.0, 6.6, 7.7, 8.8, 9.9, 10.1},
3652  {false, true, true, true});
3653  FAIL() << "Missing expected exception.";
3654  }
3655  catch (const std::exception & e)
3656  {
3657  std::string msg(e.what());
3658  ASSERT_TRUE(msg.find("Molality for species Ca++ cannot be -1: it must be non-negative") !=
3659  std::string::npos)
3660  << "Failed with unexpected error message: " << msg;
3661  }
3662 
3663  try
3664  {
3666  {"H2O", "H+", "HCO3-", "Ca++", "Calcite", "CO2(aq)", "CO3--", "CaCO3", "CaOH+", "OH-"},
3667  {1.1, 0.0, 3.3, 1.0, 1.0, 6.6, 7.7, 8.8, 9.9, 10.1},
3668  {false, true, true, true});
3669  FAIL() << "Missing expected exception.";
3670  }
3671  catch (const std::exception & e)
3672  {
3673  std::string msg(e.what());
3674  ASSERT_TRUE(msg.find("Molality for species H+ cannot be 0: it must be positive") !=
3675  std::string::npos)
3676  << "Failed with unexpected error message: " << msg;
3677  }
3678 
3679  try
3680  {
3682  {"OH-", "H+", "HCO3-", "Ca++", "Calcite", "CO2(aq)", "CO3--", "CaCO3", "CaOH+", "OH-"},
3683  {1.1, 2.2, 3.3, 4.4, 1.0, 6.6, 7.7, 8.8, 9.9, 10.1},
3684  {false, true, true, true});
3685  FAIL() << "Missing expected exception.";
3686  }
3687  catch (const std::exception & e)
3688  {
3689  std::string msg(e.what());
3690  ASSERT_TRUE(msg.find("Molality (or free mineral moles, etc - whatever is appropriate) for "
3691  "species H2O needs to be provided when setting all molalities") !=
3692  std::string::npos)
3693  << "Failed with unexpected error message: " << msg;
3694  }
3695 
3696  try
3697  {
3699  {"H2O", "H+", "HCO3-", "Ca++", "Calcite", "CO2(aq)", "CO3--", "CaCO3", "CaOH+", "OH-"},
3700  {1.1, 2.2, 3.3, 4.4, 1.0, -1.0, 7.7, 8.8, 9.9, 10.1},
3701  {false, true, true, true});
3702  FAIL() << "Missing expected exception.";
3703  }
3704  catch (const std::exception & e)
3705  {
3706  std::string msg(e.what());
3707  ASSERT_TRUE(msg.find("Molality for species CO2(aq) cannot be -1: it must be non-negative") !=
3708  std::string::npos)
3709  << "Failed with unexpected error message: " << msg;
3710  }
3711 
3712  try
3713  {
3715  {"H2O", "H+", "HCO3-", "Ca++", "Calcite", "CO2(aq)", "CO2(aq)", "CaCO3", "CaOH+", "OH-"},
3716  {1.1, 2.2, 3.3, 4.4, 1.0, 6.6, 7.7, 8.8, 9.9, 10.1},
3717  {false, true, true, true});
3718  FAIL() << "Missing expected exception.";
3719  }
3720  catch (const std::exception & e)
3721  {
3722  std::string msg(e.what());
3723  ASSERT_TRUE(
3724  msg.find("Molality for species CO3-- needs to be provided when setting all molalities") !=
3725  std::string::npos)
3726  << "Failed with unexpected error message: " << msg;
3727  }
3728 
3729  try
3730  {
3732  {"H2O", "H+", "HCO3-", "Ca++", "Calcite", "CO2(aq)", "CO3--", "CaCO3", "CaOH+", "OH-"},
3733  {1.1, 2.2, 3.3, 4.4, 1.0, 6.6, 7.7, 8.8, 9.9, 10.1},
3734  {true, true, true, true});
3735  FAIL() << "Missing expected exception.";
3736  }
3737  catch (const std::exception & e)
3738  {
3739  std::string msg(e.what());
3740  ASSERT_TRUE(
3741  msg.find("Water activity cannot be determined from molalities: constraints_from_molalities "
3742  "must be set to false for water if activity of water is fixed") !=
3743  std::string::npos)
3744  << "Failed with unexpected error message: " << msg;
3745  }
3746 }
3747 
3749 TEST_F(GeochemicalSystemTest, setMolalitiesExcept2)
3750 {
3751  const PertinentGeochemicalSystem model_gas_and_sorption(
3752  _db_calcite,
3753  {"H2O", "H+", "HCO3-", "O2(aq)", "Fe+++", ">(s)FeOH", ">(w)FeOH"},
3754  {"Goethite"},
3755  {"O2(g)"},
3756  {},
3757  {},
3758  {},
3759  "O2(aq)",
3760  "e-");
3761  ModelGeochemicalDatabase mgd = model_gas_and_sorption.modelGeochemicalDatabase();
3762  const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
3770  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
3778  GeochemicalSystem nonconst(mgd,
3779  _ac3,
3780  _is3,
3781  _swapper7,
3782  {"O2(aq)"},
3783  {"O2(g)"},
3784  "H+",
3785  {"H2O", "H+", "HCO3-", "O2(g)", "Fe+++", ">(s)FeOH", ">(w)FeOH"},
3786  {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0},
3787  cu,
3788  cm,
3789  25,
3790  0,
3791  1E-20,
3792  {},
3793  {},
3794  {});
3795 
3796  try
3797  {
3799  {"H2O",
3800  "H+",
3801  "HCO3-",
3802  "O2(g)",
3803  "Fe+++",
3804  ">(s)FeOH",
3805  ">(w)FeOH",
3806  "(O-phth)--",
3807  "CH4(aq)",
3808  "CO2(aq)",
3809  "CO3--",
3810  "OH-",
3811  ">(s)FeO-",
3812  "Goethite",
3813  "O2(aq)"},
3814  {1.1, 2.2, 3.3, 0.0, 5.5, 6.6, 7.7, 8.8, 9.9, 10, 11, 12, 14, 15, 16},
3815  {false, true, true, false, true, true, true});
3816  FAIL() << "Missing expected exception.";
3817  }
3818  catch (const std::exception & e)
3819  {
3820  std::string msg(e.what());
3821  ASSERT_TRUE(msg.find("When setting all molalities, values must be provided for every species "
3822  "and surface potentials") != std::string::npos)
3823  << "Failed with unexpected error message: " << msg;
3824  }
3825 
3826  try
3827  {
3828  nonconst.setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(
3829  {"H2O",
3830  "H+",
3831  "HCO3-",
3832  "O2(g)",
3833  "Fe+++",
3834  ">(s)FeOH",
3835  ">(w)FeOH",
3836  "(O-phth)--",
3837  "CH4(aq)",
3838  "CO2(aq)",
3839  "CO3--",
3840  "OH-",
3841  ">(s)FeO-",
3842  "Goethite",
3843  "O2(aq)",
3844  "Goethite_surface_potential_expr"},
3845  {1.1, 2.2, 3.3, 1.0, 5.5, 6.6, 7.7, 8.8, 9.9, 10, 11, 12, 14, 15, 16, 17},
3846  {false, true, true, false, true, true, true});
3847  FAIL() << "Missing expected exception.";
3848  }
3849  catch (const std::exception & e)
3850  {
3851  std::string msg(e.what());
3852  ASSERT_TRUE(msg.find("Molality for gas O2(g) cannot be 1: it must be zero") !=
3853  std::string::npos)
3854  << "Failed with unexpected error message: " << msg;
3855  }
3856 
3857  try
3858  {
3859  nonconst.setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(
3860  {"H2O",
3861  "H+",
3862  "HCO3-",
3863  "O2(g)",
3864  "Fe+++",
3865  ">(s)FeOH",
3866  ">(w)FeOH",
3867  "(O-phth)--",
3868  "CH4(aq)",
3869  "CO2(aq)",
3870  "CO3--",
3871  "OH-",
3872  ">(s)FeO-",
3873  "Goethite",
3874  "O2(aq)",
3875  "Goethite_surface_potential_expr"},
3876  {1.1, 2.2, 3.3, 0.0, 5.5, 6.6, 7.7, 8.8, 9.9, 10, 11, 12, 14, 15, -1.0, 17},
3877  {false, true, true, false, true, true, true});
3878  FAIL() << "Missing expected exception.";
3879  }
3880  catch (const std::exception & e)
3881  {
3882  std::string msg(e.what());
3883  ASSERT_TRUE(msg.find("Molality for species O2(aq) cannot be -1: it must be non-negative") !=
3884  std::string::npos)
3885  << "Failed with unexpected error message: " << msg;
3886  }
3887 
3888  try
3889  {
3890  nonconst.setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(
3891  {"H2O",
3892  "H+",
3893  "HCO3-",
3894  "O2(g)",
3895  "Fe+++",
3896  ">(s)FeOH",
3897  ">(w)FeOH",
3898  "(O-phth)--",
3899  "CH4(aq)",
3900  "CO2(aq)",
3901  "CO3--",
3902  "OH-",
3903  ">(s)FeO-",
3904  "Goethite",
3905  "O2(aq)",
3906  "Goethite"},
3907  {1.1, 2.2, 3.3, 0.0, 5.5, 6.6, 7.7, 8.8, 9.9, 10, 11, 12, 14, 15, 16, 17},
3908  {false, true, true, false, true, true, true});
3909  FAIL() << "Missing expected exception.";
3910  }
3911  catch (const std::exception & e)
3912  {
3913  std::string msg(e.what());
3914  ASSERT_TRUE(msg.find("Surface potential for mineral Goethite needs to be provided when setting "
3915  "all molalities") != std::string::npos)
3916  << "Failed with unexpected error message: " << msg;
3917  }
3918 
3919  try
3920  {
3921  nonconst.setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(
3922  {"H2O",
3923  "H+",
3924  "HCO3-",
3925  "O2(g)",
3926  "Fe+++",
3927  ">(s)FeOH",
3928  ">(w)FeOH",
3929  "(O-phth)--",
3930  "CH4(aq)",
3931  "CO2(aq)",
3932  "CO3--",
3933  "OH-",
3934  ">(s)FeO-",
3935  "Goethite",
3936  "O2(aq)",
3937  "Goethite_surface_potential_expr"},
3938  {1.1, 2.2, 3.3, 0.0, 5.5, 6.6, 7.7, 8.8, 9.9, 10, 11, 12, 14, 15, 16, 0},
3939  {false, true, true, false, true, true, true});
3940  FAIL() << "Missing expected exception.";
3941  }
3942  catch (const std::exception & e)
3943  {
3944  std::string msg(e.what());
3945  ASSERT_TRUE(msg.find("Surface-potential expression for mineral Goethite cannot be 0: it must "
3946  "be positive") != std::string::npos)
3947  << "Failed with unexpected error message: " << msg;
3948  }
3949 
3950  try
3951  {
3952  nonconst.setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(
3953  {"H2O",
3954  "H+",
3955  "HCO3-",
3956  "O2(g)",
3957  "Fe+++",
3958  ">(s)FeOH",
3959  ">(w)FeOH",
3960  "(O-phth)--",
3961  "CH4(aq)",
3962  "CO2(aq)",
3963  "CO3--",
3964  "OH-",
3965  ">(s)FeO-",
3966  "Goethite",
3967  "O2(aq)",
3968  "Goethite_surface_potential_expr"},
3969  {1.1, 2.2, 3.3, 0.0, 5.5, 6.6, 7.7, 8.8, 9.9, 10, 11, 12, 14, 15, 16, 17},
3970  {false, true, true, true, true, true, true});
3971  FAIL() << "Missing expected exception.";
3972  }
3973  catch (const std::exception & e)
3974  {
3975  std::string msg(e.what());
3976  ASSERT_TRUE(msg.find("Gas fugacity cannot be determined from molality: "
3977  "constraints_from_molalities must be set false for all gases") !=
3978  std::string::npos)
3979  << "Failed with unexpected error message: " << msg;
3980  }
3981 }
3982 
3984 TEST_F(GeochemicalSystemTest, setMolalitiesExcept3)
3985 {
3986  GeochemicalSystem nonconst = _egs_kinetic_calcite;
3987  try
3988  {
3990  {"H2O", "H+", "HCO3-", "Ca++", "CO2(aq)", "CO3--", "CaCO3", "CaOH+", "OH-", "bad"},
3991  {1.1, 2.2, 3.3, 4.4, 5.5, 6.6, 7.7, 8.8, 9.9, 10.0},
3992  {false, true, true, true});
3993  FAIL() << "Missing expected exception.";
3994  }
3995  catch (const std::exception & e)
3996  {
3997  std::string msg(e.what());
3998  ASSERT_TRUE(
3999  msg.find("Moles for species Calcite needs to be provided when setting all molalities") !=
4000  std::string::npos)
4001  << "Failed with unexpected error message: " << msg;
4002  }
4003 
4004  try
4005  {
4007  {"H2O", "H+", "HCO3-", "Ca++", "CO2(aq)", "CO3--", "CaCO3", "CaOH+", "OH-", "Calcite"},
4008  {1.1, 2.2, 3.3, 4.4, 5.5, 6.6, 7.7, 8.8, 9.9, 0.0},
4009  {false, true, true, true});
4010  FAIL() << "Missing expected exception.";
4011  }
4012  catch (const std::exception & e)
4013  {
4014  std::string msg(e.what());
4015  ASSERT_TRUE(msg.find("Mole number for kinetic species must be positive, not 0") !=
4016  std::string::npos)
4017  << "Failed with unexpected error message: " << msg;
4018  }
4019 }
4020 
4023 {
4024  const PertinentGeochemicalSystem model_gas_and_sorption(
4025  _db_calcite,
4026  {"H2O", "H+", "HCO3-", "O2(aq)", "Fe+++", ">(s)FeOH", ">(w)FeOH"},
4027  {"Goethite"},
4028  {"O2(g)"},
4029  {},
4030  {},
4031  {},
4032  "O2(aq)",
4033  "e-");
4034  ModelGeochemicalDatabase mgd = model_gas_and_sorption.modelGeochemicalDatabase();
4035  const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
4043  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
4051  GeochemicalSystem nonconst(mgd,
4052  _ac0,
4053  _is0,
4054  _swapper7,
4055  {"O2(aq)"},
4056  {"O2(g)"},
4057  "Fe+++",
4058  {"H2O", "H+", "HCO3-", "O2(g)", "Fe+++", ">(s)FeOH", ">(w)FeOH"},
4059  {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0},
4060  cu,
4061  cm,
4062  25,
4063  0,
4064  1E-20,
4065  {},
4066  {},
4067  {});
4068  const std::vector<Real> bulk_before_set = nonconst.getBulkMolesOld();
4069  const std::vector<Real> molal_before_set =
4070  nonconst.getSolventMassAndFreeMolalityAndMineralMoles();
4071 
4072  const std::vector<std::string> gold_name = {"H2O",
4073  "H+",
4074  "HCO3-",
4075  "O2(g)",
4076  "Fe+++",
4077  ">(s)FeOH",
4078  ">(w)FeOH",
4079  "(O-phth)--",
4080  "CH4(aq)",
4081  "CO2(aq)",
4082  "CO3--",
4083  "OH-",
4084  ">(s)FeO-",
4085  "Goethite",
4086  "O2(aq)",
4087  "Goethite_surface_potential_expr"};
4088  const std::vector<Real> set_molal = {
4089  1.1, 2.2, 3.3, 0.0, 5.5, 6.6, 7.7, 8.8, 9.9, 10, 11, 12, 14, 15, 16, 17};
4090  nonconst.setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(
4091  gold_name, set_molal, {true, true, true, false, false, true, false});
4092 
4093  const unsigned num_basis = nonconst.getNumInBasis();
4094  const unsigned num_eqm = nonconst.getNumInEquilibrium();
4095 
4096  // check molalities:
4097  std::vector<Real> gold_molal = set_molal;
4098  gold_molal[6] = molal_before_set[6]; // >(w)FeOH has a FREE_MOLALITY constraint and the
4099  // setSolvent... method has false
4100  gold_molal[13] = 0.0; // the setSolvent... method explicitly sets secondary mineral molality zero
4101  const std::vector<Real> molal = nonconst.getSolventMassAndFreeMolalityAndMineralMoles();
4102  for (unsigned i = 0; i < num_basis; ++i)
4103  EXPECT_EQ(molal[mgd.basis_species_index[gold_name[i]]], gold_molal[i]);
4104  for (unsigned j = num_basis; j < num_basis + num_eqm; ++j)
4105  EXPECT_EQ(nonconst.getEquilibriumMolality(mgd.eqm_species_index[gold_name[j]]), gold_molal[j]);
4106  EXPECT_NEAR(nonconst.getSurfacePotential(0),
4109  GeochemistryConstants::FARADAY * std::log(17.0),
4110  1.0E-8);
4111  EXPECT_EQ(nonconst.getSolventWaterMass(), 1.1);
4112 
4113  // check activities
4114  const std::vector<bool> act_known_gold = {false, true, false, true, false, false, false};
4115  const std::vector<bool> act_known = nonconst.getBasisActivityKnown();
4116  for (unsigned i = 0; i < num_basis; ++i)
4117  EXPECT_EQ(act_known[i], act_known_gold[i]);
4118  EXPECT_EQ(nonconst.getBasisActivity(1),
4119  2.2); // H+: activity changed by setSolvent... (activity coeffs = 1)
4120  EXPECT_EQ(nonconst.getBasisActivity(3), 4.0); // O2(g): fugacity unchanged by setSolvent...
4121 
4122  // check sorbing surface area
4123  EXPECT_EQ(nonconst.getSorbingSurfaceArea()[0], 60.0); // Goethite is not in the basis
4124 
4125  // check bulk moles
4126  const std::vector<Real> bulk = nonconst.getBulkMolesOld();
4127  for (unsigned i = 0; i < num_basis; ++i)
4128  {
4129  Real b = 0.0;
4130  if (i == 0)
4132  else
4133  b = gold_molal[i];
4134  for (unsigned j = 0; j < num_eqm; ++j)
4135  b += mgd.eqm_stoichiometry(j, i) * nonconst.getEquilibriumMolality(j);
4136  b *= gold_molal[0];
4137  if (i == 4)
4138  EXPECT_EQ(bulk[4], bulk_before_set[4]); // Fe+++ is set by original constraints
4139  else
4140  EXPECT_NEAR(bulk[mgd.basis_species_index[gold_name[i]]] / b, 1.0, 1.0E-8);
4141  }
4142 }
4143 
4146 {
4147  const PertinentGeochemicalSystem model_gas_and_sorption(
4148  _db_calcite,
4149  {"H2O", "H+", "HCO3-", "O2(aq)", "Fe+++", ">(s)FeOH", ">(w)FeOH"},
4150  {"Goethite"},
4151  {"O2(g)"},
4152  {},
4153  {},
4154  {},
4155  "O2(aq)",
4156  "e-");
4157  ModelGeochemicalDatabase mgd = model_gas_and_sorption.modelGeochemicalDatabase();
4158  const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
4166  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
4174  GeochemicalSystem nonconst(mgd,
4175  _ac3,
4176  _is3,
4177  _swapper7,
4178  {"O2(aq)"},
4179  {"O2(g)"},
4180  "Fe+++",
4181  {"H2O", "H+", "HCO3-", "O2(g)", "Fe+++", ">(s)FeOH", ">(w)FeOH"},
4182  {11.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0},
4183  cu,
4184  cm,
4185  25,
4186  0,
4187  1E-20,
4188  {},
4189  {},
4190  {});
4191  const std::vector<Real> bulk_before_set = nonconst.getBulkMolesOld();
4192  const std::vector<Real> molal_before_set =
4193  nonconst.getSolventMassAndFreeMolalityAndMineralMoles();
4194 
4195  const std::vector<std::string> gold_name = {"H2O",
4196  "H+",
4197  "HCO3-",
4198  "O2(g)",
4199  "Fe+++",
4200  ">(s)FeOH",
4201  ">(w)FeOH",
4202  "(O-phth)--",
4203  "CH4(aq)",
4204  "CO2(aq)",
4205  "CO3--",
4206  "OH-",
4207  ">(s)FeO-",
4208  "Goethite",
4209  "O2(aq)",
4210  "Goethite_surface_potential_expr"};
4211  const std::vector<Real> set_molal = {
4212  1.1, 2.2, 3.3, 0.0, 5.5, 6.6, 7.7, 8.8, 9.9, 10, 11, 12, 14, 15, 16, 17};
4213  nonconst.setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(
4214  gold_name, set_molal, {false, false, true, false, false, true, false});
4215 
4216  const unsigned num_basis = nonconst.getNumInBasis();
4217  const unsigned num_eqm = nonconst.getNumInEquilibrium();
4218 
4219  // check molalities:
4220  std::vector<Real> gold_molal = set_molal;
4221  gold_molal[1] =
4222  molal_before_set[1]; // H+ has an ACTIVITY constraint and the setSolvent... method has false
4223  gold_molal[6] = molal_before_set[6]; // >(w)FeOH has a FREE_MOLALITY constraint and the
4224  // setSolvent... method has false
4225  gold_molal[13] = 0.0; // the setSolvent... method explicitly sets secondary mineral molality zero
4226  const std::vector<Real> molal = nonconst.getSolventMassAndFreeMolalityAndMineralMoles();
4227  for (unsigned i = 0; i < num_basis; ++i)
4228  EXPECT_EQ(molal[mgd.basis_species_index[gold_name[i]]], gold_molal[i]);
4229  for (unsigned j = num_basis; j < num_basis + num_eqm; ++j)
4230  EXPECT_EQ(nonconst.getEquilibriumMolality(mgd.eqm_species_index[gold_name[j]]), gold_molal[j]);
4231  EXPECT_NEAR(nonconst.getSurfacePotential(0),
4234  GeochemistryConstants::FARADAY * std::log(17.0),
4235  1.0E-8);
4236  EXPECT_EQ(nonconst.getSolventWaterMass(), 1.1);
4237 
4238  // check activities
4239  const std::vector<bool> act_known_gold = {false, true, false, true, false, false, false};
4240  const std::vector<bool> act_known = nonconst.getBasisActivityKnown();
4241  for (unsigned i = 0; i < num_basis; ++i)
4242  EXPECT_EQ(act_known[i], act_known_gold[i]);
4243  EXPECT_EQ(nonconst.getBasisActivity(1),
4244  2.0); // H+: activity unchanged by setSolvent... (activity coeffs = 1)
4245  EXPECT_EQ(nonconst.getBasisActivity(3), 4.0); // O2(g): fugacity unchanged by setSolvent...
4246 
4247  // check sorbing surface area
4248  EXPECT_EQ(nonconst.getSorbingSurfaceArea()[0], 60.0); // Goethite is not in the basis
4249 
4250  // check bulk moles
4251  const std::vector<Real> bulk = nonconst.getBulkMolesOld();
4252  for (unsigned i = 0; i < num_basis; ++i)
4253  {
4254  Real b = 0.0;
4255  if (i == 0)
4257  else
4258  b = gold_molal[i];
4259  for (unsigned j = 0; j < num_eqm; ++j)
4260  b += mgd.eqm_stoichiometry(j, i) * nonconst.getEquilibriumMolality(j);
4261  b *= gold_molal[0];
4262  if (i == 0 || i == 4)
4263  EXPECT_EQ(bulk[i], bulk_before_set[i]); // H2O and Fe+++ are set by original constraints
4264  else
4265  EXPECT_NEAR(bulk[mgd.basis_species_index[gold_name[i]]] / b, 1.0, 1.0E-8);
4266  }
4267 }
4268 
4270 TEST_F(GeochemicalSystemTest, getConstraintMeaning)
4271 {
4272  ModelGeochemicalDatabase mgd = _model_calcite.modelGeochemicalDatabase();
4273  const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
4278  const std::vector<GeochemicalSystem::ConstraintMeaningEnum> cim = {
4283  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
4288  const GeochemicalSystem egs(mgd,
4289  _ac3,
4290  _is3,
4291  _swapper4,
4292  {"Ca++"},
4293  {"Calcite"},
4294  "H+",
4295  {"H2O", "H+", "HCO3-", "Calcite"},
4296  {1.75, 3.0, 2.0, 1.0},
4297  cu,
4298  cm,
4299  25,
4300  0,
4301  1E-20,
4302  {},
4303  {},
4304  {});
4305  for (unsigned i = 0; i < 4; ++i)
4306  EXPECT_EQ(egs.getConstraintMeaning()[i], cim[i]);
4307 }
4308 
4310 TEST_F(GeochemicalSystemTest, changeContraintToBulkExceptions)
4311 {
4312  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
4314  database, {"H2O", "H+", "HCO3-", "O2(aq)"}, {}, {"O2(g)"}, {}, {}, {}, "O2(aq)", "e-");
4316 
4317  std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
4322  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
4327 
4328  GeochemicalSystem egs(mgd,
4329  _ac3,
4330  _is3,
4331  _swapper4,
4332  {"O2(aq)"},
4333  {"O2(g)"},
4334  "H+",
4335  {"H2O", "H+", "O2(g)", "HCO3-"},
4336  {1.0, 2.0, 3.0, 4.0},
4337  cu,
4338  cm,
4339  25,
4340  0,
4341  1E-20,
4342  {},
4343  {},
4344  {});
4345 
4346  try
4347  {
4348  egs.changeConstraintToBulk(4);
4349  FAIL() << "Missing expected exception.";
4350  }
4351  catch (const std::exception & e)
4352  {
4353  std::string msg(e.what());
4354  ASSERT_TRUE(msg.find("Cannot changeConstraintToBulk for species 4 because there are only 4 "
4355  "basis species") != std::string::npos)
4356  << "Failed with unexpected error message: " << msg;
4357  }
4358 
4359  try
4360  {
4361  egs.changeConstraintToBulk(4, 1.0);
4362  FAIL() << "Missing expected exception.";
4363  }
4364  catch (const std::exception & e)
4365  {
4366  std::string msg(e.what());
4367  ASSERT_TRUE(msg.find("Cannot changeConstraintToBulk for species 4 because there are only 4 "
4368  "basis species") != std::string::npos)
4369  << "Failed with unexpected error message: " << msg;
4370  }
4371 
4372  try
4373  {
4374  egs.changeConstraintToBulk(3, 1.0);
4375  FAIL() << "Missing expected exception.";
4376  }
4377  catch (const std::exception & e)
4378  {
4379  std::string msg(e.what());
4380  ASSERT_TRUE(msg.find("Attempting to changeConstraintToBulk for a gas species. Since a swap is "
4381  "involved, you cannot specify a value for the bulk number of moles. You "
4382  "must use changeConstraintToBulk(basis_ind) method instead of "
4383  "changeConstraintToBulk(basis_ind, value)") != std::string::npos)
4384  << "Failed with unexpected error message: " << msg;
4385  }
4386 }
4387 
4390 {
4391  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
4393  database, {"H2O", "H+", "HCO3-", "O2(aq)"}, {}, {"O2(g)"}, {}, {}, {}, "O2(aq)", "e-");
4395 
4396  std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
4401  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
4406 
4407  GeochemicalSystem egs(mgd,
4408  _ac3,
4409  _is3,
4410  _swapper4,
4411  {"O2(aq)"},
4412  {"O2(g)"},
4413  "H+",
4414  {"H2O", "H+", "HCO3-", "O2(g)"},
4415  {1.0, 2.0, 3.0, 4.0},
4416  cu,
4417  cm,
4418  25,
4419  0,
4420  1E-20,
4421  {},
4422  {},
4423  {});
4424  egs.closeSystem();
4425  const std::vector<GeochemicalSystem::ConstraintMeaningEnum> cim = {
4430  for (unsigned i = 0; i < 4; ++i)
4431  EXPECT_EQ(egs.getConstraintMeaning()[i], cim[i]);
4432 }
4433 
4435 TEST_F(GeochemicalSystemTest, changeConstraintToBulk)
4436 {
4437  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
4439  database, {"H2O", "H+", "HCO3-", "O2(aq)"}, {}, {"O2(g)"}, {}, {}, {}, "O2(aq)", "e-");
4441 
4442  std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
4447  std::vector<GeochemicalSystem::ConstraintMeaningEnum> cim;
4452  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
4457 
4458  GeochemicalSystem egs(mgd,
4459  _ac3,
4460  _is3,
4461  _swapper4,
4462  {"O2(aq)"},
4463  {"O2(g)"},
4464  "H+",
4465  {"H2O", "H+", "HCO3-", "O2(g)"},
4466  {1.0, 2.0, 3.0, 4.0},
4467  cu,
4468  cm,
4469  25,
4470  0,
4471  1E-20,
4472  {},
4473  {},
4474  {});
4475  const std::vector<Real> orig_bulk = egs.getBulkMolesOld();
4476 
4477  egs.changeConstraintToBulk(0);
4479  for (unsigned i = 0; i < 4; ++i)
4480  EXPECT_EQ(egs.getConstraintMeaning()[i], cim[i]);
4481  for (unsigned i = 0; i < 4; ++i)
4482  EXPECT_EQ(egs.getBulkMolesOld()[i], orig_bulk[i]);
4483 
4484  std::vector<Real> new_bulk = orig_bulk;
4485  egs.changeConstraintToBulk(1, 1.1); // just resets the bulk constraint value
4486  new_bulk[1] = 1.1;
4487  for (unsigned i = 0; i < 4; ++i)
4488  EXPECT_EQ(egs.getConstraintMeaning()[i], cim[i]);
4489  for (unsigned i = 0; i < 4; ++i)
4490  EXPECT_EQ(egs.getBulkMolesOld()[i], new_bulk[i]);
4491 
4492  egs.changeConstraintToBulk(2); // this means charge-balance can be enforced simply: see next line
4493  new_bulk[1] = new_bulk[2];
4495  for (unsigned i = 0; i < 4; ++i)
4496  EXPECT_EQ(egs.getConstraintMeaning()[i], cim[i]);
4497  for (unsigned i = 0; i < 4; ++i)
4498  EXPECT_EQ(egs.getBulkMolesOld()[i], new_bulk[i]);
4499 
4500  egs.changeConstraintToBulk(3);
4502  for (unsigned i = 0; i < 4; ++i)
4503  EXPECT_EQ(egs.getConstraintMeaning()[i], cim[i]);
4504  EXPECT_EQ(mgd.basis_species_name[3], "(O-phth)--");
4505 }
4506 
4508 TEST_F(GeochemicalSystemTest, addToBulkMolesExceptions)
4509 {
4510  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
4512  database, {"H2O", "H+", "HCO3-", "O2(aq)"}, {}, {"O2(g)"}, {}, {}, {}, "O2(aq)", "e-");
4514 
4515  std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
4520  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
4525 
4526  GeochemicalSystem egs(mgd,
4527  _ac3,
4528  _is3,
4529  _swapper4,
4530  {"O2(aq)"},
4531  {"O2(g)"},
4532  "H+",
4533  {"H2O", "H+", "HCO3-", "O2(g)"},
4534  {1.0, 2.0, 3.0, 4.0},
4535  cu,
4536  cm,
4537  25,
4538  0,
4539  1E-20,
4540  {},
4541  {},
4542  {});
4543 
4544  try
4545  {
4546  egs.addToBulkMoles(4, 1.0);
4547  FAIL() << "Missing expected exception.";
4548  }
4549  catch (const std::exception & e)
4550  {
4551  std::string msg(e.what());
4552  ASSERT_TRUE(
4553  msg.find("Cannot addToBulkMoles for species 4 because there are only 4 basis species") !=
4554  std::string::npos)
4555  << "Failed with unexpected error message: " << msg;
4556  }
4557 }
4558 
4561 {
4562  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
4564  database, {"H2O", "H+", "O2(aq)", "HCO3-"}, {}, {"O2(g)"}, {}, {}, {}, "O2(aq)", "e-");
4566 
4567  std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
4572  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
4577 
4578  std::vector<Real> bm = {1.0, 2.0, 3.0, 4.0};
4579 
4580  GeochemicalSystem egs(mgd,
4581  _ac3,
4582  _is3,
4583  _swapper4,
4584  {},
4585  {},
4586  "H+",
4587  {"H2O", "H+", "O2(aq)", "HCO3-"},
4588  bm,
4589  cu,
4590  cm,
4591  25,
4592  0,
4593  1E-20,
4594  {},
4595  {},
4596  {});
4597 
4598  egs.addToBulkMoles(0, 1.1);
4599  bm[0] += 1.1;
4600  for (unsigned i = 0; i < 3; ++i)
4601  EXPECT_EQ(egs.getBulkMolesOld()[i], bm[i]);
4602  egs.addToBulkMoles(1, 2.2);
4603  bm[1] += 2.2;
4604  for (unsigned i = 0; i < 3; ++i)
4605  EXPECT_EQ(egs.getBulkMolesOld()[i], bm[i]);
4606  egs.addToBulkMoles(2, 3.3);
4607  bm[2] += 3.3;
4608  for (unsigned i = 0; i < 3; ++i)
4609  EXPECT_EQ(egs.getBulkMolesOld()[i], bm[i]);
4610 }
4611 
4613 TEST_F(GeochemicalSystemTest, setConstraintValueExceptions)
4614 {
4615  GeochemicalSystem nonconst = _egs_calcite;
4616  try
4617  {
4618  nonconst.setConstraintValue(4, 1.0);
4619  FAIL() << "Missing expected exception.";
4620  }
4621  catch (const std::exception & e)
4622  {
4623  std::string msg(e.what());
4624  ASSERT_TRUE(
4625  msg.find(
4626  "Cannot setConstraintValue for species 4 because there are only 4 basis species") !=
4627  std::string::npos)
4628  << "Failed with unexpected error message: " << msg;
4629  }
4630 }
4631 
4633 TEST_F(GeochemicalSystemTest, setConstraintValue)
4634 {
4635  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
4637  database, {"H2O", "H+", "O2(aq)", "HCO3-"}, {}, {"O2(g)"}, {}, {}, {}, "O2(aq)", "e-");
4639 
4640  std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
4645  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
4650 
4651  std::vector<Real> bm_in = {1.0, 2.0, 3.0, 4.0};
4652 
4653  GeochemicalSystem egs(mgd,
4654  _ac3,
4655  _is3,
4656  _swapper4,
4657  {},
4658  {},
4659  "H+",
4660  {"H2O", "H+", "O2(aq)", "HCO3-"},
4661  bm_in,
4662  cu,
4663  cm,
4664  25,
4665  0,
4666  1E-20,
4667  {},
4668  {},
4669  {});
4670 
4671  std::vector<Real> bm = egs.getBulkMolesOld();
4672  std::vector<Real> mol = egs.getSolventMassAndFreeMolalityAndMineralMoles();
4673 
4674  egs.setConstraintValue(0, 1.1);
4675  bm[0] = 1.1;
4676  for (unsigned i = 0; i < 4; ++i)
4677  EXPECT_EQ(egs.getBulkMolesOld()[i], bm[i]);
4678  for (unsigned i = 0; i < 4; ++i)
4679  EXPECT_EQ(egs.getSolventMassAndFreeMolalityAndMineralMoles()[i], mol[i]);
4680 
4681  egs.setConstraintValue(1, 2.2);
4682  bm[1] = 2.2;
4683  for (unsigned i = 0; i < 4; ++i)
4684  EXPECT_EQ(egs.getBulkMolesOld()[i], bm[i]);
4685  for (unsigned i = 0; i < 4; ++i)
4686  EXPECT_EQ(egs.getSolventMassAndFreeMolalityAndMineralMoles()[i], mol[i]);
4687 
4688  egs.setConstraintValue(3, 4.4);
4689  mol[3] = 4.4;
4690  bm = egs.getBulkMolesOld(); // setting the molality to 4.4 changes the bulk composition
4691  for (unsigned i = 0; i < 4; ++i)
4692  EXPECT_EQ(egs.getBulkMolesOld()[i], bm[i]);
4693  for (unsigned i = 0; i < 4; ++i)
4694  EXPECT_EQ(egs.getSolventMassAndFreeMolalityAndMineralMoles()[i], mol[i]);
4695 
4696  egs.setConstraintValue(2, 3.3);
4697  bm = egs.getBulkMolesOld(); // setting the activity to 3.3 changes the basis molality and hence
4698  // the bulk composition
4699  for (unsigned i = 0; i < 4; ++i)
4700  EXPECT_EQ(egs.getBulkMolesOld()[i], bm[i]);
4701  EXPECT_EQ(egs.getBasisActivity(2), 3.3);
4702 
4703  // now set the constraint on HCO3- to bulk, which allows charge-neutrality to be enforced simply
4704  egs.changeConstraintToBulk(3);
4705  // set HCO3- and see it changes H+
4706  bm = egs.getBulkMolesOld();
4707  egs.setConstraintValue(3, 7.7);
4708  bm[1] = 7.7;
4709  bm[3] = 7.7;
4710  for (unsigned i = 0; i < 4; ++i)
4711  EXPECT_NEAR(egs.getBulkMolesOld()[i], bm[i], 1.0E-8);
4712  // set H+ and see it actually makes no difference because of charge-neutrality
4713  egs.setConstraintValue(1, 1.2345);
4714  for (unsigned i = 0; i < 4; ++i)
4715  EXPECT_NEAR(egs.getBulkMolesOld()[i], bm[i], 1.0E-8);
4716 }
4717 
4719 TEST_F(GeochemicalSystemTest, getsetModelGeochemicalDatabase)
4720 {
4721  // this test involves radically changing the basis to illustrate copying, etc, works: no coder
4722  // should actually perform these radical changes otherwise a crash is virtually guaranteed!
4723 
4724  const GeochemicalDatabaseReader db("database/moose_testdb.json", true, true, false);
4725  const PertinentGeochemicalSystem pgs(db, {"H2O", "H+"}, {}, {}, {}, {}, {}, "O2(aq)", "e-");
4726  ModelGeochemicalDatabase mgd = pgs.modelGeochemicalDatabase();
4727  GeochemistrySpeciesSwapper swapper(2, 1E-6);
4728  const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
4731  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
4734  GeochemicalSystem egs(mgd,
4735  _ac3,
4736  _is3,
4737  swapper,
4738  {},
4739  {},
4740  "H+",
4741  {"H2O", "H+"},
4742  {1, 1},
4743  cu,
4744  cm,
4745  25,
4746  0,
4747  1E-20,
4748  {},
4749  {},
4750  {});
4751 
4752  const ModelGeochemicalDatabase & mgd_ref = egs.getModelGeochemicalDatabase();
4753  EXPECT_EQ(mgd_ref.basis_species_name.size(), (std::size_t)2);
4754  ModelGeochemicalDatabase copy_of_mgd = egs.getModelGeochemicalDatabase();
4755  EXPECT_EQ(copy_of_mgd.basis_species_name.size(), (std::size_t)2);
4756 
4757  // the following copies all the data of egs into copy_of_egs, but only the references (such as
4758  // _mgd) get copied: the data they refer to doesn't get copied
4759  GeochemicalSystem copy_of_egs = egs;
4760 
4761  // create a new ModelGeochemicalDatabase, with 3 basis species
4762  const PertinentGeochemicalSystem pgs3(
4763  db, {"H2O", "H+", "Ca++"}, {}, {}, {}, {}, {}, "O2(aq)", "e-");
4764  ModelGeochemicalDatabase mgd3 = pgs3.modelGeochemicalDatabase();
4765 
4766  // Copy mgd3 into the memory referred to by copy_of_egs._mgd, which is also referred to by
4767  // egs._mgd and mgd_ref
4768  copy_of_egs.setModelGeochemicalDatabase(mgd3);
4769  EXPECT_EQ(mgd_ref.basis_species_name.size(), (std::size_t)3);
4770  EXPECT_EQ(egs.getModelGeochemicalDatabase().basis_species_name.size(), (std::size_t)3);
4771  EXPECT_EQ(copy_of_egs.getModelGeochemicalDatabase().basis_species_name.size(), (std::size_t)3);
4772  EXPECT_EQ(copy_of_mgd.basis_species_name.size(),
4773  (std::size_t)2); // copy_of_mgd shouldn't be impacted
4774 
4775  // do a similar thing with 1 basis species
4776  const PertinentGeochemicalSystem pgs1(db, {"H2O"}, {}, {}, {}, {}, {}, "O2(aq)", "e-");
4777  ModelGeochemicalDatabase mgd1 = pgs1.modelGeochemicalDatabase();
4778  egs.setModelGeochemicalDatabase(mgd1);
4779  EXPECT_EQ(mgd_ref.basis_species_name.size(), (std::size_t)1);
4780  EXPECT_EQ(egs.getModelGeochemicalDatabase().basis_species_name.size(), (std::size_t)1);
4781  EXPECT_EQ(copy_of_egs.getModelGeochemicalDatabase().basis_species_name.size(), (std::size_t)1);
4782  EXPECT_EQ(copy_of_mgd.basis_species_name.size(), (std::size_t)2);
4783 }
4784 
4786 TEST_F(GeochemicalSystemTest, getsetMineralRelatedFreeMoles)
4787 {
4788  // first, create an GeochemicalSystem that has minerals and sorption-related things
4789  const PertinentGeochemicalSystem model_gas_and_sorption(
4790  _db_calcite,
4791  {"H2O", "H+", "HCO3-", "O2(aq)", "Fe+++", ">(s)FeOH", ">(w)FeOH", "Ca++"},
4792  {"Goethite", "Calcite"},
4793  {"O2(g)"},
4794  {"Calcite_asdf"},
4795  {},
4796  {},
4797  "O2(aq)",
4798  "e-");
4799  ModelGeochemicalDatabase mgd = model_gas_and_sorption.modelGeochemicalDatabase();
4800  const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
4809  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
4818  std::vector<GeochemistryUnitConverter::GeochemistryUnit> ku;
4820  GeochemicalSystem nonconst(
4821  mgd,
4822  _ac0,
4823  _is0,
4824  _swapper8,
4825  {"O2(aq)", "Ca++"},
4826  {"O2(g)", "Calcite"},
4827  "Fe+++",
4828  {"H2O", "H+", "HCO3-", "O2(g)", "Fe+++", ">(s)FeOH", ">(w)FeOH", "Calcite"},
4829  {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0},
4830  cu,
4831  cm,
4832  25,
4833  0,
4834  1E-20,
4835  {"Calcite_asdf"},
4836  {9.0},
4837  ku);
4838  const unsigned num_basis = nonconst.getNumInBasis();
4839  const unsigned num_eqm = nonconst.getNumInEquilibrium();
4840  const unsigned num_kin = nonconst.getNumKinetic();
4841 
4842  // preset some molalities
4843  const std::vector<std::string> gold_name = {
4844  "H2O", "H+", "HCO3-",
4845  "O2(g)", "Fe+++", ">(s)FeOH",
4846  ">(w)FeOH", "Calcite", "(O-phth)--",
4847  "CH4(aq)", "CO2(aq)", "CO3--",
4848  "CaCO3", "CaOH+", "OH-",
4849  ">(s)FeO-", ">(s)FeOCa+", "Goethite",
4850  "Ca++", "O2(aq)", "Goethite_surface_potential_expr",
4851  "Calcite_asdf"};
4852  const std::vector<Real> set_molal = {1.1, 2.2, 3.3, 0.0, 5.5, 6.6, 7.7, 8.8, 9.9, 10, 11,
4853  12, 13, 14, 15, 17, 18, 19, 20, 21, 22, 23};
4854  nonconst.setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(
4855  gold_name, set_molal, {true, true, true, false, true, true, true, true});
4856 
4857  // now set the moles of mineral-related things to 1234.5
4858  nonconst.setMineralRelatedFreeMoles(1234.5);
4859 
4860  std::vector<Real> gold_molal = set_molal;
4861  gold_molal[5] = 1234.5; // >(s)FeOH
4862  gold_molal[6] = 1234.5; // >(w)FeOH
4863  gold_molal[7] = 1234.5; // Calcite
4864  gold_molal[15] = 1234.5; // >(s)FeO-
4865  gold_molal[16] = 1234.5; // >(s)FeOCa+
4866  gold_molal[17] = 0.0; // Goethite is an equilibrium mineral
4867  gold_molal[21] = 1234.5; // Calcite_asdf
4868 
4869  const std::vector<Real> & basis_mol = nonconst.getSolventMassAndFreeMolalityAndMineralMoles();
4870  for (unsigned i = 0; i < num_basis; ++i)
4871  EXPECT_EQ(basis_mol[mgd.basis_species_index.at(gold_name[i])], gold_molal[i]);
4872  const std::vector<Real> & eqm_mol = nonconst.getEquilibriumMolality();
4873  for (unsigned j = num_basis; j < num_basis + num_eqm; ++j)
4874  EXPECT_EQ(eqm_mol[mgd.eqm_species_index.at(gold_name[j])], gold_molal[j]);
4875  const std::vector<Real> & kin_mol = nonconst.getKineticMoles();
4876  for (unsigned k = num_basis + num_eqm + 1; k < num_basis + num_eqm + num_kin + 1;
4877  ++k) // + 1 for surface_pot
4878  EXPECT_EQ(kin_mol[mgd.kin_species_index.at(gold_name[k])], gold_molal[k]);
4879 }
4880 
4882 TEST_F(GeochemicalSystemTest, initialKinMoleExcept)
4883 {
4884  std::vector<GeochemistryUnitConverter::GeochemistryUnit> ku;
4886  try
4887  {
4888  const GeochemicalSystem bad(_mgd_kinetic_calcite,
4889  _ac3,
4890  _is3,
4891  _swapper4,
4892  {},
4893  {},
4894  "H+",
4895  {"H2O", "Ca++", "H+", "HCO3-"},
4896  {1.75, 3.0, 2.0, 1.0},
4897  _cu_calcite,
4898  _cm_calcite,
4899  25,
4900  0,
4901  1E-20,
4902  {},
4903  {},
4904  {});
4905  FAIL() << "Missing expected exception.";
4906  }
4907  catch (const std::exception & e)
4908  {
4909  std::string msg(e.what());
4910  ASSERT_TRUE(msg.find("Initial mole number (or mass or volume) and a unit must be provided for "
4911  "each kinetic species") != std::string::npos)
4912  << "Failed with unexpected error message: " << msg;
4913  }
4914 
4915  try
4916  {
4917  const GeochemicalSystem bad(_mgd_kinetic_calcite,
4918  _ac3,
4919  _is3,
4920  _swapper4,
4921  {},
4922  {},
4923  "H+",
4924  {"H2O", "Ca++", "H+", "HCO3-"},
4925  {1.75, 3.0, 2.0, 1.0},
4926  _cu_calcite,
4927  _cm_calcite,
4928  25,
4929  0,
4930  1E-20,
4931  {"Calcite"},
4932  {},
4933  {});
4934  FAIL() << "Missing expected exception.";
4935  }
4936  catch (const std::exception & e)
4937  {
4938  std::string msg(e.what());
4939  ASSERT_TRUE(msg.find("Initial mole number (or mass or volume) and a unit must be provided for "
4940  "each kinetic species") != std::string::npos)
4941  << "Failed with unexpected error message: " << msg;
4942  }
4943 
4944  try
4945  {
4946  const GeochemicalSystem bad(_mgd_kinetic_calcite,
4947  _ac3,
4948  _is3,
4949  _swapper4,
4950  {},
4951  {},
4952  "H+",
4953  {"H2O", "Ca++", "H+", "HCO3-"},
4954  {1.75, 3.0, 2.0, 1.0},
4955  _cu_calcite,
4956  _cm_calcite,
4957  25,
4958  0,
4959  1E-20,
4960  {"Calcite"},
4961  {1.0},
4962  {});
4963  FAIL() << "Missing expected exception.";
4964  }
4965  catch (const std::exception & e)
4966  {
4967  std::string msg(e.what());
4968  ASSERT_TRUE(msg.find("Initial mole number (or mass or volume) and a unit must be provided for "
4969  "each kinetic species") != std::string::npos)
4970  << "Failed with unexpected error message: " << msg;
4971  }
4972 
4973  try
4974  {
4975  const GeochemicalSystem bad(_mgd_kinetic_calcite,
4976  _ac3,
4977  _is3,
4978  _swapper4,
4979  {},
4980  {},
4981  "H+",
4982  {"H2O", "Ca++", "H+", "HCO3-"},
4983  {1.75, 3.0, 2.0, 1.0},
4984  _cu_calcite,
4985  _cm_calcite,
4986  25,
4987  0,
4988  1E-20,
4989  {},
4990  {1.0},
4991  ku);
4992  FAIL() << "Missing expected exception.";
4993  }
4994  catch (const std::exception & e)
4995  {
4996  std::string msg(e.what());
4997  ASSERT_TRUE(msg.find("Initial mole number (or mass or volume) and a unit must be provided for "
4998  "each kinetic species") != std::string::npos)
4999  << "Failed with unexpected error message: " << msg;
5000  }
5001 
5002  try
5003  {
5004  const GeochemicalSystem bad(_mgd_kinetic_calcite,
5005  _ac3,
5006  _is3,
5007  _swapper4,
5008  {},
5009  {},
5010  "H+",
5011  {"H2O", "Ca++", "H+", "HCO3-"},
5012  {1.75, 3.0, 2.0, 1.0},
5013  _cu_calcite,
5014  _cm_calcite,
5015  25,
5016  0,
5017  1E-20,
5018  {"Ca++"},
5019  {1.0},
5020  {});
5021  FAIL() << "Missing expected exception.";
5022  }
5023  catch (const std::exception & e)
5024  {
5025  std::string msg(e.what());
5026  ASSERT_TRUE(msg.find("Initial mole number (or mass or volume) and a unit must be provided for "
5027  "each kinetic species") != std::string::npos)
5028  << "Failed with unexpected error message: " << msg;
5029  }
5030 
5031  std::vector<GeochemistryUnitConverter::GeochemistryUnit> ku_bad;
5033  try
5034  {
5035  const GeochemicalSystem bad(_mgd_kinetic_calcite,
5036  _ac3,
5037  _is3,
5038  _swapper4,
5039  {},
5040  {},
5041  "H+",
5042  {"H2O", "Ca++", "H+", "HCO3-"},
5043  {1.75, 3.0, 2.0, 1.0},
5044  _cu_calcite,
5045  _cm_calcite,
5046  25,
5047  0,
5048  1E-20,
5049  {"Calcite"},
5050  {1.0},
5051  ku_bad);
5052  FAIL() << "Missing expected exception.";
5053  }
5054  catch (const std::exception & e)
5055  {
5056  std::string msg(e.what());
5057  ASSERT_TRUE(msg.find("Kinetic species Calcite: units must be moles or mass, or volume in the "
5058  "case of minerals") != std::string::npos)
5059  << "Failed with unexpected error message: " << msg;
5060  }
5061 }
5062 
5064 TEST_F(GeochemicalSystemTest, getKineticLog10Kexcept)
5065 {
5066  try
5067  {
5068  _egs_kinetic_calcite.getKineticLog10K(1);
5069  FAIL() << "Missing expected exception.";
5070  }
5071  catch (const std::exception & e)
5072  {
5073  std::string msg(e.what());
5074  ASSERT_TRUE(msg.find("Cannot retrieve log10K for kinetic species 1 since there are only 1 "
5075  "kinetic species") != std::string::npos)
5076  << "Failed with unexpected error message: " << msg;
5077  }
5078 }
5079 
5081 TEST_F(GeochemicalSystemTest, log10KineticActivityProductExcept)
5082 {
5083  try
5084  {
5085  _egs_kinetic_calcite.log10KineticActivityProduct(1);
5086  FAIL() << "Missing expected exception.";
5087  }
5088  catch (const std::exception & e)
5089  {
5090  std::string msg(e.what());
5091  ASSERT_TRUE(
5092  msg.find("Cannot retrieve activity product for kinetic species 1 since there are only 1 "
5093  "kinetic species") != std::string::npos)
5094  << "Failed with unexpected error message: " << msg;
5095  }
5096 }
5097 
5099 TEST_F(GeochemicalSystemTest, getKineticMolesExcept)
5100 {
5101  try
5102  {
5103  _egs_kinetic_calcite.getKineticMoles(1);
5104  FAIL() << "Missing expected exception.";
5105  }
5106  catch (const std::exception & e)
5107  {
5108  std::string msg(e.what());
5109  ASSERT_TRUE(msg.find("Cannot retrieve moles for kinetic species 1 since there are only 1 "
5110  "kinetic species") != std::string::npos)
5111  << "Failed with unexpected error message: " << msg;
5112  }
5113 }
5114 
5116 TEST_F(GeochemicalSystemTest, setKineticMolesExcept)
5117 {
5118  GeochemicalSystem nonconst = _egs_kinetic_calcite;
5119 
5120  try
5121  {
5122  nonconst.setKineticMoles(1, 1.0);
5123  FAIL() << "Missing expected exception.";
5124  }
5125  catch (const std::exception & e)
5126  {
5127  std::string msg(e.what());
5128  ASSERT_TRUE(msg.find("Cannot set moles for kinetic species 1 since there are only 1 "
5129  "kinetic species") != std::string::npos)
5130  << "Failed with unexpected error message: " << msg;
5131  }
5132 
5133  try
5134  {
5135  nonconst.setKineticMoles(0, 0.0);
5136  FAIL() << "Missing expected exception.";
5137  }
5138  catch (const std::exception & e)
5139  {
5140  std::string msg(e.what());
5141  ASSERT_TRUE(msg.find("Mole number for kinetic species must be positive, not 0") !=
5142  std::string::npos)
5143  << "Failed with unexpected error message: " << msg;
5144  }
5145 }
5146 
5149 {
5150  EXPECT_EQ(_egs_kinetic_calcite.getNumKinetic(), (unsigned)1);
5151 }
5152 
5154 TEST_F(GeochemicalSystemTest, getKineticLog10K)
5155 {
5156  EXPECT_EQ(_egs_kinetic_calcite.getKineticLog10K(0), 1.7130);
5157  EXPECT_EQ(_egs_kinetic_calcite.getKineticLog10K().size(), (std::size_t)1);
5158  EXPECT_EQ(_egs_kinetic_calcite.getKineticLog10K()[0], 1.7130);
5159 }
5160 
5162 TEST_F(GeochemicalSystemTest, log10KineticActivityProduct)
5163 {
5164  const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
5169  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
5174  std::vector<GeochemistryUnitConverter::GeochemistryUnit> ku;
5176  const GeochemicalSystem gs(_mgd_kinetic_calcite,
5177  _ac3,
5178  _is3,
5179  _swapper4,
5180  {},
5181  {},
5182  "Ca++",
5183  {"H2O", "H+", "HCO3-", "Ca++"},
5184  {1.75, 3.0, 2.0, 1.0},
5185  cu,
5186  cm,
5187  25,
5188  0,
5189  1E-20,
5190  {"Calcite"},
5191  {1.1},
5192  ku);
5193  // Calcite = Ca++ + HCO3- - H+
5194  EXPECT_NEAR(
5195  gs.log10KineticActivityProduct(0), std::log10(gs.getBasisActivity(3) / 3.0 * 2.0), 1.0E-6);
5196 }
5197 
5200 {
5201  std::vector<GeochemistryUnitConverter::GeochemistryUnit> ku;
5203  GeochemicalSystem gs(_mgd_kinetic_calcite,
5204  _ac3,
5205  _is3,
5206  _swapper4,
5207  {},
5208  {},
5209  "H+",
5210  {"H2O", "H+", "HCO3-", "Ca++"},
5211  {1.75, 3.0, 2.0, 1.0},
5212  _cu_calcite,
5213  _cm_calcite,
5214  25,
5215  0,
5216  1E-20,
5217  {"Calcite"},
5218  {1.1},
5219  ku);
5220  EXPECT_EQ(gs.getKineticMoles(0), 1.1);
5221  EXPECT_EQ(gs.getKineticMoles().size(), (std::size_t)1);
5222  EXPECT_EQ(gs.getKineticMoles()[0], 1.1);
5223  gs.setKineticMoles(0, 2.2);
5224  EXPECT_EQ(gs.getKineticMoles(0), 2.2);
5225  EXPECT_EQ(gs.getKineticMoles().size(), (std::size_t)1);
5226  EXPECT_EQ(gs.getKineticMoles()[0], 2.2);
5227 
5228  std::vector<GeochemistryUnitConverter::GeochemistryUnit> kucm3;
5230  GeochemicalSystem gscm3(_mgd_kinetic_calcite,
5231  _ac3,
5232  _is3,
5233  _swapper4,
5234  {},
5235  {},
5236  "H+",
5237  {"H2O", "H+", "HCO3-", "Ca++"},
5238  {1.75, 3.0, 2.0, 1.0},
5239  _cu_calcite,
5240  _cm_calcite,
5241  25,
5242  0,
5243  1E-20,
5244  {"Calcite"},
5245  {1.1 * 36.9340},
5246  kucm3);
5247  EXPECT_NEAR(gscm3.getKineticMoles(0), 1.1, 1.0E-6);
5248  EXPECT_EQ(gscm3.getKineticMoles().size(), (std::size_t)1);
5249  EXPECT_NEAR(gscm3.getKineticMoles()[0], 1.1, 1.0E-6);
5250  gscm3.setKineticMoles(0, 2.2);
5251  EXPECT_EQ(gscm3.getKineticMoles(0), 2.2);
5252  EXPECT_EQ(gscm3.getKineticMoles().size(), (std::size_t)1);
5253  EXPECT_EQ(gscm3.getKineticMoles()[0], 2.2);
5254 }
5255 
5257 TEST_F(GeochemicalSystemTest, getBulk_with_kinetic)
5258 {
5259  const PertinentGeochemicalSystem model_without_calcite(
5260  _db_calcite, {"H2O", "H+", "HCO3-", "Ca++"}, {}, {}, {}, {}, {}, "O2(aq)", "e-");
5261  ModelGeochemicalDatabase mgd_without_calcite = model_without_calcite.modelGeochemicalDatabase();
5262  const GeochemicalSystem egs_without_kinetic(
5263  mgd_without_calcite,
5264  _ac3,
5265  _is3,
5266  _swapper4,
5267  {},
5268  {},
5269  "H+",
5270  {"H2O", "Ca++", "H+", "HCO3-"},
5271  // in _egs_kinetic_calcite, the 1.1 moles of Calcite contributes 1.1 moles of Ca++ and -1.1
5272  // moles of H+ (along with 1.1 moles of HCO3-, but this species has a free_concentration
5273  // constraint)
5274  {1.75, 3.0 + 1.1, 2.0 - 1.1, 1.0},
5275  _cu_calcite,
5276  _cm_calcite,
5277  25,
5278  0,
5279  1E-20,
5280  {},
5281  {},
5282  {});
5283  EXPECT_EQ(_egs_kinetic_calcite.getBulkMolesOld()[0],
5284  egs_without_kinetic.getBulkMolesOld()[0]); // H2O
5285  EXPECT_EQ(_egs_kinetic_calcite.getBulkMolesOld()[1],
5286  egs_without_kinetic.getBulkMolesOld()[1]); // H+ has fixed bulk
5287  EXPECT_EQ(_egs_kinetic_calcite.getBulkMolesOld()[2],
5288  egs_without_kinetic.getBulkMolesOld()[2] + 1.1); // HCO3- has fixed free molality
5289  EXPECT_EQ(_egs_kinetic_calcite.getBulkMolesOld()[3],
5290  egs_without_kinetic.getBulkMolesOld()[3]); // Ca++ has fixed bulk
5291 }
5292 
5294 TEST_F(GeochemicalSystemTest, setAlgebraicVars_kinetic)
5295 {
5296  std::vector<GeochemistryUnitConverter::GeochemistryUnit> ku;
5298  GeochemicalSystem gs(_mgd_kinetic_calcite,
5299  _ac3,
5300  _is3,
5301  _swapper4,
5302  {},
5303  {},
5304  "H+",
5305  {"H2O", "H+", "HCO3-", "Ca++"},
5306  {1.75, 3.0, 2.0, 1.0},
5307  _cu_calcite,
5308  _cm_calcite,
5309  25,
5310  0,
5311  1E-20,
5312  {"Calcite"},
5313  {1.1},
5314  ku);
5315  const DenseVector<Real> alg_vars(3, 3.3);
5316  gs.setAlgebraicVariables(alg_vars);
5317  for (unsigned i = 0; i < 3; ++i)
5318  EXPECT_EQ(gs.getAlgebraicVariableValues()[i], 3.3);
5319  EXPECT_EQ(gs.getKineticMoles()[0], 3.3);
5320 }
5321 
5323 TEST_F(GeochemicalSystemTest, free_molality_kinetic)
5324 {
5325  const PertinentGeochemicalSystem model(_db_calcite,
5326  {"H2O", "H+", "HCO3-", "Ca++"},
5327  {"Calcite"},
5328  {},
5329  {"Calcite_asdf"},
5330  {},
5331  {},
5332  "O2(aq)",
5333  "e-");
5335  const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
5340  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
5345  std::vector<GeochemistryUnitConverter::GeochemistryUnit> ku;
5348  _ac0,
5349  _is0,
5350  _swapper4,
5351  {"Ca++"},
5352  {"Calcite"},
5353  "H+",
5354  {"H2O", "H+", "HCO3-", "Calcite"},
5355  {1.75, 3.0, 2.0, 10.0},
5356  cu,
5357  cm,
5358  25,
5359  0,
5360  1E-20,
5361  {"Calcite_asdf"},
5362  {0.25},
5363  ku);
5364  const Real orig_calcite = gs.getSolventMassAndFreeMolalityAndMineralMoles()[3];
5365  gs.setKineticMoles(0, 1.25);
5366  gs.computeConsistentConfiguration();
5367  EXPECT_NEAR(orig_calcite - 2.0,
5368  gs.getSolventMassAndFreeMolalityAndMineralMoles()[3],
5369  1.0E-8); // Calcite_asdf = 2 Calcite + ...
5370 }
5371 
5374 {
5375  const PertinentGeochemicalSystem model(_db_calcite,
5376  {"H2O", "H+", "HCO3-", "Ca++"},
5377  {"Calcite"},
5378  {},
5379  {"Calcite_asdf"},
5380  {},
5381  {},
5382  "O2(aq)",
5383  "e-");
5385  const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
5390  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
5395  std::vector<GeochemistryUnitConverter::GeochemistryUnit> ku;
5398  _ac0,
5399  _is0,
5400  _swapper4,
5401  {"Ca++"},
5402  {"Calcite"},
5403  "H+",
5404  {"H2O", "H+", "HCO3-", "Calcite"},
5405  {1.75, 3.0, 2.0, 10.0},
5406  cu,
5407  cm,
5408  25,
5409  0,
5410  1E-20,
5411  {"Calcite_asdf"},
5412  {0.25},
5413  ku);
5414  const std::vector<std::string> gold_name = {"H2O",
5415  "H+",
5416  "HCO3-",
5417  "Calcite",
5418  "CO2(aq)",
5419  "CO3--",
5420  "CaCO3",
5421  "CaOH+",
5422  "OH-",
5423  "Ca++",
5424  "Calcite_asdf"};
5425  const std::vector<Real> set_molal = {1.1, 2.2, 3.3, 4.4, 5.5, 6.6, 7.7, 8.8, 9.9, 10.0, 11.1};
5427  gold_name, set_molal, {true, true, true, true});
5428 
5429  const unsigned num_basis = gs.getNumInBasis();
5430  const unsigned num_eqm = gs.getNumInEquilibrium();
5431  const unsigned num_kin = gs.getNumKinetic();
5432  const std::vector<Real> molal = gs.getSolventMassAndFreeMolalityAndMineralMoles();
5433  for (unsigned i = 0; i < num_basis; ++i)
5434  EXPECT_EQ(molal[mgd.basis_species_index.at(gold_name[i])], set_molal[i]);
5435  for (unsigned j = num_basis; j < num_basis + num_eqm; ++j)
5436  EXPECT_EQ(gs.getEquilibriumMolality(mgd.eqm_species_index.at(gold_name[j])), set_molal[j]);
5437  for (unsigned k = num_basis + num_eqm; k < num_basis + num_eqm + num_kin; ++k)
5438  EXPECT_EQ(gs.getKineticMoles()[mgd.kin_species_index.at(gold_name[k])], set_molal[k]);
5439 }
5440 
5442 TEST_F(GeochemicalSystemTest, computeJacobianExcept)
5443 {
5444  DenseMatrix<Real> jac;
5445 
5446  try
5447  {
5448  const DenseVector<Real> res(4);
5449  const DenseVector<Real> mole_additions(5);
5450  const DenseMatrix<Real> dmole_additions(5, 5);
5451  _egs_kinetic_calcite.computeJacobian(res, jac, mole_additions, dmole_additions);
5452  FAIL() << "Missing expected exception.";
5453  }
5454  catch (const std::exception & e)
5455  {
5456  std::string msg(e.what());
5457  ASSERT_TRUE(msg.find("Jacobian: residual size must be 3 but it is 4") != std::string::npos)
5458  << "Failed with unexpected error message: " << msg;
5459  }
5460 
5461  try
5462  {
5463  const DenseVector<Real> res(3);
5464  const DenseVector<Real> mole_additions(4);
5465  const DenseMatrix<Real> dmole_additions(5, 5);
5466  _egs_kinetic_calcite.computeJacobian(res, jac, mole_additions, dmole_additions);
5467  FAIL() << "Missing expected exception.";
5468  }
5469  catch (const std::exception & e)
5470  {
5471  std::string msg(e.what());
5472  ASSERT_TRUE(msg.find("Jacobian: the increment in mole numbers (mole_additions) needs to be of "
5473  "size 4 + 1 but it is of size 4") != std::string::npos)
5474  << "Failed with unexpected error message: " << msg;
5475  }
5476 
5477  try
5478  {
5479  const DenseVector<Real> res(3);
5480  const DenseVector<Real> mole_additions(5);
5481  const DenseMatrix<Real> dmole_additions(4, 5);
5482  _egs_kinetic_calcite.computeJacobian(res, jac, mole_additions, dmole_additions);
5483  FAIL() << "Missing expected exception.";
5484  }
5485  catch (const std::exception & e)
5486  {
5487  std::string msg(e.what());
5488  ASSERT_TRUE(msg.find("Jacobian: the derivative of mole additions (dmole_additions) needs to be "
5489  "of size 5x5 but it is of size 4x5") != std::string::npos)
5490  << "Failed with unexpected error message: " << msg;
5491  }
5492 
5493  try
5494  {
5495  const DenseVector<Real> res(3);
5496  const DenseVector<Real> mole_additions(5);
5497  const DenseMatrix<Real> dmole_additions(5, 4);
5498  _egs_kinetic_calcite.computeJacobian(res, jac, mole_additions, dmole_additions);
5499  FAIL() << "Missing expected exception.";
5500  }
5501  catch (const std::exception & e)
5502  {
5503  std::string msg(e.what());
5504  ASSERT_TRUE(msg.find("Jacobian: the derivative of mole additions (dmole_additions) needs to be "
5505  "of size 5x5 but it is of size 5x4") != std::string::npos)
5506  << "Failed with unexpected error message: " << msg;
5507  }
5508 }
5509 
5511 TEST_F(GeochemicalSystemTest, updateOldWithCurrent)
5512 {
5513  try
5514  {
5515  const DenseVector<Real> mole_additions(4);
5516  GeochemicalSystem nonconst = _egs_kinetic_calcite;
5517  nonconst.updateOldWithCurrent(mole_additions);
5518  FAIL() << "Missing expected exception.";
5519  }
5520  catch (const std::exception & e)
5521  {
5522  std::string msg(e.what());
5523  ASSERT_TRUE(msg.find("The increment in mole numbers (mole_additions) needs to be of size 4 + 1 "
5524  "but it is of size 4") != std::string::npos)
5525  << "Failed with unexpected error message: " << msg;
5526  }
5527 }
5528 
5530 TEST_F(GeochemicalSystemTest, addKineticRatesExcept)
5531 {
5532  GeochemicalSystem nonconst = _egs_kinetic_calcite;
5533  DenseVector<Real> mole_additions(5);
5534  DenseMatrix<Real> dmole_additions(5, 5);
5535  try
5536  {
5537  DenseVector<Real> bad(4);
5538  nonconst.addKineticRates(1.0, bad, dmole_additions);
5539  FAIL() << "Missing expected exception.";
5540  }
5541  catch (const std::exception & e)
5542  {
5543  std::string msg(e.what());
5544  ASSERT_TRUE(msg.find("addKineticRates: incorrectly sized additions: 4 5 5") !=
5545  std::string::npos)
5546  << "Failed with unexpected error message: " << msg;
5547  }
5548  try
5549  {
5550  DenseMatrix<Real> bad(4, 5);
5551  nonconst.addKineticRates(1.0, mole_additions, bad);
5552  FAIL() << "Missing expected exception.";
5553  }
5554  catch (const std::exception & e)
5555  {
5556  std::string msg(e.what());
5557  ASSERT_TRUE(msg.find("addKineticRates: incorrectly sized additions: 5 4 5") !=
5558  std::string::npos)
5559  << "Failed with unexpected error message: " << msg;
5560  }
5561  try
5562  {
5563  DenseMatrix<Real> bad(5, 4);
5564  nonconst.addKineticRates(1.0, mole_additions, bad);
5565  FAIL() << "Missing expected exception.";
5566  }
5567  catch (const std::exception & e)
5568  {
5569  std::string msg(e.what());
5570  ASSERT_TRUE(msg.find("addKineticRates: incorrectly sized additions: 5 5 4") !=
5571  std::string::npos)
5572  << "Failed with unexpected error message: " << msg;
5573  }
5574 }
5575 
5578 {
5579  PertinentGeochemicalSystem mod(_db_calcite,
5580  {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++"},
5581  {},
5582  {},
5583  {"Calcite"},
5584  {"CH4(aq)"},
5585  {},
5586  "O2(aq)",
5587  "e-");
5589  1.5,
5590  2.0,
5591  true,
5592  1.5,
5593  0.5,
5594  0.1,
5595  {"H+", "HCO3-"},
5596  {1.3, 1.2},
5597  {0.3, 0.2},
5598  {0.1, 0.2},
5599  0.8,
5600  2.5,
5601  66.0,
5602  0.003,
5604  "H2O",
5605  0.0,
5606  -1.0,
5607  0.0);
5609  7.0,
5610  6.0,
5611  false,
5612  0.0,
5613  0.0,
5614  0.0,
5615  {"H+"},
5616  {-3.0},
5617  {0.0},
5618  {0.0},
5619  2.5,
5620  0.8,
5621  55.0,
5622  0.00315,
5624  "H2O",
5625  0.0,
5626  -1.0,
5627  0.0);
5628  const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
5634  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
5640  std::vector<GeochemistryUnitConverter::GeochemistryUnit> ku;
5643  DenseVector<Real> mole_additions(7);
5644  DenseMatrix<Real> dmole_additions(7, 7);
5645  ModelGeochemicalDatabase mgd_kin = mod.modelGeochemicalDatabase();
5647  _ac3,
5648  _is3,
5649  _swapper5,
5650  {},
5651  {},
5652  "H+",
5653  {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++"},
5654  {1.75, 3.0, 2.0, 1.0, 1.5},
5655  cu,
5656  cm,
5657  25,
5658  0,
5659  1E-20,
5660  {"Calcite", "CH4(aq)"},
5661  {1.1, 2.2},
5662  ku);
5663 
5664  egs.addKineticRates(1.0, mole_additions, dmole_additions);
5665  for (unsigned i = 0; i < 7; ++i)
5666  {
5667  EXPECT_EQ(mole_additions(i), 0.0);
5668  for (unsigned j = 0; j < 7; ++j)
5669  EXPECT_EQ(dmole_additions(i, j), 0.0);
5670  }
5671 
5672  // add a rate and re-construct the geochemical system
5673  mod.addKineticRate(rate_ch4); // this adds to mole_additions(6)
5674  mgd_kin = mod.modelGeochemicalDatabase();
5675  egs.setModelGeochemicalDatabase(mgd_kin);
5676 
5677  egs.addKineticRates(1.0, mole_additions, dmole_additions);
5678  for (unsigned i = 0; i < 6; ++i)
5679  {
5680  EXPECT_EQ(mole_additions(i), 0.0);
5681  for (unsigned j = 0; j < 7; ++j)
5682  EXPECT_EQ(dmole_additions(i, j), 0.0);
5683  }
5684  const Real ch4rate = mole_additions(6);
5685  std::vector<Real> ch4deriv(7);
5686  for (unsigned j = 0; j < 7; ++j)
5687  ch4deriv[j] = dmole_additions(6, j);
5688 
5689  // add another (identical) rate and re-construct the geochemical system
5690  mod.addKineticRate(rate_ch4);
5691  mgd_kin = mod.modelGeochemicalDatabase();
5692  egs.setModelGeochemicalDatabase(mgd_kin);
5693  mole_additions.zero();
5694  dmole_additions.zero();
5695  egs.addKineticRates(1.0, mole_additions, dmole_additions);
5696  for (unsigned i = 0; i < 6; ++i)
5697  {
5698  EXPECT_EQ(mole_additions(i), 0.0);
5699  for (unsigned j = 0; j < 7; ++j)
5700  EXPECT_EQ(dmole_additions(i, j), 0.0);
5701  }
5702  EXPECT_EQ(mole_additions(6), 2 * ch4rate);
5703  for (unsigned j = 0; j < 7; ++j)
5704  EXPECT_EQ(dmole_additions(6, j), 2 * ch4deriv[j]);
5705 
5706  // check timestep size is OK
5707  mole_additions.zero();
5708  dmole_additions.zero();
5709  egs.addKineticRates(0.5, mole_additions, dmole_additions);
5710  for (unsigned i = 0; i < 6; ++i)
5711  {
5712  EXPECT_EQ(mole_additions(i), 0.0);
5713  for (unsigned j = 0; j < 7; ++j)
5714  EXPECT_EQ(dmole_additions(i, j), 0.0);
5715  }
5716  EXPECT_EQ(mole_additions(6), ch4rate);
5717  for (unsigned j = 0; j < 7; ++j)
5718  EXPECT_EQ(dmole_additions(6, j), ch4deriv[j]);
5719 
5720  // add a rate for calcite (mole_additions(5)) and re-construct the geochemical system
5721  mod.addKineticRate(rate_cal);
5722  mgd_kin = mod.modelGeochemicalDatabase();
5723  egs.setModelGeochemicalDatabase(mgd_kin);
5724  mole_additions.zero();
5725  dmole_additions.zero();
5726  egs.addKineticRates(0.5, mole_additions, dmole_additions);
5727  for (unsigned i = 0; i < 5; ++i)
5728  {
5729  EXPECT_EQ(mole_additions(i), 0.0);
5730  for (unsigned j = 0; j < 7; ++j)
5731  EXPECT_EQ(dmole_additions(i, j), 0.0);
5732  }
5733  EXPECT_NE(mole_additions(5), 0.0);
5734  const Real calciterate = mole_additions(5);
5735  EXPECT_EQ(mole_additions(6), ch4rate);
5736  for (unsigned j = 0; j < 7; ++j)
5737  EXPECT_EQ(dmole_additions(6, j), ch4deriv[j]);
5738  std::vector<Real> calcitederiv(7);
5739  Real sum_deriv = 0.0;
5740  for (unsigned j = 0; j < 7; ++j)
5741  {
5742  calcitederiv[j] = dmole_additions(5, j);
5743  sum_deriv += std::abs(dmole_additions(5, j));
5744  }
5745  EXPECT_NE(sum_deriv, 0.0);
5746 
5747  // add more rates for calcite and re-construct the geochemical system
5748  mod.addKineticRate(rate_cal);
5749  mod.addKineticRate(rate_cal);
5750  mod.addKineticRate(rate_cal);
5751  mgd_kin = mod.modelGeochemicalDatabase();
5752  egs.setModelGeochemicalDatabase(mgd_kin);
5753  mole_additions.zero();
5754  dmole_additions.zero();
5755  egs.addKineticRates(1.0, mole_additions, dmole_additions);
5756  for (unsigned i = 0; i < 5; ++i)
5757  {
5758  EXPECT_EQ(mole_additions(i), 0.0);
5759  for (unsigned j = 0; j < 7; ++j)
5760  EXPECT_EQ(dmole_additions(i, j), 0.0);
5761  }
5762  EXPECT_EQ(mole_additions(5), 8 * calciterate);
5763  EXPECT_EQ(mole_additions(6), 2 * ch4rate);
5764  for (unsigned j = 0; j < 7; ++j)
5765  {
5766  EXPECT_EQ(dmole_additions(5, j), 8 * calcitederiv[j]);
5767  EXPECT_EQ(dmole_additions(6, j), 2 * ch4deriv[j]);
5768  }
5769 }
5770 
5772 TEST_F(GeochemicalSystemTest, addKineticRates_bio_eff)
5773 {
5774  PertinentGeochemicalSystem mod(_db_calcite,
5775  {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++"},
5776  {},
5777  {},
5778  {"Calcite"},
5779  {"CH4(aq)"},
5780  {},
5781  "O2(aq)",
5782  "e-");
5784  1.5,
5785  2.0,
5786  true,
5787  1.5,
5788  0.5,
5789  0.1,
5790  {"H+", "HCO3-"},
5791  {1.3, 1.2},
5792  {0.3, 0.2},
5793  {0.1, 0.2},
5794  0.8,
5795  2.5,
5796  66.0,
5797  0.003,
5799  "H2O",
5800  0.0,
5801  -1.0,
5802  1.0);
5803  KineticRateUserDescription rate_ch4_bio("CH4(aq)",
5804  1.5,
5805  2.0,
5806  true,
5807  1.5,
5808  0.5,
5809  0.1,
5810  {"H+", "HCO3-"},
5811  {1.3, 1.2},
5812  {0.3, 0.2},
5813  {0.1, 0.2},
5814  0.8,
5815  2.5,
5816  66.0,
5817  0.003,
5819  "H2O",
5820  0.0,
5821  4.0,
5822  1.0);
5824  1.5,
5825  2.0,
5826  true,
5827  1.5,
5828  0.5,
5829  0.1,
5830  {"H+", "HCO3-"},
5831  {1.3, 1.2},
5832  {0.3, 0.2},
5833  {0.1, 0.2},
5834  0.8,
5835  2.5,
5836  66.0,
5837  0.003,
5839  "H2O",
5840  0.0,
5841  7.0,
5842  1.0);
5843  const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
5849  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
5855  std::vector<GeochemistryUnitConverter::GeochemistryUnit> ku;
5858  DenseVector<Real> mole_additions(7);
5859  DenseMatrix<Real> dmole_additions(7, 7);
5860  ModelGeochemicalDatabase mgd_kin = mod.modelGeochemicalDatabase();
5862  _ac3,
5863  _is3,
5864  _swapper5,
5865  {},
5866  {},
5867  "H+",
5868  {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++"},
5869  {1.75, 3.0, 2.0, 1.0, 1.5},
5870  cu,
5871  cm,
5872  25,
5873  0,
5874  1E-20,
5875  {"Calcite", "CH4(aq)"},
5876  {1.1, 2.2},
5877  ku);
5878 
5879  egs.addKineticRates(1.0, mole_additions, dmole_additions);
5880  for (unsigned i = 0; i < 7; ++i)
5881  {
5882  EXPECT_EQ(mole_additions(i), 0.0);
5883  for (unsigned j = 0; j < 7; ++j)
5884  EXPECT_EQ(dmole_additions(i, j), 0.0);
5885  }
5886 
5887  // add a rate and re-construct the geochemical system
5888  mod.addKineticRate(rate_ch4); // this adds to mole_additions(6) (the CH4(aq))
5889  mgd_kin = mod.modelGeochemicalDatabase();
5890  egs.setModelGeochemicalDatabase(mgd_kin);
5891 
5892  // capture the rate of CH4 with kinetic_bio_efficiency = -1. This is assumed to be correct
5893  egs.addKineticRates(2.0, mole_additions, dmole_additions);
5894  for (unsigned i = 0; i < 6; ++i)
5895  {
5896  EXPECT_EQ(mole_additions(i), 0.0);
5897  for (unsigned j = 0; j < 7; ++j)
5898  EXPECT_EQ(dmole_additions(i, j), 0.0);
5899  }
5900  const Real rate_dt = -mole_additions(6);
5901  const Real drate_dkin_dt = -dmole_additions(6, 6);
5902  std::vector<Real> drate_dmol_dt(6);
5903  for (unsigned j = 0; j < 6; ++j)
5904  drate_dmol_dt[j] = -dmole_additions(6, j);
5905 
5906  // add a rate that is identical except it has kinetic_bio_efficiency = 4, and re-construct the
5907  // geochemical system
5908  mod.addKineticRate(rate_ch4_bio);
5909  mgd_kin = mod.modelGeochemicalDatabase();
5910  egs.setModelGeochemicalDatabase(mgd_kin);
5911  mole_additions.zero();
5912  dmole_additions.zero();
5913  egs.addKineticRates(2.0, mole_additions, dmole_additions);
5914  // CH4(aq) = 1*H2O + 1*H+ + 1*HCO3- + -2*O2(aq)
5915  const std::vector<Real> stoi = {1, 1, 1, -2, 0, 0};
5916  for (unsigned i = 0; i < 6; ++i)
5917  {
5918  EXPECT_NEAR(mole_additions(i), stoi[i] * (4 + 1) * rate_dt, 1E-8);
5919  EXPECT_NEAR(dmole_additions(i, 6), stoi[i] * (4 + 1) * drate_dkin_dt, 1E-8);
5920  for (unsigned j = 0; j < 6; ++j)
5921  EXPECT_NEAR(dmole_additions(i, j), stoi[i] * (4 + 1) * drate_dmol_dt[j], 1E-8);
5922  }
5923  EXPECT_NEAR(mole_additions(6), (4 - 1) * rate_dt, 1E-8);
5924  EXPECT_NEAR(dmole_additions(6, 6), (4 - 1) * drate_dkin_dt, 1E-8);
5925  for (unsigned j = 0; j < 6; ++j)
5926  EXPECT_NEAR(dmole_additions(6, j), (4 - 1) * drate_dmol_dt[j], 1E-8);
5927 
5928  // add a rate that is identical except it has kinetic_bio_efficiency = 7 and direction = DEATH,
5929  // and re-construct the geochemical system
5930  mod.addKineticRate(rate_ch4_death);
5931  mgd_kin = mod.modelGeochemicalDatabase();
5932  egs.setModelGeochemicalDatabase(mgd_kin);
5933  mole_additions.zero();
5934  dmole_additions.zero();
5935  egs.addKineticRates(2.0, mole_additions, dmole_additions);
5936  for (unsigned i = 0; i < 6; ++i)
5937  {
5938  EXPECT_NEAR(mole_additions(i), stoi[i] * (7 + 4 + 1) * rate_dt, 1E-8);
5939  EXPECT_NEAR(dmole_additions(i, 6), stoi[i] * (7 + 4 + 1) * drate_dkin_dt, 1E-8);
5940  for (unsigned j = 0; j < 6; ++j)
5941  EXPECT_NEAR(dmole_additions(i, j), stoi[i] * (7 + 4 + 1) * drate_dmol_dt[j], 1E-8);
5942  }
5943  EXPECT_NEAR(mole_additions(6), (7 + 4 - 1) * rate_dt, 1E-8);
5944  EXPECT_NEAR(dmole_additions(6, 6), (7 + 4 - 1) * drate_dkin_dt, 1E-8);
5945  for (unsigned j = 0; j < 6; ++j)
5946  EXPECT_NEAR(dmole_additions(6, j), (7 + 4 - 1) * drate_dmol_dt[j], 1E-8);
5947 }
5948 
5950 TEST_F(GeochemicalSystemTest, addKineticRates_progeny)
5951 {
5952  PertinentGeochemicalSystem mod(_db_calcite,
5953  {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++"},
5954  {},
5955  {},
5956  {"Calcite"},
5957  {"CH4(aq)"},
5958  {},
5959  "O2(aq)",
5960  "e-");
5962  1.5,
5963  2.0,
5964  true,
5965  1.5,
5966  0.5,
5967  0.1,
5968  {"H+", "HCO3-"},
5969  {1.3, 1.2},
5970  {0.3, 0.2},
5971  {0.1, 0.2},
5972  0.8,
5973  2.5,
5974  66.0,
5975  0.003,
5977  "H2O",
5978  0.0,
5979  -1.0,
5980  0.0);
5981  KineticRateUserDescription rate_ch4_progH("CH4(aq)",
5982  1.5,
5983  2.0,
5984  true,
5985  1.5,
5986  0.5,
5987  0.1,
5988  {"H+", "HCO3-"},
5989  {1.3, 1.2},
5990  {0.3, 0.2},
5991  {0.1, 0.2},
5992  0.8,
5993  2.5,
5994  66.0,
5995  0.003,
5997  "H+",
5998  2.0,
5999  -1.0,
6000  1.0);
6001  KineticRateUserDescription rate_ch4_progOH("CH4(aq)",
6002  1.5,
6003  2.0,
6004  true,
6005  1.5,
6006  0.5,
6007  0.1,
6008  {"H+", "HCO3-"},
6009  {1.3, 1.2},
6010  {0.3, 0.2},
6011  {0.1, 0.2},
6012  0.8,
6013  2.5,
6014  66.0,
6015  0.003,
6017  "OH-",
6018  22.0,
6019  -1.0,
6020  1.0);
6021  const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
6027  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
6033  std::vector<GeochemistryUnitConverter::GeochemistryUnit> ku;
6036  DenseVector<Real> mole_additions(7);
6037  DenseMatrix<Real> dmole_additions(7, 7);
6038  ModelGeochemicalDatabase mgd_kin = mod.modelGeochemicalDatabase();
6040  _ac3,
6041  _is3,
6042  _swapper5,
6043  {},
6044  {},
6045  "H+",
6046  {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++"},
6047  {1.75, 3.0, 2.0, 1.0, 1.5},
6048  cu,
6049  cm,
6050  25,
6051  0,
6052  1E-20,
6053  {"Calcite", "CH4(aq)"},
6054  {1.1, 2.2},
6055  ku);
6056 
6057  egs.addKineticRates(1.0, mole_additions, dmole_additions);
6058  for (unsigned i = 0; i < 7; ++i)
6059  {
6060  EXPECT_EQ(mole_additions(i), 0.0);
6061  for (unsigned j = 0; j < 7; ++j)
6062  EXPECT_EQ(dmole_additions(i, j), 0.0);
6063  }
6064 
6065  // add a rate and re-construct the geochemical system
6066  mod.addKineticRate(rate_ch4); // this adds to mole_additions(6) (the CH4(aq))
6067  mgd_kin = mod.modelGeochemicalDatabase();
6068  egs.setModelGeochemicalDatabase(mgd_kin);
6069 
6070  // capture the rate of CH4 and the derivatives. These are assumed to be correct
6071  egs.addKineticRates(2.0, mole_additions, dmole_additions);
6072  for (unsigned i = 0; i < 6; ++i)
6073  {
6074  EXPECT_EQ(mole_additions(i), 0.0);
6075  for (unsigned j = 0; j < 7; ++j)
6076  EXPECT_EQ(dmole_additions(i, j), 0.0);
6077  }
6078  const Real rate_dt = -mole_additions(6);
6079  const Real drate_dkin_dt = -dmole_additions(6, 6);
6080  std::vector<Real> drate_dmol_dt(6);
6081  for (unsigned j = 0; j < 6; ++j)
6082  drate_dmol_dt[j] = -dmole_additions(6, j);
6083 
6084  // add a rate and re-construct the geochemical system
6085  mod.addKineticRate(rate_ch4_progH); // this adds 2 moles of H as 1 mole of CH4(aq) dissolves
6086  mgd_kin = mod.modelGeochemicalDatabase();
6087  egs.setModelGeochemicalDatabase(mgd_kin);
6088  mole_additions.zero();
6089  dmole_additions.zero();
6090  egs.addKineticRates(2.0, mole_additions, dmole_additions);
6091  EXPECT_NEAR(mole_additions(1), 2 * rate_dt, 1E-8);
6092  for (unsigned j = 0; j < 6; ++j)
6093  EXPECT_NEAR(dmole_additions(1, j), 2 * drate_dmol_dt[j], 1E-8);
6094  EXPECT_NEAR(dmole_additions(1, 6), 2 * drate_dkin_dt, 1E-8);
6095  for (unsigned i = 0; i < 6; ++i)
6096  {
6097  if (i == 1)
6098  continue;
6099  EXPECT_EQ(mole_additions(i), 0);
6100  for (unsigned j = 0; j < 6; ++j)
6101  EXPECT_EQ(dmole_additions(i, j), 0);
6102  }
6103  EXPECT_NEAR(mole_additions(6), -2 * rate_dt, 1E-8);
6104  EXPECT_NEAR(dmole_additions(6, 6), -2 * drate_dkin_dt, 1E-8);
6105  for (unsigned j = 0; j < 6; ++j)
6106  EXPECT_NEAR(dmole_additions(6, j), -2 * drate_dmol_dt[j], 1E-8);
6107 
6108  // add a rate and re-construct the geochemical system
6109  mod.addKineticRate(rate_ch4_progOH); // this adds 22 moles of OH-, which is 22 moles of H2O - 22
6110  // moles of H+, as 1 mole of CH4(aq) dissolves
6111  mgd_kin = mod.modelGeochemicalDatabase();
6112  egs.setModelGeochemicalDatabase(mgd_kin);
6113  mole_additions.zero();
6114  dmole_additions.zero();
6115  egs.addKineticRates(2.0, mole_additions, dmole_additions);
6116  EXPECT_NEAR(mole_additions(0), 22 * rate_dt, 1E-8);
6117  for (unsigned j = 0; j < 6; ++j)
6118  EXPECT_NEAR(dmole_additions(0, j), 22 * drate_dmol_dt[j], 1E-8);
6119  EXPECT_NEAR(dmole_additions(0, 6), 22 * drate_dkin_dt, 1E-8);
6120  EXPECT_NEAR(mole_additions(1), (-22 + 2) * rate_dt, 1E-8);
6121  for (unsigned j = 0; j < 6; ++j)
6122  EXPECT_NEAR(dmole_additions(1, j), (-22 + 2) * drate_dmol_dt[j], 1E-8);
6123  EXPECT_NEAR(dmole_additions(1, 6), (-22 + 2) * drate_dkin_dt, 1E-8);
6124  for (unsigned i = 2; i < 6; ++i)
6125  {
6126  EXPECT_EQ(mole_additions(i), 0);
6127  for (unsigned j = 0; j < 6; ++j)
6128  EXPECT_EQ(dmole_additions(i, j), 0);
6129  }
6130  EXPECT_NEAR(mole_additions(6), -3 * rate_dt, 1E-8);
6131  EXPECT_NEAR(dmole_additions(6, 6), -3 * drate_dkin_dt, 1E-8);
6132  for (unsigned j = 0; j < 6; ++j)
6133  EXPECT_NEAR(dmole_additions(6, j), -3 * drate_dmol_dt[j], 1E-8);
6134 }
6135 
6137 TEST_F(GeochemicalSystemTest, getBulkOldInOriginalBasis)
6138 {
6139  const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
6144  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
6149  std::vector<GeochemistryUnitConverter::GeochemistryUnit> ku;
6151 
6152  GeochemicalSystem egs_no_swaps(_mgd_kinetic_calcite,
6153  _ac3,
6154  _is3,
6155  _swapper4,
6156  {},
6157  {},
6158  "H+",
6159  {"H2O", "H+", "HCO3-", "Ca++"},
6160  {1.75, 1.0, 5.0, 2.0},
6161  cu,
6162  cm,
6163  25,
6164  0,
6165  1E-20,
6166  {"Calcite"},
6167  {1.1},
6168  ku);
6169  EXPECT_EQ(egs_no_swaps.getBulkOldInOriginalBasis()(0), 1.75); // H2O
6170  EXPECT_NEAR(egs_no_swaps.getBulkOldInOriginalBasis()(1), 2.1, 1E-6); // H+ from charge neutrality
6171  EXPECT_EQ(egs_no_swaps.getBulkOldInOriginalBasis()(2), 5.0 + 1.1); // HCO3-
6172  EXPECT_EQ(egs_no_swaps.getBulkOldInOriginalBasis()(3), 2.0); // Ca++
6173 
6174  ModelGeochemicalDatabase my_mgd_calcite = _mgd_calcite;
6175  GeochemicalSystem egs(my_mgd_calcite,
6176  _ac3,
6177  _is3,
6178  _swapper4,
6179  {},
6180  {},
6181  "H+",
6182  {"H2O", "H+", "HCO3-", "Calcite"},
6183  {-0.5, 2.5, 1.0, 3.5},
6184  cu,
6185  cm,
6186  25,
6187  0,
6188  1E-20,
6189  {},
6190  {},
6191  {});
6192  // In current basis, bulk(H2O) = -0.5, bulk(H+) = 1.0 (charge neutrality), bulk(HCO3-) = 1.0,
6193  // bulk(Calcite) = 3.5. So, in orig basis (Calcite = Ca++ + HCO3- - H+):
6194  EXPECT_EQ(egs.getBulkOldInOriginalBasis()(0), -0.5); // H2O
6195  EXPECT_EQ(egs.getBulkOldInOriginalBasis()(1), -2.5); // H+
6196  EXPECT_EQ(egs.getBulkOldInOriginalBasis()(2), 4.5); // HCO3-
6197  EXPECT_EQ(egs.getBulkOldInOriginalBasis()(3), 3.5); // Ca++
6198 
6199  // perform some swaps
6200  egs.performSwap(2, 0); // CO2(aq) and HCO3-
6201  egs.performSwap(3, 3); // Calcite and CaOH+
6202 
6203  // bulk in original basis should not have changed (apart from some roundoff error)
6204  EXPECT_NEAR(egs.getBulkOldInOriginalBasis()(0), -0.5, 1.0E-6); // H2O
6205  EXPECT_NEAR(egs.getBulkOldInOriginalBasis()(1), -2.5, 1.0E-6); // H+
6206  EXPECT_NEAR(egs.getBulkOldInOriginalBasis()(2), 4.5, 1.0E-6); // HCO3-
6207  EXPECT_NEAR(egs.getBulkOldInOriginalBasis()(3), 3.5, 1.0E-6); // Ca++
6208 }
6209 
6211 TEST_F(GeochemicalSystemTest, getTransportedBulkMoles)
6212 {
6213  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
6214  const std::vector<std::string> original_basis = {
6215  "H2O", "H+", ">(s)FeOH", ">(w)FeOH", "Fe+++", "HCO3-"};
6217  original_basis,
6218  {"Fe(OH)3(ppd)"},
6219  {},
6220  {"Something", "Fe(OH)3(ppd)fake"}, // non transporting
6221  {},
6222  {},
6223  "O2(aq)",
6224  "e-");
6226  const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
6233  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
6240  std::vector<GeochemistryUnitConverter::GeochemistryUnit> ku;
6243  GeochemicalSystem egs(mgd,
6244  _ac3,
6245  _is3,
6246  _swapper6,
6247  {},
6248  {},
6249  "H+",
6250  {"H2O", "H+", ">(s)FeOH", ">(w)FeOH", "Fe+++", "HCO3-"},
6251  {1.75, 1.0, 2.0, 3.0, 0.5, 1.0},
6252  cu,
6253  cm,
6254  40,
6255  0,
6256  1E-20,
6257  {"Something", "Fe(OH)3(ppd)fake"},
6258  {4.4, 5.5},
6259  ku);
6260 
6261  std::vector<Real> basis_mol = egs.getSolventMassAndFreeMolalityAndMineralMoles();
6262  // equilibrium species are (in this order):
6263  // [0] CO2(aq) = -H2O + H+ + HCO3-
6264  // [1] CO3-- = HCO3- - H+
6265  // [2] OH- = H2O - H+
6266  // [3] >(s)FeO- = >(s)FeOH - H+. This is not transporting
6267  // [4] Fe(OH)3(ppd) = 3H2O - 3H+ + Fe+++. This is not transporting
6268  std::vector<Real> eqm_mol = egs.getEquilibriumMolality();
6269 
6270  std::vector<Real> original_trans_bulk;
6271  egs.computeTransportedBulkFromMolalities(original_trans_bulk);
6272 
6273  std::vector<Real> gold_trans_bulk(6);
6274  // H2O
6275  gold_trans_bulk[0] = 1.75 * (GeochemistryConstants::MOLES_PER_KG_WATER - eqm_mol[0] + eqm_mol[2]);
6276  // "H+"
6277  gold_trans_bulk[1] = 1.75 * (basis_mol[1] + eqm_mol[0] - eqm_mol[1] - eqm_mol[2]);
6278  // ">(s)FeOH". Non transporting and all equilibrium species involving it aren't transporting
6279  gold_trans_bulk[2] = 0.0;
6280  // ">(w)FeOH". Non transporting, and there are no equilibrium species involving it
6281  gold_trans_bulk[3] = 0.0;
6282  // "Fe+++". Equilibrium species Fe(OH)3(ppd) is not transporting
6283  gold_trans_bulk[4] = 1.75 * basis_mol[4];
6284  // "HCO3-"
6285  gold_trans_bulk[5] = 1.75 * (basis_mol[5] + eqm_mol[0] + eqm_mol[1]);
6286 
6287  const Real eps = 1E-6;
6288 
6289  for (unsigned i = 0; i < 6; ++i)
6290  {
6291  if (std::abs(gold_trans_bulk[i]) <= eps)
6292  EXPECT_LE(std::abs(original_trans_bulk[i]), eps);
6293  else
6294  EXPECT_NEAR(original_trans_bulk[i] / gold_trans_bulk[i], 1.0, eps);
6295  }
6296 
6297  EXPECT_EQ(model.originalBasisNames(), original_basis);
6298 
6299  // swap out Fe+++ and replace with Fe(OH)3(ppd)
6300  egs.performSwap(4, 4);
6301 
6302  // compare the gold and computed transporting bulk
6303  basis_mol = egs.getSolventMassAndFreeMolalityAndMineralMoles();
6304  // equilibrium species are (in this order):
6305  // [0] CO2(aq) = -H2O + H+ + HCO3-
6306  // [1] CO3-- = HCO3- - H+
6307  // [2] OH- = H2O - H+
6308  // [3] >(s)FeO- = >(s)FeOH - H+. This is not transporting
6309  // [4] Fe+++ = Fe(OH)3(ppd) - 3H2O + 3H+
6310  eqm_mol = egs.getEquilibriumMolality();
6311  gold_trans_bulk[0] =
6312  1.75 * (GeochemistryConstants::MOLES_PER_KG_WATER - eqm_mol[0] + eqm_mol[2] - 3 * eqm_mol[4]);
6313  gold_trans_bulk[1] =
6314  1.75 * (basis_mol[1] + eqm_mol[0] - eqm_mol[1] - eqm_mol[2] + 3 * eqm_mol[4]);
6315  gold_trans_bulk[2] = 0.0;
6316  gold_trans_bulk[3] = 0.0;
6317  // "Fe(OH)3(ppd)". Only contribution is from Fe+++
6318  gold_trans_bulk[4] = 1.75 * eqm_mol[4];
6319  gold_trans_bulk[5] = 1.75 * (basis_mol[5] + eqm_mol[0] + eqm_mol[1]);
6320 
6321  std::vector<Real> new_trans_bulk;
6322  egs.computeTransportedBulkFromMolalities(new_trans_bulk);
6323  for (unsigned i = 0; i < 6; ++i)
6324  {
6325  if (std::abs(gold_trans_bulk[i]) <= eps)
6326  EXPECT_LE(std::abs(new_trans_bulk[i]), eps);
6327  else
6328  EXPECT_NEAR(new_trans_bulk[i] / gold_trans_bulk[i], 1.0, eps);
6329  }
6330 
6331  // now in the original basis:
6332  std::vector<Real> gold_trans_bulk_orig(6);
6333  gold_trans_bulk_orig[0] =
6334  1.75 * (GeochemistryConstants::MOLES_PER_KG_WATER - eqm_mol[0] + eqm_mol[2]);
6335  gold_trans_bulk_orig[1] = 1.75 * (basis_mol[1] + eqm_mol[0] - eqm_mol[1] - eqm_mol[2]);
6336  gold_trans_bulk_orig[2] = 0.0;
6337  gold_trans_bulk_orig[3] = 0.0;
6338  gold_trans_bulk_orig[4] = 1.75 * eqm_mol[4];
6339  gold_trans_bulk_orig[5] = 1.75 * (basis_mol[5] + eqm_mol[0] + eqm_mol[1]);
6340 
6341  DenseVector<Real> trans_bulk_original_basis = egs.getTransportedBulkInOriginalBasis();
6342  for (unsigned i = 0; i < 6; ++i)
6343  {
6344  if (std::abs(gold_trans_bulk_orig[i]) <= eps)
6345  EXPECT_LE(std::abs(trans_bulk_original_basis(i)), eps);
6346  else
6347  EXPECT_NEAR(trans_bulk_original_basis(i) / gold_trans_bulk_orig[i], 1.0, eps);
6348  }
6349 
6350  EXPECT_EQ(model.originalBasisNames(), original_basis);
6351 
6352  // now set mineral-related moles and there should be no change since these are non-transporting
6353  egs.setMineralRelatedFreeMoles(123.0);
6354  egs.computeTransportedBulkFromMolalities(new_trans_bulk);
6355  for (unsigned i = 0; i < 6; ++i)
6356  {
6357  if (std::abs(gold_trans_bulk[i]) <= eps)
6358  EXPECT_LE(std::abs(new_trans_bulk[i]), eps);
6359  else
6360  EXPECT_NEAR(new_trans_bulk[i] / gold_trans_bulk[i], 1.0, eps);
6361  }
6362  trans_bulk_original_basis = egs.getTransportedBulkInOriginalBasis();
6363  for (unsigned i = 0; i < 6; ++i)
6364  {
6365  if (std::abs(gold_trans_bulk_orig[i]) <= eps)
6366  EXPECT_LE(std::abs(trans_bulk_original_basis(i)), eps);
6367  else
6368  EXPECT_NEAR(trans_bulk_original_basis(i) / gold_trans_bulk_orig[i], 1.0, eps);
6369  }
6370 }
6371 
6373 TEST_F(GeochemicalSystemTest, getTransportedBulkMoles_kin_redox)
6374 {
6375  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
6377  database,
6378  {"H2O", "H+", "Fe++", "O2(aq)"},
6379  {},
6380  {},
6381  {},
6382  {"Fe+++"}, // is transporting, but should not be reported in transportedMoles
6383  {},
6384  "O2(aq)",
6385  "e-");
6387  const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
6392  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
6397  std::vector<GeochemistryUnitConverter::GeochemistryUnit> ku;
6399  GeochemicalSystem egs(mgd,
6400  _ac3,
6401  _is3,
6402  _swapper4,
6403  {},
6404  {},
6405  "H+",
6406  {"H2O", "H+", "Fe++", "O2(aq)"},
6407  {1.75, -2.0, 1.0, 1.0},
6408  cu,
6409  cm,
6410  40,
6411  0,
6412  1E-20,
6413  {"Fe+++"},
6414  {4.4},
6415  ku);
6416 
6417  std::vector<Real> basis_mol = egs.getSolventMassAndFreeMolalityAndMineralMoles();
6418  // equilibrium species are only OH- = H2O - H+
6419  std::vector<Real> eqm_mol = egs.getEquilibriumMolality();
6420 
6421  std::vector<Real> trans_bulk;
6422  egs.computeTransportedBulkFromMolalities(trans_bulk);
6423  std::vector<Real> gold_trans_bulk(4);
6424  // H2O
6425  gold_trans_bulk[0] = 1.75 * (GeochemistryConstants::MOLES_PER_KG_WATER + eqm_mol[0]);
6426  // "H+"
6427  gold_trans_bulk[1] = 1.75 * (basis_mol[1] - eqm_mol[0]);
6428  // "Fe++"
6429  gold_trans_bulk[2] = 1.75 * (basis_mol[2]);
6430  // "O2(aq)"
6431  gold_trans_bulk[3] = 1.75 * (basis_mol[3]);
6432 
6433  const Real eps = 1E-6;
6434  for (unsigned i = 0; i < 4; ++i)
6435  EXPECT_NEAR(trans_bulk[i] / gold_trans_bulk[i], 1.0, eps);
6436 }
6437 
6439 TEST_F(GeochemicalSystemTest, copyAssignmentExcept)
6440 {
6441  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
6443  database, {"H2O", "H+", "HCO3-"}, {}, {}, {}, {}, {}, "O2(aq)", "e-");
6445  const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
6449  std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
6453  GeochemicalSystem dest(mgd,
6454  _ac3,
6455  _is3,
6456  _swapper3,
6457  {},
6458  {},
6459  "H+",
6460  {"H2O", "H+", "HCO3-"},
6461  {1.75, 1.0, 2.0},
6462  cu,
6463  cm,
6464  40,
6465  0,
6466  1E-20,
6467  {},
6468  {},
6469  {});
6470 
6471  try
6472  {
6473  dest = _egs_calcite; // wrong number of basis species
6474  FAIL() << "Missing expected exception.";
6475  }
6476  catch (const std::exception & e)
6477  {
6478  std::string msg(e.what());
6479  ASSERT_TRUE(msg.find("GeochemicalSystem: copy assigment operator called with inconsistent "
6480  "fundamental properties") != std::string::npos)
6481  << "Failed with unexpected error message: " << msg;
6482  }
6483 
6484  GeochemicalSystem dest2(mgd,
6485  _ac3,
6486  _is3,
6487  _swapper3,
6488  {},
6489  {},
6490  "HCO3-",
6491  {"H2O", "H+", "HCO3-"},
6492  {1.75, 1.0, 2.0},
6493  cu,
6494  cm,
6495  40,
6496  0,
6497  1E-20,
6498  {},
6499  {},
6500  {});
6501  try
6502  {
6503  dest2 = dest; // incorrect original charge-balance species
6504  FAIL() << "Missing expected exception.";
6505  }
6506  catch (const std::exception & e)
6507  {
6508  std::string msg(e.what());
6509  ASSERT_TRUE(msg.find("GeochemicalSystem: copy assigment operator called with inconsistent "
6510  "fundamental properties") != std::string::npos)
6511  << "Failed with unexpected error message: " << msg;
6512  }
6513 
6514  const PertinentGeochemicalSystem model_no_kinetic_calcite(
6515  _db_calcite, {"H2O", "H+", "HCO3-", "Ca++"}, {}, {}, {}, {}, {}, "O2(aq)", "e-");
6516  ModelGeochemicalDatabase mgd_no_kinetic_calcite =
6517  model_no_kinetic_calcite.modelGeochemicalDatabase();
6518  GeochemicalSystem egs_no_kinetic_calcite(mgd_no_kinetic_calcite,
6519  _ac3,
6520  _is3,
6521  _swapper4,
6522  {},
6523  {},
6524  "H+",
6525  {"H2O", "Ca++", "H+", "HCO3-"},
6526  {1.75, 3.0, 2.0, 1.0},
6527  _cu_calcite,
6528  _cm_calcite,
6529  25,
6530  0,
6531  1E-20,
6532  {},
6533  {},
6534  {});
6535  try
6536  {
6537  egs_no_kinetic_calcite = _egs_kinetic_calcite; // wrong number of kinetic species
6538  FAIL() << "Missing expected exception.";
6539  }
6540  catch (const std::exception & e)
6541  {
6542  std::string msg(e.what());
6543  ASSERT_TRUE(msg.find("GeochemicalSystem: copy assigment operator called with inconsistent "
6544  "fundamental properties") != std::string::npos)
6545  << "Failed with unexpected error message: " << msg;
6546  }
6547 }
Real getTemperature() const
constexpr Real MOLES_PER_KG_WATER
void changeConstraintToBulk(unsigned basis_ind)
Changes the constraint to MOLES_BULK_SPECIES (or MOLES_BULK_WATER if basis_ind = 0) for the basis ind...
const std::vector< Real > & getBulkMolesOld() const
constexpr Real CELSIUS_TO_KELVIN
const Real eps
virtual void zero() override final
void setAdditionalValue(const std::string &names)
TEST_F(GeochemicalSystemTest, constructWithMultiMooseEnum)
Check MultiMooseEnum constructor.
std::unordered_map< std::string, unsigned > basis_species_index
basis_species_index[name] = index of the basis species, within all ModelGeochemicalDatabase internal ...
const ModelGeochemicalDatabase mgd
std::unordered_map< std::string, unsigned > eqm_species_index
eqm_species_index[name] = index of the equilibrium species (secondary aqueous species, redox couples in equilibrium with the aqueous solution, minerals in equilibrium with the aqueous solution, gases in equilibrium with the aqueous solution) within all ModelGeochemicalDatabase internal datastrcutres, with given name
void updateOldWithCurrent(const DenseVector< Real > &mole_additions)
Usually used at the end of a solve, to provide correct initial conditions to the next time-step...
const GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false)
const std::vector< Real > & computeAndGetEquilibriumActivity()
Computes activity for all equilibrium species (_eqm_activity) and returns a reference to the vector...
unsigned getNumInBasis() const
returns the number of species in the basis
std::unordered_map< std::string, unsigned > kin_species_index
kin_species_index[name] = index of the kinetic species, within all ModelGeochemicalDatabase internal ...
DenseMatrix< Real > eqm_stoichiometry
eqm_stoichiometry(i, j) = stoichiometric coefficient for equilibrium species "i" in terms of the basi...
KineticRateUserDescription rate_ch4("CH4(aq)", 1.5, 2.0, true, 0.0, 0.0, 0.0, {"H2O", "OH-", "O2(aq)", "CO2(aq)", "CaCO3"}, {3.0, 3.1, 3.2, 3.3, 3.4}, {0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0}, 0.8, 2.5, 66.0, 0.003, DirectionChoiceEnum::BOTH, "H2O", 0.0, -1.0, 0.0)
Real getTotalChargeOld() const
void setKineticMoles(unsigned kin, Real moles)
Sets the current AND old mole number for a kinetic species.
const GeochemicalDatabaseReader db("database/moose_testdb.json", true, true, false)
Real getBasisActivity(unsigned i) const
const ModelGeochemicalDatabase & mgd_kin
Class to swap basis species with equilibrium species.
Computes activity coefficients for non-minerals and non-gases (since these species do not have activi...
Real getEquilibriumActivity(unsigned eqm_ind) const
Returns the value of activity for the equilibrium species with index eqm_index.
const ModelGeochemicalDatabase & modelGeochemicalDatabase() const
Return a reference to the ModelGeochemicalDatabase structure.
Constructs and stores a minimal amount of information that is pertinent to the user-defined geochemic...
Calculators to compute ionic strength and stoichiometric ionic strength.
const ModelGeochemicalDatabase & getModelGeochemicalDatabase() const
void addKineticRates(Real dt, DenseVector< Real > &mole_additions, DenseMatrix< Real > &dmole_additions)
Computes the kinetic rates and their derivatives based on the current values of molality, mole number, etc, and then, using dt, calculates the appropriate mole_additions and dmole_additions.
PetscErrorCode PetscInt const PetscInt IS * is
void enforceChargeBalance()
Enforces charge balance by altering the constraint_value and bulk_moles_old of the charge-balance spe...
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-")
KineticRateUserDescription rate_cal("Calcite", 7.0, 6.0, false, 0.0, 0.0, 0.0, {"H+"}, {-3.0}, {0.0}, {0.0}, 2.5, 0.8, 55.0, 0.00315, DirectionChoiceEnum::BOTH, "H2O", 0.0, -1.0, 0.0)
void setModelGeochemicalDatabase(const ModelGeochemicalDatabase &mgd)
Copies a ModelGeochemicalDatabase into our _mgd structure.
Real getLog10K(unsigned j) const
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const std::vector< GeochemistryUnitConverter::GeochemistryUnit > cu2
std::vector< std::string > originalBasisNames() const
Real getEquilibriumMolality(unsigned j) const
void setTemperature(Real temperature)
Sets the temperature to the given quantity.
This class holds information about bulk composition, molalities, activities, activity coefficients...
Holds a user-specified description of a kinetic rate.
void setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(const std::vector< std::string > &names, const std::vector< Real > &values, const std::vector< bool > &constraints_from_molalities)
Sets:
Real getEquilibriumActivityCoefficient(unsigned j) const
Data structure to hold all relevant information from the database file.
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
std::vector< bool > eqm_species_mineral
eqm_species_mineral[i] = true iff the i^th equilibrium species is a mineral
const std::vector< Real > & getSolventMassAndFreeMolalityAndMineralMoles() const
void addToBulkMoles(unsigned basis_ind, Real value)
Add to the MOLES_BULK_SPECIES (or MOLES_BULK_WATER if basis_ind = 0) for the basis species...
KineticRateUserDescription rate_ch4_death("CH4(aq)", 1.5, 2.0, true, 0.0, 0.0, 0.0, {"H2O", "OH-", "O2(aq)", "CO2(aq)", "CaCO3"}, {3.0, 3.1, 3.2, 3.3, 3.4}, {0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0}, 0.8, 2.5, 66.0, 0.003, DirectionChoiceEnum::DEATH, "H2O", 0.0, -1.0, 0.0)
void performSwap(unsigned swap_out_of_basis, unsigned swap_into_basis)
Perform the basis swap, and ensure that the resulting system is consistent.
void setConstraintValue(unsigned basis_ind, Real value)
Set the constraint value for the basis species.
unsigned getNumInEquilibrium() const
returns the number of species in equilibrium with the basis components
void closeSystem()
Changes a KG_SOLVENT_WATER constraint to MOLES_BULK_WATER (if there is such a constraint) and all FRE...
const std::vector< GeochemicalSystem::ConstraintUserMeaningEnum > cm2
MooseUnits pow(const MooseUnits &, int)
std::vector< std::string > eqm_species_name
eqm_species_name[i] = name of the i^th eqm species
static const std::string k
Definition: NS.h:130
std::vector< std::string > basis_species_name
basis_species_name[j] = name of the j^th basis species
Class for reading geochemical reactions from a MOOSE geochemical database.
std::vector< Real > getSaturationIndices() const