16 "kg_solvent_water bulk_composition bulk_composition_with_kinetic free_concentration " 17 "free_mineral activity log10activity fugacity log10fugacity");
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");
31 {
"H2O",
"Calcite",
"H+",
"HCO3-"},
32 {1.75, 3.0, 2.0, 1.0},
34 constraint_user_meaning,
49 EXPECT_EQ(_egs_calcite.getNumInEquilibrium(), (unsigned)6);
57 {
"H2O",
"H+",
">(s)FeOH",
">(w)FeOH",
"Fe++",
"HCO3-",
"O2(aq)"},
86 FAIL() <<
"Missing expected exception.";
88 catch (
const std::exception & e)
90 std::string msg(e.what());
91 ASSERT_TRUE(msg.find(
"swap_out_of_basis must have same length as swap_into_basis") !=
93 <<
"Failed with unexpected error message: " << msg;
102 {
"H2O",
"H+",
">(s)FeOH",
">(w)FeOH",
"Fe++",
"HCO3-",
"O2(aq)"},
131 FAIL() <<
"Missing expected exception.";
133 catch (
const std::exception & e)
135 std::string msg(e.what());
136 ASSERT_TRUE(msg.find(
"Cannot swap out H+ because it is the charge-balance species") !=
138 <<
"Failed with unexpected error message: " << msg;
146 {
"H2O",
"H+",
">(s)FeOH",
">(w)FeOH",
"Fe++",
"HCO3-",
"O2(aq)"},
175 FAIL() <<
"Missing expected exception.";
177 catch (
const std::exception & e)
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;
191 database, {
"H2O",
"H+",
"HCO3-"}, {}, {}, {}, {}, {},
"O2(aq)",
"e-");
194 std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
198 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
203 std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm_poor;
206 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu_poor;
219 {
"H2O",
"H+",
"HCO3-"},
229 FAIL() <<
"Missing expected exception.";
231 catch (
const std::exception & e)
233 std::string msg(e.what());
234 ASSERT_TRUE(msg.find(
"Constrained species names must have same length as constraint values") !=
236 <<
"Failed with unexpected error message: " << msg;
248 {
"H2O",
"H+",
"HCO3-"},
258 FAIL() <<
"Missing expected exception.";
260 catch (
const std::exception & e)
262 std::string msg(e.what());
264 msg.find(
"Constrained species names must have same length as constraint meanings") !=
266 <<
"Failed with unexpected error message: " << msg;
278 {
"H2O",
"H+",
"HCO3-"},
288 FAIL() <<
"Missing expected exception.";
290 catch (
const std::exception & e)
292 std::string msg(e.what());
293 ASSERT_TRUE(msg.find(
"Constrained species names must have same length as constraint units") !=
295 <<
"Failed with unexpected error message: " << msg;
317 FAIL() <<
"Missing expected exception.";
319 catch (
const std::exception & e)
321 std::string msg(e.what());
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") !=
326 <<
"Failed with unexpected error message: " << msg;
334 database, {
"H2O",
"H+",
"HCO3-"}, {}, {}, {}, {}, {},
"O2(aq)",
"e-");
337 std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
341 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
355 {
"H2O",
"H+",
"H20"},
365 FAIL() <<
"Missing expected exception.";
367 catch (
const std::exception & e)
369 std::string msg(e.what());
370 ASSERT_TRUE(msg.find(
"The basis species HCO3- must appear in the constrained species list") !=
372 <<
"Failed with unexpected error message: " << msg;
384 {
"H2O",
"H+",
"HCO3-"},
394 FAIL() <<
"Missing expected exception.";
396 catch (
const std::exception & e)
398 std::string msg(e.what());
399 ASSERT_TRUE(msg.find(
"The basis species CO2(aq) must appear in the constrained species list") !=
401 <<
"Failed with unexpected error message: " << msg;
410 database, {
"H2O",
"H+",
"HCO3-"}, {}, {}, {}, {}, {},
"O2(aq)",
"e-");
413 std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
417 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
431 {
"H2O",
"H+",
"CO2(aq)"},
441 FAIL() <<
"Missing expected exception.";
443 catch (
const std::exception & e)
445 std::string msg(e.what());
446 ASSERT_TRUE(msg.find(
"Specified mass of solvent water must be positive: you entered -1") !=
448 <<
"Failed with unexpected error message: " << msg;
456 database, {
"H2O",
"H+",
"HCO3-"}, {}, {}, {}, {}, {},
"O2(aq)",
"e-");
459 std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
463 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
477 {
"H2O",
"H+",
"CO2(aq)"},
487 FAIL() <<
"Missing expected exception.";
489 catch (
const std::exception & e)
491 std::string msg(e.what());
492 ASSERT_TRUE(msg.find(
"Specified activity values must be positive: you entered -1") !=
494 <<
"Failed with unexpected error message: " << msg;
503 database, {
"H2O",
"H+",
"HCO3-",
"O2(aq)"}, {}, {
"O2(g)"}, {}, {}, {},
"O2(aq)",
"e-");
506 std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
511 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
526 {
"H2O",
"H+",
"O2(g)",
"HCO3-"},
527 {1.0, 2.0, -3.0, 4.0},
536 FAIL() <<
"Missing expected exception.";
538 catch (
const std::exception & e)
540 std::string msg(e.what());
541 ASSERT_TRUE(msg.find(
"Specified fugacity values must be positive: you entered -3") !=
543 <<
"Failed with unexpected error message: " << msg;
552 database, {
"H2O",
"H+",
"HCO3-"}, {}, {}, {}, {}, {},
"O2(aq)",
"e-");
555 std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
559 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
573 {
"H2O",
"H+",
"CO2(aq)"},
583 FAIL() <<
"Missing expected exception.";
585 catch (
const std::exception & e)
587 std::string msg(e.what());
588 ASSERT_TRUE(msg.find(
"Specified free concentration values must be positive: you entered -2") !=
590 <<
"Failed with unexpected error message: " << msg;
599 database, {
"H2O",
"H+",
"HCO3-",
"Ca++"}, {
"Calcite"}, {}, {}, {}, {},
"O2(aq)",
"e-");
602 std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
607 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
622 {
"H2O",
"Calcite",
"H+",
"HCO3-"},
623 {1.0, -3.0, 2.0, 1.0},
632 FAIL() <<
"Missing expected exception.";
634 catch (
const std::exception & e)
636 std::string msg(e.what());
637 ASSERT_TRUE(msg.find(
"Specified free mineral values must be positive: you entered -3") !=
639 <<
"Failed with unexpected error message: " << msg;
648 database, {
"H2O",
"H+",
"HCO3-"}, {}, {}, {}, {}, {},
"O2(aq)",
"e-");
651 std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
655 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
669 {
"H2O",
"H+",
"CO2(aq)"},
679 FAIL() <<
"Missing expected exception.";
681 catch (
const std::exception & e)
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;
695 database, {
"H2O",
"H+",
"HCO3-",
"O2(aq)"}, {}, {
"O2(g)"}, {}, {}, {},
"O2(aq)",
"e-");
698 std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
703 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
718 {
"H2O",
"H+",
"O2(g)",
"HCO3-"},
719 {1.0, 2.0, 3.0, 4.0},
728 FAIL() <<
"Missing expected exception.";
730 catch (
const std::exception & e)
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;
741 std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
746 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
761 {
"H2O",
"Calcite",
"H+",
"HCO3-"},
762 {1.0, 3.0, 2.0, 1.0},
771 FAIL() <<
"Missing expected exception.";
773 catch (
const std::exception & e)
775 std::string msg(e.what());
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") !=
780 <<
"Failed with unexpected error message: " << msg;
783 std::vector<GeochemicalSystem::ConstraintUserMeaningEnum>
cm2;
788 std::vector<GeochemistryUnitConverter::GeochemistryUnit>
cu2;
803 {
"H2O",
"Calcite",
"H+",
"HCO3-"},
804 {1.0, 3.0, 2.0, 1.0},
813 FAIL() <<
"Missing expected exception.";
815 catch (
const std::exception & e)
817 std::string msg(e.what());
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") !=
822 <<
"Failed with unexpected error message: " << msg;
829 std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
834 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
849 {
"H2O",
"Calcite",
"H+",
"HCO3-"},
850 {1.0, 3.0, 2.0, 1.0},
859 FAIL() <<
"Missing expected exception.";
861 catch (
const std::exception & e)
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;
869 std::vector<GeochemicalSystem::ConstraintUserMeaningEnum>
cm2;
874 std::vector<GeochemistryUnitConverter::GeochemistryUnit>
cu2;
889 {
"H2O",
"Calcite",
"H+",
"HCO3-"},
890 {1.0, 3.0, 2.0, 1.0},
899 FAIL() <<
"Missing expected exception.";
901 catch (
const std::exception & e)
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;
915 database, {
"H2O",
"H+",
"HCO3-"}, {}, {}, {}, {}, {},
"O2(aq)",
"e-");
918 std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
922 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
936 {
"H2O",
"H+",
"HCO3-"},
946 FAIL() <<
"Missing expected exception.";
948 catch (
const std::exception & e)
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;
961 database, {
"H2O",
"H+",
"HCO3-"}, {}, {}, {}, {}, {},
"O2(aq)",
"e-");
964 std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
968 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
982 {
"H2O",
"H+",
"HCO3-"},
992 FAIL() <<
"Missing expected exception.";
994 catch (
const std::exception & e)
996 std::string msg(e.what());
997 ASSERT_TRUE(msg.find(
"Species H+: units for bulk composition must be moles or mass") !=
999 <<
"Failed with unexpected error message: " << msg;
1008 database, {
"H2O",
"H+",
"HCO3-"}, {}, {}, {}, {}, {},
"O2(aq)",
"e-");
1011 std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
1015 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
1029 {
"H2O",
"H+",
"HCO3-"},
1039 FAIL() <<
"Missing expected exception.";
1041 catch (
const std::exception & e)
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;
1055 database, {
"H2O",
"H+",
"HCO3-",
"Ca++"}, {
"Calcite"}, {}, {}, {}, {},
"O2(aq)",
"e-");
1058 std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
1063 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
1078 {
"H2O",
"Calcite",
"H+",
"HCO3-"},
1079 {1.0, 3.0, 2.0, 1.0},
1088 FAIL() <<
"Missing expected exception.";
1090 catch (
const std::exception & e)
1092 std::string msg(e.what());
1095 "Species Calcite: units for free mineral quantities must be moles, mass or volume") !=
1097 <<
"Failed with unexpected error message: " << msg;
1106 database, {
"H2O",
"H+",
"HCO3-"}, {}, {}, {}, {}, {},
"O2(aq)",
"e-");
1109 std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
1113 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
1127 {
"H2O",
"H+",
"HCO3-"},
1137 FAIL() <<
"Missing expected exception.";
1139 catch (
const std::exception & e)
1141 std::string msg(e.what());
1143 msg.find(
"Species HCO3-: dimensionless units must be used when specifying activity") !=
1145 <<
"Failed with unexpected error message: " << msg;
1158 {
"H2O",
"H+",
"HCO3-"},
1168 FAIL() <<
"Missing expected exception.";
1170 catch (
const std::exception & e)
1172 std::string msg(e.what());
1174 msg.find(
"Species HCO3-: dimensionless units must be used when specifying log10activity") !=
1176 <<
"Failed with unexpected error message: " << msg;
1185 database, {
"H2O",
"H+",
"HCO3-",
"O2(aq)"}, {}, {
"O2(g)"}, {}, {}, {},
"O2(aq)",
"e-");
1188 std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
1193 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
1208 {
"H2O",
"H+",
"O2(g)",
"HCO3-"},
1209 {1.0, 2.0, 3.0, 4.0},
1218 FAIL() <<
"Missing expected exception.";
1220 catch (
const std::exception & e)
1222 std::string msg(e.what());
1224 msg.find(
"Species O2(g): dimensionless units must be used when specifying fugacity") !=
1226 <<
"Failed with unexpected error message: " << msg;
1240 {
"H2O",
"H+",
"O2(g)",
"HCO3-"},
1241 {1.0, 2.0, 3.0, 4.0},
1250 FAIL() <<
"Missing expected exception.";
1252 catch (
const std::exception & e)
1254 std::string msg(e.what());
1256 msg.find(
"Species O2(g): dimensionless units must be used when specifying log10fugacity") !=
1258 <<
"Failed with unexpected error message: " << msg;
1265 std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
1270 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
1283 {
"H2O",
"H+",
"HCO3-",
"Calcite"},
1284 {100.0, 3.0, 1.0, 2.0},
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);
1301 EXPECT_NEAR(egs.getBulkMolesOld()[1], 3.0, 1.0E-6);
1303 egs.getSolventMassAndFreeMolalityAndMineralMoles()[2], 1.0 / 61.0171, 1.0E-6);
1305 egs.getSolventMassAndFreeMolalityAndMineralMoles()[3], 2.0 / 36.9340, 1.0E-6);
1313 database, {
"H2O",
"H+",
"HCO3-",
"O2(aq)"}, {}, {
"O2(g)"}, {}, {}, {},
"O2(aq)",
"e-");
1316 std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
1321 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
1334 {
"H2O",
"H+",
"HCO3-",
"O2(g)"},
1335 {1.0, 2.0, 3.0, 4.0},
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);
1352 EXPECT_NEAR(egs.getBasisActivity(1), 1E2, 1.0E-6);
1353 EXPECT_NEAR(egs.getBulkMolesOld()[2], 2.0E-6 / 61.0171, 1.0E-6);
1354 EXPECT_NEAR(egs.getBasisActivity(3), 1E4, 1.0E-6);
1360 std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
1365 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
1380 {
"H2O",
"Calcite",
"H+",
"HCO3-"},
1381 {1.0, 3.0, 2.0, 1.0},
1390 FAIL() <<
"Missing expected exception.";
1392 catch (
const std::exception & e)
1394 std::string msg(e.what());
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") !=
1401 <<
"Failed with unexpected error message: " << msg;
1404 std::vector<GeochemicalSystem::ConstraintUserMeaningEnum>
cm2;
1409 std::vector<GeochemistryUnitConverter::GeochemistryUnit>
cu2;
1424 {
"H2O",
"Calcite",
"H+",
"HCO3-"},
1425 {1.0, 3.0, 2.0, 1.0},
1434 FAIL() <<
"Missing expected exception.";
1436 catch (
const std::exception & e)
1438 std::string msg(e.what());
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") !=
1445 <<
"Failed with unexpected error message: " << msg;
1461 {
"H2O",
"Calcite",
"H+",
"HCO3-"},
1462 {1.0, 3.0, 2.0, 1.0},
1471 FAIL() <<
"Missing expected exception.";
1473 catch (
const std::exception & e)
1475 std::string msg(e.what());
1477 msg.find(
"Cannot enforce charge balance using OH- because it is not in the basis") !=
1479 <<
"Failed with unexpected error message: " << msg;
1495 {
"H2O",
"Calcite",
"H+",
"HCO3-"},
1496 {1.0, 3.0, 2.0, 1.0},
1505 FAIL() <<
"Missing expected exception.";
1507 catch (
const std::exception & e)
1509 std::string msg(e.what());
1511 msg.find(
"Cannot enforce charge balance using Calcite because it has zero charge") !=
1513 <<
"Failed with unexpected error message: " << msg;
1520 EXPECT_EQ(_egs_calcite.getChargeBalanceBasisIndex(), (unsigned)1);
1528 EXPECT_EQ(_egs_calcite.getLog10K(123), 1);
1529 FAIL() <<
"Missing expected exception.";
1531 catch (
const std::exception & e)
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;
1538 EXPECT_EQ(_egs_calcite.getLog10K(0), -6.3660);
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]);
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);
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);
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);
1580 EXPECT_EQ(_egs_calcite.getSolventWaterMass(), 1.0);
1582 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
1587 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
1599 {
"H2O",
"Calcite",
"H+",
"HCO3-"},
1600 {2.5, 3.0, 2.0, 1.0},
1609 EXPECT_EQ(egs.getSolventWaterMass(), 2.5);
1615 EXPECT_EQ(_egs_calcite.getBulkMolesOld()[3], 3.0);
1616 EXPECT_EQ(_egs_calcite.getBulkMolesOld()[1], 2.0);
1618 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
1623 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
1635 {
"H2O",
"Calcite",
"H+",
"HCO3-"},
1636 {-0.5, 3.5, 2.5, 1.0},
1645 EXPECT_EQ(egs.getBulkMolesOld()[0], -0.5);
1646 EXPECT_EQ(egs.getBulkMolesOld()[1], 1.0);
1647 EXPECT_EQ(egs.getBulkMolesOld()[2], 1.0);
1648 EXPECT_EQ(egs.getBulkMolesOld()[3], 3.5);
1654 EXPECT_EQ(_egs_calcite.getSolventMassAndFreeMolalityAndMineralMoles()[2], 1.0);
1655 EXPECT_EQ(_egs_calcite.getAlgebraicVariableValues()[0],
1656 _egs_calcite.getSolventMassAndFreeMolalityAndMineralMoles()[1]);
1657 EXPECT_EQ(_egs_calcite.getAlgebraicBasisValues()[0],
1658 _egs_calcite.getSolventMassAndFreeMolalityAndMineralMoles()[1]);
1659 EXPECT_EQ(_egs_calcite.getAlgebraicVariableDenseValues()(0),
1660 _egs_calcite.getSolventMassAndFreeMolalityAndMineralMoles()[1]);
1662 EXPECT_EQ(_egs_kinetic_calcite.getSolventMassAndFreeMolalityAndMineralMoles()[2], 1.0);
1663 EXPECT_EQ(_egs_kinetic_calcite.getAlgebraicVariableValues()[0],
1664 _egs_kinetic_calcite.getSolventMassAndFreeMolalityAndMineralMoles()[1]);
1665 EXPECT_EQ(_egs_kinetic_calcite.getAlgebraicBasisValues()[0],
1666 _egs_kinetic_calcite.getSolventMassAndFreeMolalityAndMineralMoles()[1]);
1667 EXPECT_EQ(_egs_kinetic_calcite.getAlgebraicVariableDenseValues()(0),
1668 _egs_kinetic_calcite.getSolventMassAndFreeMolalityAndMineralMoles()[1]);
1669 EXPECT_EQ(_egs_kinetic_calcite.getAlgebraicVariableValues()[1],
1670 _egs_kinetic_calcite.getSolventMassAndFreeMolalityAndMineralMoles()[3]);
1671 EXPECT_EQ(_egs_kinetic_calcite.getAlgebraicBasisValues()[1],
1672 _egs_kinetic_calcite.getSolventMassAndFreeMolalityAndMineralMoles()[3]);
1673 EXPECT_EQ(_egs_kinetic_calcite.getAlgebraicVariableDenseValues()(1),
1674 _egs_kinetic_calcite.getSolventMassAndFreeMolalityAndMineralMoles()[3]);
1675 EXPECT_EQ(_egs_kinetic_calcite.getAlgebraicVariableValues()[2], 1.1);
1676 EXPECT_EQ(_egs_kinetic_calcite.getAlgebraicVariableDenseValues()(2), 1.1);
1678 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
1683 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
1695 {
"H2O",
"Calcite",
"H+",
"HCO3-"},
1696 {0.5, 3.5, 2.5, 1.0},
1705 EXPECT_EQ(egs.getSolventMassAndFreeMolalityAndMineralMoles()[0], 0.5);
1706 EXPECT_EQ(egs.getSolventMassAndFreeMolalityAndMineralMoles()[2], 1.0);
1707 EXPECT_EQ(egs.getSolventMassAndFreeMolalityAndMineralMoles()[3], 3.5);
1708 EXPECT_EQ(_egs_calcite.getAlgebraicVariableValues()[0],
1709 _egs_calcite.getSolventMassAndFreeMolalityAndMineralMoles()[1]);
1710 EXPECT_EQ(_egs_calcite.getAlgebraicBasisValues()[0],
1711 _egs_calcite.getSolventMassAndFreeMolalityAndMineralMoles()[1]);
1712 EXPECT_EQ(_egs_calcite.getAlgebraicVariableDenseValues()(0),
1713 _egs_calcite.getSolventMassAndFreeMolalityAndMineralMoles()[1]);
1719 EXPECT_EQ(_egs_calcite.getBasisActivityKnown()[0],
true);
1720 EXPECT_EQ(_egs_calcite.getBasisActivityKnown()[1],
false);
1721 EXPECT_EQ(_egs_calcite.getBasisActivityKnown()[2],
false);
1722 EXPECT_EQ(_egs_calcite.getBasisActivityKnown()[3],
true);
1726 database, {
"H2O",
"H+",
"HCO3-",
"O2(aq)"}, {}, {
"O2(g)"}, {}, {}, {},
"O2(aq)",
"e-");
1729 std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
1734 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
1747 {
"H2O",
"H+",
"O2(g)",
"HCO3-"},
1748 {1.0, 2.0, 3.0, 4.0},
1758 EXPECT_EQ(egs.getBasisActivityKnown()[0],
false);
1759 EXPECT_EQ(egs.getBasisActivityKnown()[1],
false);
1760 EXPECT_EQ(egs.getBasisActivityKnown()[2],
false);
1761 EXPECT_EQ(egs.getBasisActivityKnown()[3],
true);
1767 EXPECT_EQ(_egs_calcite.getBasisActivity(0), 1.75);
1768 EXPECT_EQ(_egs_calcite.getBasisActivity(3), 1.0);
1769 EXPECT_EQ(_egs_calcite.getBasisActivity()[0], 1.75);
1770 EXPECT_EQ(_egs_calcite.getBasisActivity()[3], 1.0);
1773 EXPECT_EQ(_egs_calcite.getBasisActivity(4), 1);
1774 FAIL() <<
"Missing expected exception.";
1776 catch (
const std::exception & e)
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;
1786 database, {
"H2O",
"H+",
"HCO3-",
"O2(aq)"}, {}, {
"O2(g)"}, {}, {}, {},
"O2(aq)",
"e-");
1789 std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
1794 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
1807 {
"H2O",
"H+",
"O2(g)",
"HCO3-"},
1808 {1.0, 1.5, 3.0, 2.5},
1818 EXPECT_EQ(egs.getBasisActivity(1), egs.getSolventMassAndFreeMolalityAndMineralMoles()[1]);
1819 EXPECT_EQ(egs.getBasisActivity(2), 2.5);
1820 EXPECT_EQ(egs.getBasisActivity(3), 3.0);
1821 EXPECT_EQ(egs.getBasisActivity()[2], 2.5);
1822 EXPECT_EQ(egs.getBasisActivity()[3], 3.0);
1830 EXPECT_EQ(_egs_calcite.getBasisActivityCoefficient(4), 1);
1831 FAIL() <<
"Missing expected exception.";
1833 catch (
const std::exception & e)
1835 std::string msg(e.what());
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;
1844 database, {
"H2O",
"H+",
"HCO3-",
"O2(aq)"}, {}, {
"O2(g)"}, {}, {}, {},
"O2(aq)",
"e-");
1847 std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
1852 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
1865 {
"H2O",
"H+",
"O2(g)",
"HCO3-"},
1866 {1.0, 1.5, 3.0, 2.5},
1876 EXPECT_EQ(egs.getBasisActivityCoefficient(1), 1.0);
1877 EXPECT_EQ(egs.getBasisActivityCoefficient(2), 1.0);
1878 EXPECT_EQ(egs.getBasisActivityCoefficient()[1], 1.0);
1879 EXPECT_EQ(egs.getBasisActivityCoefficient()[2], 1.0);
1887 EXPECT_EQ(_egs_calcite.getEquilibriumActivityCoefficient(44), 1);
1888 FAIL() <<
"Missing expected exception.";
1890 catch (
const std::exception & e)
1892 std::string msg(e.what());
1893 ASSERT_TRUE(msg.find(
"Cannot retrieve activity coefficient for equilibrium species 44 since " 1895 "equilibrium species") != std::string::npos)
1896 <<
"Failed with unexpected error message: " << msg;
1905 EXPECT_EQ(_egs_calcite.getEquilibriumMolality(44), 1);
1906 FAIL() <<
"Missing expected exception.";
1908 catch (
const std::exception & e)
1910 std::string msg(e.what());
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;
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))),
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))),
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))),
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))),
1946 EXPECT_NEAR(_egs_calcite.getEquilibriumMolality(2) /
1947 (1.0 / _egs_calcite.getEquilibriumActivityCoefficient(2) /
1948 std::pow(10.0, _egs_calcite.getLog10K(2))),
1951 EXPECT_NEAR(_egs_calcite.getEquilibriumMolality()[2] /
1952 (1.0 / _egs_calcite.getEquilibriumActivityCoefficient()[2] /
1953 std::pow(10.0, _egs_calcite.getLog10K(2))),
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))),
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))),
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))),
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))),
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))),
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))),
2009 EXPECT_EQ(eqm_act[
j], 1.0);
2011 EXPECT_NEAR(eqm_act[
j],
2024 EXPECT_EQ(eqm_act[
j], eqm_m[
j] * eqm_ga[
j]);
2033 EXPECT_EQ(_egs_calcite.getEquilibriumActivity(6), 1);
2034 FAIL() <<
"Missing expected exception.";
2036 catch (
const std::exception & e)
2038 std::string msg(e.what());
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;
2051 database, {
"H2O",
"StoiCheckBasis",
"Fe+++"}, {}, {}, {}, {}, {},
"O2(aq)",
"e-");
2053 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
2057 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
2068 {
"H2O",
"Fe+++",
"StoiCheckBasis"},
2078 EXPECT_EQ(egs.getBulkMolesOld()[2], 5.0);
2079 EXPECT_NEAR(egs.getBulkMolesOld()[1],
2082 EXPECT_NEAR(egs.getTotalChargeOld(), 0.0, 1.0E-12);
2088 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
2093 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
2105 {
"H2O",
"Calcite",
"HCO3-",
"H+"},
2106 {175.0, 3.0, 2.0, 1.0},
2116 EXPECT_EQ(egs.getNumInAlgebraicSystem(), (unsigned)2);
2117 EXPECT_EQ(egs.getNumBasisInAlgebraicSystem(), (unsigned)2);
2118 EXPECT_EQ(egs.getNumSurfacePotentials(), (unsigned)0);
2120 const DenseVector<Real> mole_additions(4);
2124 EXPECT_EQ(egs.getResidualComponent(3, mole_additions), 2);
2125 FAIL() <<
"Missing expected exception.";
2127 catch (
const std::exception & e)
2129 std::string msg(e.what());
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") !=
2135 <<
"Failed with unexpected error message: " << msg;
2140 const DenseVector<Real> bad(5);
2141 EXPECT_EQ(egs.getResidualComponent(0, bad), 2);
2142 FAIL() <<
"Missing expected exception.";
2144 catch (
const std::exception & e)
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;
2152 const Real nw = egs.getSolventMassAndFreeMolalityAndMineralMoles()[0];
2153 const Real mHCO3 = egs.getSolventMassAndFreeMolalityAndMineralMoles()[2];
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);
2160 res = -egs.getBulkMolesOld()[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);
2169 _db_calcite, {
"H2O",
"H+",
"StoiCheckBasis"}, {}, {}, {}, {}, {},
"O2(aq)",
"e-");
2171 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum>
cm2 = {
2175 std::vector<GeochemistryUnitConverter::GeochemistryUnit>
cu2;
2186 {
"H2O",
"H+",
"StoiCheckBasis"},
2196 DenseVector<Real> mole_add3(3);
2198 EXPECT_EQ(egs2.getNumInAlgebraicSystem(), (unsigned)2);
2199 EXPECT_EQ(egs2.getNumBasisInAlgebraicSystem(), (unsigned)2);
2200 EXPECT_EQ(egs2.getNumSurfacePotentials(), (unsigned)0);
2202 const Real nw2 = egs2.getSolventMassAndFreeMolalityAndMineralMoles()[0];
2203 const Real mscb = egs2.getSolventMassAndFreeMolalityAndMineralMoles()[2];
2204 const Real mscr = egs2.getEquilibriumMolality(0);
2207 EXPECT_NEAR(egs2.getResidualComponent(0, mole_add3) / resh2o, 1.0, 1.0E-12);
2209 const Real resscb = (1.0 / 2.5) * egs2.getBulkMolesOld()[1] +
2210 nw2 * (mscb + 1.5 * mscr);
2211 EXPECT_NEAR(egs2.getResidualComponent(1, mole_add3) / resscb, 1.0, 1.0E-12);
2213 for (
unsigned i = 0; i < 3; ++i)
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);
2240 {
"H2O",
"Calcite",
"H+",
"HCO3-"},
2241 {1.75, 3.0, 2.0, 1.0},
2251 DenseVector<Real> var_incorrect_size(123);
2254 egs.setAlgebraicVariables(var_incorrect_size);
2255 FAIL() <<
"Missing expected exception.";
2257 catch (
const std::exception & e)
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;
2264 DenseVector<Real> var_neg(egs.getNumInAlgebraicSystem());
2268 egs.setAlgebraicVariables(var_neg);
2269 FAIL() <<
"Missing expected exception.";
2271 catch (
const std::exception & e)
2273 std::string msg(e.what());
2274 ASSERT_TRUE(msg.find(
"Cannot set algebraic variables to non-positive values such as -1.25") !=
2276 <<
"Failed with unexpected error message: " << msg;
2287 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
2292 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
2305 {
"H2O",
"Calcite",
"HCO3-",
"H+"},
2306 {175.0, 3.0, 3.2E-2, 1E-8},
2316 DenseVector<Real> mole_add(4);
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();
2330 egs.computeJacobian(res_orig, jac, mole_add, dmole_add);
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]);
2339 for (
unsigned var_num = 0; var_num < 2; ++var_num)
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)
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);
2353 for (
unsigned i = 0; i < 4; ++i)
2354 mole_add(i) = i + 1;
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);
2362 egs.computeJacobian(res_orig, jac, mole_add, dmole_add);
2363 for (
unsigned var_num = 0; var_num < 2; ++var_num)
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)
2371 mole_add(i) = i + 1 + dmole_add(i, 0) *
eps;
2373 mole_add(i) = i + 1 + dmole_add(i, 2) *
eps;
2375 for (
unsigned res_comp = 0; res_comp < 2; ++res_comp)
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);
2391 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
2396 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
2408 {
"H2O",
"Calcite",
"HCO3-",
"H+"},
2409 {1.0, 3.0, 3.2E-4, 1E-4},
2418 DenseVector<Real> mole_add(4);
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();
2430 egs.computeJacobian(res_orig, jac, mole_add, dmole_add);
2433 for (
unsigned var_num = 0; var_num < 2; ++var_num)
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)
2440 const Real expected =
2441 (egs.getResidualComponent(res_comp, mole_add) - res_orig(res_comp)) /
eps;
2443 EXPECT_NEAR(jac(res_comp, var_num) / expected,
2447 EXPECT_NEAR(jac(res_comp, var_num) / expected, 1.0, 1.0E-4);
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;
2460 egs.computeJacobian(res_orig, jac, mole_add, dmole_add);
2461 for (
unsigned var_num = 0; var_num < 2; ++var_num)
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)
2469 mole_add(i) = (i + 1) * 1E4 + dmole_add(i, 1) *
eps;
2471 mole_add(i) = (i + 1) * 1E4 + dmole_add(i, 2) *
eps;
2473 for (
unsigned res_comp = 0; res_comp < 2; ++res_comp)
2475 const Real expected =
2476 (egs.getResidualComponent(res_comp, mole_add) - res_orig(res_comp)) /
eps;
2478 EXPECT_NEAR(jac(res_comp, var_num) / expected,
2482 EXPECT_NEAR(jac(res_comp, var_num) / expected, 1.0, 1.0E-4);
2494 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
2499 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
2511 {
"H2O",
"Calcite",
"H+",
"HCO3-"},
2512 {1.0, 3.0, 3.2E-4, 1E-4},
2521 DenseVector<Real> mole_add(4);
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();
2532 egs.computeJacobian(res_orig, jac, mole_add, dmole_add);
2535 std::vector<Real> var_new = var_orig;
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);
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;
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;
2554 expected = (egs.getResidualComponent(0, mole_add) - res_orig(0)) /
eps;
2555 EXPECT_NEAR(jac(0, 0) / expected, 1.0, 1.0E-7);
2566 _db_calcite, {
"H2O",
"H+",
"Ca++",
"HCO3-"}, {
"Calcite"}, {}, {}, {}, {},
"O2(aq)",
"e-");
2568 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
2573 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
2585 {
"H2O",
"H+",
"Ca++",
"HCO3-"},
2586 {175.0, 1E1, 4E1, 1E-4},
2595 DenseVector<Real> mole_add(4);
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();
2607 egs.computeJacobian(res_orig, jac, mole_add, dmole_add);
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]);
2615 for (
unsigned var_num = 0; var_num < 3; ++var_num)
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)
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);
2627 EXPECT_NEAR(expected / jac(res_comp, var_num), 1.0, 1.0E-6);
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);
2641 egs.computeJacobian(res_orig, jac, mole_add, dmole_add);
2642 for (
unsigned var_num = 0; var_num < 3; ++var_num)
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)
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);
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;
2681 {
"H2O",
"Ca++",
"H+",
"HCO3-"},
2682 {1.11, 3.0, 2.0, 1.5},
2693 if (name_index.first ==
"Calcite")
2695 si[name_index.second],
2696 std::log10(egs.getBasisActivity(3) * egs.getBasisActivity(2) / egs.getBasisActivity(1)) -
2700 EXPECT_EQ(si[name_index.second], 0.0);
2710 FAIL() <<
"Missing expected exception.";
2712 catch (
const std::exception & e)
2714 std::string msg(e.what());
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;
2726 FAIL() <<
"Missing expected exception.";
2728 catch (
const std::exception & e)
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;
2738 database, {
"H2O",
"H+",
"HCO3-",
"O2(aq)"}, {}, {
"O2(g)"}, {}, {}, {},
"O2(aq)",
"e-");
2741 std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
2746 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
2759 {
"H2O",
"H+",
"O2(aq)",
"HCO3-"},
2760 {1.0, 1.5, 3.0, 2.5},
2773 FAIL() <<
"Missing expected exception.";
2775 catch (
const std::exception & e)
2777 std::string msg(e.what());
2778 ASSERT_TRUE(msg.find(
"GeochemicalSystem: attempting to swap a gas into the basis") !=
2780 <<
"Failed with unexpected error message: " << msg;
2783 std::vector<GeochemicalSystem::ConstraintUserMeaningEnum>
cm2;
2788 std::vector<GeochemistryUnitConverter::GeochemistryUnit>
cu2;
2800 {
"H2O",
"H+",
"O2(g)",
"HCO3-"},
2801 {1.0, 1.5, 3.0, 2.5},
2814 FAIL() <<
"Missing expected exception.";
2816 catch (
const std::exception & e)
2818 std::string msg(e.what());
2819 ASSERT_TRUE(msg.find(
"GeochemicalSystem: attempting to swap a gas out of the basis") !=
2821 <<
"Failed with unexpected error message: " << msg;
2828 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
2833 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
2846 {
"H2O",
"Calcite",
"H+",
"HCO3-"},
2847 {1.75, 3.0, 2.0, 1.0},
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();
2867 egs.performSwap(3, 5);
2869 EXPECT_EQ(egs.getNumInAlgebraicSystem(), (unsigned)2);
2870 EXPECT_EQ(egs.getNumBasisInAlgebraicSystem(),
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);
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);
2887 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
2892 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
2907 {
"H2O",
"Ca++",
"H+",
"HCO3-"},
2908 {1.75, 3.0E-1, 2.0E-1, 1.0E-1},
2917 EXPECT_EQ(egs_small_IS.getIonicStrength(), 0.0078125);
2918 EXPECT_EQ(egs_small_IS.getStoichiometricIonicStrength(), 0.0078125);
2930 {
"H2O",
"Ca++",
"H+",
"HCO3-"},
2931 {1.75, 1E-10, 1E-10, 1.0},
2940 EXPECT_NEAR(egs.getIonicStrength(), 0.5, 1.0E-8);
2941 EXPECT_NEAR(egs.getStoichiometricIonicStrength(), 0.5, 1.0E-8);
2949 database, {
"H2O",
"StoiCheckBasis",
"HCO3-"}, {}, {}, {}, {}, {},
"O2(aq)",
"e-");
2951 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
2955 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
2966 {
"H2O",
"StoiCheckBasis",
"HCO3-"},
2977 EXPECT_NEAR(egs.getTotalChargeOld(), 0.0, 1.0E-12);
2978 EXPECT_NEAR(egs.getBulkMolesOld()[1],
2981 EXPECT_EQ(egs.getBulkMolesOld()[2], 1.0);
2982 EXPECT_EQ(egs.getChargeBalanceBasisIndex(), (unsigned)1);
2986 egs.alterChargeBalanceSpecies(0.3);
2987 EXPECT_EQ(egs.getChargeBalanceBasisIndex(), (unsigned)2);
2988 EXPECT_EQ(egs.getBulkMolesOld()[1], 5.0);
2989 EXPECT_EQ(egs.getBulkMolesOld()[2],
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);
3005 EXPECT_EQ(_egs_redox.getOriginalRedoxLHS(),
"e-");
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;
3019 _egs_redox.getRedoxLog10K(ophth_slot), -boa * 542.8292 + 20.7757 - 0.25 * (-2.8990), 1E-8);
3022 _egs_redox.getRedoxLog10K(ch4_slot), -boa * 144.1080 + 20.7757 - 0.25 * (-2.8990), 1E-8);
3025 _egs_redox.getRedoxLog10K(3);
3026 FAIL() <<
"Missing expected exception.";
3028 catch (
const std::exception & e)
3030 std::string msg(e.what());
3033 "Cannot retrieve log10K for redox species 3 since there are only 3 redox species") !=
3035 <<
"Failed with unexpected error message: " << msg;
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);
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);
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);
3060 _egs_redox.log10RedoxActivityProduct(3);
3061 FAIL() <<
"Missing expected exception.";
3063 catch (
const std::exception & e)
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;
3077 _egs_calcite.getSurfacePotential(0);
3078 FAIL() <<
"Missing expected exception.";
3080 catch (
const std::exception & e)
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;
3094 _egs_calcite.getSurfaceCharge(0);
3095 FAIL() <<
"Missing expected exception.";
3097 catch (
const std::exception & e)
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;
3111 {
"H2O",
"H+",
">(s)FeOH",
">(w)FeOH",
"Fe+++"},
3120 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
3126 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
3139 {
"H2O",
"H+",
">(s)FeOH",
">(w)FeOH",
"Fe(OH)3(ppd)"},
3140 {1.75, 1.0, 2.0, 1.0, 1.5},
3149 EXPECT_EQ(egs.getSorbingSurfaceArea().size(), (std::size_t)1);
3150 EXPECT_EQ(egs.getSorbingSurfaceArea()[0], 1.5 * 106.8689 * 600.0);
3158 {
"H2O",
"H+",
">(s)FeOH",
">(w)FeOH",
"Fe+++",
"HCO3-"},
3167 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
3174 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
3181 const Real temp = 45.0;
3189 {
"H2O",
"H+",
">(s)FeOH",
">(w)FeOH",
"Fe+++",
"HCO3-"},
3190 {1.75, 1.0, 2.0, 3.0, 1.0, 1.0},
3199 const DenseVector<Real> mole_add(6);
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],
3208 EXPECT_EQ(egs.getBasisIndexOfAlgebraicSystem()[1],
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);
3224 egs.setAlgebraicVariables(alg_vars);
3225 FAIL() <<
"Missing expected exception.";
3227 catch (
const std::exception & e)
3229 std::string msg(e.what());
3230 ASSERT_TRUE(msg.find(
"Cannot set algebraic variables to non-positive values such as -1") !=
3232 <<
"Failed with unexpected error message: " << msg;
3237 egs.setAlgebraicVariables(alg_vars);
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));
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);
3250 egs.getSurfacePotential(1);
3251 FAIL() <<
"Missing expected exception.";
3253 catch (
const std::exception & e)
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;
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);
3270 egs.getSurfaceCharge(1);
3271 FAIL() <<
"Missing expected exception.";
3273 catch (
const std::exception & e)
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;
3282 const Real act_h = egs.getBasisActivity(1);
3283 const Real act_bicarb = egs.getBasisActivity(2);
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)) *
3293 egs.getResidualComponent(2, mole_add), -pref * 600.0 + 1.75 * (-1.0) * mol_surf, 1.0E-9);
3296 egs.getResidualComponent(3, mole_add);
3297 FAIL() <<
"Missing expected exception.";
3299 catch (
const std::exception & e)
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") !=
3305 <<
"Failed with unexpected error message: " << msg;
3314 {
"H2O",
"H+",
">(s)FeOH",
">(w)FeOH",
"Fe+++",
"HCO3-"},
3323 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
3330 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
3337 const Real temp = 45.0;
3347 {
"H2O",
"H+",
">(s)FeOH",
">(w)FeOH",
"Fe+++",
"HCO3-"},
3348 {1.75, 1.0, 2.0, 3.0, 0.5, 1.0},
3357 const DenseVector<Real> mole_add(6);
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();
3366 egs.computeJacobian(res_orig, jac, mole_add, dmole_add);
3369 for (
unsigned var_num = 0; var_num < num_alg; ++var_num)
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)
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);
3381 EXPECT_NEAR(expected / jac(res_comp, var_num), 1.0, 1.0E-2);
3391 {
"H2O",
"H+",
">(s)FeOH",
">(w)FeOH",
"Fe+++",
"HCO3-"},
3394 {
"Something",
"Fe(OH)3(ppd)fake"},
3400 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
3407 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
3414 std::vector<GeochemistryUnitConverter::GeochemistryUnit> ku;
3418 const Real temp = 45.0;
3428 {
"H2O",
"H+",
">(s)FeOH",
">(w)FeOH",
"Fe+++",
"HCO3-"},
3429 {1.75, 1.0, 2.0, 3.0, 0.5, 1.0},
3435 {
"Something",
"Fe(OH)3(ppd)fake"},
3438 const DenseVector<Real> mole_add(8);
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();
3447 egs.computeJacobian(res_orig, jac, mole_add, dmole_add);
3450 for (
unsigned var_num = 0; var_num < num_alg; ++var_num)
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)
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);
3462 EXPECT_NEAR(expected / jac(res_comp, var_num), 1.0, 1.0E-2);
3472 {
"H2O",
"H+",
">(s)FeOH",
">(w)FeOH",
"Fe+++",
"HCO3-"},
3475 {
"Something",
"Fe(OH)3(ppd)fake"},
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;
3508 {
"H2O",
"H+",
">(s)FeOH",
">(w)FeOH",
"Fe+++",
"HCO3-"},
3509 {1.75, 1.0, 2.0, 3.0, 0.5, 1.0},
3515 {
"Something",
"Fe(OH)3(ppd)fake"},
3518 DenseVector<Real> mole_add(8);
3519 for (
unsigned i = 0; i < 8; ++i)
3520 mole_add(i) = 1.0 + i;
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);
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();
3533 egs.computeJacobian(res_orig, jac, mole_add, dmole_add);
3536 for (
unsigned var_num = 0; var_num < num_alg; ++var_num)
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)
3545 mole_add(i) = 1.0 + i + dmole_add(i, var_num) * del;
3546 else if (var_num == 3)
3547 mole_add(i) = 1.0 + i + dmole_add(i, 5) * del;
3548 else if (var_num == 4)
3549 mole_add(i) = 1.0 + i;
3550 else if (var_num == 5 || var_num == 6)
3551 mole_add(i) = 1.0 + i + dmole_add(i, var_num + 1) * del;
3554 for (
unsigned res_comp = 0; res_comp < num_alg; ++res_comp)
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);
3561 EXPECT_NEAR(expected / jac(res_comp, var_num), 1.0, 1.0E-2);
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.";
3588 catch (
const std::exception & e)
3590 std::string msg(e.what());
3591 ASSERT_TRUE(msg.find(
"When setting all molalities, names and values must have same size") !=
3593 <<
"Failed with unexpected error message: " << msg;
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.";
3604 catch (
const std::exception & e)
3606 std::string msg(e.what());
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;
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.";
3622 catch (
const std::exception & e)
3624 std::string msg(e.what());
3626 msg.find(
"constraints_from_molalities has size 3 while number of basis species is 4") !=
3628 <<
"Failed with unexpected error message: " << msg;
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.";
3639 catch (
const std::exception & e)
3641 std::string msg(e.what());
3642 ASSERT_TRUE(msg.find(
"Molality for mineral Calcite cannot be -1: it must be non-negative") !=
3644 <<
"Failed with unexpected error message: " << msg;
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.";
3655 catch (
const std::exception & e)
3657 std::string msg(e.what());
3658 ASSERT_TRUE(msg.find(
"Molality for species Ca++ cannot be -1: it must be non-negative") !=
3660 <<
"Failed with unexpected error message: " << msg;
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.";
3671 catch (
const std::exception & e)
3673 std::string msg(e.what());
3674 ASSERT_TRUE(msg.find(
"Molality for species H+ cannot be 0: it must be positive") !=
3676 <<
"Failed with unexpected error message: " << msg;
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.";
3687 catch (
const std::exception & e)
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") !=
3693 <<
"Failed with unexpected error message: " << msg;
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.";
3704 catch (
const std::exception & e)
3706 std::string msg(e.what());
3707 ASSERT_TRUE(msg.find(
"Molality for species CO2(aq) cannot be -1: it must be non-negative") !=
3709 <<
"Failed with unexpected error message: " << msg;
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.";
3720 catch (
const std::exception & e)
3722 std::string msg(e.what());
3724 msg.find(
"Molality for species CO3-- needs to be provided when setting all molalities") !=
3726 <<
"Failed with unexpected error message: " << msg;
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.";
3737 catch (
const std::exception & e)
3739 std::string msg(e.what());
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") !=
3744 <<
"Failed with unexpected error message: " << msg;
3753 {
"H2O",
"H+",
"HCO3-",
"O2(aq)",
"Fe+++",
">(s)FeOH",
">(w)FeOH"},
3762 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
3770 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
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},
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.";
3818 catch (
const std::exception & e)
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;
3828 nonconst.setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(
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.";
3849 catch (
const std::exception & e)
3851 std::string msg(e.what());
3852 ASSERT_TRUE(msg.find(
"Molality for gas O2(g) cannot be 1: it must be zero") !=
3854 <<
"Failed with unexpected error message: " << msg;
3859 nonconst.setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(
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.";
3880 catch (
const std::exception & e)
3882 std::string msg(e.what());
3883 ASSERT_TRUE(msg.find(
"Molality for species O2(aq) cannot be -1: it must be non-negative") !=
3885 <<
"Failed with unexpected error message: " << msg;
3890 nonconst.setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(
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.";
3911 catch (
const std::exception & e)
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;
3921 nonconst.setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(
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.";
3942 catch (
const std::exception & e)
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;
3952 nonconst.setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(
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.";
3973 catch (
const std::exception & e)
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") !=
3979 <<
"Failed with unexpected error message: " << msg;
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.";
3995 catch (
const std::exception & e)
3997 std::string msg(e.what());
3999 msg.find(
"Moles for species Calcite needs to be provided when setting all molalities") !=
4001 <<
"Failed with unexpected error message: " << msg;
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.";
4012 catch (
const std::exception & e)
4014 std::string msg(e.what());
4015 ASSERT_TRUE(msg.find(
"Mole number for kinetic species must be positive, not 0") !=
4017 <<
"Failed with unexpected error message: " << msg;
4026 {
"H2O",
"H+",
"HCO3-",
"O2(aq)",
"Fe+++",
">(s)FeOH",
">(w)FeOH"},
4035 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
4043 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
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},
4069 const std::vector<Real> molal_before_set =
4070 nonconst.getSolventMassAndFreeMolalityAndMineralMoles();
4072 const std::vector<std::string> gold_name = {
"H2O",
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});
4093 const unsigned num_basis = nonconst.getNumInBasis();
4094 const unsigned num_eqm = nonconst.getNumInEquilibrium();
4097 std::vector<Real> gold_molal = set_molal;
4098 gold_molal[6] = molal_before_set[6];
4100 gold_molal[13] = 0.0;
4101 const std::vector<Real> molal = nonconst.getSolventMassAndFreeMolalityAndMineralMoles();
4102 for (
unsigned i = 0; i < num_basis; ++i)
4104 for (
unsigned j = num_basis;
j < num_basis + num_eqm; ++
j)
4106 EXPECT_NEAR(nonconst.getSurfacePotential(0),
4111 EXPECT_EQ(nonconst.getSolventWaterMass(), 1.1);
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),
4120 EXPECT_EQ(nonconst.getBasisActivity(3), 4.0);
4123 EXPECT_EQ(nonconst.getSorbingSurfaceArea()[0], 60.0);
4126 const std::vector<Real> bulk = nonconst.getBulkMolesOld();
4127 for (
unsigned i = 0; i < num_basis; ++i)
4134 for (
unsigned j = 0;
j < num_eqm; ++
j)
4138 EXPECT_EQ(bulk[4], bulk_before_set[4]);
4149 {
"H2O",
"H+",
"HCO3-",
"O2(aq)",
"Fe+++",
">(s)FeOH",
">(w)FeOH"},
4158 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
4166 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
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},
4192 const std::vector<Real> molal_before_set =
4193 nonconst.getSolventMassAndFreeMolalityAndMineralMoles();
4195 const std::vector<std::string> gold_name = {
"H2O",
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});
4216 const unsigned num_basis = nonconst.getNumInBasis();
4217 const unsigned num_eqm = nonconst.getNumInEquilibrium();
4220 std::vector<Real> gold_molal = set_molal;
4222 molal_before_set[1];
4223 gold_molal[6] = molal_before_set[6];
4225 gold_molal[13] = 0.0;
4226 const std::vector<Real> molal = nonconst.getSolventMassAndFreeMolalityAndMineralMoles();
4227 for (
unsigned i = 0; i < num_basis; ++i)
4229 for (
unsigned j = num_basis;
j < num_basis + num_eqm; ++
j)
4231 EXPECT_NEAR(nonconst.getSurfacePotential(0),
4236 EXPECT_EQ(nonconst.getSolventWaterMass(), 1.1);
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),
4245 EXPECT_EQ(nonconst.getBasisActivity(3), 4.0);
4248 EXPECT_EQ(nonconst.getSorbingSurfaceArea()[0], 60.0);
4251 const std::vector<Real> bulk = nonconst.getBulkMolesOld();
4252 for (
unsigned i = 0; i < num_basis; ++i)
4259 for (
unsigned j = 0;
j < num_eqm; ++
j)
4262 if (i == 0 || i == 4)
4263 EXPECT_EQ(bulk[i], bulk_before_set[i]);
4273 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
4278 const std::vector<GeochemicalSystem::ConstraintMeaningEnum> cim = {
4283 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
4295 {
"H2O",
"H+",
"HCO3-",
"Calcite"},
4296 {1.75, 3.0, 2.0, 1.0},
4305 for (
unsigned i = 0; i < 4; ++i)
4306 EXPECT_EQ(egs.getConstraintMeaning()[i], cim[i]);
4314 database, {
"H2O",
"H+",
"HCO3-",
"O2(aq)"}, {}, {
"O2(g)"}, {}, {}, {},
"O2(aq)",
"e-");
4317 std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
4322 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
4335 {
"H2O",
"H+",
"O2(g)",
"HCO3-"},
4336 {1.0, 2.0, 3.0, 4.0},
4349 FAIL() <<
"Missing expected exception.";
4351 catch (
const std::exception & e)
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;
4361 egs.changeConstraintToBulk(4, 1.0);
4362 FAIL() <<
"Missing expected exception.";
4364 catch (
const std::exception & e)
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;
4374 egs.changeConstraintToBulk(3, 1.0);
4375 FAIL() <<
"Missing expected exception.";
4377 catch (
const std::exception & e)
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;
4393 database, {
"H2O",
"H+",
"HCO3-",
"O2(aq)"}, {}, {
"O2(g)"}, {}, {}, {},
"O2(aq)",
"e-");
4396 std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
4401 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
4414 {
"H2O",
"H+",
"HCO3-",
"O2(g)"},
4415 {1.0, 2.0, 3.0, 4.0},
4425 const std::vector<GeochemicalSystem::ConstraintMeaningEnum> cim = {
4430 for (
unsigned i = 0; i < 4; ++i)
4431 EXPECT_EQ(egs.getConstraintMeaning()[i], cim[i]);
4439 database, {
"H2O",
"H+",
"HCO3-",
"O2(aq)"}, {}, {
"O2(g)"}, {}, {}, {},
"O2(aq)",
"e-");
4442 std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
4447 std::vector<GeochemicalSystem::ConstraintMeaningEnum> cim;
4452 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
4465 {
"H2O",
"H+",
"HCO3-",
"O2(g)"},
4466 {1.0, 2.0, 3.0, 4.0},
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]);
4484 std::vector<Real> new_bulk = orig_bulk;
4485 egs.changeConstraintToBulk(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]);
4492 egs.changeConstraintToBulk(2);
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]);
4500 egs.changeConstraintToBulk(3);
4502 for (
unsigned i = 0; i < 4; ++i)
4503 EXPECT_EQ(egs.getConstraintMeaning()[i], cim[i]);
4512 database, {
"H2O",
"H+",
"HCO3-",
"O2(aq)"}, {}, {
"O2(g)"}, {}, {}, {},
"O2(aq)",
"e-");
4515 std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
4520 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
4533 {
"H2O",
"H+",
"HCO3-",
"O2(g)"},
4534 {1.0, 2.0, 3.0, 4.0},
4547 FAIL() <<
"Missing expected exception.";
4549 catch (
const std::exception & e)
4551 std::string msg(e.what());
4553 msg.find(
"Cannot addToBulkMoles for species 4 because there are only 4 basis species") !=
4555 <<
"Failed with unexpected error message: " << msg;
4564 database, {
"H2O",
"H+",
"O2(aq)",
"HCO3-"}, {}, {
"O2(g)"}, {}, {}, {},
"O2(aq)",
"e-");
4567 std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
4572 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
4578 std::vector<Real> bm = {1.0, 2.0, 3.0, 4.0};
4587 {
"H2O",
"H+",
"O2(aq)",
"HCO3-"},
4600 for (
unsigned i = 0; i < 3; ++i)
4601 EXPECT_EQ(egs.getBulkMolesOld()[i], bm[i]);
4602 egs.addToBulkMoles(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);
4608 for (
unsigned i = 0; i < 3; ++i)
4609 EXPECT_EQ(egs.getBulkMolesOld()[i], bm[i]);
4619 FAIL() <<
"Missing expected exception.";
4621 catch (
const std::exception & e)
4623 std::string msg(e.what());
4626 "Cannot setConstraintValue for species 4 because there are only 4 basis species") !=
4628 <<
"Failed with unexpected error message: " << msg;
4637 database, {
"H2O",
"H+",
"O2(aq)",
"HCO3-"}, {}, {
"O2(g)"}, {}, {}, {},
"O2(aq)",
"e-");
4640 std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm;
4645 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
4651 std::vector<Real> bm_in = {1.0, 2.0, 3.0, 4.0};
4660 {
"H2O",
"H+",
"O2(aq)",
"HCO3-"},
4672 std::vector<Real> mol = egs.getSolventMassAndFreeMolalityAndMineralMoles();
4674 egs.setConstraintValue(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]);
4681 egs.setConstraintValue(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]);
4688 egs.setConstraintValue(3, 4.4);
4690 bm = egs.getBulkMolesOld();
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]);
4696 egs.setConstraintValue(2, 3.3);
4697 bm = egs.getBulkMolesOld();
4699 for (
unsigned i = 0; i < 4; ++i)
4700 EXPECT_EQ(egs.getBulkMolesOld()[i], bm[i]);
4701 EXPECT_EQ(egs.getBasisActivity(2), 3.3);
4704 egs.changeConstraintToBulk(3);
4706 bm = egs.getBulkMolesOld();
4707 egs.setConstraintValue(3, 7.7);
4710 for (
unsigned i = 0; i < 4; ++i)
4711 EXPECT_NEAR(egs.getBulkMolesOld()[i], bm[i], 1.0E-8);
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);
4728 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
4731 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
4753 EXPECT_EQ(mgd_ref.basis_species_name.size(), (std::size_t)2);
4763 db, {
"H2O",
"H+",
"Ca++"}, {}, {}, {}, {}, {},
"O2(aq)",
"e-");
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);
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);
4791 {
"H2O",
"H+",
"HCO3-",
"O2(aq)",
"Fe+++",
">(s)FeOH",
">(w)FeOH",
"Ca++"},
4792 {
"Goethite",
"Calcite"},
4800 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
4809 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
4818 std::vector<GeochemistryUnitConverter::GeochemistryUnit> ku;
4826 {
"O2(g)",
"Calcite"},
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},
4839 const unsigned num_eqm = nonconst.getNumInEquilibrium();
4840 const unsigned num_kin = nonconst.getNumKinetic();
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",
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});
4858 nonconst.setMineralRelatedFreeMoles(1234.5);
4860 std::vector<Real> gold_molal = set_molal;
4861 gold_molal[5] = 1234.5;
4862 gold_molal[6] = 1234.5;
4863 gold_molal[7] = 1234.5;
4864 gold_molal[15] = 1234.5;
4865 gold_molal[16] = 1234.5;
4866 gold_molal[17] = 0.0;
4867 gold_molal[21] = 1234.5;
4869 const std::vector<Real> & basis_mol = nonconst.getSolventMassAndFreeMolalityAndMineralMoles();
4870 for (
unsigned i = 0; i < num_basis; ++i)
4872 const std::vector<Real> & eqm_mol = nonconst.getEquilibriumMolality();
4873 for (
unsigned j = num_basis;
j < num_basis + num_eqm; ++
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;
4884 std::vector<GeochemistryUnitConverter::GeochemistryUnit> ku;
4895 {
"H2O",
"Ca++",
"H+",
"HCO3-"},
4896 {1.75, 3.0, 2.0, 1.0},
4905 FAIL() <<
"Missing expected exception.";
4907 catch (
const std::exception & e)
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;
4924 {
"H2O",
"Ca++",
"H+",
"HCO3-"},
4925 {1.75, 3.0, 2.0, 1.0},
4934 FAIL() <<
"Missing expected exception.";
4936 catch (
const std::exception & e)
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;
4953 {
"H2O",
"Ca++",
"H+",
"HCO3-"},
4954 {1.75, 3.0, 2.0, 1.0},
4963 FAIL() <<
"Missing expected exception.";
4965 catch (
const std::exception & e)
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;
4982 {
"H2O",
"Ca++",
"H+",
"HCO3-"},
4983 {1.75, 3.0, 2.0, 1.0},
4992 FAIL() <<
"Missing expected exception.";
4994 catch (
const std::exception & e)
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;
5011 {
"H2O",
"Ca++",
"H+",
"HCO3-"},
5012 {1.75, 3.0, 2.0, 1.0},
5021 FAIL() <<
"Missing expected exception.";
5023 catch (
const std::exception & e)
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;
5031 std::vector<GeochemistryUnitConverter::GeochemistryUnit> ku_bad;
5042 {
"H2O",
"Ca++",
"H+",
"HCO3-"},
5043 {1.75, 3.0, 2.0, 1.0},
5052 FAIL() <<
"Missing expected exception.";
5054 catch (
const std::exception & e)
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;
5068 _egs_kinetic_calcite.getKineticLog10K(1);
5069 FAIL() <<
"Missing expected exception.";
5071 catch (
const std::exception & e)
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;
5085 _egs_kinetic_calcite.log10KineticActivityProduct(1);
5086 FAIL() <<
"Missing expected exception.";
5088 catch (
const std::exception & e)
5090 std::string msg(e.what());
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;
5103 _egs_kinetic_calcite.getKineticMoles(1);
5104 FAIL() <<
"Missing expected exception.";
5106 catch (
const std::exception & e)
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;
5123 FAIL() <<
"Missing expected exception.";
5125 catch (
const std::exception & e)
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;
5136 FAIL() <<
"Missing expected exception.";
5138 catch (
const std::exception & e)
5140 std::string msg(e.what());
5141 ASSERT_TRUE(msg.find(
"Mole number for kinetic species must be positive, not 0") !=
5143 <<
"Failed with unexpected error message: " << msg;
5150 EXPECT_EQ(_egs_kinetic_calcite.getNumKinetic(), (unsigned)1);
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);
5164 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
5169 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
5174 std::vector<GeochemistryUnitConverter::GeochemistryUnit> ku;
5183 {
"H2O",
"H+",
"HCO3-",
"Ca++"},
5184 {1.75, 3.0, 2.0, 1.0},
5195 gs.log10KineticActivityProduct(0), std::log10(gs.getBasisActivity(3) / 3.0 * 2.0), 1.0E-6);
5201 std::vector<GeochemistryUnitConverter::GeochemistryUnit> ku;
5210 {
"H2O",
"H+",
"HCO3-",
"Ca++"},
5211 {1.75, 3.0, 2.0, 1.0},
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);
5228 std::vector<GeochemistryUnitConverter::GeochemistryUnit> kucm3;
5237 {
"H2O",
"H+",
"HCO3-",
"Ca++"},
5238 {1.75, 3.0, 2.0, 1.0},
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);
5260 _db_calcite, {
"H2O",
"H+",
"HCO3-",
"Ca++"}, {}, {}, {}, {}, {},
"O2(aq)",
"e-");
5263 mgd_without_calcite,
5270 {
"H2O",
"Ca++",
"H+",
"HCO3-"},
5274 {1.75, 3.0 + 1.1, 2.0 - 1.1, 1.0},
5283 EXPECT_EQ(_egs_kinetic_calcite.getBulkMolesOld()[0],
5284 egs_without_kinetic.getBulkMolesOld()[0]);
5285 EXPECT_EQ(_egs_kinetic_calcite.getBulkMolesOld()[1],
5286 egs_without_kinetic.getBulkMolesOld()[1]);
5287 EXPECT_EQ(_egs_kinetic_calcite.getBulkMolesOld()[2],
5288 egs_without_kinetic.getBulkMolesOld()[2] + 1.1);
5289 EXPECT_EQ(_egs_kinetic_calcite.getBulkMolesOld()[3],
5290 egs_without_kinetic.getBulkMolesOld()[3]);
5296 std::vector<GeochemistryUnitConverter::GeochemistryUnit> ku;
5305 {
"H2O",
"H+",
"HCO3-",
"Ca++"},
5306 {1.75, 3.0, 2.0, 1.0},
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);
5326 {
"H2O",
"H+",
"HCO3-",
"Ca++"},
5335 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
5340 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
5345 std::vector<GeochemistryUnitConverter::GeochemistryUnit> ku;
5354 {
"H2O",
"H+",
"HCO3-",
"Calcite"},
5355 {1.75, 3.0, 2.0, 10.0},
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],
5376 {
"H2O",
"H+",
"HCO3-",
"Ca++"},
5385 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
5390 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
5395 std::vector<GeochemistryUnitConverter::GeochemistryUnit> ku;
5404 {
"H2O",
"H+",
"HCO3-",
"Calcite"},
5405 {1.75, 3.0, 2.0, 10.0},
5414 const std::vector<std::string> gold_name = {
"H2O",
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});
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)
5435 for (
unsigned j = num_basis;
j < num_basis + num_eqm; ++
j)
5437 for (
unsigned k = num_basis + num_eqm;
k < num_basis + num_eqm + num_kin; ++
k)
5448 const DenseVector<Real> res(4);
5449 const DenseVector<Real> mole_additions(5);
5451 _egs_kinetic_calcite.computeJacobian(res, jac, mole_additions, dmole_additions);
5452 FAIL() <<
"Missing expected exception.";
5454 catch (
const std::exception & e)
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;
5463 const DenseVector<Real> res(3);
5464 const DenseVector<Real> mole_additions(4);
5466 _egs_kinetic_calcite.computeJacobian(res, jac, mole_additions, dmole_additions);
5467 FAIL() <<
"Missing expected exception.";
5469 catch (
const std::exception & e)
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;
5479 const DenseVector<Real> res(3);
5480 const DenseVector<Real> mole_additions(5);
5482 _egs_kinetic_calcite.computeJacobian(res, jac, mole_additions, dmole_additions);
5483 FAIL() <<
"Missing expected exception.";
5485 catch (
const std::exception & e)
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;
5495 const DenseVector<Real> res(3);
5496 const DenseVector<Real> mole_additions(5);
5498 _egs_kinetic_calcite.computeJacobian(res, jac, mole_additions, dmole_additions);
5499 FAIL() <<
"Missing expected exception.";
5501 catch (
const std::exception & e)
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;
5515 const DenseVector<Real> mole_additions(4);
5518 FAIL() <<
"Missing expected exception.";
5520 catch (
const std::exception & e)
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;
5533 DenseVector<Real> mole_additions(5);
5537 DenseVector<Real> bad(4);
5539 FAIL() <<
"Missing expected exception.";
5541 catch (
const std::exception & e)
5543 std::string msg(e.what());
5544 ASSERT_TRUE(msg.find(
"addKineticRates: incorrectly sized additions: 4 5 5") !=
5546 <<
"Failed with unexpected error message: " << msg;
5552 FAIL() <<
"Missing expected exception.";
5554 catch (
const std::exception & e)
5556 std::string msg(e.what());
5557 ASSERT_TRUE(msg.find(
"addKineticRates: incorrectly sized additions: 5 4 5") !=
5559 <<
"Failed with unexpected error message: " << msg;
5565 FAIL() <<
"Missing expected exception.";
5567 catch (
const std::exception & e)
5569 std::string msg(e.what());
5570 ASSERT_TRUE(msg.find(
"addKineticRates: incorrectly sized additions: 5 5 4") !=
5572 <<
"Failed with unexpected error message: " << msg;
5580 {
"H2O",
"H+",
"HCO3-",
"O2(aq)",
"Ca++"},
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);
5653 {
"H2O",
"H+",
"HCO3-",
"O2(aq)",
"Ca++"},
5654 {1.75, 3.0, 2.0, 1.0, 1.5},
5660 {
"Calcite",
"CH4(aq)"},
5665 for (
unsigned i = 0; i < 7; ++i)
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);
5674 mgd_kin = mod.modelGeochemicalDatabase();
5675 egs.setModelGeochemicalDatabase(
mgd_kin);
5677 egs.addKineticRates(1.0, mole_additions, dmole_additions);
5678 for (
unsigned i = 0; i < 6; ++i)
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);
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);
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)
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);
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]);
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)
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);
5716 EXPECT_EQ(mole_additions(6), ch4rate);
5717 for (
unsigned j = 0;
j < 7; ++
j)
5718 EXPECT_EQ(dmole_additions(6,
j), ch4deriv[
j]);
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)
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);
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)
5742 calcitederiv[
j] = dmole_additions(5,
j);
5743 sum_deriv += std::abs(dmole_additions(5,
j));
5745 EXPECT_NE(sum_deriv, 0.0);
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)
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);
5762 EXPECT_EQ(mole_additions(5), 8 * calciterate);
5763 EXPECT_EQ(mole_additions(6), 2 * ch4rate);
5764 for (
unsigned j = 0;
j < 7; ++
j)
5766 EXPECT_EQ(dmole_additions(5,
j), 8 * calcitederiv[
j]);
5767 EXPECT_EQ(dmole_additions(6,
j), 2 * ch4deriv[
j]);
5775 {
"H2O",
"H+",
"HCO3-",
"O2(aq)",
"Ca++"},
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);
5868 {
"H2O",
"H+",
"HCO3-",
"O2(aq)",
"Ca++"},
5869 {1.75, 3.0, 2.0, 1.0, 1.5},
5875 {
"Calcite",
"CH4(aq)"},
5880 for (
unsigned i = 0; i < 7; ++i)
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);
5889 mgd_kin = mod.modelGeochemicalDatabase();
5890 egs.setModelGeochemicalDatabase(
mgd_kin);
5893 egs.addKineticRates(2.0, mole_additions, dmole_additions);
5894 for (
unsigned i = 0; i < 6; ++i)
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);
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);
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);
5915 const std::vector<Real> stoi = {1, 1, 1, -2, 0, 0};
5916 for (
unsigned i = 0; i < 6; ++i)
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);
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);
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)
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);
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);
5953 {
"H2O",
"H+",
"HCO3-",
"O2(aq)",
"Ca++"},
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);
6046 {
"H2O",
"H+",
"HCO3-",
"O2(aq)",
"Ca++"},
6047 {1.75, 3.0, 2.0, 1.0, 1.5},
6053 {
"Calcite",
"CH4(aq)"},
6058 for (
unsigned i = 0; i < 7; ++i)
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);
6067 mgd_kin = mod.modelGeochemicalDatabase();
6068 egs.setModelGeochemicalDatabase(
mgd_kin);
6071 egs.addKineticRates(2.0, mole_additions, dmole_additions);
6072 for (
unsigned i = 0; i < 6; ++i)
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);
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);
6085 mod.addKineticRate(rate_ch4_progH);
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)
6099 EXPECT_EQ(mole_additions(i), 0);
6100 for (
unsigned j = 0;
j < 6; ++
j)
6101 EXPECT_EQ(dmole_additions(i,
j), 0);
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);
6109 mod.addKineticRate(rate_ch4_progOH);
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)
6126 EXPECT_EQ(mole_additions(i), 0);
6127 for (
unsigned j = 0;
j < 6; ++
j)
6128 EXPECT_EQ(dmole_additions(i,
j), 0);
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);
6139 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
6144 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
6149 std::vector<GeochemistryUnitConverter::GeochemistryUnit> ku;
6159 {
"H2O",
"H+",
"HCO3-",
"Ca++"},
6160 {1.75, 1.0, 5.0, 2.0},
6169 EXPECT_EQ(egs_no_swaps.getBulkOldInOriginalBasis()(0), 1.75);
6170 EXPECT_NEAR(egs_no_swaps.getBulkOldInOriginalBasis()(1), 2.1, 1E-6);
6171 EXPECT_EQ(egs_no_swaps.getBulkOldInOriginalBasis()(2), 5.0 + 1.1);
6172 EXPECT_EQ(egs_no_swaps.getBulkOldInOriginalBasis()(3), 2.0);
6182 {
"H2O",
"H+",
"HCO3-",
"Calcite"},
6183 {-0.5, 2.5, 1.0, 3.5},
6194 EXPECT_EQ(egs.getBulkOldInOriginalBasis()(0), -0.5);
6195 EXPECT_EQ(egs.getBulkOldInOriginalBasis()(1), -2.5);
6196 EXPECT_EQ(egs.getBulkOldInOriginalBasis()(2), 4.5);
6197 EXPECT_EQ(egs.getBulkOldInOriginalBasis()(3), 3.5);
6200 egs.performSwap(2, 0);
6201 egs.performSwap(3, 3);
6204 EXPECT_NEAR(egs.getBulkOldInOriginalBasis()(0), -0.5, 1.0E-6);
6205 EXPECT_NEAR(egs.getBulkOldInOriginalBasis()(1), -2.5, 1.0E-6);
6206 EXPECT_NEAR(egs.getBulkOldInOriginalBasis()(2), 4.5, 1.0E-6);
6207 EXPECT_NEAR(egs.getBulkOldInOriginalBasis()(3), 3.5, 1.0E-6);
6214 const std::vector<std::string> original_basis = {
6215 "H2O",
"H+",
">(s)FeOH",
">(w)FeOH",
"Fe+++",
"HCO3-"};
6220 {
"Something",
"Fe(OH)3(ppd)fake"},
6226 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
6233 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
6240 std::vector<GeochemistryUnitConverter::GeochemistryUnit> ku;
6250 {
"H2O",
"H+",
">(s)FeOH",
">(w)FeOH",
"Fe+++",
"HCO3-"},
6251 {1.75, 1.0, 2.0, 3.0, 0.5, 1.0},
6257 {
"Something",
"Fe(OH)3(ppd)fake"},
6268 std::vector<Real> eqm_mol = egs.getEquilibriumMolality();
6270 std::vector<Real> original_trans_bulk;
6271 egs.computeTransportedBulkFromMolalities(original_trans_bulk);
6273 std::vector<Real> gold_trans_bulk(6);
6277 gold_trans_bulk[1] = 1.75 * (basis_mol[1] + eqm_mol[0] - eqm_mol[1] - eqm_mol[2]);
6279 gold_trans_bulk[2] = 0.0;
6281 gold_trans_bulk[3] = 0.0;
6283 gold_trans_bulk[4] = 1.75 * basis_mol[4];
6285 gold_trans_bulk[5] = 1.75 * (basis_mol[5] + eqm_mol[0] + eqm_mol[1]);
6289 for (
unsigned i = 0; i < 6; ++i)
6291 if (std::abs(gold_trans_bulk[i]) <=
eps)
6292 EXPECT_LE(std::abs(original_trans_bulk[i]),
eps);
6294 EXPECT_NEAR(original_trans_bulk[i] / gold_trans_bulk[i], 1.0,
eps);
6300 egs.performSwap(4, 4);
6303 basis_mol = egs.getSolventMassAndFreeMolalityAndMineralMoles();
6310 eqm_mol = egs.getEquilibriumMolality();
6311 gold_trans_bulk[0] =
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;
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]);
6321 std::vector<Real> new_trans_bulk;
6322 egs.computeTransportedBulkFromMolalities(new_trans_bulk);
6323 for (
unsigned i = 0; i < 6; ++i)
6325 if (std::abs(gold_trans_bulk[i]) <=
eps)
6326 EXPECT_LE(std::abs(new_trans_bulk[i]),
eps);
6328 EXPECT_NEAR(new_trans_bulk[i] / gold_trans_bulk[i], 1.0,
eps);
6332 std::vector<Real> gold_trans_bulk_orig(6);
6333 gold_trans_bulk_orig[0] =
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]);
6341 DenseVector<Real> trans_bulk_original_basis = egs.getTransportedBulkInOriginalBasis();
6342 for (
unsigned i = 0; i < 6; ++i)
6344 if (std::abs(gold_trans_bulk_orig[i]) <=
eps)
6345 EXPECT_LE(std::abs(trans_bulk_original_basis(i)),
eps);
6347 EXPECT_NEAR(trans_bulk_original_basis(i) / gold_trans_bulk_orig[i], 1.0,
eps);
6353 egs.setMineralRelatedFreeMoles(123.0);
6354 egs.computeTransportedBulkFromMolalities(new_trans_bulk);
6355 for (
unsigned i = 0; i < 6; ++i)
6357 if (std::abs(gold_trans_bulk[i]) <=
eps)
6358 EXPECT_LE(std::abs(new_trans_bulk[i]),
eps);
6360 EXPECT_NEAR(new_trans_bulk[i] / gold_trans_bulk[i], 1.0,
eps);
6362 trans_bulk_original_basis = egs.getTransportedBulkInOriginalBasis();
6363 for (
unsigned i = 0; i < 6; ++i)
6365 if (std::abs(gold_trans_bulk_orig[i]) <=
eps)
6366 EXPECT_LE(std::abs(trans_bulk_original_basis(i)),
eps);
6368 EXPECT_NEAR(trans_bulk_original_basis(i) / gold_trans_bulk_orig[i], 1.0,
eps);
6378 {
"H2O",
"H+",
"Fe++",
"O2(aq)"},
6387 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
6392 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
6397 std::vector<GeochemistryUnitConverter::GeochemistryUnit> ku;
6406 {
"H2O",
"H+",
"Fe++",
"O2(aq)"},
6407 {1.75, -2.0, 1.0, 1.0},
6419 std::vector<Real> eqm_mol = egs.getEquilibriumMolality();
6421 std::vector<Real> trans_bulk;
6422 egs.computeTransportedBulkFromMolalities(trans_bulk);
6423 std::vector<Real> gold_trans_bulk(4);
6427 gold_trans_bulk[1] = 1.75 * (basis_mol[1] - eqm_mol[0]);
6429 gold_trans_bulk[2] = 1.75 * (basis_mol[2]);
6431 gold_trans_bulk[3] = 1.75 * (basis_mol[3]);
6434 for (
unsigned i = 0; i < 4; ++i)
6435 EXPECT_NEAR(trans_bulk[i] / gold_trans_bulk[i], 1.0,
eps);
6443 database, {
"H2O",
"H+",
"HCO3-"}, {}, {}, {}, {}, {},
"O2(aq)",
"e-");
6445 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
6449 std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu;
6460 {
"H2O",
"H+",
"HCO3-"},
6473 dest = _egs_calcite;
6474 FAIL() <<
"Missing expected exception.";
6476 catch (
const std::exception & e)
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;
6491 {
"H2O",
"H+",
"HCO3-"},
6504 FAIL() <<
"Missing expected exception.";
6506 catch (
const std::exception & e)
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;
6515 _db_calcite, {
"H2O",
"H+",
"HCO3-",
"Ca++"}, {}, {}, {}, {}, {},
"O2(aq)",
"e-");
6517 model_no_kinetic_calcite.modelGeochemicalDatabase();
6525 {
"H2O",
"Ca++",
"H+",
"HCO3-"},
6526 {1.75, 3.0, 2.0, 1.0},
6537 egs_no_kinetic_calcite = _egs_kinetic_calcite;
6538 FAIL() <<
"Missing expected exception.";
6540 catch (
const std::exception & e)
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;
Real getTemperature() const
constexpr Real MOLES_PER_KG_WATER
void changeConstraintToBulk(unsigned basis_ind)
Changes the constraint to MOLES_BULK_SPECIES (or MOLES_BULK_WATER if basis_ind = 0) for the basis ind...
const std::vector< Real > & getBulkMolesOld() const
constexpr Real CELSIUS_TO_KELVIN
virtual void zero() override final
void setAdditionalValue(const std::string &names)
TEST_F(GeochemicalSystemTest, constructWithMultiMooseEnum)
Check MultiMooseEnum constructor.
std::unordered_map< std::string, unsigned > basis_species_index
basis_species_index[name] = index of the basis species, within all ModelGeochemicalDatabase internal ...
const ModelGeochemicalDatabase mgd
std::unordered_map< std::string, unsigned > eqm_species_index
eqm_species_index[name] = index of the equilibrium species (secondary aqueous species, redox couples in equilibrium with the aqueous solution, minerals in equilibrium with the aqueous solution, gases in equilibrium with the aqueous solution) within all ModelGeochemicalDatabase internal datastrcutres, with given name
void updateOldWithCurrent(const DenseVector< Real > &mole_additions)
Usually used at the end of a solve, to provide correct initial conditions to the next time-step...
const GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false)
const std::vector< Real > & computeAndGetEquilibriumActivity()
Computes activity for all equilibrium species (_eqm_activity) and returns a reference to the vector...
unsigned getNumInBasis() const
returns the number of species in the basis
std::unordered_map< std::string, unsigned > kin_species_index
kin_species_index[name] = index of the kinetic species, within all ModelGeochemicalDatabase internal ...
DenseMatrix< Real > eqm_stoichiometry
eqm_stoichiometry(i, j) = stoichiometric coefficient for equilibrium species "i" in terms of the basi...
KineticRateUserDescription rate_ch4("CH4(aq)", 1.5, 2.0, true, 0.0, 0.0, 0.0, {"H2O", "OH-", "O2(aq)", "CO2(aq)", "CaCO3"}, {3.0, 3.1, 3.2, 3.3, 3.4}, {0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0}, 0.8, 2.5, 66.0, 0.003, DirectionChoiceEnum::BOTH, "H2O", 0.0, -1.0, 0.0)
Real getTotalChargeOld() const
void setKineticMoles(unsigned kin, Real moles)
Sets the current AND old mole number for a kinetic species.
const GeochemicalDatabaseReader db("database/moose_testdb.json", true, true, false)
Real getBasisActivity(unsigned i) const
const ModelGeochemicalDatabase & mgd_kin
Class to swap basis species with equilibrium species.
Computes activity coefficients for non-minerals and non-gases (since these species do not have activi...
Real getEquilibriumActivity(unsigned eqm_ind) const
Returns the value of activity for the equilibrium species with index eqm_index.
const ModelGeochemicalDatabase & modelGeochemicalDatabase() const
Return a reference to the ModelGeochemicalDatabase structure.
Constructs and stores a minimal amount of information that is pertinent to the user-defined geochemic...
Calculators to compute ionic strength and stoichiometric ionic strength.
const ModelGeochemicalDatabase & getModelGeochemicalDatabase() const
void addKineticRates(Real dt, DenseVector< Real > &mole_additions, DenseMatrix< Real > &dmole_additions)
Computes the kinetic rates and their derivatives based on the current values of molality, mole number, etc, and then, using dt, calculates the appropriate mole_additions and dmole_additions.
PetscErrorCode PetscInt const PetscInt IS * is
void enforceChargeBalance()
Enforces charge balance by altering the constraint_value and bulk_moles_old of the charge-balance spe...
const PertinentGeochemicalSystem model(database, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++", ">(s)FeOH", "radius_neg1", "radius_neg1.5"}, {"Calcite"}, {}, {"Calcite_asdf"}, {"CH4(aq)"}, {">(s)FeOCa+"}, "O2(aq)", "e-")
KineticRateUserDescription rate_cal("Calcite", 7.0, 6.0, false, 0.0, 0.0, 0.0, {"H+"}, {-3.0}, {0.0}, {0.0}, 2.5, 0.8, 55.0, 0.00315, DirectionChoiceEnum::BOTH, "H2O", 0.0, -1.0, 0.0)
void setModelGeochemicalDatabase(const ModelGeochemicalDatabase &mgd)
Copies a ModelGeochemicalDatabase into our _mgd structure.
Real getLog10K(unsigned j) const
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const std::vector< GeochemistryUnitConverter::GeochemistryUnit > cu2
std::vector< std::string > originalBasisNames() const
Real getEquilibriumMolality(unsigned j) const
void setTemperature(Real temperature)
Sets the temperature to the given quantity.
constexpr Real GAS_CONSTANT
This class holds information about bulk composition, molalities, activities, activity coefficients...
Holds a user-specified description of a kinetic rate.
void setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(const std::vector< std::string > &names, const std::vector< Real > &values, const std::vector< bool > &constraints_from_molalities)
Sets:
Real getEquilibriumActivityCoefficient(unsigned j) const
Data structure to hold all relevant information from the database file.
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
std::vector< bool > eqm_species_mineral
eqm_species_mineral[i] = true iff the i^th equilibrium species is a mineral
const std::vector< Real > & getSolventMassAndFreeMolalityAndMineralMoles() const
void addToBulkMoles(unsigned basis_ind, Real value)
Add to the MOLES_BULK_SPECIES (or MOLES_BULK_WATER if basis_ind = 0) for the basis species...
KineticRateUserDescription rate_ch4_death("CH4(aq)", 1.5, 2.0, true, 0.0, 0.0, 0.0, {"H2O", "OH-", "O2(aq)", "CO2(aq)", "CaCO3"}, {3.0, 3.1, 3.2, 3.3, 3.4}, {0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0}, 0.8, 2.5, 66.0, 0.003, DirectionChoiceEnum::DEATH, "H2O", 0.0, -1.0, 0.0)
void performSwap(unsigned swap_out_of_basis, unsigned swap_into_basis)
Perform the basis swap, and ensure that the resulting system is consistent.
void setConstraintValue(unsigned basis_ind, Real value)
Set the constraint value for the basis species.
unsigned getNumInEquilibrium() const
returns the number of species in equilibrium with the basis components
void closeSystem()
Changes a KG_SOLVENT_WATER constraint to MOLES_BULK_WATER (if there is such a constraint) and all FRE...
const std::vector< GeochemicalSystem::ConstraintUserMeaningEnum > cm2
MooseUnits pow(const MooseUnits &, int)
std::vector< std::string > eqm_species_name
eqm_species_name[i] = name of the i^th eqm species
static const std::string k
std::vector< std::string > basis_species_name
basis_species_name[j] = name of the j^th basis species
Class for reading geochemical reactions from a MOOSE geochemical database.
std::vector< Real > getSaturationIndices() const