https://mooseframework.inl.gov
Loading...
Searching...
No Matches
GeochemistryIonicStrengthTest.C
Go to the documentation of this file.
1//* This file is part of the MOOSE framework
2//* https://mooseframework.inl.gov
3//*
4//* All rights reserved, see COPYRIGHT for full restrictions
5//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6//*
7//* Licensed under LGPL 2.1, please see LICENSE for details
8//* https://www.gnu.org/licenses/lgpl-2.1.html
9
10#include "gtest/gtest.h"
11
12#include "GeochemistryIonicStrength.h"
13#include "GeochemistrySpeciesSwapper.h"
14
15const GeochemicalDatabaseReader database("database/moose_testdb.json");
16// The following system has secondary species: CO2(aq), CO3--, CaCO3, CaOH+, OH-, (O-phth)--,
17// >(s)FeO-
18const PertinentGeochemicalSystem model(database,
19 {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++", ">(s)FeOH"},
20 {"Calcite"},
21 {},
22 {"Calcite_asdf"},
23 {"CH4(aq)"},
24 {">(s)FeOCa+"},
25 "O2(aq)",
26 "e-");
27const ModelGeochemicalDatabase mgd = model.modelGeochemicalDatabase();
28const GeochemistryIonicStrength is_calc(1E9, 1E9, false, false);
29const GeochemistryIonicStrength is_basis_only(1E9, 1E9, true, false);
30const GeochemistryIonicStrength is_cl_only(1E9, 1E9, false, true);
31
32TEST(GeochemistryIonicStrengthTest, sizeExceptions)
33{
34 const std::vector<Real> basis_m_good(6);
35 const std::vector<Real> eqm_m_good(8);
36 const std::vector<Real> kin_m_good(3);
37
38 try
39 {
40 is_calc.ionicStrength(mgd, {}, eqm_m_good, kin_m_good);
41 FAIL() << "Missing expected exception.";
42 }
43 catch (const std::exception & e)
44 {
45 std::string msg(e.what());
46 ASSERT_TRUE(
47 msg.find(
48 "Ionic strength calculation: Number of basis species in mgd not equal to the size of "
49 "basis_species_molality") != std::string::npos)
50 << "Failed with unexpected error message: " << msg;
51 }
52
53 try
54 {
55 is_calc.ionicStrength(mgd, basis_m_good, {}, kin_m_good);
56 FAIL() << "Missing expected exception.";
57 }
58 catch (const std::exception & e)
59 {
60 std::string msg(e.what());
61 ASSERT_TRUE(msg.find("Ionic strength calculation: Number of equilibrium species in mgd not "
62 "equal to the size of "
63 "eqm_species_molality") != std::string::npos)
64 << "Failed with unexpected error message: " << msg;
65 }
66
67 try
68 {
69 is_calc.ionicStrength(mgd, basis_m_good, eqm_m_good, {});
70 FAIL() << "Missing expected exception.";
71 }
72 catch (const std::exception & e)
73 {
74 std::string msg(e.what());
75 ASSERT_TRUE(msg.find("Ionic strength calculation: Number of kinetic species in mgd not "
76 "equal to the size of "
77 "kin_species_molality") != std::string::npos)
78 << "Failed with unexpected error message: " << msg;
79 }
80
81 try
82 {
83 is_calc.stoichiometricIonicStrength(mgd, {}, eqm_m_good, kin_m_good);
84 FAIL() << "Missing expected exception.";
85 }
86 catch (const std::exception & e)
87 {
88 std::string msg(e.what());
89 ASSERT_TRUE(msg.find("Stoichiometric ionic strength calculation: Number of basis species in "
90 "mgd not equal to the size of "
91 "basis_species_molality") != std::string::npos)
92 << "Failed with unexpected error message: " << msg;
93 }
94
95 try
96 {
97 is_calc.stoichiometricIonicStrength(mgd, basis_m_good, {}, kin_m_good);
98 FAIL() << "Missing expected exception.";
99 }
100 catch (const std::exception & e)
101 {
102 std::string msg(e.what());
103 ASSERT_TRUE(
104 msg.find(
105 "Stoichiometric ionic strength calculation: Number of equilibrium species in mgd not "
106 "equal to the size of "
107 "eqm_species_molality") != std::string::npos)
108 << "Failed with unexpected error message: " << msg;
109 }
110
111 try
112 {
113 is_calc.stoichiometricIonicStrength(mgd, basis_m_good, eqm_m_good, {});
114 FAIL() << "Missing expected exception.";
115 }
116 catch (const std::exception & e)
117 {
118 std::string msg(e.what());
119 ASSERT_TRUE(
120 msg.find("Stoichiometric ionic strength calculation: Number of kinetic species in mgd not "
121 "equal to the size of "
122 "kin_species_molality") != std::string::npos)
123 << "Failed with unexpected error message: " << msg;
124 }
125}
126
128TEST(GeochemistryIonicStrengthTest, getsetMax)
129{
130 GeochemistryIonicStrength is(1.0, 2.0, false, false);
131 EXPECT_EQ(is.getMaxIonicStrength(), 1.0);
132 EXPECT_EQ(is.getMaxStoichiometricIonicStrength(), 2.0);
133
134 is.setMaxIonicStrength(3.25);
135 is.setMaxStoichiometricIonicStrength(4.5);
136
137 EXPECT_EQ(is.getMaxIonicStrength(), 3.25);
138 EXPECT_EQ(is.getMaxStoichiometricIonicStrength(), 4.5);
139}
140
142TEST(GeochemistryIonicStrengthTest, getsetOnly)
143{
144 EXPECT_FALSE(is_calc.getUseOnlyBasisMolality());
145 EXPECT_TRUE(is_basis_only.getUseOnlyBasisMolality());
146 EXPECT_FALSE(is_cl_only.getUseOnlyBasisMolality());
147 EXPECT_FALSE(is_calc.getUseOnlyClMolality());
148 EXPECT_FALSE(is_basis_only.getUseOnlyClMolality());
149 EXPECT_TRUE(is_cl_only.getUseOnlyClMolality());
150
151 GeochemistryIonicStrength is_tmp(0.125, 1E9, false, false);
152 EXPECT_FALSE(is_tmp.getUseOnlyBasisMolality());
153 EXPECT_FALSE(is_tmp.getUseOnlyClMolality());
154 is_tmp.setUseOnlyBasisMolality(true);
155 EXPECT_TRUE(is_tmp.getUseOnlyBasisMolality());
156 EXPECT_FALSE(is_tmp.getUseOnlyClMolality());
157 is_tmp.setUseOnlyClMolality(true);
158 EXPECT_TRUE(is_tmp.getUseOnlyBasisMolality());
159 EXPECT_TRUE(is_tmp.getUseOnlyClMolality());
160 is_tmp.setUseOnlyBasisMolality(false);
161 EXPECT_FALSE(is_tmp.getUseOnlyBasisMolality());
162 EXPECT_TRUE(is_tmp.getUseOnlyClMolality());
163}
164
166TEST(GeochemistryIonicStrengthTest, ionicStrength)
167{
168 Real gold_ionic_str = 0.0;
169
170 std::vector<Real> basis_m(6);
171 basis_m[mgd.basis_species_index.at("H2O")] = 1.0;
172 basis_m[mgd.basis_species_index.at("H+")] = 2.0;
173 gold_ionic_str += 2.0;
174 basis_m[mgd.basis_species_index.at("HCO3-")] = 3.0;
175 gold_ionic_str += 3.0;
176 basis_m[mgd.basis_species_index.at("O2(aq)")] = 4.0;
177 basis_m[mgd.basis_species_index.at("Ca++")] = 5.0;
178 gold_ionic_str += 4 * 5.0;
179 basis_m[mgd.basis_species_index.at(">(s)FeOH")] = 6.0;
180
181 const Real gold_basis_only = 0.5 * gold_ionic_str;
182
183 std::vector<Real> eqm_m(8);
184 eqm_m[mgd.eqm_species_index.at("CO2(aq)")] = 1.1;
185 eqm_m[mgd.eqm_species_index.at("CO3--")] = 2.2;
186 gold_ionic_str += 4 * 2.2;
187 eqm_m[mgd.eqm_species_index.at("CaCO3")] = 3.3;
188 eqm_m[mgd.eqm_species_index.at("CaOH+")] = 4.4;
189 gold_ionic_str += 4.4;
190 eqm_m[mgd.eqm_species_index.at("OH-")] = 5.5;
191 gold_ionic_str += 5.5;
192 eqm_m[mgd.eqm_species_index.at("(O-phth)--")] = 6.6;
193 gold_ionic_str += 4 * 6.6;
194 eqm_m[mgd.eqm_species_index.at(">(s)FeO-")] = 7.7;
195 gold_ionic_str += 7.7;
196 eqm_m[mgd.eqm_species_index.at("Calcite")] = 9.9;
197
198 std::vector<Real> kin_m(3);
199 kin_m[mgd.kin_species_index.at("Calcite_asdf")] = -1.1;
200 kin_m[mgd.kin_species_index.at("CH4(aq)")] = -2.2;
201 kin_m[mgd.kin_species_index.at(">(s)FeOCa+")] = -3.3;
202 gold_ionic_str += -3.3;
203
204 gold_ionic_str *= 0.5;
205
206 const Real ionic_str = is_calc.ionicStrength(mgd, basis_m, eqm_m, kin_m);
207 const Real ionic_str_basis_only = is_basis_only.ionicStrength(mgd, basis_m, eqm_m, kin_m);
208
209 EXPECT_NEAR(ionic_str, gold_ionic_str, 1E-9);
210 EXPECT_NEAR(ionic_str_basis_only, gold_basis_only, 1E-9);
211
212 GeochemistryIonicStrength is_max(0.125, 1E9, false, false);
213 EXPECT_EQ(is_max.ionicStrength(mgd, basis_m, eqm_m, kin_m), 0.125);
214}
215
217TEST(GeochemistryIonicStrengthTest, stoichiometricIonicStrength)
218{
219 Real gold_ionic_str = 0.0;
220
221 std::vector<Real> basis_m(6);
222 basis_m[mgd.basis_species_index.at("H2O")] = 1.0;
223 basis_m[mgd.basis_species_index.at("H+")] = 2.0;
224 gold_ionic_str += 2.0;
225 basis_m[mgd.basis_species_index.at("HCO3-")] = 3.0;
226 gold_ionic_str += 3.0;
227 basis_m[mgd.basis_species_index.at("O2(aq)")] = 4.0;
228 basis_m[mgd.basis_species_index.at("Ca++")] = 5.0;
229 gold_ionic_str += 5.0 * 4;
230 basis_m[mgd.basis_species_index.at(">(s)FeOH")] = 6.0;
231
232 const Real gold_basis_only = 0.5 * gold_ionic_str;
233
234 std::vector<Real> eqm_m(8);
235 eqm_m[mgd.eqm_species_index.at("CO2(aq)")] = 1.1;
236 gold_ionic_str += 1.1 * (1 + 1);
237 eqm_m[mgd.eqm_species_index.at("CO3--")] = 2.2;
238 gold_ionic_str += 2.2 * 4;
239 eqm_m[mgd.eqm_species_index.at("CaCO3")] = 3.3;
240 gold_ionic_str += 3.3 * (4 + 1 - 1);
241 eqm_m[mgd.eqm_species_index.at("CaOH+")] = 4.4;
242 gold_ionic_str += 4.4 * (1);
243 eqm_m[mgd.eqm_species_index.at("OH-")] = 5.5;
244 gold_ionic_str += 5.5 * (1);
245 eqm_m[mgd.eqm_species_index.at("(O-phth)--")] = 6.6;
246 gold_ionic_str += 6.6 * (4);
247 eqm_m[mgd.eqm_species_index.at(">(s)FeO-")] = 7.7;
248 gold_ionic_str += 7.7 * (1);
249 eqm_m[mgd.eqm_species_index.at("Calcite")] = 9.9;
250 gold_ionic_str += 9.9 * (4 + 1 - 1);
251
252 std::vector<Real> kin_m(3);
253 kin_m[mgd.kin_species_index.at("Calcite_asdf")] = -1.1;
254 gold_ionic_str += 0.0;
255 kin_m[mgd.kin_species_index.at("CH4(aq)")] = -2.2;
256 gold_ionic_str += 0.0;
257 kin_m[mgd.kin_species_index.at(">(s)FeOCa+")] = -3.3;
258 gold_ionic_str += -3.3 * (1);
259
260 gold_ionic_str *= 0.5;
261
262 const Real ionic_str = is_calc.stoichiometricIonicStrength(mgd, basis_m, eqm_m, kin_m);
263 const Real ionic_str_basis_only =
264 is_basis_only.stoichiometricIonicStrength(mgd, basis_m, eqm_m, kin_m);
265
266 EXPECT_NEAR(ionic_str, gold_ionic_str, 1E-9);
267 EXPECT_NEAR(ionic_str_basis_only, gold_basis_only, 1E-9);
268
269 GeochemistryIonicStrength is_max(1E9, 0.125, false, false);
270 EXPECT_EQ(is_max.stoichiometricIonicStrength(mgd, basis_m, eqm_m, kin_m), 0.125);
271}
272
274TEST(GeochemistryIonicStrengthTest, stoichiometricIonicStrength_except)
275{
276 std::vector<Real> basis_m(6);
277 std::vector<Real> eqm_m(8);
278 std::vector<Real> kin_m(3);
279 try
280 {
281 is_cl_only.stoichiometricIonicStrength(mgd, basis_m, eqm_m, kin_m);
282 FAIL() << "Missing expected exception.";
283 }
284 catch (const std::exception & e)
285 {
286 std::string msg(e.what());
287 ASSERT_TRUE(
288 msg.find(
289 "GeochemistryIonicStrength: attempting to compute stoichiometric ionic strength "
290 "using only the Cl- molality, but Cl- does not appear in the geochemical system") !=
291 std::string::npos)
292 << "Failed with unexpected error message: " << msg;
293 }
294}
295
297TEST(GeochemistryIonicStrengthTest, stoichiometricIonicStrengthOnlyCl)
298{
299 const GeochemicalDatabaseReader db_cl("../database/moose_geochemdb.json", true, true, false);
300
301 // Cl is a basis species
302 const PertinentGeochemicalSystem model_cl(
303 db_cl, {"H2O", "H+", "Cl-"}, {}, {}, {}, {}, {}, "O2(aq)", "e-");
304 const ModelGeochemicalDatabase mgd_cl = model_cl.modelGeochemicalDatabase();
305 const GeochemistryIonicStrength is_cl(1E9, 1E9, false, true);
306 const std::vector<Real> basis_m{0.75, 2.5, 3.5};
307 const std::vector<Real> eqm_m{1.5, 0.25};
308 EXPECT_EQ(is_cl.stoichiometricIonicStrength(mgd_cl, basis_m, eqm_m, {}), 3.5);
309
310 // Cl is an equilibrium species
311 const PertinentGeochemicalSystem model_cl_eqm(
312 db_cl, {"H2O", "H+", "Cl-"}, {}, {}, {}, {}, {}, "O2(aq)", "e-");
313 ModelGeochemicalDatabase mgd_cl_eqm = model_cl_eqm.modelGeochemicalDatabase();
314 GeochemistrySpeciesSwapper swapper(3, 1E-6);
315 swapper.performSwap(mgd_cl_eqm, "Cl-", "HCl");
316 const GeochemistryIonicStrength is_cl_eqm(1E9, 1E9, false, true);
317 const std::vector<Real> basis_m_eqm{0.75, 2.5, 3.5};
318 const std::vector<Real> eqm_m_eqm{1.5, 0.25};
319 const unsigned pos = mgd_cl_eqm.eqm_species_name[0] == "Cl-" ? 0 : 1;
320 EXPECT_EQ(is_cl_eqm.stoichiometricIonicStrength(mgd_cl_eqm, basis_m_eqm, eqm_m_eqm, {}),
321 eqm_m_eqm[pos]);
322
323 // Cl is a kinetic species
324 const GeochemicalDatabaseReader db_cl_kin(
325 "database/moose_testdb_cl_kinetic.json", true, true, false);
326 const PertinentGeochemicalSystem model_cl_kin(
327 db_cl_kin, {"H2O"}, {}, {}, {}, {"Cl-"}, {}, "O2(aq)", "e-");
328 ModelGeochemicalDatabase mgd_cl_kin = model_cl_kin.modelGeochemicalDatabase();
329 const GeochemistryIonicStrength is_cl_kin(1E9, 1E9, false, true);
330 const std::vector<Real> basis_m_kin{0.75};
331 const std::vector<Real> kin_m_kin{7.75};
332 EXPECT_EQ(is_cl_kin.stoichiometricIonicStrength(mgd_cl_kin, basis_m_kin, {}, kin_m_kin), 7.75);
333}
model
const PertinentGeochemicalSystem model(database, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++", ">(s)FeOH"}, {"Calcite"}, {}, {"Calcite_asdf"}, {"CH4(aq)"}, {">(s)FeOCa+"}, "O2(aq)", "e-")
TEST
TEST(GeochemistryIonicStrengthTest, sizeExceptions)
Definition GeochemistryIonicStrengthTest.C:32
is_calc
const GeochemistryIonicStrength is_calc(1E9, 1E9, false, false)
mgd
const ModelGeochemicalDatabase mgd
Definition GeochemistryIonicStrengthTest.C:27
is_cl_only
const GeochemistryIonicStrength is_cl_only(1E9, 1E9, false, true)
is_basis_only
const GeochemistryIonicStrength is_basis_only(1E9, 1E9, true, false)
database
const GeochemicalDatabaseReader database("database/moose_testdb.json")
GeochemicalDatabaseReader
Class for reading geochemical reactions from a MOOSE geochemical database.
Definition GeochemicalDatabaseReader.h:216
GeochemistryIonicStrength
Calculators to compute ionic strength and stoichiometric ionic strength.
Definition GeochemistryIonicStrength.h:18
GeochemistryIonicStrength::stoichiometricIonicStrength
Real stoichiometricIonicStrength(const ModelGeochemicalDatabase &mgd, const std::vector< Real > &basis_species_molality, const std::vector< Real > &eqm_species_molality, const std::vector< Real > &kin_species_molality) const
Compute stoichiometric ionic strength.
Definition GeochemistryIonicStrength.C:61
GeochemistryIonicStrength::setUseOnlyBasisMolality
void setUseOnlyBasisMolality(bool use_only_basis_molality)
Set the value of use_only_basis_molality.
Definition GeochemistryIonicStrength.C:147
GeochemistryIonicStrength::ionicStrength
Real ionicStrength(const ModelGeochemicalDatabase &mgd, const std::vector< Real > &basis_species_molality, const std::vector< Real > &eqm_species_molality, const std::vector< Real > &kin_species_molality) const
Compute ionic strength.
Definition GeochemistryIonicStrength.C:25
GeochemistryIonicStrength::getUseOnlyClMolality
Real getUseOnlyClMolality() const
Return the value of use_only_Cl_molality.
Definition GeochemistryIonicStrength.C:165
GeochemistryIonicStrength::setUseOnlyClMolality
void setUseOnlyClMolality(bool use_only_Cl_molality)
Set the value of use_only_Cl_molality.
Definition GeochemistryIonicStrength.C:159
GeochemistryIonicStrength::getUseOnlyBasisMolality
Real getUseOnlyBasisMolality() const
Return the value of use_only_basis_molality.
Definition GeochemistryIonicStrength.C:153
GeochemistrySpeciesSwapper
Class to swap basis species with equilibrium species.
Definition GeochemistrySpeciesSwapper.h:20
GeochemistrySpeciesSwapper::performSwap
void performSwap(ModelGeochemicalDatabase &mgd, const std::string &replace_this, const std::string &with_this)
Check that replacing the named basis species with the named equilibrium species is valid,...
Definition GeochemistrySpeciesSwapper.C:114
PertinentGeochemicalSystem
Constructs and stores a minimal amount of information that is pertinent to the user-defined geochemic...
Definition PertinentGeochemicalSystem.h:387
PertinentGeochemicalSystem::modelGeochemicalDatabase
const ModelGeochemicalDatabase & modelGeochemicalDatabase() const
Return a reference to the ModelGeochemicalDatabase structure.
Definition PertinentGeochemicalSystem.C:850
ModelGeochemicalDatabase
Data structure to hold all relevant information from the database file.
Definition PertinentGeochemicalSystem.h:102
ModelGeochemicalDatabase::eqm_species_index
std::unordered_map< std::string, unsigned > eqm_species_index
eqm_species_index[name] = index of the equilibrium species (secondary aqueous species,...
Definition PertinentGeochemicalSystem.h:194
ModelGeochemicalDatabase::basis_species_index
std::unordered_map< std::string, unsigned > basis_species_index
basis_species_index[name] = index of the basis species, within all ModelGeochemicalDatabase internal ...
Definition PertinentGeochemicalSystem.h:162
ModelGeochemicalDatabase::kin_species_index
std::unordered_map< std::string, unsigned > kin_species_index
kin_species_index[name] = index of the kinetic species, within all ModelGeochemicalDatabase internal ...
Definition PertinentGeochemicalSystem.h:304
Generated on Fri Sep 11 2026 20:36:06 for https://mooseframework.inl.gov by doxygen 1.9.8