https://mooseframework.inl.gov
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 
13 
14 const Real eps =
15  1E-12; // accounts for precision loss when substituting reactions and inverting the swap matrix
16 
18 TEST(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-
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-");
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 
53 TEST(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-
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-");
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 
209 TEST(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
215  database, {"H2O", "Ca++", "HCO3-", "H+"}, {"Calcite"}, {}, {}, {}, {}, "O2(aq)", "e-");
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 
548 TEST(GeochemistrySpeciesSwapperTest, swap2)
549 {
550  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
551 
552  // eqm species are: OH-, CO2(aq), CO3--, StoiCheckRedox, StoiCheckGas
554  {"H2O", "StoiCheckBasis", "H+", "HCO3-"},
555  {},
556  {"StoiCheckGas"},
557  {},
558  {},
559  {},
560  "O2(aq)",
561  "e-");
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 
850 TEST(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+++
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,
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,
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,
922  "Fe+++",
923  2.0,
924  -1.0,
925  0.0);
926  model.addKineticRate(rate3);
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 
1426 TEST(GeochemistrySpeciesSwapperTest, swap_redox)
1427 {
1428  GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false);
1429 
1431  {"H2O", "H+", "Fe++", "O2(aq)", "HCO3-", "(O-phth)--", "Fe+++"},
1432  {},
1433  {},
1434  {},
1435  {},
1436  {},
1437  "O2(aq)",
1438  "e-");
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
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 
1533 TEST(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
1539  database, {"H2O", "Ca++", "HCO3-", "H+"}, {"Calcite"}, {}, {}, {}, {}, "O2(aq)", "e-");
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 
1573 TEST(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
1579  database, {"H2O", "Ca++", "HCO3-", "H+"}, {"Calcite"}, {}, {}, {}, {}, "O2(aq)", "e-");
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)
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 }
DenseMatrix< Real > redox_stoichiometry
redox_stoichiometry(i, j) = stoichiometric coefficients for i^th redox species that is in disequilibr...
std::vector< Real > eqm_species_molecular_weight
all quantities have a molecular weight (g)
std::string redox_lhs
the name of the species on the left-hand side of the redox equations.
std::vector< bool > kin_species_mineral
kin_species_mineral[j] = true iff the j^th kinetic species is a mineral
std::vector< bool > basis_species_transported
basis_species_transported[j] = true iff the j^th basis species is transported in reactive-transport s...
std::vector< Real > kin_species_charge
all kinetic quantities have a charge (mineral charge = 0)
const Real eps
std::vector< bool > eqm_species_gas
eqm_species_gas[i] = true iff the i^th equilibrium species is a gas
std::unordered_map< std::string, unsigned > basis_species_index
basis_species_index[name] = index of the basis species, within all ModelGeochemicalDatabase internal ...
const ModelGeochemicalDatabase mgd
DenseMatrix< Real > kin_stoichiometry
kin_stoichiometry(i, j) = stoichiometric coefficient for kinetic species "i" in terms of the basis sp...
std::vector< Real > kin_species_molecular_weight
all quantities have a molecular weight (g/mol)
std::vector< Real > kin_species_molecular_volume
all quantities have a molecular volume (cm^3/mol) (only nonzero for minerals, however) ...
unsigned int m() const
std::unordered_map< std::string, unsigned > eqm_species_index
eqm_species_index[name] = index of the equilibrium species (secondary aqueous species, redox couples in equilibrium with the aqueous solution, minerals in equilibrium with the aqueous solution, gases in equilibrium with the aqueous solution) within all ModelGeochemicalDatabase internal datastrcutres, with given name
DenseMatrix< Real > redox_log10K
redox_log10K(i, j) = log10(equilibrium constant) for i^th redox species at the j^th temperature point...
DenseMatrix< Real > swap_to_original_basis
Swap matrix that allows expression in terms of the original basis.
const GeochemicalDatabaseReader database("database/moose_testdb.json", true, true, false)
std::unordered_map< std::string, unsigned > kin_species_index
kin_species_index[name] = index of the kinetic species, within all ModelGeochemicalDatabase internal ...
DenseMatrix< Real > eqm_stoichiometry
eqm_stoichiometry(i, j) = stoichiometric coefficient for equilibrium species "i" in terms of the basi...
std::vector< Real > basis_species_radius
all quantities have an ionic radius (Angstrom) for computing activity (mineral radius = 0...
Class to swap basis species with equilibrium species.
bool findBestEqmSwap(unsigned basis_ind, const ModelGeochemicalDatabase &mgd, const std::vector< Real > &eqm_molality, bool minerals_allowed, bool gas_allowed, bool sorption_allowed, unsigned &best_eqm_species) const
For the the given basis index, find the equilibrium index that it should be swapped with...
const ModelGeochemicalDatabase & modelGeochemicalDatabase() const
Return a reference to the ModelGeochemicalDatabase structure.
void addKineticRate(const KineticRateUserDescription &description)
Adds a rate description for kinetic_species.
Constructs and stores a minimal amount of information that is pertinent to the user-defined geochemic...
DenseMatrix< Real > kin_log10K
kin_log10K(i, j) = log10(equilibrium constant for the i^th kinetic species at the j^th temperature po...
std::vector< bool > kin_species_transported
kin_species_transported[j] = true iff the j^th kinetic species is transported in reactive-transport s...
const PertinentGeochemicalSystem model(database, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++", ">(s)FeOH", "radius_neg1", "radius_neg1.5"}, {"Calcite"}, {}, {"Calcite_asdf"}, {"CH4(aq)"}, {">(s)FeOCa+"}, "O2(aq)", "e-")
std::vector< Real > eqm_species_radius
all quantities have an ionic radius (Angstrom) for computing activity (mineral radius = 0...
DenseMatrix< Real > eqm_log10K
eqm_log10K(i, j) = log10(equilibrium constant) for i^th equilibrium species at the j^th temperature p...
std::vector< Real > basis_species_molecular_weight
all quantities have a molecular weight (g)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
std::vector< bool > eqm_species_transported
eqm_species_transported[i] = true iff the i^th eqm species is transported in reactive-transport sims ...
std::vector< Real > eqm_species_molecular_volume
all quantities have a molecular volume (cm^3) (only nonzero for minerals, however) ...
std::vector< Real > basis_species_charge
all quantities have a charge (mineral charge = 0, gas charge = 0, oxide charge = 0) ...
Holds a user-specified description of a kinetic rate.
TEST(GeochemistrySpeciesSwapperTest, bulkCompositionException)
Check bulk_composition exception.
std::vector< unsigned > have_swapped_into_basis
Species that have been swapped into the basis.
std::vector< bool > basis_species_gas
basis_species_gas[j] = true iff the j^th basis species is a gas
std::vector< bool > basis_species_mineral
basis_species_mineral[j] = true iff the j^th basis species is a mineral
std::vector< Real > eqm_species_charge
all quantities have a charge (mineral charge = 0, gas charge = 0, oxide charge = 0) ...
std::vector< std::string > kin_species_name
kin_species_name[j] = name of the j^th kinetic species
Data structure to hold all relevant information from the database file.
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
std::vector< bool > eqm_species_mineral
eqm_species_mineral[i] = true iff the i^th equilibrium species is a mineral
std::vector< Real > basis_species_molecular_volume
all quantities have a molecular volume (cm^3) (only nonzero for minerals, however) ...
unsigned int n() const
DenseMatrix sub_matrix(unsigned int row_id, unsigned int row_size, unsigned int col_id, unsigned int col_size) const
std::vector< std::string > eqm_species_name
eqm_species_name[i] = name of the i^th eqm species
std::vector< std::string > basis_species_name
basis_species_name[j] = name of the j^th basis species
Class for reading geochemical reactions from a MOOSE geochemical database.
std::vector< unsigned > have_swapped_out_of_basis
Species that have been swapped out of the basis.
std::vector< KineticRateDefinition > kin_rate
rates given to kinetic species.
void checkSwap(const ModelGeochemicalDatabase &mgd, const std::string &replace_this, const std::string &with_this)
Check that replacing the named basis species with the named equilibrium species is valid...