https://mooseframework.inl.gov
Loading...
Searching...
No Matches
GeochemistrySpeciesSwapperTest.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 "GeochemistrySpeciesSwapper.h"
13
14const Real eps =
15 1E-12; // accounts for precision loss when substituting reactions and inverting the swap matrix
16
18TEST(GeochemistrySpeciesSwapperTest, bulkCompositionException)
19{
20 GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
21
22 // The following system has secondary species: CO2(aq), CO3--, OH-,
23 // (O-phth)--, CH4(aq), Fe+++,
24 // >(s)FeO-
25 PertinentGeochemicalSystem model(database,
26 {"H2O", "H+", ">(s)FeOH", ">(w)FeOH", "Fe++", "HCO3-", "O2(aq)"},
27 {"Fe(OH)3(ppd)fake"},
28 {"CH4(g)fake"},
29 {},
30 {},
31 {},
32 "O2(aq)",
33 "e-");
34 ModelGeochemicalDatabase mgd = model.modelGeochemicalDatabase();
35 GeochemistrySpeciesSwapper swapper(mgd.basis_species_index.size(), 1E-6);
36
37 try
38 {
39 DenseVector<Real> bulk_composition(8);
40 swapper.performSwap(mgd, bulk_composition, "H+", "(O-phth)--");
41 FAIL() << "Missing expected exception.";
42 }
43 catch (const std::exception & e)
44 {
45 std::string msg(e.what());
46 ASSERT_TRUE(msg.find("GeochemistrySpeciesSwapper: bulk_composition has size 8 which differs "
47 "from the basis size") != std::string::npos)
48 << "Failed with unexpected error message: " << msg;
49 }
50}
51
53TEST(GeochemistrySpeciesSwapperTest, swapExceptions)
54{
55 GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
56
57 // The following system has secondary species: CO2(aq), CO3--, OH-,
58 // (O-phth)--, CH4(aq), Fe+++,
59 // >(s)FeO-
60 PertinentGeochemicalSystem model(database,
61 {"H2O", "H+", ">(s)FeOH", ">(w)FeOH", "Fe++", "HCO3-", "O2(aq)"},
62 {"Fe(OH)3(ppd)"},
63 {"CH4(g)fake"},
64 {},
65 {},
66 {},
67 "O2(aq)",
68 "e-");
69 ModelGeochemicalDatabase mgd = model.modelGeochemicalDatabase();
70 GeochemistrySpeciesSwapper swapper(mgd.basis_species_index.size(), 1E-6);
71
72 try
73 {
74 swapper.checkSwap(mgd, "CO2(aq)", "CO2(aq)");
75 FAIL() << "Missing expected exception.";
76 }
77 catch (const std::exception & e)
78 {
79 std::string msg(e.what());
80 ASSERT_TRUE(msg.find("CO2(aq) is not in the basis, so cannot be removed "
81 "from the basis") != std::string::npos)
82 << "Failed with unexpected error message: " << msg;
83 }
84
85 try
86 {
87 swapper.checkSwap(mgd, "CH4(g)fake", "CO2(aq)");
88 FAIL() << "Missing expected exception.";
89 }
90 catch (const std::exception & e)
91 {
92 std::string msg(e.what());
93 ASSERT_TRUE(msg.find("CH4(g)fake is not in the basis, so cannot be removed "
94 "from the basis") != std::string::npos)
95 << "Failed with unexpected error message: " << msg;
96 }
97
98 try
99 {
100 swapper.checkSwap(mgd, "H+", ">(s)FeOH");
101 FAIL() << "Missing expected exception.";
102 }
103 catch (const std::exception & e)
104 {
105 std::string msg(e.what());
106 ASSERT_TRUE(msg.find(">(s)FeOH is not an equilibrium species, so cannot be "
107 "removed from the "
108 "equilibrium species list") != std::string::npos)
109 << "Failed with unexpected error message: " << msg;
110 }
111
112 try
113 {
114 swapper.checkSwap(mgd, "H2O", "CO2(aq)");
115 FAIL() << "Missing expected exception.";
116 }
117 catch (const std::exception & e)
118 {
119 std::string msg(e.what());
120 ASSERT_TRUE(msg.find("Cannot remove H2O from the basis") != std::string::npos)
121 << "Failed with unexpected error message: " << msg;
122 }
123
124 try
125 {
126 swapper.checkSwap(mgd, 123, 0);
127 FAIL() << "Missing expected exception.";
128 }
129 catch (const std::exception & e)
130 {
131 std::string msg(e.what());
132 ASSERT_TRUE(msg.find("123 exceeds the number of basis species in the problem") !=
133 std::string::npos)
134 << "Failed with unexpected error message: " << msg;
135 }
136
137 try
138 {
139 swapper.checkSwap(mgd, 0, 0);
140 FAIL() << "Missing expected exception.";
141 }
142 catch (const std::exception & e)
143 {
144 std::string msg(e.what());
145 ASSERT_TRUE(msg.find("Cannot remove H2O from the basis") != std::string::npos)
146 << "Failed with unexpected error message: " << msg;
147 }
148
149 try
150 {
151 swapper.checkSwap(mgd, 1, 123);
152 FAIL() << "Missing expected exception.";
153 }
154 catch (const std::exception & e)
155 {
156 std::string msg(e.what());
157 ASSERT_TRUE(msg.find("123 exceeds the number of equilibrium species in the problem") !=
158 std::string::npos)
159 << "Failed with unexpected error message: " << msg;
160 }
161
162 try
163 {
164 swapper.performSwap(mgd, ">(s)FeOH", "CO2(aq)");
165 FAIL() << "Missing expected exception.";
166 }
167 catch (const std::exception & e)
168 {
169 std::string msg(e.what());
170 ASSERT_TRUE(msg.find("Matrix is not invertible, which signals an invalid basis swap") !=
171 std::string::npos)
172 << "Failed with unexpected error message: " << msg;
173 }
174
175 try
176 {
177 swapper = GeochemistrySpeciesSwapper(8, 1E-6);
178 swapper.checkSwap(mgd, "Fe++", "Fe+++");
179 FAIL() << "Missing expected exception.";
180 }
181 catch (const std::exception & e)
182 {
183 std::string msg(e.what());
184 ASSERT_TRUE(
185 msg.find("GeochemistrySpeciesSwapper constructed with incorrect basis_species size") !=
186 std::string::npos)
187 << "Failed with unexpected error message: " << msg;
188 }
189
190 try
191 {
192 swapper.checkSwap(mgd, ">(s)FeOH", ">(s)FeO-");
193 FAIL() << "Missing expected exception.";
194 }
195 catch (const std::exception & e)
196 {
197 std::string msg(e.what());
198 ASSERT_TRUE(
199 msg.find("Equilibrium species >(s)FeO- is involved in surface sorption so cannot be "
200 "swapped into the basis. If this is truly necessary, code enhancements will need "
201 "to be made including: recording whether basis species are involved in surface "
202 "sorption, including them in the surface-potential calculations, and carefully "
203 "swapping surface-potential-modified equilibrium constants") != std::string::npos)
204 << "Failed with unexpected error message: " << msg;
205 }
206}
207
209TEST(GeochemistrySpeciesSwapperTest, swap1)
210{
211 GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
212
213 // eqm species are: CO2(aq), CO3--, CaCO3, CaOH+, OH-, Calcite
214 PertinentGeochemicalSystem model(
215 database, {"H2O", "Ca++", "HCO3-", "H+"}, {"Calcite"}, {}, {}, {}, {}, "O2(aq)", "e-");
216 ModelGeochemicalDatabase mgd = model.modelGeochemicalDatabase();
217 GeochemistrySpeciesSwapper swapper(mgd.basis_species_index.size(), 1E-6);
218 DenseVector<Real> bulk_composition(4);
219 bulk_composition(0) = 0.5;
220 bulk_composition(1) = 1.0;
221 bulk_composition(2) = 2.5;
222 bulk_composition(3) = 3.0;
223
224 ASSERT_EQ(mgd.basis_species_index.size(), (std::size_t)4);
225 for (const auto & sp : mgd.basis_species_index)
226 ASSERT_EQ(mgd.basis_species_name[sp.second], sp.first);
227
228 ASSERT_EQ(mgd.eqm_species_index.size(), (std::size_t)6);
229 for (const auto & sp : mgd.eqm_species_index)
230 ASSERT_EQ(mgd.eqm_species_name[sp.second], sp.first);
231
232 ASSERT_EQ(mgd.swap_to_original_basis.n(), (unsigned)0);
233 ASSERT_EQ(mgd.have_swapped_out_of_basis.size(), (std::size_t)0);
234 ASSERT_EQ(mgd.have_swapped_into_basis.size(), (std::size_t)0);
235
236 for (const auto & species : mgd.basis_species_index)
237 ASSERT_EQ(mgd.basis_species_mineral[species.second], false);
238 for (const auto & species : mgd.eqm_species_index)
239 if (species.first == "Calcite")
240 ASSERT_EQ(mgd.eqm_species_mineral[species.second], true);
241 else
242 ASSERT_EQ(mgd.eqm_species_mineral[species.second], false);
243
244 for (const auto & species : mgd.basis_species_index)
245 ASSERT_EQ(mgd.basis_species_gas[species.second], false);
246 for (const auto & species : mgd.eqm_species_index)
247 ASSERT_EQ(mgd.eqm_species_gas[species.second], false);
248
249 for (const auto & species : mgd.basis_species_index)
250 ASSERT_EQ(mgd.basis_species_transported[species.second], true);
251 for (const auto & species : mgd.eqm_species_index)
252 if (species.first == "Calcite")
253 ASSERT_EQ(mgd.eqm_species_transported[species.second], false);
254 else
255 ASSERT_EQ(mgd.eqm_species_transported[species.second], true);
256
257 ASSERT_EQ(mgd.eqm_species_index.size(), (std::size_t)6);
258 for (const auto & sp : {"CO2(aq)", "CO3--", "CaCO3", "CaOH+", "OH-", "Calcite"})
259 ASSERT_EQ(mgd.eqm_species_index.count(sp), (std::size_t)1);
260
261 std::map<std::string, Real> charge_gold;
262 charge_gold["H2O"] = 0.0;
263 charge_gold["H+"] = 1.0;
264 charge_gold["HCO3-"] = -1.0;
265 charge_gold["Ca++"] = 2.0;
266 charge_gold["CO2(aq)"] = 0.0;
267 charge_gold["CO3--"] = -2.0;
268 charge_gold["CaCO3"] = 0.0;
269 charge_gold["CaOH+"] = 1.0;
270 charge_gold["OH-"] = -1.0;
271 charge_gold["Calcite"] = 0.0;
272 for (const auto & sp : mgd.basis_species_index)
273 ASSERT_EQ(mgd.basis_species_charge[sp.second], charge_gold[sp.first]);
274 for (const auto & sp : mgd.eqm_species_index)
275 ASSERT_EQ(mgd.eqm_species_charge[sp.second], charge_gold[sp.first]);
276
277 std::map<std::string, Real> radius_gold;
278 radius_gold["H2O"] = 0.0;
279 radius_gold["H+"] = 9.0;
280 radius_gold["HCO3-"] = 4.5;
281 radius_gold["Ca++"] = 6.0;
282 radius_gold["CO2(aq)"] = 4.0;
283 radius_gold["CO3--"] = 4.5;
284 radius_gold["CaCO3"] = 4.0;
285 radius_gold["CaOH+"] = 4.0;
286 radius_gold["OH-"] = 3.5;
287 radius_gold["Calcite"] = 0.0;
288 for (const auto & sp : mgd.basis_species_index)
289 ASSERT_EQ(mgd.basis_species_radius[sp.second], radius_gold[sp.first]);
290 for (const auto & sp : mgd.eqm_species_index)
291 ASSERT_EQ(mgd.eqm_species_radius[sp.second], radius_gold[sp.first]);
292
293 std::map<std::string, Real> molecular_weight_gold;
294 molecular_weight_gold["H2O"] = 18.0152;
295 molecular_weight_gold["H+"] = 1.0079;
296 molecular_weight_gold["HCO3-"] = 61.0171;
297 molecular_weight_gold["Ca++"] = 40.0800;
298 molecular_weight_gold["CO2(aq)"] = 44.0098;
299 molecular_weight_gold["CO3--"] = 60.0092;
300 molecular_weight_gold["CaCO3"] = 100.0892;
301 molecular_weight_gold["CaOH+"] = 57.0873;
302 molecular_weight_gold["OH-"] = 17.0073;
303 molecular_weight_gold["Calcite"] = 100.0892;
304 for (const auto & sp : mgd.basis_species_index)
305 ASSERT_EQ(mgd.basis_species_molecular_weight[sp.second], molecular_weight_gold[sp.first]);
306 for (const auto & sp : mgd.eqm_species_index)
307 ASSERT_EQ(mgd.eqm_species_molecular_weight[sp.second], molecular_weight_gold[sp.first]);
308
309 std::map<std::string, Real> molecular_volume_gold;
310 molecular_volume_gold["H2O"] = 0.0;
311 molecular_volume_gold["H+"] = 0.0;
312 molecular_volume_gold["HCO3-"] = 0.0;
313 molecular_volume_gold["O2(aq)"] = 0.0;
314 molecular_volume_gold["Ca++"] = 0.0;
315 molecular_volume_gold["CO2(aq)"] = 0.0;
316 molecular_volume_gold["CO3--"] = 0.0;
317 molecular_volume_gold["CaCO3"] = 0.0;
318 molecular_volume_gold["CaOH+"] = 0.0;
319 molecular_volume_gold["OH-"] = 0.0;
320 molecular_volume_gold["Calcite"] = 36.9340;
321 for (const auto & sp : mgd.basis_species_index)
322 ASSERT_EQ(mgd.basis_species_molecular_volume[sp.second], molecular_volume_gold[sp.first]);
323 for (const auto & sp : mgd.eqm_species_index)
324 ASSERT_EQ(mgd.eqm_species_molecular_volume[sp.second], molecular_volume_gold[sp.first]);
325
326 std::map<std::string, DenseMatrix<Real>> stoi_gold;
327 for (const auto & sp : {"CO2(aq)", "CO3--", "CaCO3", "CaOH+", "OH-", "Calcite"})
328 stoi_gold[sp] = DenseMatrix<Real>(1, 4);
329 // remember the order of primaries: {"H2O", "Ca++", "HCO3-", "H+"}
330 stoi_gold["CO2(aq)"](0, 0) = -1;
331 stoi_gold["CO2(aq)"](0, 3) = 1;
332 stoi_gold["CO2(aq)"](0, 2) = 1;
333 stoi_gold["CO3--"](0, 2) = 1;
334 stoi_gold["CO3--"](0, 3) = -1;
335 stoi_gold["CaCO3"](0, 1) = 1;
336 stoi_gold["CaCO3"](0, 2) = 1;
337 stoi_gold["CaCO3"](0, 3) = -1;
338 stoi_gold["CaOH+"](0, 1) = 1;
339 stoi_gold["CaOH+"](0, 0) = 1;
340 stoi_gold["CaOH+"](0, 3) = -1;
341 stoi_gold["OH-"](0, 0) = 1;
342 stoi_gold["OH-"](0, 3) = -1;
343 stoi_gold["Calcite"](0, 1) = 1;
344 stoi_gold["Calcite"](0, 2) = 1;
345 stoi_gold["Calcite"](0, 3) = -1;
346
347 for (const auto & sp : {"CO2(aq)", "CO3--", "CaCO3", "CaOH+", "OH-", "Calcite"})
348 {
349 const unsigned row = mgd.eqm_species_index[sp];
350 ASSERT_EQ(mgd.eqm_stoichiometry.sub_matrix(row, 1, 0, 4), stoi_gold[sp]);
351 }
352
353 ASSERT_EQ(mgd.eqm_log10K(mgd.eqm_species_index["CO2(aq)"], 0), -6.5570);
354 ASSERT_EQ(mgd.eqm_log10K(mgd.eqm_species_index["CO2(aq)"], 1), -6.3660);
355 ASSERT_EQ(mgd.eqm_log10K(mgd.eqm_species_index["CO2(aq)"], 2), -6.3325);
356 ASSERT_EQ(mgd.eqm_log10K(mgd.eqm_species_index["CO2(aq)"], 3), -6.4330);
357 ASSERT_EQ(mgd.eqm_log10K(mgd.eqm_species_index["CO2(aq)"], 4), -6.7420);
358 ASSERT_EQ(mgd.eqm_log10K(mgd.eqm_species_index["CO2(aq)"], 5), -7.1880);
359 ASSERT_EQ(mgd.eqm_log10K(mgd.eqm_species_index["CO2(aq)"], 6), -7.7630);
360 ASSERT_EQ(mgd.eqm_log10K(mgd.eqm_species_index["CO2(aq)"], 7), -8.4650);
361 ASSERT_EQ(mgd.eqm_log10K(mgd.eqm_species_index["CO3--"], 0), 10.6169);
362 ASSERT_EQ(mgd.eqm_log10K(mgd.eqm_species_index["CaCO3"], 0), 7.5520);
363 ASSERT_EQ(mgd.eqm_log10K(mgd.eqm_species_index["CaOH+"], 0), 13.7095);
364 ASSERT_EQ(mgd.eqm_log10K(mgd.eqm_species_index["OH-"], 0), 14.9325);
365 ASSERT_EQ(mgd.eqm_log10K(mgd.eqm_species_index["Calcite"], 0), 2.0683);
366
367 const unsigned ca_posn = mgd.basis_species_index["Ca++"];
368 const unsigned calcite_posn = mgd.eqm_species_index["Calcite"];
369
370 swapper.performSwap(mgd, bulk_composition, "Ca++", "Calcite");
371
372 // check names are swapped correctly
373 ASSERT_EQ(mgd.basis_species_index.size(), (std::size_t)4);
374 for (const auto & sp : {"Calcite", "H2O", "HCO3-", "H+"})
375 ASSERT_EQ(mgd.basis_species_index.count(sp), (std::size_t)1);
376 for (const auto & sp : mgd.basis_species_index)
377 ASSERT_EQ(mgd.basis_species_name[sp.second], sp.first);
378 ASSERT_EQ(mgd.basis_species_index["Calcite"], ca_posn);
379
380 // check names are swapped correctly
381 ASSERT_EQ(mgd.eqm_species_index.size(), (std::size_t)6);
382 for (const auto & sp : {"CO2(aq)", "CO3--", "CaCO3", "CaOH+", "OH-", "Ca++"})
383 ASSERT_EQ(mgd.eqm_species_index.count(sp), (std::size_t)1);
384 for (const auto & sp : mgd.eqm_species_index)
385 ASSERT_EQ(mgd.eqm_species_name[sp.second], sp.first);
386 ASSERT_EQ(mgd.eqm_species_index["Ca++"], calcite_posn);
387
388 // check the swap is recorded correctly
389 ASSERT_EQ(mgd.have_swapped_out_of_basis.size(), (std::size_t)1);
390 ASSERT_EQ(mgd.have_swapped_into_basis.size(), (std::size_t)1);
391 ASSERT_EQ(mgd.have_swapped_out_of_basis[0], ca_posn);
392 ASSERT_EQ(mgd.have_swapped_into_basis[0], calcite_posn);
393 ASSERT_EQ(mgd.swap_to_original_basis.n(), (unsigned)4);
394 DenseMatrix<Real> gold_swap_matrix(4, 4);
395 for (unsigned i = 0; i < 4; ++i)
396 gold_swap_matrix(i, i) = 1.0;
397 gold_swap_matrix(ca_posn, 0) = 0.0; // H2O
398 gold_swap_matrix(ca_posn, 1) = 1.0; // Ca++
399 gold_swap_matrix(ca_posn, 2) = 1.0; // HCO3-
400 gold_swap_matrix(ca_posn, 3) = -1.0; // H+
401 for (unsigned i = 0; i < 4; ++i)
402 for (unsigned j = 0; j < 4; ++j)
403 EXPECT_EQ(mgd.swap_to_original_basis(i, j), gold_swap_matrix(i, j));
404
405 // check charges swapped correctly
406 for (const auto & sp : mgd.basis_species_index)
407 ASSERT_EQ(mgd.basis_species_charge[sp.second], charge_gold[sp.first]);
408 for (const auto & sp : mgd.eqm_species_index)
409 ASSERT_EQ(mgd.eqm_species_charge[sp.second], charge_gold[sp.first]);
410
411 // check radii swapped correctly
412 for (const auto & sp : mgd.basis_species_index)
413 ASSERT_EQ(mgd.basis_species_radius[sp.second], radius_gold[sp.first]);
414 for (const auto & sp : mgd.eqm_species_index)
415 ASSERT_EQ(mgd.eqm_species_radius[sp.second], radius_gold[sp.first]);
416
417 // check molecular weights swapped correctly
418 for (const auto & sp : mgd.basis_species_index)
419 ASSERT_EQ(mgd.basis_species_molecular_weight[sp.second], molecular_weight_gold[sp.first]);
420 for (const auto & sp : mgd.eqm_species_index)
421 ASSERT_EQ(mgd.eqm_species_molecular_weight[sp.second], molecular_weight_gold[sp.first]);
422
423 // check molecular volumes swapped correctly
424 for (const auto & sp : mgd.basis_species_index)
425 ASSERT_EQ(mgd.basis_species_molecular_volume[sp.second], molecular_volume_gold[sp.first]);
426 for (const auto & sp : mgd.eqm_species_index)
427 ASSERT_EQ(mgd.eqm_species_molecular_volume[sp.second], molecular_volume_gold[sp.first]);
428
429 // check isMineral information has swapped correctly
430 for (const auto & species : mgd.basis_species_index)
431 if (species.first == "Calcite")
432 ASSERT_EQ(mgd.basis_species_mineral[species.second], true);
433 else
434 ASSERT_EQ(mgd.basis_species_mineral[species.second], false);
435 for (const auto & species : mgd.eqm_species_index)
436 ASSERT_EQ(mgd.eqm_species_mineral[species.second], false);
437
438 // check isGas information has swapped correctly
439 for (const auto & species : mgd.basis_species_index)
440 ASSERT_EQ(mgd.basis_species_gas[species.second], false);
441 for (const auto & species : mgd.eqm_species_index)
442 ASSERT_EQ(mgd.eqm_species_gas[species.second], false);
443
444 // check isTransported information has swapped correctly
445 for (const auto & species : mgd.basis_species_index)
446 if (species.first == "Calcite")
447 ASSERT_EQ(mgd.basis_species_transported[species.second], false);
448 else
449 ASSERT_EQ(mgd.basis_species_transported[species.second], true);
450 for (const auto & species : mgd.eqm_species_index)
451 ASSERT_EQ(mgd.eqm_species_transported[species.second], true);
452
453 // check stoichiometry
454 stoi_gold = std::map<std::string, DenseMatrix<Real>>();
455 for (const auto & sp : {"CO2(aq)", "CO3--", "CaCO3", "CaOH+", "OH-", "Ca++"})
456 stoi_gold[sp] = DenseMatrix<Real>(1, 4);
457 // remember the order of primaries: {"H2O", "Calcite", "HCO3-", "H+"}
458 stoi_gold["CO2(aq)"](0, 0) = -1;
459 stoi_gold["CO2(aq)"](0, 3) = 1;
460 stoi_gold["CO2(aq)"](0, 2) = 1;
461 stoi_gold["CO3--"](0, 2) = 1;
462 stoi_gold["CO3--"](0, 3) = -1;
463 stoi_gold["CaCO3"](0, 1) = 1;
464 stoi_gold["CaOH+"](0, 1) = 1;
465 stoi_gold["CaOH+"](0, 0) = 1;
466 stoi_gold["CaOH+"](0, 2) = -1;
467 stoi_gold["OH-"](0, 0) = 1;
468 stoi_gold["OH-"](0, 3) = -1;
469 stoi_gold["Ca++"](0, 1) = 1;
470 stoi_gold["Ca++"](0, 2) = -1;
471 stoi_gold["Ca++"](0, 3) = 1;
472
473 for (const auto & sp : {"CO2(aq)", "CO3--", "CaCO3", "CaOH+", "OH-", "Ca++"})
474 {
475 const unsigned row = mgd.eqm_species_index[sp];
476 for (unsigned i = 0; i < 4; ++i)
477 ASSERT_NEAR(mgd.eqm_stoichiometry(row, i), stoi_gold[sp](0, i), eps);
478 }
479
480 // check eqm constants
481 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CO2(aq)"], 0), -6.5570, eps);
482 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CO2(aq)"], 1), -6.3660, eps);
483 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CO2(aq)"], 2), -6.3325, eps);
484 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CO2(aq)"], 3), -6.4330, eps);
485 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CO2(aq)"], 4), -6.7420, eps);
486 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CO2(aq)"], 5), -7.1880, eps);
487 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CO2(aq)"], 6), -7.7630, eps);
488 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CO2(aq)"], 7), -8.4650, eps);
489 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CO3--"], 0), 10.6169 - 0 * 2.0683, eps);
490 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CO3--"], 1), 10.3439 - 0 * 1.7130, eps);
491 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CO3--"], 2), 10.2092 - 0 * 1.2133, eps);
492 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CO3--"], 3), 10.2793 - 0 * 0.6871, eps);
493 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CO3--"], 4), 10.5131 - 0 * 0.0762, eps);
494 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CO3--"], 5), 10.8637 - 0 * (-0.5349), eps);
495 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CO3--"], 6), 11.2860 - 0 * (-1.2301), eps);
496 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CO3--"], 7), 11.6319 - 0 * (-2.2107), eps);
497 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CaCO3"], 0), 7.5520 - 1 * 2.0683, eps);
498 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CaCO3"], 1), 7.1280 - 1 * 1.7130, eps);
499 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CaCO3"], 2), 6.7340 - 1 * 1.2133, eps);
500 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CaCO3"], 3), 6.4350 - 1 * 0.6871, eps);
501 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CaCO3"], 4), 6.1810 - 1 * 0.0762, eps);
502 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CaCO3"], 5), 5.9320 - 1 * (-0.5349), eps);
503 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CaCO3"], 6), 5.5640 - 1 * (-1.2301), eps);
504 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CaCO3"], 7), 4.7890 - 1 * (-2.2107), eps);
505 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CaOH+"], 0), 13.7095 - 1 * 2.0683, eps);
506 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CaOH+"], 1), 12.6887 - 1 * 1.7130, eps);
507 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CaOH+"], 2), 11.5069 - 1 * 1.2133, eps);
508 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CaOH+"], 3), 10.4366 - 1 * 0.6871, eps);
509 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CaOH+"], 4), 9.3958 - 1 * 0.0762, eps);
510 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CaOH+"], 5), 8.5583 - 1 * (-0.5349), eps);
511 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CaOH+"], 6), 7.8155 - 1 * (-1.2301), eps);
512 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CaOH+"], 7), 7.0306 - 1 * (-2.2107), eps);
513 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["OH-"], 0), 14.9325 - 0 * 2.0683, eps);
514 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["OH-"], 1), 13.9868 - 0 * 1.7130, eps);
515 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["OH-"], 2), 13.0199 - 0 * 1.2133, eps);
516 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["OH-"], 3), 12.2403 - 0 * 0.6871, eps);
517 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["OH-"], 4), 11.5940 - 0 * 0.0762, eps);
518 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["OH-"], 5), 11.2191 - 0 * (-0.5349), eps);
519 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["OH-"], 6), 11.0880 - 0 * (-1.2301), eps);
520 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["OH-"], 7), 1001.2844 - 0 * (-2.2107), eps);
521 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["Ca++"], 0), 0 - 1 * 2.0683, eps);
522 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["Ca++"], 1), 0 - 1 * 1.7130, eps);
523 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["Ca++"], 2), 0 - 1 * 1.2133, eps);
524 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["Ca++"], 3), 0 - 1 * 0.6871, eps);
525 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["Ca++"], 4), 0 - 1 * 0.0762, eps);
526 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["Ca++"], 5), 0 - 1 * (-0.5349), eps);
527 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["Ca++"], 6), 0 - 1 * (-1.2301), eps);
528 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["Ca++"], 7), 0 - 1 * (-2.2107), eps);
529
530 // check bulk composition
531 ASSERT_NEAR(bulk_composition(0), 0.5, eps);
532 ASSERT_NEAR(bulk_composition(1), 1.0, eps);
533 ASSERT_NEAR(bulk_composition(2), 1.5, eps);
534 ASSERT_NEAR(bulk_composition(3), 4.0, eps);
535
536 // swap back, and check that swap_to_original_basis is the identity
537 swapper.performSwap(mgd, bulk_composition, "Calcite", "Ca++");
538 ASSERT_EQ(mgd.swap_to_original_basis.n(), (unsigned)4);
539 for (unsigned i = 0; i < 4; ++i)
540 for (unsigned j = 0; j < 4; ++j)
541 if (i == j)
542 EXPECT_NEAR(mgd.swap_to_original_basis(i, j), 1.0, 1.0E-6);
543 else
544 EXPECT_NEAR(mgd.swap_to_original_basis(i, j), 0.0, 1.0E-6);
545}
546
548TEST(GeochemistrySpeciesSwapperTest, swap2)
549{
550 GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
551
552 // eqm species are: OH-, CO2(aq), CO3--, StoiCheckRedox, StoiCheckGas
553 PertinentGeochemicalSystem model(database,
554 {"H2O", "StoiCheckBasis", "H+", "HCO3-"},
555 {},
556 {"StoiCheckGas"},
557 {},
558 {},
559 {},
560 "O2(aq)",
561 "e-");
562 ModelGeochemicalDatabase mgd = model.modelGeochemicalDatabase();
563 GeochemistrySpeciesSwapper swapper(mgd.basis_species_index.size(), 1E-6);
564
565 ASSERT_EQ(mgd.basis_species_index.size(), (std::size_t)4);
566 for (const auto & sp : mgd.basis_species_index)
567 ASSERT_EQ(mgd.basis_species_name[sp.second], sp.first);
568
569 ASSERT_EQ(mgd.eqm_species_index.size(), (std::size_t)5);
570 for (const auto & sp : mgd.eqm_species_index)
571 ASSERT_EQ(mgd.eqm_species_name[sp.second], sp.first);
572
573 ASSERT_EQ(mgd.swap_to_original_basis.n(), (unsigned)0);
574 ASSERT_EQ(mgd.have_swapped_out_of_basis.size(), (std::size_t)0);
575 ASSERT_EQ(mgd.have_swapped_into_basis.size(), (std::size_t)0);
576
577 for (const auto & species : mgd.basis_species_index)
578 ASSERT_EQ(mgd.basis_species_mineral[species.second], false);
579 for (const auto & species : mgd.eqm_species_index)
580 ASSERT_EQ(mgd.eqm_species_mineral[species.second], false);
581
582 for (const auto & species : mgd.basis_species_index)
583 ASSERT_EQ(mgd.basis_species_transported[species.second], true);
584 for (const auto & species : mgd.eqm_species_index)
585 ASSERT_EQ(mgd.eqm_species_transported[species.second], true);
586
587 ASSERT_EQ(mgd.eqm_species_index.size(), (std::size_t)5);
588 for (const auto & sp : {"OH-", "CO2(aq)", "CO3--", "StoiCheckRedox", "StoiCheckGas"})
589 ASSERT_EQ(mgd.eqm_species_index.count(sp), (std::size_t)1);
590
591 for (const auto & species : mgd.basis_species_index)
592 ASSERT_EQ(mgd.basis_species_gas[species.second], false);
593 for (const auto & species : mgd.eqm_species_index)
594 if (species.first == "StoiCheckGas")
595 ASSERT_EQ(mgd.eqm_species_gas[species.second], true);
596 else
597 ASSERT_EQ(mgd.eqm_species_gas[species.second], false);
598
599 std::map<std::string, Real> charge_gold;
600 charge_gold["H2O"] = 0.0;
601 charge_gold["StoiCheckBasis"] = 2.5;
602 charge_gold["H+"] = 1.0;
603 charge_gold["HCO3-"] = -1.0;
604 charge_gold["CO2(aq)"] = 0.0;
605 charge_gold["CO3--"] = -2.0;
606 charge_gold["OH-"] = -1.0;
607 charge_gold["StoiCheckRedox"] = 3.3;
608 charge_gold["StoiCheckGas"] = 0.0;
609 for (const auto & sp : mgd.basis_species_index)
610 ASSERT_EQ(mgd.basis_species_charge[sp.second], charge_gold[sp.first]);
611 for (const auto & sp : mgd.eqm_species_index)
612 ASSERT_EQ(mgd.eqm_species_charge[sp.second], charge_gold[sp.first]);
613
614 std::map<std::string, Real> radius_gold;
615 radius_gold["H2O"] = 0.0;
616 radius_gold["StoiCheckBasis"] = 6.54;
617 radius_gold["H+"] = 9.0;
618 radius_gold["HCO3-"] = 4.5;
619 radius_gold["CO2(aq)"] = 4.0;
620 radius_gold["CO3--"] = 4.5;
621 radius_gold["OH-"] = 3.5;
622 radius_gold["StoiCheckRedox"] = 9.9;
623 radius_gold["StoiCheckGas"] = 0.0;
624 for (const auto & sp : mgd.basis_species_index)
625 ASSERT_EQ(mgd.basis_species_radius[sp.second], radius_gold[sp.first]);
626 for (const auto & sp : mgd.eqm_species_index)
627 ASSERT_EQ(mgd.eqm_species_radius[sp.second], radius_gold[sp.first]);
628
629 std::map<std::string, Real> molecular_weight_gold;
630 molecular_weight_gold["H2O"] = 18.0152;
631 molecular_weight_gold["StoiCheckBasis"] = 55.8470;
632 molecular_weight_gold["H+"] = 1.0079;
633 molecular_weight_gold["HCO3-"] = 61.0171;
634 molecular_weight_gold["CO2(aq)"] = 44.0098;
635 molecular_weight_gold["CO3--"] = 60.0092;
636 molecular_weight_gold["OH-"] = 17.0073;
637 molecular_weight_gold["StoiCheckRedox"] = 55.8470;
638 molecular_weight_gold["StoiCheckGas"] = 28.0134;
639 for (const auto & sp : mgd.basis_species_index)
640 ASSERT_EQ(mgd.basis_species_molecular_weight[sp.second], molecular_weight_gold[sp.first]);
641 for (const auto & sp : mgd.eqm_species_index)
642 ASSERT_EQ(mgd.eqm_species_molecular_weight[sp.second], molecular_weight_gold[sp.first]);
643
644 std::map<std::string, Real> molecular_volume_gold;
645 molecular_volume_gold["H2O"] = 0.0;
646 molecular_volume_gold["StoiCheckBasis"] = 0.0;
647 molecular_volume_gold["H+"] = 0.0;
648 molecular_volume_gold["HCO3-"] = 0.0;
649 molecular_volume_gold["CO2(aq)"] = 0.0;
650 molecular_volume_gold["CO3--"] = 0.0;
651 molecular_volume_gold["OH-"] = 0.0;
652 molecular_volume_gold["StoiCheckRedox"] = 0.0;
653 molecular_volume_gold["StoiCheckGas"] = 0.0;
654 for (const auto & sp : mgd.basis_species_index)
655 ASSERT_EQ(mgd.basis_species_molecular_volume[sp.second], molecular_volume_gold[sp.first]);
656 for (const auto & sp : mgd.eqm_species_index)
657 ASSERT_EQ(mgd.eqm_species_molecular_volume[sp.second], molecular_volume_gold[sp.first]);
658
659 std::map<std::string, DenseMatrix<Real>> stoi_gold;
660 for (const auto & sp : {"OH-", "CO2(aq)", "CO3--", "StoiCheckRedox", "StoiCheckGas"})
661 stoi_gold[sp] = DenseMatrix<Real>(1, 4);
662 // remember the order of primaries: {"H2O", "StoiCheckBasis", "H+", "HCO3-"}
663 stoi_gold["CO2(aq)"](0, 0) = -1;
664 stoi_gold["CO2(aq)"](0, 2) = 1;
665 stoi_gold["CO2(aq)"](0, 3) = 1;
666 stoi_gold["CO3--"](0, 2) = -1;
667 stoi_gold["CO3--"](0, 3) = 1;
668 stoi_gold["OH-"](0, 0) = 1;
669 stoi_gold["OH-"](0, 2) = -1;
670 stoi_gold["StoiCheckRedox"](0, 0) = -0.5;
671 stoi_gold["StoiCheckRedox"](0, 1) = 1.5;
672 stoi_gold["StoiCheckRedox"](0, 2) = -1;
673 stoi_gold["StoiCheckGas"](0, 0) = 2.0;
674 stoi_gold["StoiCheckGas"](0, 1) = 3.0;
675 stoi_gold["StoiCheckGas"](0, 2) = -5.0;
676
677 for (const auto & sp : {"OH-", "CO2(aq)", "CO3--", "StoiCheckRedox", "StoiCheckGas"})
678 {
679 const unsigned row = mgd.eqm_species_index[sp];
680 ASSERT_EQ(mgd.eqm_stoichiometry.sub_matrix(row, 1, 0, 4), stoi_gold[sp]);
681 }
682
683 ASSERT_EQ(mgd.eqm_log10K(mgd.eqm_species_index["CO2(aq)"], 0), -6.5570);
684 ASSERT_EQ(mgd.eqm_log10K(mgd.eqm_species_index["CO3--"], 0), 10.6169);
685 ASSERT_EQ(mgd.eqm_log10K(mgd.eqm_species_index["OH-"], 0), 14.9325);
686 ASSERT_EQ(mgd.eqm_log10K(mgd.eqm_species_index["StoiCheckRedox"], 0), -10.0553);
687 ASSERT_EQ(mgd.eqm_log10K(mgd.eqm_species_index["StoiCheckGas"], 0), -2.9620 + 2 * (-10.0553));
688
689 const unsigned sc_basis_posn = mgd.basis_species_index["StoiCheckBasis"];
690 const unsigned sc_gas_posn = mgd.eqm_species_index["StoiCheckGas"];
691
692 swapper.performSwap(mgd, "StoiCheckBasis", "StoiCheckGas");
693
694 // check names are swapped correctly
695 ASSERT_EQ(mgd.basis_species_index.size(), (std::size_t)4);
696 for (const auto & sp : {"H2O", "StoiCheckGas", "H+", "HCO3-"})
697 ASSERT_EQ(mgd.basis_species_index.count(sp), (std::size_t)1);
698 for (const auto & sp : mgd.basis_species_index)
699 ASSERT_EQ(mgd.basis_species_name[sp.second], sp.first);
700 ASSERT_EQ(mgd.basis_species_index["StoiCheckGas"], sc_basis_posn);
701
702 // check names are swapped correctly
703 ASSERT_EQ(mgd.eqm_species_index.size(), (std::size_t)5);
704 for (const auto & sp : {"OH-", "CO2(aq)", "CO3--", "StoiCheckRedox", "StoiCheckBasis"})
705 ASSERT_EQ(mgd.eqm_species_index.count(sp), (std::size_t)1);
706 for (const auto & sp : mgd.eqm_species_index)
707 ASSERT_EQ(mgd.eqm_species_name[sp.second], sp.first);
708 ASSERT_EQ(mgd.eqm_species_index["StoiCheckBasis"], sc_gas_posn);
709
710 // check the swap is recorded correctly
711 ASSERT_EQ(mgd.swap_to_original_basis.n(), (unsigned)4);
712 ASSERT_EQ(mgd.have_swapped_out_of_basis.size(), (std::size_t)1);
713 ASSERT_EQ(mgd.have_swapped_into_basis.size(), (std::size_t)1);
714 ASSERT_EQ(mgd.have_swapped_out_of_basis[0], sc_basis_posn);
715 ASSERT_EQ(mgd.have_swapped_into_basis[0], sc_gas_posn);
716
717 // check charges swapped correctly
718 for (const auto & sp : mgd.basis_species_index)
719 ASSERT_EQ(mgd.basis_species_charge[sp.second], charge_gold[sp.first]);
720 for (const auto & sp : mgd.eqm_species_index)
721 ASSERT_EQ(mgd.eqm_species_charge[sp.second], charge_gold[sp.first]);
722
723 // check radii swapped correctly
724 for (const auto & sp : mgd.basis_species_index)
725 ASSERT_EQ(mgd.basis_species_radius[sp.second], radius_gold[sp.first]);
726 for (const auto & sp : mgd.eqm_species_index)
727 ASSERT_EQ(mgd.eqm_species_radius[sp.second], radius_gold[sp.first]);
728
729 // check molecular weights swapped correctly
730 for (const auto & sp : mgd.basis_species_index)
731 ASSERT_EQ(mgd.basis_species_molecular_weight[sp.second], molecular_weight_gold[sp.first]);
732 for (const auto & sp : mgd.eqm_species_index)
733 ASSERT_EQ(mgd.eqm_species_molecular_weight[sp.second], molecular_weight_gold[sp.first]);
734
735 // check molecular volumes swapped correctly
736 for (const auto & sp : mgd.basis_species_index)
737 ASSERT_EQ(mgd.basis_species_molecular_volume[sp.second], molecular_volume_gold[sp.first]);
738 for (const auto & sp : mgd.eqm_species_index)
739 ASSERT_EQ(mgd.eqm_species_molecular_volume[sp.second], molecular_volume_gold[sp.first]);
740
741 // check isMineral information has swapped correctly
742 for (const auto & species : mgd.basis_species_index)
743 ASSERT_EQ(mgd.basis_species_mineral[species.second], false);
744 for (const auto & species : mgd.eqm_species_index)
745 ASSERT_EQ(mgd.eqm_species_mineral[species.second], false);
746
747 // check isGas information has swapped correctly
748 for (const auto & species : mgd.basis_species_index)
749 if (species.first == "StoiCheckGas")
750 ASSERT_EQ(mgd.basis_species_gas[species.second], true);
751 else
752 ASSERT_EQ(mgd.basis_species_gas[species.second], false);
753 for (const auto & species : mgd.eqm_species_index)
754 ASSERT_EQ(mgd.eqm_species_gas[species.second], false);
755
756 // check isTransported information has swapped correctly
757 for (const auto & species : mgd.basis_species_index)
758 ASSERT_EQ(mgd.basis_species_transported[species.second], true);
759 for (const auto & species : mgd.eqm_species_index)
760 ASSERT_EQ(mgd.eqm_species_transported[species.second], true);
761
762 // check stoichiometry
763 stoi_gold = std::map<std::string, DenseMatrix<Real>>();
764 for (const auto & sp : {"OH-", "CO2(aq)", "CO3--", "StoiCheckRedox", "StoiCheckBasis"})
765 stoi_gold[sp] = DenseMatrix<Real>(1, 4);
766 // remember the order of primaries: {"H2O", "StoiCheckGas", "H+", "HCO3-"}
767 stoi_gold["CO2(aq)"](0, 0) = -1;
768 stoi_gold["CO2(aq)"](0, 2) = 1;
769 stoi_gold["CO2(aq)"](0, 3) = 1;
770 stoi_gold["CO3--"](0, 2) = -1;
771 stoi_gold["CO3--"](0, 3) = 1;
772 stoi_gold["OH-"](0, 0) = 1;
773 stoi_gold["OH-"](0, 2) = -1;
774 stoi_gold["StoiCheckRedox"](0, 0) = -1.5;
775 stoi_gold["StoiCheckRedox"](0, 1) = 0.5;
776 stoi_gold["StoiCheckRedox"](0, 2) = 1.5;
777 stoi_gold["StoiCheckBasis"](0, 0) = -2.0 / 3.0;
778 stoi_gold["StoiCheckBasis"](0, 1) = 1.0 / 3.0;
779 stoi_gold["StoiCheckBasis"](0, 2) = 5.0 / 3.0;
780
781 for (const auto & sp : {"OH-", "CO2(aq)", "CO3--", "StoiCheckRedox", "StoiCheckBasis"})
782 {
783 const unsigned row = mgd.eqm_species_index[sp];
784 for (unsigned i = 0; i < 4; ++i)
785 ASSERT_NEAR(mgd.eqm_stoichiometry(row, i), stoi_gold[sp](0, i), eps);
786 }
787
788 // check eqm constants
789 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CO2(aq)"], 0), -6.5570, eps);
790 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CO2(aq)"], 1), -6.3660, eps);
791 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CO2(aq)"], 2), -6.3325, eps);
792 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CO2(aq)"], 3), -6.4330, eps);
793 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CO2(aq)"], 4), -6.7420, eps);
794 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CO2(aq)"], 5), -7.1880, eps);
795 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CO2(aq)"], 6), -7.7630, eps);
796 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CO2(aq)"], 7), -8.4650, eps);
797 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CO3--"], 0), 10.6169, eps);
798 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CO3--"], 1), 10.3439, eps);
799 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CO3--"], 2), 10.2092, eps);
800 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CO3--"], 3), 10.2793, eps);
801 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CO3--"], 4), 10.5131, eps);
802 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CO3--"], 5), 10.8637, eps);
803 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CO3--"], 6), 11.2860, eps);
804 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["CO3--"], 7), 11.6319, eps);
805 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["OH-"], 0), 14.9325, eps);
806 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["OH-"], 1), 13.9868, eps);
807 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["OH-"], 2), 13.0199, eps);
808 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["OH-"], 3), 12.2403, eps);
809 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["OH-"], 4), 11.5940, eps);
810 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["OH-"], 5), 11.2191, eps);
811 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["OH-"], 6), 11.0880, eps);
812 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["OH-"], 7), 1001.2844, eps);
813 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["StoiCheckRedox"], 0), -0.5 * (-2.9620), eps);
814 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["StoiCheckRedox"], 1), -0.5 * (-3.1848), eps);
815 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["StoiCheckRedox"], 2), -0.5 * (-3.3320), eps);
816 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["StoiCheckRedox"], 3), -0.5 * (-3.2902), eps);
817 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["StoiCheckRedox"], 4), -0.5 * (-3.1631), eps);
818 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["StoiCheckRedox"], 5), -0.5 * (-2.9499), eps);
819 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["StoiCheckRedox"], 6), -0.5 * (-2.7827), eps);
820 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["StoiCheckRedox"], 7), -0.5 * (-2.3699), eps);
821 const Real ot = -1.0 / 3.0;
822 const Real tt = -2.0 / 3.0;
823 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["StoiCheckBasis"], 0),
824 tt * (-10.0553) + ot * (-2.9620),
825 eps);
826 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["StoiCheckBasis"], 1),
827 tt * (-8.4878) + ot * (-3.1848),
828 eps);
829 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["StoiCheckBasis"], 2),
830 tt * (-6.6954) + ot * (-3.3320),
831 eps);
832 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["StoiCheckBasis"], 3),
833 tt * (-5.0568) + ot * (-3.2902),
834 eps);
835 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["StoiCheckBasis"], 4),
836 tt * (-3.4154) + ot * (-3.1631),
837 eps);
838 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["StoiCheckBasis"], 5),
839 tt * (-2.0747) + ot * (-2.9499),
840 eps);
841 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["StoiCheckBasis"], 6),
842 tt * (-.8908) + ot * (-2.7827),
843 eps);
844 ASSERT_NEAR(mgd.eqm_log10K(mgd.eqm_species_index["StoiCheckBasis"], 7),
845 tt * (.2679) + ot * (-2.3699),
846 eps);
847}
848
850TEST(GeochemistrySpeciesSwapperTest, swap3)
851{
852 GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
853
854 // The following system has secondary species: CO2(aq), CO3--, OH-, CH4(aq), Fe+++
855 PertinentGeochemicalSystem model(database,
856 {"H2O", "H+", ">(s)FeOH", ">(w)FeOH", "Fe++", "HCO3-", "O2(aq)"},
857 {},
858 {"CH4(g)fake"},
859 {"Fe(OH)3(ppd)", "Fe(OH)3(ppd)fake"},
860 {"(O-phth)--"},
861 {">(s)FeO-"},
862 "O2(aq)",
863 "e-");
864 KineticRateUserDescription rate1("Fe(OH)3(ppd)",
865 1.0,
866 2.0,
867 true,
868 0.0,
869 0.0,
870 0.0,
871 {"H+", "O2(aq)", "CO3--"},
872 {1.1, 2.2, 3.3},
873 {-0.1, 0.2, -0.3},
874 {0.5, 0.25, 0.125},
875 5.0,
876 6.0,
877 7.0,
878 8.0,
879 DirectionChoiceEnum::BOTH,
880 "O2(aq)",
881 1.0,
882 -1.0,
883 0.0);
884 model.addKineticRate(rate1);
885 KineticRateUserDescription rate2("Fe(OH)3(ppd)",
886 -1.0,
887 -2.0,
888 false,
889 0.0,
890 0.0,
891 0.0,
892 {"O2(aq)", "Fe+++", "H2O"},
893 {-1.1, -2.2, -3.3},
894 {0.75, 0.875, 1.0},
895 {0.2, 0.3, 0.4},
896 -5.0,
897 -6.0,
898 -7.0,
899 -8.0,
900 DirectionChoiceEnum::BOTH,
901 "Fe++",
902 -2.0,
903 -1.0,
904 0.0);
905 model.addKineticRate(rate2);
906 KineticRateUserDescription rate3(">(s)FeO-",
907 1.1,
908 2.2,
909 false,
910 0.0,
911 0.0,
912 0.0,
913 {"CO2(aq)", "Fe+++", "Fe++"},
914 {1.25, 2.25, 3.25},
915 {0.1, -0.12, 1.23},
916 {0.5, 0.6, 0.7},
917 5.5,
918 -6.6,
919 -7.7,
920 -8.8,
921 DirectionChoiceEnum::BOTH,
922 "Fe+++",
923 2.0,
924 -1.0,
925 0.0);
926 model.addKineticRate(rate3);
927 ModelGeochemicalDatabase mgd = model.modelGeochemicalDatabase();
928 GeochemistrySpeciesSwapper swapper(mgd.basis_species_index.size(), 1E-6);
929
930 ASSERT_EQ(mgd.basis_species_index.size(), (std::size_t)7);
931 for (const auto & sp : mgd.basis_species_index)
932 ASSERT_EQ(mgd.basis_species_name[sp.second], sp.first);
933
934 ASSERT_EQ(mgd.eqm_species_index.size(), (std::size_t)6);
935 for (const auto & sp : mgd.eqm_species_index)
936 ASSERT_EQ(mgd.eqm_species_name[sp.second], sp.first);
937
938 ASSERT_EQ(mgd.kin_species_index.size(), (std::size_t)4);
939 for (const auto & sp : mgd.kin_species_index)
940 ASSERT_EQ(mgd.kin_species_name[sp.second], sp.first);
941
942 ASSERT_EQ(mgd.swap_to_original_basis.n(), (unsigned)0);
943 ASSERT_EQ(mgd.have_swapped_out_of_basis.size(), (std::size_t)0);
944 ASSERT_EQ(mgd.have_swapped_into_basis.size(), (std::size_t)0);
945
946 for (const auto & species : mgd.basis_species_index)
947 ASSERT_EQ(mgd.basis_species_mineral[species.second], false);
948 for (const auto & species : mgd.eqm_species_index)
949 ASSERT_EQ(mgd.eqm_species_mineral[species.second], false);
950 for (const auto & species : mgd.kin_species_index)
951 if (species.first == "Fe(OH)3(ppd)" || species.first == "Fe(OH)3(ppd)fake")
952 ASSERT_EQ(mgd.kin_species_mineral[species.second], true);
953 else
954 ASSERT_EQ(mgd.kin_species_mineral[species.second], false);
955
956 for (const auto & species : mgd.basis_species_index)
957 ASSERT_EQ(mgd.basis_species_gas[species.second], false);
958 for (const auto & species : mgd.eqm_species_index)
959 if (species.first == "CH4(g)fake")
960 ASSERT_EQ(mgd.eqm_species_gas[species.second], true);
961 else
962 ASSERT_EQ(mgd.eqm_species_gas[species.second], false);
963
964 for (const auto & species : mgd.basis_species_index)
965 if (species.first == ">(s)FeOH" || species.first == ">(w)FeOH")
966 ASSERT_EQ(mgd.basis_species_transported[species.second], false);
967 else
968 ASSERT_EQ(mgd.basis_species_transported[species.second], true);
969 for (const auto & species : mgd.eqm_species_index)
970 ASSERT_EQ(mgd.eqm_species_transported[species.second], true);
971 for (const auto & species : mgd.kin_species_index)
972 if (species.first == "Fe(OH)3(ppd)" || species.first == "Fe(OH)3(ppd)fake" ||
973 species.first == ">(s)FeO-")
974 ASSERT_EQ(mgd.kin_species_transported[species.second], false);
975 else
976 ASSERT_EQ(mgd.kin_species_transported[species.second], true);
977
978 ASSERT_EQ(mgd.eqm_species_index.size(), (std::size_t)6);
979 for (const auto & sp : {"CO2(aq)", "CO3--", "OH-", "CH4(aq)", "Fe+++", "CH4(g)fake"})
980 ASSERT_EQ(mgd.eqm_species_index.count(sp), (std::size_t)1);
981
982 ASSERT_EQ(mgd.kin_species_index.size(), (std::size_t)4);
983 for (const auto & sp : {"Fe(OH)3(ppd)", "Fe(OH)3(ppd)fake", "(O-phth)--", ">(s)FeO-"})
984 ASSERT_EQ(mgd.kin_species_index.count(sp), (std::size_t)1);
985
986 std::map<std::string, Real> charge_gold;
987 charge_gold["H2O"] = 0.0;
988 charge_gold["H+"] = 1.0;
989 charge_gold[">(s)FeOH"] = 0.0;
990 charge_gold[">(w)FeOH"] = 0.0;
991 charge_gold["Fe++"] = 2.0;
992 charge_gold["HCO3-"] = -1.0;
993 charge_gold["O2(aq)"] = 0.0;
994 charge_gold["CO2(aq)"] = 0.0;
995 charge_gold["CO3--"] = -2.0;
996 charge_gold["OH-"] = -1.0;
997 charge_gold["CH4(aq)"] = 0.0;
998 charge_gold["Fe+++"] = 3.0;
999 charge_gold["CH4(g)fake"] = 0.0;
1000 charge_gold["Fe(OH)3(ppd)"] = 0.0;
1001 charge_gold["Fe(OH)3(ppd)fake"] = 0.0;
1002 charge_gold["(O-phth)--"] = -2.0;
1003 charge_gold[">(s)FeO-"] = -1.0;
1004 for (const auto & sp : mgd.basis_species_index)
1005 ASSERT_EQ(mgd.basis_species_charge[sp.second], charge_gold[sp.first]);
1006 for (const auto & sp : mgd.eqm_species_index)
1007 ASSERT_EQ(mgd.eqm_species_charge[sp.second], charge_gold[sp.first]);
1008 for (const auto & sp : mgd.kin_species_index)
1009 ASSERT_EQ(mgd.kin_species_charge[sp.second], charge_gold[sp.first]);
1010
1011 std::map<std::string, Real> radius_gold;
1012 radius_gold["H2O"] = 0.0;
1013 radius_gold["H+"] = 9.0;
1014 radius_gold[">(s)FeOH"] = 0.0;
1015 radius_gold[">(w)FeOH"] = 0.0;
1016 radius_gold["Fe++"] = 6.0;
1017 radius_gold["HCO3-"] = 4.5;
1018 radius_gold["O2(aq)"] = -0.5;
1019 radius_gold["CO2(aq)"] = 4.0;
1020 radius_gold["CO3--"] = 4.5;
1021 radius_gold["OH-"] = 3.5;
1022 radius_gold["CH4(aq)"] = -0.5;
1023 radius_gold["Fe+++"] = 9.0;
1024 radius_gold["CH4(g)fake"] = 0.0;
1025 for (const auto & sp : mgd.basis_species_index)
1026 ASSERT_EQ(mgd.basis_species_radius[sp.second], radius_gold[sp.first]);
1027 for (const auto & sp : mgd.eqm_species_index)
1028 ASSERT_EQ(mgd.eqm_species_radius[sp.second], radius_gold[sp.first]);
1029
1030 std::map<std::string, Real> molecular_weight_gold;
1031 molecular_weight_gold["H2O"] = 18.0152;
1032 molecular_weight_gold["H+"] = 1.0079;
1033 molecular_weight_gold[">(s)FeOH"] = 72.8543;
1034 molecular_weight_gold[">(w)FeOH"] = 1234.567;
1035 molecular_weight_gold["Fe++"] = 55.8470;
1036 molecular_weight_gold["HCO3-"] = 61.0171;
1037 molecular_weight_gold["O2(aq)"] = 31.9988;
1038 molecular_weight_gold["CO2(aq)"] = 44.0098;
1039 molecular_weight_gold["CO3--"] = 60.0092;
1040 molecular_weight_gold["OH-"] = 17.0073;
1041 molecular_weight_gold["CH4(aq)"] = 16.0426;
1042 molecular_weight_gold["Fe+++"] = 55.8470;
1043 molecular_weight_gold["CH4(g)fake"] = 16.0426;
1044 molecular_weight_gold["Fe(OH)3(ppd)"] = 106.8689;
1045 molecular_weight_gold["Fe(OH)3(ppd)fake"] = 106.8689;
1046 molecular_weight_gold["(O-phth)--"] = 164.1172;
1047 molecular_weight_gold[">(s)FeO-"] = 71.8464;
1048 for (const auto & sp : mgd.basis_species_index)
1049 ASSERT_EQ(mgd.basis_species_molecular_weight[sp.second], molecular_weight_gold[sp.first]);
1050 for (const auto & sp : mgd.eqm_species_index)
1051 ASSERT_EQ(mgd.eqm_species_molecular_weight[sp.second], molecular_weight_gold[sp.first]);
1052 for (const auto & sp : mgd.kin_species_index)
1053 ASSERT_EQ(mgd.kin_species_molecular_weight[sp.second], molecular_weight_gold[sp.first]);
1054
1055 std::map<std::string, Real> molecular_volume_gold;
1056 molecular_volume_gold["H2O"] = 0.0;
1057 molecular_volume_gold["H+"] = 0.0;
1058 molecular_volume_gold[">(s)FeOH"] = 0.0;
1059 molecular_volume_gold[">(w)FeOH"] = 0.0;
1060 molecular_volume_gold["Fe++"] = 0.0;
1061 molecular_volume_gold["HCO3-"] = 0.0;
1062 molecular_volume_gold["O2(aq)"] = 0.0;
1063 molecular_volume_gold["CO2(aq)"] = 0.0;
1064 molecular_volume_gold["CO3--"] = 0.0;
1065 molecular_volume_gold["OH-"] = 0.0;
1066 molecular_volume_gold["CH4(aq)"] = 0.0;
1067 molecular_volume_gold["Fe+++"] = 0.0;
1068 molecular_volume_gold["CH4(g)fake"] = 0.0;
1069 molecular_volume_gold["Fe(OH)3(ppd)"] = 34.3200;
1070 molecular_volume_gold["Fe(OH)3(ppd)fake"] = 34.3200;
1071 molecular_volume_gold["(O-phth)--"] = 0.0;
1072 molecular_volume_gold[">(s)FeO-"] = 0.0;
1073 for (const auto & sp : mgd.basis_species_index)
1074 ASSERT_EQ(mgd.basis_species_molecular_volume[sp.second], molecular_volume_gold[sp.first]);
1075 for (const auto & sp : mgd.eqm_species_index)
1076 ASSERT_EQ(mgd.eqm_species_molecular_volume[sp.second], molecular_volume_gold[sp.first]);
1077 for (const auto & sp : mgd.kin_species_index)
1078 ASSERT_EQ(mgd.kin_species_molecular_volume[sp.second], molecular_volume_gold[sp.first]);
1079
1080 std::map<std::string, DenseMatrix<Real>> stoi_gold;
1081 for (const auto & sp : {"CO2(aq)",
1082 "CO3--",
1083 "OH-",
1084 "CH4(aq)",
1085 "Fe+++",
1086 "CH4(g)fake",
1087 "Fe(OH)3(ppd)",
1088 "Fe(OH)3(ppd)fake",
1089 "(O-phth)--",
1090 ">(s)FeO-"})
1091 stoi_gold[sp] = DenseMatrix<Real>(1, 7);
1092 // remember the order of primaries:
1093 // {"H2O", "H+", ">(s)FeOH", ">(w)FeOH", "Fe++", "HCO3-", "O2(aq)"}
1094 stoi_gold["CO2(aq)"](0, 0) = -1;
1095 stoi_gold["CO2(aq)"](0, 1) = 1;
1096 stoi_gold["CO2(aq)"](0, 5) = 1;
1097 stoi_gold["CO3--"](0, 5) = 1;
1098 stoi_gold["CO3--"](0, 1) = -1;
1099 stoi_gold["OH-"](0, 0) = 1;
1100 stoi_gold["OH-"](0, 1) = -1;
1101 stoi_gold["CH4(aq)"](0, 0) = 1;
1102 stoi_gold["CH4(aq)"](0, 1) = 1;
1103 stoi_gold["CH4(aq)"](0, 5) = 1;
1104 stoi_gold["CH4(aq)"](0, 6) = -2;
1105 stoi_gold["Fe+++"](0, 0) = -0.5;
1106 stoi_gold["Fe+++"](0, 4) = 1;
1107 stoi_gold["Fe+++"](0, 1) = 1;
1108 stoi_gold["Fe+++"](0, 6) = 0.25;
1109 stoi_gold["CH4(g)fake"](0, 0) = 3;
1110 stoi_gold["CH4(g)fake"](0, 4) = -2;
1111 stoi_gold["CH4(g)fake"](0, 5) = 3.5;
1112 stoi_gold["CH4(g)fake"](0, 6) = -4.5;
1113 stoi_gold["Fe(OH)3(ppd)"](0, 1) = -2;
1114 stoi_gold["Fe(OH)3(ppd)"](0, 4) = 1;
1115 stoi_gold["Fe(OH)3(ppd)"](0, 0) = 2.5;
1116 stoi_gold["Fe(OH)3(ppd)"](0, 6) = 0.25;
1117 stoi_gold["Fe(OH)3(ppd)fake"](0, 1) = -1;
1118 stoi_gold["Fe(OH)3(ppd)fake"](0, 0) = 2;
1119 stoi_gold["Fe(OH)3(ppd)fake"](0, 4) = 2;
1120 stoi_gold["Fe(OH)3(ppd)fake"](0, 6) = 0.5;
1121 stoi_gold["(O-phth)--"](0, 0) = -5;
1122 stoi_gold["(O-phth)--"](0, 5) = 8;
1123 stoi_gold["(O-phth)--"](0, 1) = 6;
1124 stoi_gold["(O-phth)--"](0, 6) = -7.5;
1125 stoi_gold[">(s)FeO-"](0, 2) = 1;
1126 stoi_gold[">(s)FeO-"](0, 1) = -1;
1127
1128 for (const auto & sp : {"CO2(aq)", "CO3--", "OH-", "CH4(aq)", "Fe+++", "CH4(g)fake"})
1129 {
1130 const unsigned row = mgd.eqm_species_index[sp];
1131 ASSERT_EQ(mgd.eqm_stoichiometry.sub_matrix(row, 1, 0, 7), stoi_gold[sp]);
1132 }
1133
1134 for (const auto & sp : {"Fe(OH)3(ppd)", "Fe(OH)3(ppd)fake", "(O-phth)--", ">(s)FeO-"})
1135 {
1136 const unsigned row = mgd.kin_species_index[sp];
1137 ASSERT_EQ(mgd.kin_stoichiometry.sub_matrix(row, 1, 0, 7), stoi_gold[sp]);
1138 }
1139
1140 const std::vector<Real> feoh3ppd = {
1141 6.1946, 4.8890, 3.4608, 2.2392, 1.1150, 0.2446, -0.5504, -1.5398};
1142 const std::vector<Real> fe3 = {10.0553, 8.4878, 6.6954, 5.0568, 3.4154, 2.0747, 0.8908, -0.2679};
1143 const std::vector<Real> ophth = {
1144 594.3211, 542.8292, 482.3612, 425.9738, 368.7004, 321.8658, 281.8216, 246.4849};
1145 const std::vector<Real> feom = {8.93,
1146 8.93 - 0.3 * 25,
1147 8.93 - 0.3 * 60,
1148 8.93 - 0.3 * 100,
1149 8.93 - 0.3 * 150,
1150 8.93 - 0.3 * 200,
1151 8.93 - 0.3 * 250,
1152 8.93 - 0.3 * 300};
1153
1154 for (unsigned t = 0; t < 8; ++t)
1155 {
1156 ASSERT_NEAR(
1157 mgd.kin_log10K(mgd.kin_species_index.at("Fe(OH)3(ppd)"), t), feoh3ppd[t] - fe3[t], eps);
1158 ASSERT_NEAR(mgd.kin_log10K(mgd.kin_species_index.at("Fe(OH)3(ppd)fake"), t),
1159 feoh3ppd[t] - 2 * fe3[t],
1160 eps);
1161 ASSERT_NEAR(mgd.kin_log10K(mgd.kin_species_index.at("(O-phth)--"), t), ophth[t], eps);
1162 ASSERT_NEAR(mgd.kin_log10K(mgd.kin_species_index.at(">(s)FeO-"), t), feom[t], eps);
1163 }
1164
1165 ASSERT_EQ(mgd.kin_rate.size(), (std::size_t)3);
1166 std::vector<std::vector<Real>> kin_rate_gold(3, std::vector<Real>(7 + 6, 0.0));
1167 std::vector<std::vector<Real>> kin_monod_gold(3, std::vector<Real>(7 + 6, 0.0));
1168 std::vector<std::vector<Real>> kin_k_gold(3, std::vector<Real>(7 + 6, 0.0));
1169 kin_rate_gold[0][mgd.basis_species_index.at("H+")] = 1.1;
1170 kin_rate_gold[0][mgd.basis_species_index.at("O2(aq)")] = 2.2;
1171 kin_rate_gold[0][7 + mgd.eqm_species_index.at("CO3--")] = 3.3;
1172 kin_monod_gold[0][mgd.basis_species_index.at("H+")] = -0.1;
1173 kin_monod_gold[0][mgd.basis_species_index.at("O2(aq)")] = 0.2;
1174 kin_monod_gold[0][7 + mgd.eqm_species_index.at("CO3--")] = -0.3;
1175 kin_k_gold[0][mgd.basis_species_index.at("H+")] = 0.5;
1176 kin_k_gold[0][mgd.basis_species_index.at("O2(aq)")] = 0.25;
1177 kin_k_gold[0][7 + mgd.eqm_species_index.at("CO3--")] = 0.125;
1178 kin_rate_gold[1][mgd.basis_species_index.at("O2(aq)")] = -1.1;
1179 kin_rate_gold[1][7 + mgd.eqm_species_index.at("Fe+++")] = -2.2;
1180 kin_rate_gold[1][mgd.basis_species_index.at("H2O")] = -3.3;
1181 kin_monod_gold[1][mgd.basis_species_index.at("O2(aq)")] = 0.75;
1182 kin_monod_gold[1][7 + mgd.eqm_species_index.at("Fe+++")] = 0.875;
1183 kin_monod_gold[1][mgd.basis_species_index.at("H2O")] = 1.0;
1184 kin_k_gold[1][mgd.basis_species_index.at("O2(aq)")] = 0.2;
1185 kin_k_gold[1][7 + mgd.eqm_species_index.at("Fe+++")] = 0.3;
1186 kin_k_gold[1][mgd.basis_species_index.at("H2O")] = 0.4;
1187 kin_rate_gold[2][7 + mgd.eqm_species_index.at("CO2(aq)")] = 1.25;
1188 kin_rate_gold[2][7 + mgd.eqm_species_index.at("Fe+++")] = 2.25;
1189 kin_rate_gold[2][mgd.basis_species_index.at("Fe++")] = 3.25;
1190 kin_monod_gold[2][7 + mgd.eqm_species_index.at("CO2(aq)")] = 0.1;
1191 kin_monod_gold[2][7 + mgd.eqm_species_index.at("Fe+++")] = -0.12;
1192 kin_monod_gold[2][mgd.basis_species_index.at("Fe++")] = 1.23;
1193 kin_k_gold[2][7 + mgd.eqm_species_index.at("CO2(aq)")] = 0.5;
1194 kin_k_gold[2][7 + mgd.eqm_species_index.at("Fe+++")] = 0.6;
1195 kin_k_gold[2][mgd.basis_species_index.at("Fe++")] = 0.7;
1196 for (unsigned i = 0; i < 3; ++i)
1197 for (unsigned j = 0; j < 7 + 6; ++j)
1198 {
1199 EXPECT_EQ(mgd.kin_rate[i].promoting_indices[j], kin_rate_gold[i][j]);
1200 EXPECT_EQ(mgd.kin_rate[i].promoting_monod_indices[j], kin_monod_gold[i][j]);
1201 EXPECT_EQ(mgd.kin_rate[i].promoting_half_saturation[j], kin_k_gold[i][j]);
1202 }
1203
1204 EXPECT_EQ(mgd.kin_rate[0].progeny_index, mgd.basis_species_index.at("O2(aq)"));
1205 EXPECT_EQ(mgd.kin_rate[1].progeny_index, mgd.basis_species_index.at("Fe++"));
1206 EXPECT_EQ(mgd.kin_rate[2].progeny_index, 7 + mgd.eqm_species_index.at("Fe+++"));
1207
1208 const unsigned o2aq_posn = mgd.basis_species_index["O2(aq)"];
1209 const unsigned fe3_posn = mgd.eqm_species_index["Fe+++"];
1210
1211 swapper.performSwap(mgd, "O2(aq)", "Fe+++");
1212
1213 // check names are swapped correctly
1214 ASSERT_EQ(mgd.basis_species_index.size(), (std::size_t)7);
1215 for (const auto & sp : {"H2O", "H+", ">(s)FeOH", ">(w)FeOH", "Fe++", "HCO3-", "Fe+++"})
1216 ASSERT_EQ(mgd.basis_species_index.count(sp), (std::size_t)1);
1217 for (const auto & sp : mgd.basis_species_index)
1218 ASSERT_EQ(mgd.basis_species_name[sp.second], sp.first);
1219 ASSERT_EQ(mgd.basis_species_index["Fe+++"], o2aq_posn);
1220
1221 // check names are swapped correctly
1222 ASSERT_EQ(mgd.eqm_species_index.size(), (std::size_t)6);
1223 for (const auto & sp : {"CO2(aq)", "CO3--", "OH-", "CH4(aq)", "O2(aq)", "CH4(g)fake"})
1224 ASSERT_EQ(mgd.eqm_species_index.count(sp), (std::size_t)1);
1225 for (const auto & sp : mgd.eqm_species_index)
1226 ASSERT_EQ(mgd.eqm_species_name[sp.second], sp.first);
1227 ASSERT_EQ(mgd.eqm_species_index["O2(aq)"], fe3_posn);
1228
1229 ASSERT_EQ(mgd.kin_species_index.size(), (std::size_t)4);
1230 for (const auto & sp : mgd.kin_species_index)
1231 ASSERT_EQ(mgd.kin_species_name[sp.second], sp.first);
1232 for (const auto & sp : {"Fe(OH)3(ppd)", "Fe(OH)3(ppd)fake", "(O-phth)--", ">(s)FeO-"})
1233 ASSERT_EQ(mgd.kin_species_index.count(sp), (std::size_t)1);
1234
1235 // check the swap is recorded correctly
1236 ASSERT_EQ(mgd.swap_to_original_basis.n(), (unsigned)7);
1237 ASSERT_EQ(mgd.have_swapped_out_of_basis.size(), (std::size_t)1);
1238 ASSERT_EQ(mgd.have_swapped_into_basis.size(), (std::size_t)1);
1239 ASSERT_EQ(mgd.have_swapped_out_of_basis[0], o2aq_posn);
1240 ASSERT_EQ(mgd.have_swapped_into_basis[0], fe3_posn);
1241
1242 // check charges swapped correctly
1243 for (const auto & sp : mgd.basis_species_index)
1244 ASSERT_EQ(mgd.basis_species_charge[sp.second], charge_gold[sp.first]);
1245 for (const auto & sp : mgd.eqm_species_index)
1246 ASSERT_EQ(mgd.eqm_species_charge[sp.second], charge_gold[sp.first]);
1247 for (const auto & sp : mgd.kin_species_index)
1248 ASSERT_EQ(mgd.kin_species_charge[sp.second], charge_gold[sp.first]);
1249
1250 // check radii swapped correctly
1251 for (const auto & sp : mgd.basis_species_index)
1252 ASSERT_EQ(mgd.basis_species_radius[sp.second], radius_gold[sp.first]);
1253 for (const auto & sp : mgd.eqm_species_index)
1254 ASSERT_EQ(mgd.eqm_species_radius[sp.second], radius_gold[sp.first]);
1255
1256 // check molecular weights swapped correctly
1257 for (const auto & sp : mgd.basis_species_index)
1258 ASSERT_EQ(mgd.basis_species_molecular_weight[sp.second], molecular_weight_gold[sp.first]);
1259 for (const auto & sp : mgd.eqm_species_index)
1260 ASSERT_EQ(mgd.eqm_species_molecular_weight[sp.second], molecular_weight_gold[sp.first]);
1261 for (const auto & sp : mgd.kin_species_index)
1262 ASSERT_EQ(mgd.kin_species_molecular_weight[sp.second], molecular_weight_gold[sp.first]);
1263
1264 // check molecular volumes swapped correctly
1265 for (const auto & sp : mgd.basis_species_index)
1266 ASSERT_EQ(mgd.basis_species_molecular_volume[sp.second], molecular_volume_gold[sp.first]);
1267 for (const auto & sp : mgd.eqm_species_index)
1268 ASSERT_EQ(mgd.eqm_species_molecular_volume[sp.second], molecular_volume_gold[sp.first]);
1269 for (const auto & sp : mgd.kin_species_index)
1270 ASSERT_EQ(mgd.kin_species_molecular_volume[sp.second], molecular_volume_gold[sp.first]);
1271
1272 // check isMineral information has swapped correctly
1273 for (const auto & species : mgd.basis_species_index)
1274 ASSERT_EQ(mgd.basis_species_mineral[species.second], false);
1275 for (const auto & species : mgd.eqm_species_index)
1276 ASSERT_EQ(mgd.eqm_species_mineral[species.second], false);
1277 for (const auto & species : mgd.kin_species_index)
1278 if (species.first == "Fe(OH)3(ppd)" || species.first == "Fe(OH)3(ppd)fake")
1279 ASSERT_EQ(mgd.kin_species_mineral[species.second], true);
1280 else
1281 ASSERT_EQ(mgd.kin_species_mineral[species.second], false);
1282
1283 // check isGas information has swapped correctly
1284 for (const auto & species : mgd.basis_species_index)
1285 ASSERT_EQ(mgd.basis_species_gas[species.second], false);
1286 for (const auto & species : mgd.eqm_species_index)
1287 if (species.first == "CH4(g)fake")
1288 ASSERT_EQ(mgd.eqm_species_gas[species.second], true);
1289 else
1290 ASSERT_EQ(mgd.eqm_species_gas[species.second], false);
1291
1292 // check isTransported information has swapped correctly
1293 for (const auto & species : mgd.basis_species_index)
1294 if (species.first == ">(s)FeOH" || species.first == ">(w)FeOH")
1295 ASSERT_EQ(mgd.basis_species_transported[species.second], false);
1296 else
1297 ASSERT_EQ(mgd.basis_species_transported[species.second], true);
1298 for (const auto & species : mgd.eqm_species_index)
1299 ASSERT_EQ(mgd.eqm_species_transported[species.second], true);
1300 for (const auto & species : mgd.kin_species_index)
1301 if (species.first == "Fe(OH)3(ppd)" || species.first == "Fe(OH)3(ppd)fake" ||
1302 species.first == ">(s)FeO-")
1303 ASSERT_EQ(mgd.kin_species_transported[species.second], false);
1304 else
1305 ASSERT_EQ(mgd.kin_species_transported[species.second], true);
1306
1307 // check stoichiometry
1308 for (const auto & sp : {"CO2(aq)",
1309 "CO3--",
1310 "OH-",
1311 "CH4(aq)",
1312 "O2(aq)",
1313 "CH4(g)fake",
1314 "Fe(OH)3(ppd)",
1315 "Fe(OH)3(ppd)fake",
1316 "(O-phth)--",
1317 ">(s)FeO-"})
1318 stoi_gold[sp] = DenseMatrix<Real>(1, 7);
1319 // remember the order of primaries:
1320 // {"H2O", "H+", ">(s)FeOH", ">(w)FeOH", "Fe++", "HCO3-", "Fe+++"}
1321 stoi_gold["CO2(aq)"](0, 0) = -1;
1322 stoi_gold["CO2(aq)"](0, 1) = 1;
1323 stoi_gold["CO2(aq)"](0, 5) = 1;
1324 stoi_gold["CO3--"](0, 5) = 1;
1325 stoi_gold["CO3--"](0, 1) = -1;
1326 stoi_gold["OH-"](0, 0) = 1;
1327 stoi_gold["OH-"](0, 1) = -1;
1328 stoi_gold["CH4(aq)"](0, 0) = -3;
1329 stoi_gold["CH4(aq)"](0, 1) = 9;
1330 stoi_gold["CH4(aq)"](0, 4) = 8;
1331 stoi_gold["CH4(aq)"](0, 5) = 1;
1332 stoi_gold["CH4(aq)"](0, 6) = -8;
1333 stoi_gold["O2(aq)"](0, 0) = 2;
1334 stoi_gold["O2(aq)"](0, 1) = -4;
1335 stoi_gold["O2(aq)"](0, 4) = -4;
1336 stoi_gold["O2(aq)"](0, 6) = 4;
1337 stoi_gold["CH4(g)fake"](0, 0) = -6;
1338 stoi_gold["CH4(g)fake"](0, 1) = 18;
1339 stoi_gold["CH4(g)fake"](0, 4) = 16;
1340 stoi_gold["CH4(g)fake"](0, 5) = 3.5;
1341 stoi_gold["CH4(g)fake"](0, 6) = -18;
1342 stoi_gold["Fe(OH)3(ppd)"](0, 0) = 3;
1343 stoi_gold["Fe(OH)3(ppd)"](0, 1) = -3;
1344 stoi_gold["Fe(OH)3(ppd)"](0, 6) = 1;
1345 stoi_gold["Fe(OH)3(ppd)fake"](0, 0) = 3;
1346 stoi_gold["Fe(OH)3(ppd)fake"](0, 1) = -3;
1347 stoi_gold["Fe(OH)3(ppd)fake"](0, 6) = 2;
1348 stoi_gold["(O-phth)--"](0, 0) = -5 - 7.5 * 2;
1349 stoi_gold["(O-phth)--"](0, 1) = 6 + 7.5 * 4;
1350 stoi_gold["(O-phth)--"](0, 4) = 7.5 * 4;
1351 stoi_gold["(O-phth)--"](0, 5) = 8;
1352 stoi_gold["(O-phth)--"](0, 6) = -7.5 * 4;
1353 stoi_gold[">(s)FeO-"](0, 2) = 1;
1354 stoi_gold[">(s)FeO-"](0, 1) = -1;
1355
1356 for (const auto & sp : {"CO2(aq)", "CO3--", "OH-", "CH4(aq)", "O2(aq)", "CH4(g)fake"})
1357 {
1358 const unsigned row = mgd.eqm_species_index[sp];
1359 for (unsigned i = 0; i < 7; ++i)
1360 ASSERT_NEAR(mgd.eqm_stoichiometry(row, i), stoi_gold[sp](0, i), eps);
1361 }
1362
1363 for (const auto & sp : {"Fe(OH)3(ppd)", "Fe(OH)3(ppd)fake", "(O-phth)--", ">(s)FeO-"})
1364 {
1365 const unsigned row = mgd.kin_species_index[sp];
1366 for (unsigned i = 0; i < 7; ++i)
1367 ASSERT_NEAR(mgd.kin_stoichiometry(row, i), stoi_gold[sp](0, i), eps);
1368 }
1369
1370 for (unsigned t = 0; t < 8; ++t)
1371 {
1372 ASSERT_NEAR(mgd.kin_log10K(mgd.kin_species_index.at("Fe(OH)3(ppd)"), t), feoh3ppd[t], eps);
1373 ASSERT_NEAR(mgd.kin_log10K(mgd.kin_species_index.at("Fe(OH)3(ppd)fake"), t), feoh3ppd[t], eps);
1374 ASSERT_NEAR(mgd.kin_log10K(mgd.kin_species_index.at("(O-phth)--"), t),
1375 ophth[t] - (7.5 / 0.25) * fe3[t],
1376 eps);
1377 ASSERT_NEAR(mgd.kin_log10K(mgd.kin_species_index.at(">(s)FeO-"), t), feom[t], eps);
1378 }
1379
1380 ASSERT_EQ(mgd.kin_rate.size(), (std::size_t)3);
1381 std::vector<std::vector<Real>> new_kin_rate_gold(3, std::vector<Real>(7 + 6, 0.0));
1382 std::vector<std::vector<Real>> new_kin_monod_gold(3, std::vector<Real>(7 + 6, 0.0));
1383 std::vector<std::vector<Real>> new_kin_k_gold(3, std::vector<Real>(7 + 6, 0.0));
1384 new_kin_rate_gold[0][mgd.basis_species_index.at("H+")] = 1.1;
1385 new_kin_rate_gold[0][7 + mgd.eqm_species_index.at("O2(aq)")] = 2.2;
1386 new_kin_rate_gold[0][7 + mgd.eqm_species_index.at("CO3--")] = 3.3;
1387 new_kin_monod_gold[0][mgd.basis_species_index.at("H+")] = -0.1;
1388 new_kin_monod_gold[0][7 + mgd.eqm_species_index.at("O2(aq)")] = 0.2;
1389 new_kin_monod_gold[0][7 + mgd.eqm_species_index.at("CO3--")] = -0.3;
1390 new_kin_k_gold[0][mgd.basis_species_index.at("H+")] = 0.5;
1391 new_kin_k_gold[0][7 + mgd.eqm_species_index.at("O2(aq)")] = 0.25;
1392 new_kin_k_gold[0][7 + mgd.eqm_species_index.at("CO3--")] = 0.125;
1393 new_kin_rate_gold[1][7 + mgd.eqm_species_index.at("O2(aq)")] = -1.1;
1394 new_kin_rate_gold[1][mgd.basis_species_index.at("Fe+++")] = -2.2;
1395 new_kin_rate_gold[1][mgd.basis_species_index.at("H2O")] = -3.3;
1396 new_kin_monod_gold[1][7 + mgd.eqm_species_index.at("O2(aq)")] = 0.75;
1397 new_kin_monod_gold[1][mgd.basis_species_index.at("Fe+++")] = 0.875;
1398 new_kin_monod_gold[1][mgd.basis_species_index.at("H2O")] = 1.0;
1399 new_kin_k_gold[1][7 + mgd.eqm_species_index.at("O2(aq)")] = 0.2;
1400 new_kin_k_gold[1][mgd.basis_species_index.at("Fe+++")] = 0.3;
1401 new_kin_k_gold[1][mgd.basis_species_index.at("H2O")] = 0.4;
1402 new_kin_rate_gold[2][7 + mgd.eqm_species_index.at("CO2(aq)")] = 1.25;
1403 new_kin_rate_gold[2][mgd.basis_species_index.at("Fe+++")] = 2.25;
1404 new_kin_rate_gold[2][mgd.basis_species_index.at("Fe++")] = 3.25;
1405 new_kin_monod_gold[2][7 + mgd.eqm_species_index.at("CO2(aq)")] = 0.1;
1406 new_kin_monod_gold[2][mgd.basis_species_index.at("Fe+++")] = -0.12;
1407 new_kin_monod_gold[2][mgd.basis_species_index.at("Fe++")] = 1.23;
1408 new_kin_k_gold[2][7 + mgd.eqm_species_index.at("CO2(aq)")] = 0.5;
1409 new_kin_k_gold[2][mgd.basis_species_index.at("Fe+++")] = 0.6;
1410 new_kin_k_gold[2][mgd.basis_species_index.at("Fe++")] = 0.7;
1411 for (unsigned i = 0; i < 3; ++i)
1412 for (unsigned j = 0; j < 7 + 6; ++j)
1413 {
1414 EXPECT_EQ(mgd.kin_rate[i].promoting_indices[j], new_kin_rate_gold[i][j]);
1415 EXPECT_EQ(mgd.kin_rate[i].promoting_monod_indices[j], new_kin_monod_gold[i][j]);
1416 EXPECT_EQ(mgd.kin_rate[i].promoting_half_saturation[j], new_kin_k_gold[i][j]);
1417 }
1418
1419 // check progeny_index is swapped correctly
1420 EXPECT_EQ(mgd.kin_rate[0].progeny_index, 7 + fe3_posn);
1421 EXPECT_EQ(mgd.kin_rate[1].progeny_index, mgd.basis_species_index.at("Fe++"));
1422 EXPECT_EQ(mgd.kin_rate[2].progeny_index, o2aq_posn);
1423}
1424
1426TEST(GeochemistrySpeciesSwapperTest, swap_redox)
1427{
1428 GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
1429
1430 PertinentGeochemicalSystem model(database,
1431 {"H2O", "H+", "Fe++", "O2(aq)", "HCO3-", "(O-phth)--", "Fe+++"},
1432 {},
1433 {},
1434 {},
1435 {},
1436 {},
1437 "O2(aq)",
1438 "e-");
1439 ModelGeochemicalDatabase mgd = model.modelGeochemicalDatabase();
1440 GeochemistrySpeciesSwapper swapper(mgd.basis_species_index.size(), 1E-6);
1441
1442 EXPECT_EQ(mgd.redox_lhs, "e-");
1443 EXPECT_EQ(mgd.redox_stoichiometry.m(), (unsigned)3);
1444 const Real p = 1.0 / 4.0 / 7.5;
1445 // not sure of the order of the redox_stoichiometry stuff: this tells us
1446 const bool fe3_is_slot_one =
1447 (std::abs(mgd.redox_stoichiometry(1, 2)) > 1.0E-8); // note: involves Fe++
1448 const unsigned fe3_slot = (fe3_is_slot_one ? 1 : 2);
1449 const unsigned ophth_slot = (fe3_is_slot_one ? 2 : 1);
1450 // e- = p*(O-phth)-- + (5p+0.5)*H20 - 8p*HCO3- + (-6p-1)*H+
1451 EXPECT_NEAR(mgd.redox_stoichiometry(ophth_slot, 0), 5.0 * p + 0.5, 1.0E-8);
1452 EXPECT_NEAR(mgd.redox_stoichiometry(ophth_slot, 1), -6.0 * p - 1.0, 1.0E-8);
1453 EXPECT_NEAR(mgd.redox_stoichiometry(ophth_slot, 2), 0.0, 1.0E-8);
1454 EXPECT_NEAR(mgd.redox_stoichiometry(ophth_slot, 3), 0.0, 1.0E-8);
1455 EXPECT_NEAR(mgd.redox_stoichiometry(ophth_slot, 4), -8.0 * p, 1.0E-8);
1456 EXPECT_NEAR(mgd.redox_stoichiometry(ophth_slot, 5), p, 1.0E-8);
1457 EXPECT_NEAR(mgd.redox_stoichiometry(ophth_slot, 6), 0.0, 1.0E-8);
1458 // e- = Fe++ - Fe+++
1459 EXPECT_NEAR(mgd.redox_stoichiometry(fe3_slot, 0), 0.0, 1.0E-8);
1460 EXPECT_NEAR(mgd.redox_stoichiometry(fe3_slot, 1), 0.0, 1.0E-8);
1461 EXPECT_NEAR(mgd.redox_stoichiometry(fe3_slot, 2), 1.0, 1.0E-8);
1462 EXPECT_NEAR(mgd.redox_stoichiometry(fe3_slot, 3), 0.0, 1.0E-8);
1463 EXPECT_NEAR(mgd.redox_stoichiometry(fe3_slot, 4), 0.0, 1.0E-8);
1464 EXPECT_NEAR(mgd.redox_stoichiometry(fe3_slot, 5), 0.0, 1.0E-8);
1465 EXPECT_NEAR(mgd.redox_stoichiometry(fe3_slot, 6), -1.0, 1.0E-8);
1466 // record stoichiometry and log10K to compare with below
1467 DenseMatrix<Real> orig_stoi = mgd.redox_stoichiometry;
1468 DenseMatrix<Real> orig_log10K = mgd.redox_log10K;
1469
1470 ASSERT_EQ(mgd.swap_to_original_basis.n(), (unsigned)0);
1471 ASSERT_EQ(mgd.have_swapped_out_of_basis.size(), (std::size_t)0);
1472 ASSERT_EQ(mgd.have_swapped_into_basis.size(), (std::size_t)0);
1473
1474 const unsigned hco3_posn = mgd.basis_species_index["HCO3-"];
1475 const unsigned co3_posn = mgd.eqm_species_index["CO3--"];
1476
1477 // swap CO3-- into the basis in place of HCO3-
1478 swapper.performSwap(mgd, "HCO3-", "CO3--");
1479
1480 ASSERT_EQ(mgd.have_swapped_out_of_basis.size(), (std::size_t)1);
1481 ASSERT_EQ(mgd.have_swapped_into_basis.size(), (std::size_t)1);
1482 ASSERT_EQ(mgd.have_swapped_out_of_basis[0], hco3_posn);
1483 ASSERT_EQ(mgd.have_swapped_into_basis[0], co3_posn);
1484 ASSERT_EQ(mgd.swap_to_original_basis.n(), (unsigned)7);
1485 DenseMatrix<Real> gold_swap_matrix(7, 7);
1486 for (unsigned i = 0; i < 7; ++i)
1487 {
1488 gold_swap_matrix(i, i) = 1.0;
1489 gold_swap_matrix(hco3_posn, i) = 0.0;
1490 }
1491 gold_swap_matrix(hco3_posn, 1) = -1.0; // H+
1492 gold_swap_matrix(hco3_posn, 4) = 1.0; // HCO3-
1493 for (unsigned i = 0; i < 7; ++i)
1494 for (unsigned j = 0; j < 7; ++j)
1495 EXPECT_EQ(mgd.swap_to_original_basis(i, j), gold_swap_matrix(i, j));
1496
1497 EXPECT_EQ(mgd.redox_lhs, "e-");
1498 EXPECT_EQ(mgd.redox_stoichiometry.m(), (unsigned)3);
1499
1500 // Since HCO3- = CO3-- + H+
1501 // e- = p*(O-phth)-- + (5p+0.5)*H20 - 8p*CO3-- + (-6p-1-8p)*H+
1502 EXPECT_NEAR(mgd.redox_stoichiometry(ophth_slot, 0), 5.0 * p + 0.5, 1.0E-8);
1503 EXPECT_NEAR(mgd.redox_stoichiometry(ophth_slot, 1), -6.0 * p - 1.0 - 8.0 * p, 1.0E-8);
1504 EXPECT_NEAR(mgd.redox_stoichiometry(ophth_slot, 2), 0.0, 1.0E-8);
1505 EXPECT_NEAR(mgd.redox_stoichiometry(ophth_slot, 3), 0.0, 1.0E-8);
1506 EXPECT_NEAR(mgd.redox_stoichiometry(ophth_slot, 4), -8.0 * p, 1.0E-8);
1507 EXPECT_NEAR(mgd.redox_stoichiometry(ophth_slot, 5), p, 1.0E-8);
1508 EXPECT_NEAR(mgd.redox_stoichiometry(ophth_slot, 6), 0.0, 1.0E-8);
1509 EXPECT_NEAR(
1510 mgd.redox_log10K(ophth_slot, 0), orig_log10K(ophth_slot, 0) + 8.0 * p * 10.6169, 1.0E-8);
1511 EXPECT_NEAR(
1512 mgd.redox_log10K(ophth_slot, 1), orig_log10K(ophth_slot, 1) + 8.0 * p * 10.3439, 1.0E-8);
1513 EXPECT_NEAR(
1514 mgd.redox_log10K(ophth_slot, 2), orig_log10K(ophth_slot, 2) + 8.0 * p * 10.2092, 1.0E-8);
1515 EXPECT_NEAR(
1516 mgd.redox_log10K(ophth_slot, 3), orig_log10K(ophth_slot, 3) + 8.0 * p * 10.2793, 1.0E-8);
1517 EXPECT_NEAR(
1518 mgd.redox_log10K(ophth_slot, 4), orig_log10K(ophth_slot, 4) + 8.0 * p * 10.5131, 1.0E-8);
1519 EXPECT_NEAR(
1520 mgd.redox_log10K(ophth_slot, 5), orig_log10K(ophth_slot, 5) + 8.0 * p * 10.8637, 1.0E-8);
1521 EXPECT_NEAR(
1522 mgd.redox_log10K(ophth_slot, 6), orig_log10K(ophth_slot, 6) + 8.0 * p * 11.2860, 1.0E-8);
1523 EXPECT_NEAR(
1524 mgd.redox_log10K(ophth_slot, 7), orig_log10K(ophth_slot, 7) + 8.0 * p * 11.6319, 1.0E-8);
1525 // e- = Fe++ - Fe+++ is unimpacted by the swap
1526 for (unsigned basis_i = 0; basis_i < 7; ++basis_i)
1527 EXPECT_NEAR(mgd.redox_stoichiometry(fe3_slot, basis_i), orig_stoi(fe3_slot, basis_i), 1.0E-8);
1528 for (unsigned temp = 0; temp < 8; ++temp)
1529 EXPECT_NEAR(mgd.redox_log10K(fe3_slot, temp), orig_log10K(fe3_slot, temp), 1.0E-8);
1530}
1531
1533TEST(GeochemistrySpeciesSwapperTest, findBestEqmSwapException)
1534{
1535 GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
1536
1537 // eqm species are: CO2(aq), CO3--, CaCO3, CaOH+, OH-, Calcite
1538 PertinentGeochemicalSystem model(
1539 database, {"H2O", "Ca++", "HCO3-", "H+"}, {"Calcite"}, {}, {}, {}, {}, "O2(aq)", "e-");
1540 ModelGeochemicalDatabase mgd = model.modelGeochemicalDatabase();
1541 GeochemistrySpeciesSwapper swapper(mgd.basis_species_index.size(), 1E-6);
1542 unsigned best;
1543
1544 try
1545 {
1546 const std::vector<Real> eqm_molality(6, 1.0);
1547 swapper.findBestEqmSwap(4, mgd, eqm_molality, false, false, false, best);
1548 FAIL() << "Missing expected exception.";
1549 }
1550 catch (const std::exception & e)
1551 {
1552 std::string msg(e.what());
1553 ASSERT_TRUE(msg.find("basis index 4 must be less than 4") != std::string::npos)
1554 << "Failed with unexpected error message: " << msg;
1555 }
1556
1557 try
1558 {
1559 const std::vector<Real> eqm_molality(5, 1.0);
1560 swapper.findBestEqmSwap(1, mgd, eqm_molality, false, false, false, best);
1561 FAIL() << "Missing expected exception.";
1562 }
1563 catch (const std::exception & e)
1564 {
1565 std::string msg(e.what());
1566 ASSERT_TRUE(msg.find("Size of eqm_molality is 5 which is not equal to the number of "
1567 "equilibrium species 6") != std::string::npos)
1568 << "Failed with unexpected error message: " << msg;
1569 }
1570}
1571
1573TEST(GeochemistrySpeciesSwapperTest, findBestEqmSwap)
1574{
1575 GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
1576
1577 // eqm species are: CO2(aq), CO3--, CaCO3, CaOH+, OH-, Calcite
1578 PertinentGeochemicalSystem model(
1579 database, {"H2O", "Ca++", "HCO3-", "H+"}, {"Calcite"}, {}, {}, {}, {}, "O2(aq)", "e-");
1580 ModelGeochemicalDatabase mgd = model.modelGeochemicalDatabase();
1581 GeochemistrySpeciesSwapper swapper(mgd.basis_species_index.size(), 1E-6);
1582 unsigned best;
1583 // the following equilibrium molality has molality=5.0 for equilibrium mineral Calcite, which
1584 // violates the assumption in the Geochemistry module that all equilibrium minerals have zero
1585 // molality, but i'm setting this for testing only
1586 const std::vector<Real> eqm_molality = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0};
1587 bool legit = swapper.findBestEqmSwap(1, mgd, eqm_molality, false, false, false, best);
1588 EXPECT_TRUE(legit);
1589 EXPECT_EQ(best, (unsigned)3);
1590 legit = swapper.findBestEqmSwap(1, mgd, eqm_molality, true, false, false, best);
1591 EXPECT_TRUE(legit);
1592 EXPECT_EQ(best, (unsigned)5);
1593 legit = swapper.findBestEqmSwap(1, mgd, eqm_molality, true, true, true, best);
1594 EXPECT_TRUE(legit);
1595 EXPECT_EQ(best, (unsigned)5);
1596
1597 // The following system has mineral, gas and sorption, with secondary species: CO2(aq), CO3--,
1598 // OH-, (O-phth)--, CH4(aq), Fe+++, CO2(aq), CO3--, OH-, >(s)FeO- (surface-sorption related),
1599 // Fe(OH)3(ppd) (mineral), CH4(g)fake (gas)
1600 PertinentGeochemicalSystem model_s(
1601 database,
1602 {"H2O", "H+", ">(s)FeOH", ">(w)FeOH", "Fe++", "HCO3-", "O2(aq)"},
1603 {"Fe(OH)3(ppd)"},
1604 {"CH4(g)fake"},
1605 {},
1606 {},
1607 {},
1608 "O2(aq)",
1609 "e-");
1610 ModelGeochemicalDatabase mgd_s = model_s.modelGeochemicalDatabase();
1611 GeochemistrySpeciesSwapper swapper_s(mgd_s.basis_species_index.size(), 1E-6);
1612
1613 const std::vector<Real> eqm_molality_s = {0.0, 1.0, 2.0, 3.0, 4.0, 50.0, 6.0E3, 7.0E3, 8.0E3};
1614 legit = swapper_s.findBestEqmSwap(1, mgd_s, eqm_molality_s, false, false, false, best);
1615 EXPECT_TRUE(legit);
1616 EXPECT_EQ(best, (unsigned)5);
1617 legit = swapper_s.findBestEqmSwap(1, mgd_s, eqm_molality_s, true, false, false, best);
1618 EXPECT_TRUE(legit);
1619 EXPECT_EQ(best, (unsigned)7);
1620 legit = swapper_s.findBestEqmSwap(5, mgd_s, eqm_molality_s, false, true, false, best);
1621 EXPECT_TRUE(legit);
1622 EXPECT_EQ(best, (unsigned)8);
1623 legit = swapper_s.findBestEqmSwap(1, mgd_s, eqm_molality_s, false, false, true, best);
1624 EXPECT_TRUE(legit);
1625 EXPECT_EQ(best, (unsigned)6);
1626}
p
const Real p
Definition GasLiquidMassTransferTest.C:18
model
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-")
database
const GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false)
mgd
const ModelGeochemicalDatabase mgd
Definition GeochemistryActivityCoefficientsDebyeHuckelTest.C:28
TEST
TEST(GeochemistrySpeciesSwapperTest, bulkCompositionException)
Check bulk_composition exception.
Definition GeochemistrySpeciesSwapperTest.C:18
GeochemicalDatabaseReader
Class for reading geochemical reactions from a MOOSE geochemical database.
Definition GeochemicalDatabaseReader.h:216
GeochemistrySpeciesSwapper
Class to swap basis species with equilibrium species.
Definition GeochemistrySpeciesSwapper.h:20
GeochemistrySpeciesSwapper::checkSwap
void checkSwap(const ModelGeochemicalDatabase &mgd, const std::string &replace_this, const std::string &with_this)
Check that replacing the named basis species with the named equilibrium species is valid.
Definition GeochemistrySpeciesSwapper.C:24
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
GeochemistrySpeciesSwapper::findBestEqmSwap
bool findBestEqmSwap(unsigned basis_ind, const ModelGeochemicalDatabase &mgd, const std::vector< Real > &eqm_molality, bool minerals_allowed, bool gas_allowed, bool sorption_allowed, unsigned &best_eqm_species) const
For the the given basis index, find the equilibrium index that it should be swapped with.
Definition GeochemistrySpeciesSwapper.C:368
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
PertinentGeochemicalSystem::addKineticRate
void addKineticRate(const KineticRateUserDescription &description)
Adds a rate description for kinetic_species.
Definition PertinentGeochemicalSystem.C:856
libMesh::DenseMatrix::n
unsigned int n() const
libMesh::DenseMatrix::m
unsigned int m() const
libMesh::DenseMatrix::sub_matrix
DenseMatrix sub_matrix(unsigned int row_id, unsigned int row_size, unsigned int col_id, unsigned int col_size) const
KineticRateUserDescription
Holds a user-specified description of a kinetic rate.
Definition GeochemistryKineticRateCalculator.h:76
ModelGeochemicalDatabase
Data structure to hold all relevant information from the database file.
Definition PertinentGeochemicalSystem.h:102
ModelGeochemicalDatabase::eqm_species_mineral
std::vector< bool > eqm_species_mineral
eqm_species_mineral[i] = true iff the i^th equilibrium species is a mineral
Definition PertinentGeochemicalSystem.h:200
ModelGeochemicalDatabase::eqm_species_radius
std::vector< Real > eqm_species_radius
all quantities have an ionic radius (Angstrom) for computing activity (mineral radius = 0,...
Definition PertinentGeochemicalSystem.h:212
ModelGeochemicalDatabase::eqm_stoichiometry
DenseMatrix< Real > eqm_stoichiometry
eqm_stoichiometry(i, j) = stoichiometric coefficient for equilibrium species "i" in terms of the basi...
Definition PertinentGeochemicalSystem.h:224
ModelGeochemicalDatabase::basis_species_molecular_weight
std::vector< Real > basis_species_molecular_weight
all quantities have a molecular weight (g)
Definition PertinentGeochemicalSystem.h:183
ModelGeochemicalDatabase::kin_stoichiometry
DenseMatrix< Real > kin_stoichiometry
kin_stoichiometry(i, j) = stoichiometric coefficient for kinetic species "i" in terms of the basis sp...
Definition PertinentGeochemicalSystem.h:334
ModelGeochemicalDatabase::basis_species_radius
std::vector< Real > basis_species_radius
all quantities have an ionic radius (Angstrom) for computing activity (mineral radius = 0,...
Definition PertinentGeochemicalSystem.h:180
ModelGeochemicalDatabase::basis_species_charge
std::vector< Real > basis_species_charge
all quantities have a charge (mineral charge = 0, gas charge = 0, oxide charge = 0)
Definition PertinentGeochemicalSystem.h:177
ModelGeochemicalDatabase::eqm_species_molecular_volume
std::vector< Real > eqm_species_molecular_volume
all quantities have a molecular volume (cm^3) (only nonzero for minerals, however)
Definition PertinentGeochemicalSystem.h:218
ModelGeochemicalDatabase::redox_lhs
std::string redox_lhs
the name of the species on the left-hand side of the redox equations.
Definition PertinentGeochemicalSystem.h:261
ModelGeochemicalDatabase::kin_species_mineral
std::vector< bool > kin_species_mineral
kin_species_mineral[j] = true iff the j^th kinetic species is a mineral
Definition PertinentGeochemicalSystem.h:310
ModelGeochemicalDatabase::eqm_species_molecular_weight
std::vector< Real > eqm_species_molecular_weight
all quantities have a molecular weight (g)
Definition PertinentGeochemicalSystem.h:215
ModelGeochemicalDatabase::have_swapped_into_basis
std::vector< unsigned > have_swapped_into_basis
Species that have been swapped into the basis.
Definition PertinentGeochemicalSystem.h:355
ModelGeochemicalDatabase::basis_species_transported
std::vector< bool > basis_species_transported
basis_species_transported[j] = true iff the j^th basis species is transported in reactive-transport s...
Definition PertinentGeochemicalSystem.h:174
ModelGeochemicalDatabase::kin_species_name
std::vector< std::string > kin_species_name
kin_species_name[j] = name of the j^th kinetic species
Definition PertinentGeochemicalSystem.h:307
ModelGeochemicalDatabase::redox_stoichiometry
DenseMatrix< Real > redox_stoichiometry
redox_stoichiometry(i, j) = stoichiometric coefficients for i^th redox species that is in disequilibr...
Definition PertinentGeochemicalSystem.h:276
ModelGeochemicalDatabase::swap_to_original_basis
DenseMatrix< Real > swap_to_original_basis
Swap matrix that allows expression in terms of the original basis.
Definition PertinentGeochemicalSystem.h:371
ModelGeochemicalDatabase::basis_species_mineral
std::vector< bool > basis_species_mineral
basis_species_mineral[j] = true iff the j^th basis species is a mineral
Definition PertinentGeochemicalSystem.h:168
ModelGeochemicalDatabase::kin_rate
std::vector< KineticRateDefinition > kin_rate
rates given to kinetic species.
Definition PertinentGeochemicalSystem.h:341
ModelGeochemicalDatabase::kin_species_molecular_volume
std::vector< Real > kin_species_molecular_volume
all quantities have a molecular volume (cm^3/mol) (only nonzero for minerals, however)
Definition PertinentGeochemicalSystem.h:322
ModelGeochemicalDatabase::redox_log10K
DenseMatrix< Real > redox_log10K
redox_log10K(i, j) = log10(equilibrium constant) for i^th redox species at the j^th temperature point
Definition PertinentGeochemicalSystem.h:282
ModelGeochemicalDatabase::kin_log10K
DenseMatrix< Real > kin_log10K
kin_log10K(i, j) = log10(equilibrium constant for the i^th kinetic species at the j^th temperature po...
Definition PertinentGeochemicalSystem.h:328
ModelGeochemicalDatabase::have_swapped_out_of_basis
std::vector< unsigned > have_swapped_out_of_basis
Species that have been swapped out of the basis.
Definition PertinentGeochemicalSystem.h:348
ModelGeochemicalDatabase::basis_species_molecular_volume
std::vector< Real > basis_species_molecular_volume
all quantities have a molecular volume (cm^3) (only nonzero for minerals, however)
Definition PertinentGeochemicalSystem.h:186
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_name
std::vector< std::string > basis_species_name
basis_species_name[j] = name of the j^th basis species
Definition PertinentGeochemicalSystem.h:165
ModelGeochemicalDatabase::kin_species_transported
std::vector< bool > kin_species_transported
kin_species_transported[j] = true iff the j^th kinetic species is transported in reactive-transport s...
Definition PertinentGeochemicalSystem.h:313
ModelGeochemicalDatabase::kin_species_charge
std::vector< Real > kin_species_charge
all kinetic quantities have a charge (mineral charge = 0)
Definition PertinentGeochemicalSystem.h:316
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::eqm_log10K
DenseMatrix< Real > eqm_log10K
eqm_log10K(i, j) = log10(equilibrium constant) for i^th equilibrium species at the j^th temperature p...
Definition PertinentGeochemicalSystem.h:230
ModelGeochemicalDatabase::basis_species_gas
std::vector< bool > basis_species_gas
basis_species_gas[j] = true iff the j^th basis species is a gas
Definition PertinentGeochemicalSystem.h:171
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
ModelGeochemicalDatabase::eqm_species_name
std::vector< std::string > eqm_species_name
eqm_species_name[i] = name of the i^th eqm species
Definition PertinentGeochemicalSystem.h:197
ModelGeochemicalDatabase::eqm_species_charge
std::vector< Real > eqm_species_charge
all quantities have a charge (mineral charge = 0, gas charge = 0, oxide charge = 0)
Definition PertinentGeochemicalSystem.h:209
ModelGeochemicalDatabase::eqm_species_gas
std::vector< bool > eqm_species_gas
eqm_species_gas[i] = true iff the i^th equilibrium species is a gas
Definition PertinentGeochemicalSystem.h:203
ModelGeochemicalDatabase::kin_species_molecular_weight
std::vector< Real > kin_species_molecular_weight
all quantities have a molecular weight (g/mol)
Definition PertinentGeochemicalSystem.h:319
ModelGeochemicalDatabase::eqm_species_transported
std::vector< bool > eqm_species_transported
eqm_species_transported[i] = true iff the i^th eqm species is transported in reactive-transport sims
Definition PertinentGeochemicalSystem.h:206
Generated on Fri Aug 21 2026 13:41:03 for https://mooseframework.inl.gov by doxygen 1.9.8