https://mooseframework.inl.gov
Loading...
Searching...
No Matches
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
11
13TEST_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
44TEST_F(GeochemicalSystemTest, numInBasis) { EXPECT_EQ(_egs_calcite.getNumInBasis(), (unsigned)4); }
45
47TEST_F(GeochemicalSystemTest, numInEquilibrium)
48{
49 EXPECT_EQ(_egs_calcite.getNumInEquilibrium(), (unsigned)6);
50}
51
53TEST_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
98TEST_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
187TEST_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}
330TEST_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
406TEST_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}
452TEST_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
499TEST_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
548TEST_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
595TEST_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
691TEST_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
739TEST_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
827TEST_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
911TEST_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
957TEST_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
1004TEST_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 {
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
1051TEST_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 {
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
1102TEST_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 {
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 {
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
1181TEST_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 {
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
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
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
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
1358TEST_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
1450TEST_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
1484TEST_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
1518TEST_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
1542TEST_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
1551TEST_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
1562TEST_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
1570TEST_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
1578TEST_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
1652TEST_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
1717TEST_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
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
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
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
1826TEST_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
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
1883TEST_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
1998TEST_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
2029TEST_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
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
2222TEST_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;
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;
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
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
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;
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();
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
2945TEST_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;
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
3000TEST_F(GeochemicalSystemTest, getNumRedox) { EXPECT_EQ(_egs_redox.getNumRedox(), (unsigned)3); }
3001
3003TEST_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
3040TEST_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
3073TEST_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
3090TEST_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
3107TEST_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;
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;
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);
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);
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);
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
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
3576TEST_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
3749TEST_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
3984TEST_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
4270TEST_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
4310TEST_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
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 {
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
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
4435TEST_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
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
4508TEST_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
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
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
4613TEST_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
4633TEST_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
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
4719TEST_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;
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
4786TEST_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
4882TEST_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
5064TEST_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
5081TEST_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
5099TEST_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
5116TEST_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
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
5162TEST_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
5257TEST_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
5294TEST_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
5323TEST_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
5442TEST_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
5511TEST_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
5530TEST_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
5772TEST_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
5950TEST_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
6137TEST_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
6211TEST_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;
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
6373TEST_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;
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
6439TEST_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;
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}
const std::vector< GeochemicalSystem::ConstraintUserMeaningEnum > cm2
const std::vector< GeochemistryUnitConverter::GeochemistryUnit > cu2
TEST_F(GeochemicalSystemTest, constructWithMultiMooseEnum)
Check MultiMooseEnum constructor.
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-")
const GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false)
const ModelGeochemicalDatabase mgd
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)
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)
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)
const ModelGeochemicalDatabase & mgd_kin
const GeochemicalDatabaseReader db("database/moose_testdb.json", true, true, false)
Class for reading geochemical reactions from a MOOSE geochemical database.
This class holds information about bulk composition, molalities, activities, activity coefficients,...
std::vector< Real > getSaturationIndices() const
Real getEquilibriumActivity(unsigned eqm_ind) const
Returns the value of activity for the equilibrium species with index eqm_index.
void enforceChargeBalance()
Enforces charge balance by altering the constraint_value and bulk_moles_old of the charge-balance spe...
Real getLog10K(unsigned j) 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,...
Real getEquilibriumActivityCoefficient(unsigned j) const
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,...
void setKineticMoles(unsigned kin, Real moles)
Sets the current AND old mole number for a kinetic species.
const std::vector< Real > & computeAndGetEquilibriumActivity()
Computes activity for all equilibrium species (_eqm_activity) and returns a reference to the vector.
void changeConstraintToBulk(unsigned basis_ind)
Changes the constraint to MOLES_BULK_SPECIES (or MOLES_BULK_WATER if basis_ind = 0) for the basis ind...
void performSwap(unsigned swap_out_of_basis, unsigned swap_into_basis)
Perform the basis swap, and ensure that the resulting system is consistent.
unsigned getNumInBasis() const
returns the number of species in the basis
unsigned getNumInEquilibrium() const
returns the number of species in equilibrium with the basis components
Real getTotalChargeOld() const
void setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(const std::vector< std::string > &names, const std::vector< Real > &values, const std::vector< bool > &constraints_from_molalities)
Sets:
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.
void setConstraintValue(unsigned basis_ind, Real value)
Set the constraint value for the basis species.
const std::vector< Real > & getBulkMolesOld() const
const ModelGeochemicalDatabase & getModelGeochemicalDatabase() const
const std::vector< Real > & getSolventMassAndFreeMolalityAndMineralMoles() const
Real getBasisActivity(unsigned i) const
void closeSystem()
Changes a KG_SOLVENT_WATER constraint to MOLES_BULK_WATER (if there is such a constraint) and all FRE...
Real getEquilibriumMolality(unsigned j) const
void setTemperature(Real temperature)
Sets the temperature to the given quantity.
void setModelGeochemicalDatabase(const ModelGeochemicalDatabase &mgd)
Copies a ModelGeochemicalDatabase into our _mgd structure.
Real getTemperature() const
Computes activity coefficients for non-minerals and non-gases (since these species do not have activi...
Calculators to compute ionic strength and stoichiometric ionic strength.
Class to swap basis species with equilibrium species.
void setAdditionalValue(const std::string &names)
Constructs and stores a minimal amount of information that is pertinent to the user-defined geochemic...
const ModelGeochemicalDatabase & modelGeochemicalDatabase() const
Return a reference to the ModelGeochemicalDatabase structure.
std::vector< std::string > originalBasisNames() const
Holds a user-specified description of a kinetic rate.
Data structure to hold all relevant information from the database file.
std::vector< bool > eqm_species_mineral
eqm_species_mineral[i] = true iff the i^th equilibrium species is a mineral
DenseMatrix< Real > eqm_stoichiometry
eqm_stoichiometry(i, j) = stoichiometric coefficient for equilibrium species "i" in terms of the basi...
std::unordered_map< std::string, unsigned > eqm_species_index
eqm_species_index[name] = index of the equilibrium species (secondary aqueous species,...
std::vector< std::string > basis_species_name
basis_species_name[j] = name of the j^th basis species
std::unordered_map< std::string, unsigned > basis_species_index
basis_species_index[name] = index of the basis species, within all ModelGeochemicalDatabase internal ...
std::unordered_map< std::string, unsigned > kin_species_index
kin_species_index[name] = index of the kinetic species, within all ModelGeochemicalDatabase internal ...
std::vector< std::string > eqm_species_name
eqm_species_name[i] = name of the i^th eqm species