10 #include "gtest/gtest.h" 18 TEST(GeochemistrySpeciesSwapperTest, bulkCompositionException)
26 {
"H2O",
"H+",
">(s)FeOH",
">(w)FeOH",
"Fe++",
"HCO3-",
"O2(aq)"},
39 DenseVector<Real> bulk_composition(8);
40 swapper.performSwap(
mgd, bulk_composition,
"H+",
"(O-phth)--");
41 FAIL() <<
"Missing expected exception.";
43 catch (
const std::exception & e)
45 std::string msg(e.what());
46 ASSERT_TRUE(msg.find(
"GeochemistrySpeciesSwapper: bulk_composition has size 8 which differs " 47 "from the basis size") != std::string::npos)
48 <<
"Failed with unexpected error message: " << msg;
53 TEST(GeochemistrySpeciesSwapperTest, swapExceptions)
61 {
"H2O",
"H+",
">(s)FeOH",
">(w)FeOH",
"Fe++",
"HCO3-",
"O2(aq)"},
75 FAIL() <<
"Missing expected exception.";
77 catch (
const std::exception & e)
79 std::string msg(e.what());
80 ASSERT_TRUE(msg.find(
"CO2(aq) is not in the basis, so cannot be removed " 81 "from the basis") != std::string::npos)
82 <<
"Failed with unexpected error message: " << msg;
87 swapper.checkSwap(
mgd,
"CH4(g)fake",
"CO2(aq)");
88 FAIL() <<
"Missing expected exception.";
90 catch (
const std::exception & e)
92 std::string msg(e.what());
93 ASSERT_TRUE(msg.find(
"CH4(g)fake is not in the basis, so cannot be removed " 94 "from the basis") != std::string::npos)
95 <<
"Failed with unexpected error message: " << msg;
100 swapper.checkSwap(
mgd,
"H+",
">(s)FeOH");
101 FAIL() <<
"Missing expected exception.";
103 catch (
const std::exception & e)
105 std::string msg(e.what());
106 ASSERT_TRUE(msg.find(
">(s)FeOH is not an equilibrium species, so cannot be " 108 "equilibrium species list") != std::string::npos)
109 <<
"Failed with unexpected error message: " << msg;
114 swapper.checkSwap(
mgd,
"H2O",
"CO2(aq)");
115 FAIL() <<
"Missing expected exception.";
117 catch (
const std::exception & e)
119 std::string msg(e.what());
120 ASSERT_TRUE(msg.find(
"Cannot remove H2O from the basis") != std::string::npos)
121 <<
"Failed with unexpected error message: " << msg;
126 swapper.checkSwap(
mgd, 123, 0);
127 FAIL() <<
"Missing expected exception.";
129 catch (
const std::exception & e)
131 std::string msg(e.what());
132 ASSERT_TRUE(msg.find(
"123 exceeds the number of basis species in the problem") !=
134 <<
"Failed with unexpected error message: " << msg;
139 swapper.checkSwap(
mgd, 0, 0);
140 FAIL() <<
"Missing expected exception.";
142 catch (
const std::exception & e)
144 std::string msg(e.what());
145 ASSERT_TRUE(msg.find(
"Cannot remove H2O from the basis") != std::string::npos)
146 <<
"Failed with unexpected error message: " << msg;
151 swapper.checkSwap(
mgd, 1, 123);
152 FAIL() <<
"Missing expected exception.";
154 catch (
const std::exception & e)
156 std::string msg(e.what());
157 ASSERT_TRUE(msg.find(
"123 exceeds the number of equilibrium species in the problem") !=
159 <<
"Failed with unexpected error message: " << msg;
164 swapper.performSwap(
mgd,
">(s)FeOH",
"CO2(aq)");
165 FAIL() <<
"Missing expected exception.";
167 catch (
const std::exception & e)
169 std::string msg(e.what());
170 ASSERT_TRUE(msg.find(
"Matrix is not invertible, which signals an invalid basis swap") !=
172 <<
"Failed with unexpected error message: " << msg;
178 swapper.checkSwap(
mgd,
"Fe++",
"Fe+++");
179 FAIL() <<
"Missing expected exception.";
181 catch (
const std::exception & e)
183 std::string msg(e.what());
185 msg.find(
"GeochemistrySpeciesSwapper constructed with incorrect basis_species size") !=
187 <<
"Failed with unexpected error message: " << msg;
192 swapper.checkSwap(
mgd,
">(s)FeOH",
">(s)FeO-");
193 FAIL() <<
"Missing expected exception.";
195 catch (
const std::exception & e)
197 std::string msg(e.what());
199 msg.find(
"Equilibrium species >(s)FeO- is involved in surface sorption so cannot be " 200 "swapped into the basis. If this is truly necessary, code enhancements will need " 201 "to be made including: recording whether basis species are involved in surface " 202 "sorption, including them in the surface-potential calculations, and carefully " 203 "swapping surface-potential-modified equilibrium constants") != std::string::npos)
204 <<
"Failed with unexpected error message: " << msg;
209 TEST(GeochemistrySpeciesSwapperTest, swap1)
215 database, {
"H2O",
"Ca++",
"HCO3-",
"H+"}, {
"Calcite"}, {}, {}, {}, {},
"O2(aq)",
"e-");
218 DenseVector<Real> bulk_composition(4);
219 bulk_composition(0) = 0.5;
220 bulk_composition(1) = 1.0;
221 bulk_composition(2) = 2.5;
222 bulk_composition(3) = 3.0;
239 if (species.first ==
"Calcite")
252 if (species.first ==
"Calcite")
258 for (
const auto & sp : {
"CO2(aq)",
"CO3--",
"CaCO3",
"CaOH+",
"OH-",
"Calcite"})
261 std::map<std::string, Real> charge_gold;
262 charge_gold[
"H2O"] = 0.0;
263 charge_gold[
"H+"] = 1.0;
264 charge_gold[
"HCO3-"] = -1.0;
265 charge_gold[
"Ca++"] = 2.0;
266 charge_gold[
"CO2(aq)"] = 0.0;
267 charge_gold[
"CO3--"] = -2.0;
268 charge_gold[
"CaCO3"] = 0.0;
269 charge_gold[
"CaOH+"] = 1.0;
270 charge_gold[
"OH-"] = -1.0;
271 charge_gold[
"Calcite"] = 0.0;
277 std::map<std::string, Real> radius_gold;
278 radius_gold[
"H2O"] = 0.0;
279 radius_gold[
"H+"] = 9.0;
280 radius_gold[
"HCO3-"] = 4.5;
281 radius_gold[
"Ca++"] = 6.0;
282 radius_gold[
"CO2(aq)"] = 4.0;
283 radius_gold[
"CO3--"] = 4.5;
284 radius_gold[
"CaCO3"] = 4.0;
285 radius_gold[
"CaOH+"] = 4.0;
286 radius_gold[
"OH-"] = 3.5;
287 radius_gold[
"Calcite"] = 0.0;
293 std::map<std::string, Real> molecular_weight_gold;
294 molecular_weight_gold[
"H2O"] = 18.0152;
295 molecular_weight_gold[
"H+"] = 1.0079;
296 molecular_weight_gold[
"HCO3-"] = 61.0171;
297 molecular_weight_gold[
"Ca++"] = 40.0800;
298 molecular_weight_gold[
"CO2(aq)"] = 44.0098;
299 molecular_weight_gold[
"CO3--"] = 60.0092;
300 molecular_weight_gold[
"CaCO3"] = 100.0892;
301 molecular_weight_gold[
"CaOH+"] = 57.0873;
302 molecular_weight_gold[
"OH-"] = 17.0073;
303 molecular_weight_gold[
"Calcite"] = 100.0892;
309 std::map<std::string, Real> molecular_volume_gold;
310 molecular_volume_gold[
"H2O"] = 0.0;
311 molecular_volume_gold[
"H+"] = 0.0;
312 molecular_volume_gold[
"HCO3-"] = 0.0;
313 molecular_volume_gold[
"O2(aq)"] = 0.0;
314 molecular_volume_gold[
"Ca++"] = 0.0;
315 molecular_volume_gold[
"CO2(aq)"] = 0.0;
316 molecular_volume_gold[
"CO3--"] = 0.0;
317 molecular_volume_gold[
"CaCO3"] = 0.0;
318 molecular_volume_gold[
"CaOH+"] = 0.0;
319 molecular_volume_gold[
"OH-"] = 0.0;
320 molecular_volume_gold[
"Calcite"] = 36.9340;
326 std::map<std::string, DenseMatrix<Real>> stoi_gold;
327 for (
const auto & sp : {
"CO2(aq)",
"CO3--",
"CaCO3",
"CaOH+",
"OH-",
"Calcite"})
330 stoi_gold[
"CO2(aq)"](0, 0) = -1;
331 stoi_gold[
"CO2(aq)"](0, 3) = 1;
332 stoi_gold[
"CO2(aq)"](0, 2) = 1;
333 stoi_gold[
"CO3--"](0, 2) = 1;
334 stoi_gold[
"CO3--"](0, 3) = -1;
335 stoi_gold[
"CaCO3"](0, 1) = 1;
336 stoi_gold[
"CaCO3"](0, 2) = 1;
337 stoi_gold[
"CaCO3"](0, 3) = -1;
338 stoi_gold[
"CaOH+"](0, 1) = 1;
339 stoi_gold[
"CaOH+"](0, 0) = 1;
340 stoi_gold[
"CaOH+"](0, 3) = -1;
341 stoi_gold[
"OH-"](0, 0) = 1;
342 stoi_gold[
"OH-"](0, 3) = -1;
343 stoi_gold[
"Calcite"](0, 1) = 1;
344 stoi_gold[
"Calcite"](0, 2) = 1;
345 stoi_gold[
"Calcite"](0, 3) = -1;
347 for (
const auto & sp : {
"CO2(aq)",
"CO3--",
"CaCO3",
"CaOH+",
"OH-",
"Calcite"})
370 swapper.performSwap(
mgd, bulk_composition,
"Ca++",
"Calcite");
374 for (
const auto & sp : {
"Calcite",
"H2O",
"HCO3-",
"H+"})
382 for (
const auto & sp : {
"CO2(aq)",
"CO3--",
"CaCO3",
"CaOH+",
"OH-",
"Ca++"})
395 for (
unsigned i = 0; i < 4; ++i)
396 gold_swap_matrix(i, i) = 1.0;
397 gold_swap_matrix(ca_posn, 0) = 0.0;
398 gold_swap_matrix(ca_posn, 1) = 1.0;
399 gold_swap_matrix(ca_posn, 2) = 1.0;
400 gold_swap_matrix(ca_posn, 3) = -1.0;
401 for (
unsigned i = 0; i < 4; ++i)
402 for (
unsigned j = 0;
j < 4; ++
j)
431 if (species.first ==
"Calcite")
446 if (species.first ==
"Calcite")
454 stoi_gold = std::map<std::string, DenseMatrix<Real>>();
455 for (
const auto & sp : {
"CO2(aq)",
"CO3--",
"CaCO3",
"CaOH+",
"OH-",
"Ca++"})
458 stoi_gold[
"CO2(aq)"](0, 0) = -1;
459 stoi_gold[
"CO2(aq)"](0, 3) = 1;
460 stoi_gold[
"CO2(aq)"](0, 2) = 1;
461 stoi_gold[
"CO3--"](0, 2) = 1;
462 stoi_gold[
"CO3--"](0, 3) = -1;
463 stoi_gold[
"CaCO3"](0, 1) = 1;
464 stoi_gold[
"CaOH+"](0, 1) = 1;
465 stoi_gold[
"CaOH+"](0, 0) = 1;
466 stoi_gold[
"CaOH+"](0, 2) = -1;
467 stoi_gold[
"OH-"](0, 0) = 1;
468 stoi_gold[
"OH-"](0, 3) = -1;
469 stoi_gold[
"Ca++"](0, 1) = 1;
470 stoi_gold[
"Ca++"](0, 2) = -1;
471 stoi_gold[
"Ca++"](0, 3) = 1;
473 for (
const auto & sp : {
"CO2(aq)",
"CO3--",
"CaCO3",
"CaOH+",
"OH-",
"Ca++"})
476 for (
unsigned i = 0; i < 4; ++i)
531 ASSERT_NEAR(bulk_composition(0), 0.5,
eps);
532 ASSERT_NEAR(bulk_composition(1), 1.0,
eps);
533 ASSERT_NEAR(bulk_composition(2), 1.5,
eps);
534 ASSERT_NEAR(bulk_composition(3), 4.0,
eps);
537 swapper.performSwap(
mgd, bulk_composition,
"Calcite",
"Ca++");
539 for (
unsigned i = 0; i < 4; ++i)
540 for (
unsigned j = 0;
j < 4; ++
j)
548 TEST(GeochemistrySpeciesSwapperTest, swap2)
554 {
"H2O",
"StoiCheckBasis",
"H+",
"HCO3-"},
588 for (
const auto & sp : {
"OH-",
"CO2(aq)",
"CO3--",
"StoiCheckRedox",
"StoiCheckGas"})
594 if (species.first ==
"StoiCheckGas")
599 std::map<std::string, Real> charge_gold;
600 charge_gold[
"H2O"] = 0.0;
601 charge_gold[
"StoiCheckBasis"] = 2.5;
602 charge_gold[
"H+"] = 1.0;
603 charge_gold[
"HCO3-"] = -1.0;
604 charge_gold[
"CO2(aq)"] = 0.0;
605 charge_gold[
"CO3--"] = -2.0;
606 charge_gold[
"OH-"] = -1.0;
607 charge_gold[
"StoiCheckRedox"] = 3.3;
608 charge_gold[
"StoiCheckGas"] = 0.0;
614 std::map<std::string, Real> radius_gold;
615 radius_gold[
"H2O"] = 0.0;
616 radius_gold[
"StoiCheckBasis"] = 6.54;
617 radius_gold[
"H+"] = 9.0;
618 radius_gold[
"HCO3-"] = 4.5;
619 radius_gold[
"CO2(aq)"] = 4.0;
620 radius_gold[
"CO3--"] = 4.5;
621 radius_gold[
"OH-"] = 3.5;
622 radius_gold[
"StoiCheckRedox"] = 9.9;
623 radius_gold[
"StoiCheckGas"] = 0.0;
629 std::map<std::string, Real> molecular_weight_gold;
630 molecular_weight_gold[
"H2O"] = 18.0152;
631 molecular_weight_gold[
"StoiCheckBasis"] = 55.8470;
632 molecular_weight_gold[
"H+"] = 1.0079;
633 molecular_weight_gold[
"HCO3-"] = 61.0171;
634 molecular_weight_gold[
"CO2(aq)"] = 44.0098;
635 molecular_weight_gold[
"CO3--"] = 60.0092;
636 molecular_weight_gold[
"OH-"] = 17.0073;
637 molecular_weight_gold[
"StoiCheckRedox"] = 55.8470;
638 molecular_weight_gold[
"StoiCheckGas"] = 28.0134;
644 std::map<std::string, Real> molecular_volume_gold;
645 molecular_volume_gold[
"H2O"] = 0.0;
646 molecular_volume_gold[
"StoiCheckBasis"] = 0.0;
647 molecular_volume_gold[
"H+"] = 0.0;
648 molecular_volume_gold[
"HCO3-"] = 0.0;
649 molecular_volume_gold[
"CO2(aq)"] = 0.0;
650 molecular_volume_gold[
"CO3--"] = 0.0;
651 molecular_volume_gold[
"OH-"] = 0.0;
652 molecular_volume_gold[
"StoiCheckRedox"] = 0.0;
653 molecular_volume_gold[
"StoiCheckGas"] = 0.0;
659 std::map<std::string, DenseMatrix<Real>> stoi_gold;
660 for (
const auto & sp : {
"OH-",
"CO2(aq)",
"CO3--",
"StoiCheckRedox",
"StoiCheckGas"})
663 stoi_gold[
"CO2(aq)"](0, 0) = -1;
664 stoi_gold[
"CO2(aq)"](0, 2) = 1;
665 stoi_gold[
"CO2(aq)"](0, 3) = 1;
666 stoi_gold[
"CO3--"](0, 2) = -1;
667 stoi_gold[
"CO3--"](0, 3) = 1;
668 stoi_gold[
"OH-"](0, 0) = 1;
669 stoi_gold[
"OH-"](0, 2) = -1;
670 stoi_gold[
"StoiCheckRedox"](0, 0) = -0.5;
671 stoi_gold[
"StoiCheckRedox"](0, 1) = 1.5;
672 stoi_gold[
"StoiCheckRedox"](0, 2) = -1;
673 stoi_gold[
"StoiCheckGas"](0, 0) = 2.0;
674 stoi_gold[
"StoiCheckGas"](0, 1) = 3.0;
675 stoi_gold[
"StoiCheckGas"](0, 2) = -5.0;
677 for (
const auto & sp : {
"OH-",
"CO2(aq)",
"CO3--",
"StoiCheckRedox",
"StoiCheckGas"})
692 swapper.performSwap(
mgd,
"StoiCheckBasis",
"StoiCheckGas");
696 for (
const auto & sp : {
"H2O",
"StoiCheckGas",
"H+",
"HCO3-"})
704 for (
const auto & sp : {
"OH-",
"CO2(aq)",
"CO3--",
"StoiCheckRedox",
"StoiCheckBasis"})
749 if (species.first ==
"StoiCheckGas")
763 stoi_gold = std::map<std::string, DenseMatrix<Real>>();
764 for (
const auto & sp : {
"OH-",
"CO2(aq)",
"CO3--",
"StoiCheckRedox",
"StoiCheckBasis"})
767 stoi_gold[
"CO2(aq)"](0, 0) = -1;
768 stoi_gold[
"CO2(aq)"](0, 2) = 1;
769 stoi_gold[
"CO2(aq)"](0, 3) = 1;
770 stoi_gold[
"CO3--"](0, 2) = -1;
771 stoi_gold[
"CO3--"](0, 3) = 1;
772 stoi_gold[
"OH-"](0, 0) = 1;
773 stoi_gold[
"OH-"](0, 2) = -1;
774 stoi_gold[
"StoiCheckRedox"](0, 0) = -1.5;
775 stoi_gold[
"StoiCheckRedox"](0, 1) = 0.5;
776 stoi_gold[
"StoiCheckRedox"](0, 2) = 1.5;
777 stoi_gold[
"StoiCheckBasis"](0, 0) = -2.0 / 3.0;
778 stoi_gold[
"StoiCheckBasis"](0, 1) = 1.0 / 3.0;
779 stoi_gold[
"StoiCheckBasis"](0, 2) = 5.0 / 3.0;
781 for (
const auto & sp : {
"OH-",
"CO2(aq)",
"CO3--",
"StoiCheckRedox",
"StoiCheckBasis"})
784 for (
unsigned i = 0; i < 4; ++i)
821 const Real ot = -1.0 / 3.0;
822 const Real tt = -2.0 / 3.0;
824 tt * (-10.0553) + ot * (-2.9620),
827 tt * (-8.4878) + ot * (-3.1848),
830 tt * (-6.6954) + ot * (-3.3320),
833 tt * (-5.0568) + ot * (-3.2902),
836 tt * (-3.4154) + ot * (-3.1631),
839 tt * (-2.0747) + ot * (-2.9499),
842 tt * (-.8908) + ot * (-2.7827),
845 tt * (.2679) + ot * (-2.3699),
850 TEST(GeochemistrySpeciesSwapperTest, swap3)
856 {
"H2O",
"H+",
">(s)FeOH",
">(w)FeOH",
"Fe++",
"HCO3-",
"O2(aq)"},
859 {
"Fe(OH)3(ppd)",
"Fe(OH)3(ppd)fake"},
871 {
"H+",
"O2(aq)",
"CO3--"},
892 {
"O2(aq)",
"Fe+++",
"H2O"},
913 {
"CO2(aq)",
"Fe+++",
"Fe++"},
951 if (species.first ==
"Fe(OH)3(ppd)" || species.first ==
"Fe(OH)3(ppd)fake")
959 if (species.first ==
"CH4(g)fake")
965 if (species.first ==
">(s)FeOH" || species.first ==
">(w)FeOH")
972 if (species.first ==
"Fe(OH)3(ppd)" || species.first ==
"Fe(OH)3(ppd)fake" ||
973 species.first ==
">(s)FeO-")
979 for (
const auto & sp : {
"CO2(aq)",
"CO3--",
"OH-",
"CH4(aq)",
"Fe+++",
"CH4(g)fake"})
983 for (
const auto & sp : {
"Fe(OH)3(ppd)",
"Fe(OH)3(ppd)fake",
"(O-phth)--",
">(s)FeO-"})
986 std::map<std::string, Real> charge_gold;
987 charge_gold[
"H2O"] = 0.0;
988 charge_gold[
"H+"] = 1.0;
989 charge_gold[
">(s)FeOH"] = 0.0;
990 charge_gold[
">(w)FeOH"] = 0.0;
991 charge_gold[
"Fe++"] = 2.0;
992 charge_gold[
"HCO3-"] = -1.0;
993 charge_gold[
"O2(aq)"] = 0.0;
994 charge_gold[
"CO2(aq)"] = 0.0;
995 charge_gold[
"CO3--"] = -2.0;
996 charge_gold[
"OH-"] = -1.0;
997 charge_gold[
"CH4(aq)"] = 0.0;
998 charge_gold[
"Fe+++"] = 3.0;
999 charge_gold[
"CH4(g)fake"] = 0.0;
1000 charge_gold[
"Fe(OH)3(ppd)"] = 0.0;
1001 charge_gold[
"Fe(OH)3(ppd)fake"] = 0.0;
1002 charge_gold[
"(O-phth)--"] = -2.0;
1003 charge_gold[
">(s)FeO-"] = -1.0;
1011 std::map<std::string, Real> radius_gold;
1012 radius_gold[
"H2O"] = 0.0;
1013 radius_gold[
"H+"] = 9.0;
1014 radius_gold[
">(s)FeOH"] = 0.0;
1015 radius_gold[
">(w)FeOH"] = 0.0;
1016 radius_gold[
"Fe++"] = 6.0;
1017 radius_gold[
"HCO3-"] = 4.5;
1018 radius_gold[
"O2(aq)"] = -0.5;
1019 radius_gold[
"CO2(aq)"] = 4.0;
1020 radius_gold[
"CO3--"] = 4.5;
1021 radius_gold[
"OH-"] = 3.5;
1022 radius_gold[
"CH4(aq)"] = -0.5;
1023 radius_gold[
"Fe+++"] = 9.0;
1024 radius_gold[
"CH4(g)fake"] = 0.0;
1030 std::map<std::string, Real> molecular_weight_gold;
1031 molecular_weight_gold[
"H2O"] = 18.0152;
1032 molecular_weight_gold[
"H+"] = 1.0079;
1033 molecular_weight_gold[
">(s)FeOH"] = 72.8543;
1034 molecular_weight_gold[
">(w)FeOH"] = 1234.567;
1035 molecular_weight_gold[
"Fe++"] = 55.8470;
1036 molecular_weight_gold[
"HCO3-"] = 61.0171;
1037 molecular_weight_gold[
"O2(aq)"] = 31.9988;
1038 molecular_weight_gold[
"CO2(aq)"] = 44.0098;
1039 molecular_weight_gold[
"CO3--"] = 60.0092;
1040 molecular_weight_gold[
"OH-"] = 17.0073;
1041 molecular_weight_gold[
"CH4(aq)"] = 16.0426;
1042 molecular_weight_gold[
"Fe+++"] = 55.8470;
1043 molecular_weight_gold[
"CH4(g)fake"] = 16.0426;
1044 molecular_weight_gold[
"Fe(OH)3(ppd)"] = 106.8689;
1045 molecular_weight_gold[
"Fe(OH)3(ppd)fake"] = 106.8689;
1046 molecular_weight_gold[
"(O-phth)--"] = 164.1172;
1047 molecular_weight_gold[
">(s)FeO-"] = 71.8464;
1055 std::map<std::string, Real> molecular_volume_gold;
1056 molecular_volume_gold[
"H2O"] = 0.0;
1057 molecular_volume_gold[
"H+"] = 0.0;
1058 molecular_volume_gold[
">(s)FeOH"] = 0.0;
1059 molecular_volume_gold[
">(w)FeOH"] = 0.0;
1060 molecular_volume_gold[
"Fe++"] = 0.0;
1061 molecular_volume_gold[
"HCO3-"] = 0.0;
1062 molecular_volume_gold[
"O2(aq)"] = 0.0;
1063 molecular_volume_gold[
"CO2(aq)"] = 0.0;
1064 molecular_volume_gold[
"CO3--"] = 0.0;
1065 molecular_volume_gold[
"OH-"] = 0.0;
1066 molecular_volume_gold[
"CH4(aq)"] = 0.0;
1067 molecular_volume_gold[
"Fe+++"] = 0.0;
1068 molecular_volume_gold[
"CH4(g)fake"] = 0.0;
1069 molecular_volume_gold[
"Fe(OH)3(ppd)"] = 34.3200;
1070 molecular_volume_gold[
"Fe(OH)3(ppd)fake"] = 34.3200;
1071 molecular_volume_gold[
"(O-phth)--"] = 0.0;
1072 molecular_volume_gold[
">(s)FeO-"] = 0.0;
1080 std::map<std::string, DenseMatrix<Real>> stoi_gold;
1081 for (
const auto & sp : {
"CO2(aq)",
1094 stoi_gold[
"CO2(aq)"](0, 0) = -1;
1095 stoi_gold[
"CO2(aq)"](0, 1) = 1;
1096 stoi_gold[
"CO2(aq)"](0, 5) = 1;
1097 stoi_gold[
"CO3--"](0, 5) = 1;
1098 stoi_gold[
"CO3--"](0, 1) = -1;
1099 stoi_gold[
"OH-"](0, 0) = 1;
1100 stoi_gold[
"OH-"](0, 1) = -1;
1101 stoi_gold[
"CH4(aq)"](0, 0) = 1;
1102 stoi_gold[
"CH4(aq)"](0, 1) = 1;
1103 stoi_gold[
"CH4(aq)"](0, 5) = 1;
1104 stoi_gold[
"CH4(aq)"](0, 6) = -2;
1105 stoi_gold[
"Fe+++"](0, 0) = -0.5;
1106 stoi_gold[
"Fe+++"](0, 4) = 1;
1107 stoi_gold[
"Fe+++"](0, 1) = 1;
1108 stoi_gold[
"Fe+++"](0, 6) = 0.25;
1109 stoi_gold[
"CH4(g)fake"](0, 0) = 3;
1110 stoi_gold[
"CH4(g)fake"](0, 4) = -2;
1111 stoi_gold[
"CH4(g)fake"](0, 5) = 3.5;
1112 stoi_gold[
"CH4(g)fake"](0, 6) = -4.5;
1113 stoi_gold[
"Fe(OH)3(ppd)"](0, 1) = -2;
1114 stoi_gold[
"Fe(OH)3(ppd)"](0, 4) = 1;
1115 stoi_gold[
"Fe(OH)3(ppd)"](0, 0) = 2.5;
1116 stoi_gold[
"Fe(OH)3(ppd)"](0, 6) = 0.25;
1117 stoi_gold[
"Fe(OH)3(ppd)fake"](0, 1) = -1;
1118 stoi_gold[
"Fe(OH)3(ppd)fake"](0, 0) = 2;
1119 stoi_gold[
"Fe(OH)3(ppd)fake"](0, 4) = 2;
1120 stoi_gold[
"Fe(OH)3(ppd)fake"](0, 6) = 0.5;
1121 stoi_gold[
"(O-phth)--"](0, 0) = -5;
1122 stoi_gold[
"(O-phth)--"](0, 5) = 8;
1123 stoi_gold[
"(O-phth)--"](0, 1) = 6;
1124 stoi_gold[
"(O-phth)--"](0, 6) = -7.5;
1125 stoi_gold[
">(s)FeO-"](0, 2) = 1;
1126 stoi_gold[
">(s)FeO-"](0, 1) = -1;
1128 for (
const auto & sp : {
"CO2(aq)",
"CO3--",
"OH-",
"CH4(aq)",
"Fe+++",
"CH4(g)fake"})
1134 for (
const auto & sp : {
"Fe(OH)3(ppd)",
"Fe(OH)3(ppd)fake",
"(O-phth)--",
">(s)FeO-"})
1140 const std::vector<Real> feoh3ppd = {
1141 6.1946, 4.8890, 3.4608, 2.2392, 1.1150, 0.2446, -0.5504, -1.5398};
1142 const std::vector<Real> fe3 = {10.0553, 8.4878, 6.6954, 5.0568, 3.4154, 2.0747, 0.8908, -0.2679};
1143 const std::vector<Real> ophth = {
1144 594.3211, 542.8292, 482.3612, 425.9738, 368.7004, 321.8658, 281.8216, 246.4849};
1145 const std::vector<Real> feom = {8.93,
1154 for (
unsigned t = 0; t < 8; ++t)
1159 feoh3ppd[t] - 2 * fe3[t],
1166 std::vector<std::vector<Real>> kin_rate_gold(3, std::vector<Real>(7 + 6, 0.0));
1167 std::vector<std::vector<Real>> kin_monod_gold(3, std::vector<Real>(7 + 6, 0.0));
1168 std::vector<std::vector<Real>> kin_k_gold(3, std::vector<Real>(7 + 6, 0.0));
1196 for (
unsigned i = 0; i < 3; ++i)
1197 for (
unsigned j = 0;
j < 7 + 6; ++
j)
1199 EXPECT_EQ(
mgd.
kin_rate[i].promoting_indices[
j], kin_rate_gold[i][
j]);
1200 EXPECT_EQ(
mgd.
kin_rate[i].promoting_monod_indices[
j], kin_monod_gold[i][
j]);
1201 EXPECT_EQ(
mgd.
kin_rate[i].promoting_half_saturation[
j], kin_k_gold[i][
j]);
1211 swapper.performSwap(
mgd,
"O2(aq)",
"Fe+++");
1215 for (
const auto & sp : {
"H2O",
"H+",
">(s)FeOH",
">(w)FeOH",
"Fe++",
"HCO3-",
"Fe+++"})
1223 for (
const auto & sp : {
"CO2(aq)",
"CO3--",
"OH-",
"CH4(aq)",
"O2(aq)",
"CH4(g)fake"})
1232 for (
const auto & sp : {
"Fe(OH)3(ppd)",
"Fe(OH)3(ppd)fake",
"(O-phth)--",
">(s)FeO-"})
1278 if (species.first ==
"Fe(OH)3(ppd)" || species.first ==
"Fe(OH)3(ppd)fake")
1287 if (species.first ==
"CH4(g)fake")
1294 if (species.first ==
">(s)FeOH" || species.first ==
">(w)FeOH")
1301 if (species.first ==
"Fe(OH)3(ppd)" || species.first ==
"Fe(OH)3(ppd)fake" ||
1302 species.first ==
">(s)FeO-")
1308 for (
const auto & sp : {
"CO2(aq)",
1321 stoi_gold[
"CO2(aq)"](0, 0) = -1;
1322 stoi_gold[
"CO2(aq)"](0, 1) = 1;
1323 stoi_gold[
"CO2(aq)"](0, 5) = 1;
1324 stoi_gold[
"CO3--"](0, 5) = 1;
1325 stoi_gold[
"CO3--"](0, 1) = -1;
1326 stoi_gold[
"OH-"](0, 0) = 1;
1327 stoi_gold[
"OH-"](0, 1) = -1;
1328 stoi_gold[
"CH4(aq)"](0, 0) = -3;
1329 stoi_gold[
"CH4(aq)"](0, 1) = 9;
1330 stoi_gold[
"CH4(aq)"](0, 4) = 8;
1331 stoi_gold[
"CH4(aq)"](0, 5) = 1;
1332 stoi_gold[
"CH4(aq)"](0, 6) = -8;
1333 stoi_gold[
"O2(aq)"](0, 0) = 2;
1334 stoi_gold[
"O2(aq)"](0, 1) = -4;
1335 stoi_gold[
"O2(aq)"](0, 4) = -4;
1336 stoi_gold[
"O2(aq)"](0, 6) = 4;
1337 stoi_gold[
"CH4(g)fake"](0, 0) = -6;
1338 stoi_gold[
"CH4(g)fake"](0, 1) = 18;
1339 stoi_gold[
"CH4(g)fake"](0, 4) = 16;
1340 stoi_gold[
"CH4(g)fake"](0, 5) = 3.5;
1341 stoi_gold[
"CH4(g)fake"](0, 6) = -18;
1342 stoi_gold[
"Fe(OH)3(ppd)"](0, 0) = 3;
1343 stoi_gold[
"Fe(OH)3(ppd)"](0, 1) = -3;
1344 stoi_gold[
"Fe(OH)3(ppd)"](0, 6) = 1;
1345 stoi_gold[
"Fe(OH)3(ppd)fake"](0, 0) = 3;
1346 stoi_gold[
"Fe(OH)3(ppd)fake"](0, 1) = -3;
1347 stoi_gold[
"Fe(OH)3(ppd)fake"](0, 6) = 2;
1348 stoi_gold[
"(O-phth)--"](0, 0) = -5 - 7.5 * 2;
1349 stoi_gold[
"(O-phth)--"](0, 1) = 6 + 7.5 * 4;
1350 stoi_gold[
"(O-phth)--"](0, 4) = 7.5 * 4;
1351 stoi_gold[
"(O-phth)--"](0, 5) = 8;
1352 stoi_gold[
"(O-phth)--"](0, 6) = -7.5 * 4;
1353 stoi_gold[
">(s)FeO-"](0, 2) = 1;
1354 stoi_gold[
">(s)FeO-"](0, 1) = -1;
1356 for (
const auto & sp : {
"CO2(aq)",
"CO3--",
"OH-",
"CH4(aq)",
"O2(aq)",
"CH4(g)fake"})
1359 for (
unsigned i = 0; i < 7; ++i)
1363 for (
const auto & sp : {
"Fe(OH)3(ppd)",
"Fe(OH)3(ppd)fake",
"(O-phth)--",
">(s)FeO-"})
1366 for (
unsigned i = 0; i < 7; ++i)
1370 for (
unsigned t = 0; t < 8; ++t)
1375 ophth[t] - (7.5 / 0.25) * fe3[t],
1381 std::vector<std::vector<Real>> new_kin_rate_gold(3, std::vector<Real>(7 + 6, 0.0));
1382 std::vector<std::vector<Real>> new_kin_monod_gold(3, std::vector<Real>(7 + 6, 0.0));
1383 std::vector<std::vector<Real>> new_kin_k_gold(3, std::vector<Real>(7 + 6, 0.0));
1411 for (
unsigned i = 0; i < 3; ++i)
1412 for (
unsigned j = 0;
j < 7 + 6; ++
j)
1414 EXPECT_EQ(
mgd.
kin_rate[i].promoting_indices[
j], new_kin_rate_gold[i][
j]);
1415 EXPECT_EQ(
mgd.
kin_rate[i].promoting_monod_indices[
j], new_kin_monod_gold[i][
j]);
1416 EXPECT_EQ(
mgd.
kin_rate[i].promoting_half_saturation[
j], new_kin_k_gold[i][
j]);
1420 EXPECT_EQ(
mgd.
kin_rate[0].progeny_index, 7 + fe3_posn);
1422 EXPECT_EQ(
mgd.
kin_rate[2].progeny_index, o2aq_posn);
1426 TEST(GeochemistrySpeciesSwapperTest, swap_redox)
1431 {
"H2O",
"H+",
"Fe++",
"O2(aq)",
"HCO3-",
"(O-phth)--",
"Fe+++"},
1444 const Real p = 1.0 / 4.0 / 7.5;
1446 const bool fe3_is_slot_one =
1448 const unsigned fe3_slot = (fe3_is_slot_one ? 1 : 2);
1449 const unsigned ophth_slot = (fe3_is_slot_one ? 2 : 1);
1478 swapper.performSwap(
mgd,
"HCO3-",
"CO3--");
1486 for (
unsigned i = 0; i < 7; ++i)
1488 gold_swap_matrix(i, i) = 1.0;
1489 gold_swap_matrix(hco3_posn, i) = 0.0;
1491 gold_swap_matrix(hco3_posn, 1) = -1.0;
1492 gold_swap_matrix(hco3_posn, 4) = 1.0;
1493 for (
unsigned i = 0; i < 7; ++i)
1494 for (
unsigned j = 0;
j < 7; ++
j)
1510 mgd.
redox_log10K(ophth_slot, 0), orig_log10K(ophth_slot, 0) + 8.0 * p * 10.6169, 1.0E-8);
1512 mgd.
redox_log10K(ophth_slot, 1), orig_log10K(ophth_slot, 1) + 8.0 * p * 10.3439, 1.0E-8);
1514 mgd.
redox_log10K(ophth_slot, 2), orig_log10K(ophth_slot, 2) + 8.0 * p * 10.2092, 1.0E-8);
1516 mgd.
redox_log10K(ophth_slot, 3), orig_log10K(ophth_slot, 3) + 8.0 * p * 10.2793, 1.0E-8);
1518 mgd.
redox_log10K(ophth_slot, 4), orig_log10K(ophth_slot, 4) + 8.0 * p * 10.5131, 1.0E-8);
1520 mgd.
redox_log10K(ophth_slot, 5), orig_log10K(ophth_slot, 5) + 8.0 * p * 10.8637, 1.0E-8);
1522 mgd.
redox_log10K(ophth_slot, 6), orig_log10K(ophth_slot, 6) + 8.0 * p * 11.2860, 1.0E-8);
1524 mgd.
redox_log10K(ophth_slot, 7), orig_log10K(ophth_slot, 7) + 8.0 * p * 11.6319, 1.0E-8);
1526 for (
unsigned basis_i = 0; basis_i < 7; ++basis_i)
1528 for (
unsigned temp = 0; temp < 8; ++temp)
1529 EXPECT_NEAR(
mgd.
redox_log10K(fe3_slot, temp), orig_log10K(fe3_slot, temp), 1.0E-8);
1533 TEST(GeochemistrySpeciesSwapperTest, findBestEqmSwapException)
1539 database, {
"H2O",
"Ca++",
"HCO3-",
"H+"}, {
"Calcite"}, {}, {}, {}, {},
"O2(aq)",
"e-");
1546 const std::vector<Real> eqm_molality(6, 1.0);
1547 swapper.findBestEqmSwap(4,
mgd, eqm_molality,
false,
false,
false, best);
1548 FAIL() <<
"Missing expected exception.";
1550 catch (
const std::exception & e)
1552 std::string msg(e.what());
1553 ASSERT_TRUE(msg.find(
"basis index 4 must be less than 4") != std::string::npos)
1554 <<
"Failed with unexpected error message: " << msg;
1559 const std::vector<Real> eqm_molality(5, 1.0);
1560 swapper.findBestEqmSwap(1,
mgd, eqm_molality,
false,
false,
false, best);
1561 FAIL() <<
"Missing expected exception.";
1563 catch (
const std::exception & e)
1565 std::string msg(e.what());
1566 ASSERT_TRUE(msg.find(
"Size of eqm_molality is 5 which is not equal to the number of " 1567 "equilibrium species 6") != std::string::npos)
1568 <<
"Failed with unexpected error message: " << msg;
1573 TEST(GeochemistrySpeciesSwapperTest, findBestEqmSwap)
1579 database, {
"H2O",
"Ca++",
"HCO3-",
"H+"}, {
"Calcite"}, {}, {}, {}, {},
"O2(aq)",
"e-");
1586 const std::vector<Real> eqm_molality = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0};
1587 bool legit = swapper.
findBestEqmSwap(1,
mgd, eqm_molality,
false,
false,
false, best);
1589 EXPECT_EQ(best, (
unsigned)3);
1590 legit = swapper.findBestEqmSwap(1,
mgd, eqm_molality,
true,
false,
false, best);
1592 EXPECT_EQ(best, (
unsigned)5);
1593 legit = swapper.findBestEqmSwap(1,
mgd, eqm_molality,
true,
true,
true, best);
1595 EXPECT_EQ(best, (
unsigned)5);
1602 {
"H2O",
"H+",
">(s)FeOH",
">(w)FeOH",
"Fe++",
"HCO3-",
"O2(aq)"},
1613 const std::vector<Real> eqm_molality_s = {0.0, 1.0, 2.0, 3.0, 4.0, 50.0, 6.0E3, 7.0E3, 8.0E3};
1614 legit = swapper_s.
findBestEqmSwap(1, mgd_s, eqm_molality_s,
false,
false,
false, best);
1616 EXPECT_EQ(best, (
unsigned)5);
1617 legit = swapper_s.findBestEqmSwap(1, mgd_s, eqm_molality_s,
true,
false,
false, best);
1619 EXPECT_EQ(best, (
unsigned)7);
1620 legit = swapper_s.findBestEqmSwap(5, mgd_s, eqm_molality_s,
false,
true,
false, best);
1622 EXPECT_EQ(best, (
unsigned)8);
1623 legit = swapper_s.findBestEqmSwap(1, mgd_s, eqm_molality_s,
false,
false,
true, best);
1625 EXPECT_EQ(best, (
unsigned)6);
DenseMatrix< Real > redox_stoichiometry
redox_stoichiometry(i, j) = stoichiometric coefficients for i^th redox species that is in disequilibr...
std::vector< Real > eqm_species_molecular_weight
all quantities have a molecular weight (g)
std::string redox_lhs
the name of the species on the left-hand side of the redox equations.
std::vector< bool > kin_species_mineral
kin_species_mineral[j] = true iff the j^th kinetic species is a mineral
std::vector< bool > basis_species_transported
basis_species_transported[j] = true iff the j^th basis species is transported in reactive-transport s...
std::vector< Real > kin_species_charge
all kinetic quantities have a charge (mineral charge = 0)
std::vector< bool > eqm_species_gas
eqm_species_gas[i] = true iff the i^th equilibrium species is a gas
std::unordered_map< std::string, unsigned > basis_species_index
basis_species_index[name] = index of the basis species, within all ModelGeochemicalDatabase internal ...
const ModelGeochemicalDatabase mgd
DenseMatrix< Real > kin_stoichiometry
kin_stoichiometry(i, j) = stoichiometric coefficient for kinetic species "i" in terms of the basis sp...
std::vector< Real > kin_species_molecular_weight
all quantities have a molecular weight (g/mol)
std::vector< Real > kin_species_molecular_volume
all quantities have a molecular volume (cm^3/mol) (only nonzero for minerals, however) ...
std::unordered_map< std::string, unsigned > eqm_species_index
eqm_species_index[name] = index of the equilibrium species (secondary aqueous species, redox couples in equilibrium with the aqueous solution, minerals in equilibrium with the aqueous solution, gases in equilibrium with the aqueous solution) within all ModelGeochemicalDatabase internal datastrcutres, with given name
DenseMatrix< Real > redox_log10K
redox_log10K(i, j) = log10(equilibrium constant) for i^th redox species at the j^th temperature point...
DenseMatrix< Real > swap_to_original_basis
Swap matrix that allows expression in terms of the original basis.
const GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false)
std::unordered_map< std::string, unsigned > kin_species_index
kin_species_index[name] = index of the kinetic species, within all ModelGeochemicalDatabase internal ...
DenseMatrix< Real > eqm_stoichiometry
eqm_stoichiometry(i, j) = stoichiometric coefficient for equilibrium species "i" in terms of the basi...
std::vector< Real > basis_species_radius
all quantities have an ionic radius (Angstrom) for computing activity (mineral radius = 0...
Class to swap basis species with equilibrium species.
bool findBestEqmSwap(unsigned basis_ind, const ModelGeochemicalDatabase &mgd, const std::vector< Real > &eqm_molality, bool minerals_allowed, bool gas_allowed, bool sorption_allowed, unsigned &best_eqm_species) const
For the the given basis index, find the equilibrium index that it should be swapped with...
const ModelGeochemicalDatabase & modelGeochemicalDatabase() const
Return a reference to the ModelGeochemicalDatabase structure.
void addKineticRate(const KineticRateUserDescription &description)
Adds a rate description for kinetic_species.
Constructs and stores a minimal amount of information that is pertinent to the user-defined geochemic...
DenseMatrix< Real > kin_log10K
kin_log10K(i, j) = log10(equilibrium constant for the i^th kinetic species at the j^th temperature po...
std::vector< bool > kin_species_transported
kin_species_transported[j] = true iff the j^th kinetic species is transported in reactive-transport s...
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-")
std::vector< Real > eqm_species_radius
all quantities have an ionic radius (Angstrom) for computing activity (mineral radius = 0...
DenseMatrix< Real > eqm_log10K
eqm_log10K(i, j) = log10(equilibrium constant) for i^th equilibrium species at the j^th temperature p...
std::vector< Real > basis_species_molecular_weight
all quantities have a molecular weight (g)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
std::vector< bool > eqm_species_transported
eqm_species_transported[i] = true iff the i^th eqm species is transported in reactive-transport sims ...
std::vector< Real > eqm_species_molecular_volume
all quantities have a molecular volume (cm^3) (only nonzero for minerals, however) ...
std::vector< Real > basis_species_charge
all quantities have a charge (mineral charge = 0, gas charge = 0, oxide charge = 0) ...
Holds a user-specified description of a kinetic rate.
TEST(GeochemistrySpeciesSwapperTest, bulkCompositionException)
Check bulk_composition exception.
std::vector< unsigned > have_swapped_into_basis
Species that have been swapped into the basis.
std::vector< bool > basis_species_gas
basis_species_gas[j] = true iff the j^th basis species is a gas
std::vector< bool > basis_species_mineral
basis_species_mineral[j] = true iff the j^th basis species is a mineral
std::vector< Real > eqm_species_charge
all quantities have a charge (mineral charge = 0, gas charge = 0, oxide charge = 0) ...
std::vector< std::string > kin_species_name
kin_species_name[j] = name of the j^th kinetic species
Data structure to hold all relevant information from the database file.
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
std::vector< bool > eqm_species_mineral
eqm_species_mineral[i] = true iff the i^th equilibrium species is a mineral
std::vector< Real > basis_species_molecular_volume
all quantities have a molecular volume (cm^3) (only nonzero for minerals, however) ...
DenseMatrix sub_matrix(unsigned int row_id, unsigned int row_size, unsigned int col_id, unsigned int col_size) const
std::vector< std::string > eqm_species_name
eqm_species_name[i] = name of the i^th eqm species
std::vector< std::string > basis_species_name
basis_species_name[j] = name of the j^th basis species
Class for reading geochemical reactions from a MOOSE geochemical database.
std::vector< unsigned > have_swapped_out_of_basis
Species that have been swapped out of the basis.
std::vector< KineticRateDefinition > kin_rate
rates given to kinetic species.
void checkSwap(const ModelGeochemicalDatabase &mgd, const std::string &replace_this, const std::string &with_this)
Check that replacing the named basis species with the named equilibrium species is valid...