https://mooseframework.inl.gov
Loading...
Searching...
No Matches
Functions | Variables
GeochemicalSolverTest.C File Reference

Go to the source code of this file.

Functions

 TEST (GeochemicalSolverTest, exception)
 Check exception.
 
 TEST (GeochemicalSolverTest, setgetMaxInitialResidual)
 Check set/getMaxInitialResidual.
 
 TEST (GeochemicalSolverTest, solve1)
 Solve super-simple case.
 
 TEST (GeochemicalSolverTest, solve2)
 Solve realistic case with minerals and gases, but no precipitation, no redox, no sorption.
 
 TEST (GeochemicalSolverTest, solve3)
 Solve realistic case with minerals that precipitate, no redox, no sorption.
 
 TEST (GeochemicalSolverTest, solve3_restore)
 Solve realistic case with minerals that precipitate, no redox, no sorption, and then "restore" to check it works correctly, and then check copy-assignment of GeochemicalSystem works properly.
 
 TEST (GeochemicalSolverTest, solve4)
 Solve realistic case with redox disequilibrium, no minerals, no sorption.
 
 TEST (GeochemicalSolverTest, solve4_restore)
 Solve realistic case with redox disequilibrium, no minerals, no sorption, then "restore" the solution and check.
 
 TEST (GeochemicalSolverTest, solve5)
 Solve realistic case with sorption and minerals, no redox.
 
 TEST (GeochemicalSolverTest, solve5_restore)
 Repeat the realistic case with sorption and minerals, no redox, this time with a "restore", and then use the copy assignment operator of GeochemicalSystem to create a new geochemical system.
 
 TEST (GeochemicalSolverTest, solve_addH)
 Test progressively adding H+ to a system.
 
 TEST (GeochemicalSolverTest, maxSwapsException)
 Test that the max swaps allowed works OK.
 
 TEST (GeochemicalSolverTest, setRampMaxIonicStrength)
 Check setRampMaxIonicStrength.
 
 TEST (GeochemicalSolverTest, solve_kinetic1)
 Solve case that involves kinetic species with zero rates (so kinetic species should have no impact except to modify the bulk composition)
 
 TEST (GeochemicalSolverTest, solve_kinetic2)
 Solve case that involves kinetic species with constant rates.
 
 TEST (GeochemicalSolverTest, solve_kinetic3)
 Solve case that involves kinetic species with promoting indices and implicit solve.
 

Variables

const GeochemicalDatabaseReader db_solver ("database/moose_testdb.json", true, true, false)
 
const GeochemicalDatabaseReader db_full ("../database/moose_geochemdb.json", true, true, false)
 
const GeochemicalDatabaseReader db_ferric ("../test/database/ferric_hydroxide_sorption.json", true, true)
 
const PertinentGeochemicalSystem model_simplest (db_solver, {"H2O", "H+"}, {}, {}, {}, {}, {}, "O2(aq)", "e-")
 
GeochemistrySpeciesSwapper swapper2 (2, 1E-6)
 
GeochemistrySpeciesSwapper swapper_kin (4, 1E-6)
 
const std::vector< GeochemicalSystem::ConstraintUserMeaningEnum > cm2
 
const std::vector< GeochemistryUnitConverter::GeochemistryUnit > cu2
 
const std::vector< GeochemicalSystem::ConstraintUserMeaningEnum > cm4
 
const std::vector< GeochemistryUnitConverter::GeochemistryUnit > cu4
 
GeochemistryIonicStrength is_solver (3.0, 3.0, false, false)
 
GeochemistryActivityCoefficientsDebyeHuckel ac_solver (is_solver, db_solver)
 

Function Documentation

◆ TEST() [1/16]

TEST ( GeochemicalSolverTest  ,
exception   
)

Check exception.

Definition at line 44 of file GeochemicalSolverTest.C.

45{
46 ModelGeochemicalDatabase mgd = model_simplest.modelGeochemicalDatabase();
47 GeochemicalSystem egs(mgd,
48 ac_solver,
49 is_solver,
50 swapper2,
51 {},
52 {},
53 "H+",
54 {"H2O", "H+"},
55 {1.75, 3.0},
56 cu2,
57 cm2,
58 25,
59 0,
60 1E-20,
61 {},
62 {},
63 {});
64 try
65 {
66 const GeochemicalSolver solver(mgd.basis_species_name.size(),
67 mgd.kin_species_name.size(),
68 is_solver,
69 1.0,
70 0.1,
71 1,
72 1E100,
73 0.1,
74 1,
75 {},
76 3.0,
77 10,
78 true);
79 FAIL() << "Missing expected exception.";
80 }
81 catch (const std::exception & e)
82 {
83 std::string msg(e.what());
84 ASSERT_TRUE(msg.find("GeochemicalSolver: ramp_max_ionic_strength must be less than max_iter") !=
85 std::string::npos)
86 << "Failed with unexpected error message: " << msg;
87 }
88
89 try
90 {
91 const GeochemicalSolver solver(mgd.basis_species_name.size(),
92 mgd.kin_species_name.size(),
93 is_solver,
94 1.0,
95 0.1,
96 100,
97 0.0,
98 0.1,
99 1,
100 {},
101 3.0,
102 10,
103 true);
104 FAIL() << "Missing expected exception.";
105 }
106 catch (const std::exception & e)
107 {
108 std::string msg(e.what());
109 ASSERT_TRUE(msg.find("GeochemicalSolver: max_initial_residual must be positive") !=
110 std::string::npos)
111 << "Failed with unexpected error message: " << msg;
112 }
113
114 try
115 {
116 const GeochemicalSolver solver(mgd.basis_species_name.size(),
117 mgd.kin_species_name.size(),
118 is_solver,
119 1.0,
120 0.1,
121 100,
122 1E100,
123 0.1,
124 1,
125 {},
126 -1.0,
127 10,
128 true);
129 FAIL() << "Missing expected exception.";
130 }
131 catch (const std::exception & e)
132 {
133 std::string msg(e.what());
134 ASSERT_TRUE(msg.find("GeochemicalSolver: max_ionic_strength must not be negative") !=
135 std::string::npos)
136 << "Failed with unexpected error message: " << msg;
137 }
138
139 try
140 {
141 const GeochemicalSolver solver(mgd.basis_species_name.size(),
142 mgd.kin_species_name.size(),
143 is_solver,
144 1.0,
145 -0.1,
146 100,
147 1E100,
148 0.1,
149 1,
150 {},
151 1.0,
152 10,
153 true);
154 FAIL() << "Missing expected exception.";
155 }
156 catch (const std::exception & e)
157 {
158 std::string msg(e.what());
159 ASSERT_TRUE(msg.find("GeochemicalSolver: rel_tol must not be negative") != std::string::npos)
160 << "Failed with unexpected error message: " << msg;
161 }
162
163 try
164 {
165 const GeochemicalSolver solver(mgd.basis_species_name.size(),
166 mgd.kin_species_name.size(),
167 is_solver,
168 -1.0,
169 0.1,
170 100,
171 1E100,
172 0.1,
173 1,
174 {},
175 1.0,
176 10,
177 true);
178 FAIL() << "Missing expected exception.";
179 }
180 catch (const std::exception & e)
181 {
182 std::string msg(e.what());
183 ASSERT_TRUE(msg.find("GeochemicalSolver: abs_tol must not be negative") != std::string::npos)
184 << "Failed with unexpected error message: " << msg;
185 }
186
187 try
188 {
189 const GeochemicalSolver solver(mgd.basis_species_name.size(),
190 mgd.kin_species_name.size(),
191 is_solver,
192 0.0,
193 0.0,
194 100,
195 1E100,
196 0.1,
197 1,
198 {},
199 1.0,
200 10,
201 true);
202 FAIL() << "Missing expected exception.";
203 }
204 catch (const std::exception & e)
205 {
206 std::string msg(e.what());
207 ASSERT_TRUE(msg.find("GeochemicalSolver: either rel_tol or abs_tol must be positive") !=
208 std::string::npos)
209 << "Failed with unexpected error message: " << msg;
210 }
211
212 try
213 {
214 GeochemicalSolver solver(mgd.basis_species_name.size(),
215 mgd.kin_species_name.size(),
216 is_solver,
217 1.0,
218 0.1,
219 100,
220 1E100,
221 0.1,
222 1,
223 {},
224 1.0,
225 10,
226 true);
227 solver.setMaxInitialResidual(0.0);
228 FAIL() << "Missing expected exception.";
229 }
230 catch (const std::exception & e)
231 {
232 std::string msg(e.what());
233 ASSERT_TRUE(msg.find("GeochemicalSolver: max_initial_residual must be positive") !=
234 std::string::npos)
235 << "Failed with unexpected error message: " << msg;
236 }
237}
model_simplest
const PertinentGeochemicalSystem model_simplest(db_solver, {"H2O", "H+"}, {}, {}, {}, {}, {}, "O2(aq)", "e-")
swapper2
GeochemistrySpeciesSwapper swapper2(2, 1E-6)
is_solver
GeochemistryIonicStrength is_solver(3.0, 3.0, false, false)
cm2
const std::vector< GeochemicalSystem::ConstraintUserMeaningEnum > cm2
Definition GeochemicalSolverTest.C:24
ac_solver
GeochemistryActivityCoefficientsDebyeHuckel ac_solver(is_solver, db_solver)
cu2
const std::vector< GeochemistryUnitConverter::GeochemistryUnit > cu2
Definition GeochemicalSolverTest.C:27
mgd
const ModelGeochemicalDatabase mgd
Definition GeochemistryActivityCoefficientsDebyeHuckelTest.C:28
GeochemicalSolver
This class contains methods to solve the algebraic system in GeochemicalSystem.
Definition GeochemicalSolver.h:19
GeochemicalSystem
This class holds information about bulk composition, molalities, activities, activity coefficients,...
Definition GeochemicalSystem.h:40
PertinentGeochemicalSystem::modelGeochemicalDatabase
const ModelGeochemicalDatabase & modelGeochemicalDatabase() const
Return a reference to the ModelGeochemicalDatabase structure.
Definition PertinentGeochemicalSystem.C:850
ModelGeochemicalDatabase
Data structure to hold all relevant information from the database file.
Definition PertinentGeochemicalSystem.h:102
ModelGeochemicalDatabase::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::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

◆ TEST() [2/16]

TEST ( GeochemicalSolverTest  ,
maxSwapsException   
)

Test that the max swaps allowed works OK.

Definition at line 2218 of file GeochemicalSolverTest.C.

2219{
2220 // build the model
2221 const PertinentGeochemicalSystem model(
2222 db_full,
2223 {"H2O", "H+", "Cl-", "Na+", "SO4--", "Mg++", "Ca++", "K+", "HCO3-", "SiO2(aq)", "O2(aq)"},
2224 {"Antigorite",
2225 "Tremolite",
2226 "Talc",
2227 "Chrysotile",
2228 "Sepiolite",
2229 "Anthophyllite",
2230 "Dolomite",
2231 "Dolomite-ord",
2232 "Huntite",
2233 "Dolomite-dis",
2234 "Magnesite",
2235 "Calcite",
2236 "Aragonite",
2237 "Quartz"},
2238 {},
2239 {},
2240 {},
2241 {},
2242 "O2(aq)",
2243 "e-");
2244 ModelGeochemicalDatabase mgd = model.modelGeochemicalDatabase();
2245
2246 // build the equilibrium system
2247 GeochemistrySpeciesSwapper swapper(11, 1E-6);
2248 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
2249 GeochemicalSystem::ConstraintUserMeaningEnum::KG_SOLVENT_WATER,
2250 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
2251 GeochemicalSystem::ConstraintUserMeaningEnum::FREE_CONCENTRATION,
2252 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
2253 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
2254 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
2255 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
2256 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
2257 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
2258 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
2259 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION};
2260 const std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu = {
2261 GeochemistryUnitConverter::GeochemistryUnit::KG,
2262 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
2263 GeochemistryUnitConverter::GeochemistryUnit::MOLAL,
2264 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
2265 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
2266 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
2267 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
2268 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
2269 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
2270 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
2271 GeochemistryUnitConverter::GeochemistryUnit::MOLES};
2272 GeochemicalSystem egs(
2273 mgd,
2274 ac_solver,
2275 is_solver,
2276 swapper,
2277 {"H+"},
2278 {"MgCO3"},
2279 "Cl-",
2280 {"H2O", "MgCO3", "O2(aq)", "Cl-", "Na+", "SO4--", "Mg++", "Ca++", "K+", "HCO3-", "SiO2(aq)"},
2281 {1.0,
2282 0.196E-3,
2283 0.2151E-3,
2284 0.5656,
2285 0.4850,
2286 0.02924,
2287 0.05501,
2288 0.01063,
2289 0.010576055,
2290 0.002412,
2291 0.00010349},
2292 cu,
2293 cm,
2294 25,
2295 0,
2296 1E-20,
2297 {},
2298 {},
2299 {});
2300
2301 // build solver (4 swaps are needed to solve this system)
2302 GeochemicalSolver solver(mgd.basis_species_name.size(),
2303 mgd.kin_species_name.size(),
2304 is_solver,
2305 1.0E-15,
2306 1.0E-200,
2307 100,
2308 1E3,
2309 0.1,
2310 3,
2311 {"Dolomite-ord", "Dolomite-dis"},
2312 3.0,
2313 10,
2314 false);
2315
2316 // solve
2317 std::stringstream ss;
2318 unsigned tot_iter;
2319 Real abs_residual;
2320 try
2321 {
2322 DenseVector<Real> mole_additions(egs.getNumInBasis() + egs.getNumKinetic());
2323 DenseMatrix<Real> dmole_additions(egs.getNumInBasis() + egs.getNumKinetic(),
2324 egs.getNumInBasis() + egs.getNumKinetic());
2325 solver.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
2326 FAIL() << "Missing expected exception.";
2327 }
2328 catch (const std::exception & e)
2329 {
2330 std::string msg(e.what());
2331 ASSERT_TRUE(msg.find("Maximum number of swaps performed during solve") != std::string::npos)
2332 << "Failed with unexpected error message: " << msg;
2333 }
2334}
db_full
const GeochemicalDatabaseReader db_full("../database/moose_geochemdb.json", true, true, false)
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-")
GeochemistrySpeciesSwapper
Class to swap basis species with equilibrium species.
Definition GeochemistrySpeciesSwapper.h:20
PertinentGeochemicalSystem
Constructs and stores a minimal amount of information that is pertinent to the user-defined geochemic...
Definition PertinentGeochemicalSystem.h:387
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

◆ TEST() [3/16]

TEST ( GeochemicalSolverTest  ,
setgetMaxInitialResidual   
)

Check set/getMaxInitialResidual.

Definition at line 240 of file GeochemicalSolverTest.C.

241{
242 ModelGeochemicalDatabase mgd = model_simplest.modelGeochemicalDatabase();
243 GeochemicalSystem egs(mgd,
244 ac_solver,
245 is_solver,
246 swapper2,
247 {},
248 {},
249 "H+",
250 {"H2O", "H+"},
251 {1.75, 3.0},
252 cu2,
253 cm2,
254 25,
255 0,
256 1E-20,
257 {},
258 {},
259 {});
260 GeochemicalSolver solver(mgd.basis_species_name.size(),
261 mgd.kin_species_name.size(),
262 is_solver,
263 1.0,
264 0.1,
265 100,
266 99.0,
267 0.1,
268 1,
269 {},
270 1.0,
271 10,
272 true);
273 EXPECT_EQ(solver.getMaxInitialResidual(), 99.0);
274 solver.setMaxInitialResidual(123.0);
275 EXPECT_EQ(solver.getMaxInitialResidual(), 123.0);
276}

◆ TEST() [4/16]

TEST ( GeochemicalSolverTest  ,
setRampMaxIonicStrength   
)

Check setRampMaxIonicStrength.

Definition at line 2337 of file GeochemicalSolverTest.C.

2338{
2339 ModelGeochemicalDatabase mgd = model_simplest.modelGeochemicalDatabase();
2340 GeochemicalSystem egs(mgd,
2341 ac_solver,
2342 is_solver,
2343 swapper2,
2344 {},
2345 {},
2346 "H+",
2347 {"H2O", "H+"},
2348 {1.75, 3.0},
2349 cu2,
2350 cm2,
2351 25,
2352 0,
2353 1E-20,
2354 {},
2355 {},
2356 {});
2357 GeochemicalSolver solver(mgd.basis_species_name.size(),
2358 mgd.kin_species_name.size(),
2359 is_solver,
2360 1.0,
2361 0.1,
2362 100,
2363 1E100,
2364 0.1,
2365 1,
2366 {},
2367 3.0,
2368 10,
2369 true);
2370
2371 ASSERT_EQ(solver.getRampMaxIonicStrength(), (unsigned)10);
2372 try
2373 {
2374 solver.setRampMaxIonicStrength(101);
2375 FAIL() << "Missing expected exception.";
2376 }
2377 catch (const std::exception & e)
2378 {
2379 std::string msg(e.what());
2380 ASSERT_TRUE(msg.find("GeochemicalSolver: ramp_max_ionic_strength must be less than max_iter") !=
2381 std::string::npos)
2382 << "Failed with unexpected error message: " << msg;
2383 }
2384 solver.setRampMaxIonicStrength(21);
2385 ASSERT_EQ(solver.getRampMaxIonicStrength(), (unsigned)21);
2386}

◆ TEST() [5/16]

TEST ( GeochemicalSolverTest  ,
solve1   
)

Solve super-simple case.

Definition at line 279 of file GeochemicalSolverTest.C.

280{
281 ModelGeochemicalDatabase mgd = model_simplest.modelGeochemicalDatabase();
282 GeochemicalSystem egs(mgd,
283 ac_solver,
284 is_solver,
285 swapper2,
286 {},
287 {},
288 "H+",
289 {"H2O", "H+"},
290 {1.75, 1},
291 cu2,
292 cm2,
293 25,
294 0,
295 1E-20,
296 {},
297 {},
298 {});
299 GeochemicalSolver solver(mgd.basis_species_name.size(),
300 mgd.kin_species_name.size(),
301 is_solver,
302 1.0E-15,
303 1.0E-200,
304 100,
305 10.0,
306 0.1,
307 1,
308 {},
309 3.0,
310 0,
311 false);
312
313 std::stringstream ss;
314 unsigned tot_iter;
315 Real abs_residual;
316 DenseVector<Real> mole_additions(egs.getNumInBasis() + egs.getNumKinetic());
317 DenseMatrix<Real> dmole_additions(egs.getNumInBasis() + egs.getNumKinetic(),
318 egs.getNumInBasis() + egs.getNumKinetic());
319 solver.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
320
321 // check Newton has converged
322 EXPECT_LE(tot_iter, (unsigned)100);
323 EXPECT_LE(abs_residual, 1.0E-15);
324
325 // check that equilibrium molality is set correctly
326 EXPECT_NEAR(egs.getEquilibriumMolality(0) /
327 (egs.getBasisActivity(0) / egs.getBasisActivity(1) /
328 std::pow(10.0, egs.getLog10K(0)) / egs.getEquilibriumActivityCoefficient(0)),
329 1.0,
330 1.0E-9);
331
332 // check that basis molality is correct
333 EXPECT_EQ(egs.getSolventWaterMass(), 1.75);
334 EXPECT_NEAR(egs.getSolventMassAndFreeMolalityAndMineralMoles()[1] / egs.getEquilibriumMolality(0),
335 1.0,
336 1.0E-9);
337
338 // check that basis bulk composition is correct
339 EXPECT_NEAR(egs.getBulkMolesOld()[0],
340 1.75 * (GeochemistryConstants::MOLES_PER_KG_WATER + egs.getEquilibriumMolality(0)),
341 1.0E-9);
342 EXPECT_NEAR(egs.getBulkMolesOld()[1], 0.0, 1.0E-9);
343}

◆ TEST() [6/16]

TEST ( GeochemicalSolverTest  ,
solve2   
)

Solve realistic case with minerals and gases, but no precipitation, no redox, no sorption.

Definition at line 346 of file GeochemicalSolverTest.C.

347{
348 // build the model
349 const PertinentGeochemicalSystem model(
350 db_full,
351 {"H2O", "H+", "Cl-", "Na+", "SO4--", "Mg++", "Ca++", "K+", "HCO3-", "SiO2(aq)", "O2(aq)"},
352 {"Antigorite",
353 "Tremolite",
354 "Talc",
355 "Chrysotile",
356 "Sepiolite",
357 "Anthophyllite",
358 "Dolomite",
359 "Dolomite-ord",
360 "Huntite",
361 "Dolomite-dis",
362 "Magnesite",
363 "Calcite",
364 "Aragonite",
365 "Quartz"},
366 {"O2(g)", "CO2(g)"},
367 {},
368 {},
369 {},
370 "O2(aq)",
371 "e-");
372 ModelGeochemicalDatabase mgd = model.modelGeochemicalDatabase();
373
374 // build the equilibrium system
375 GeochemistrySpeciesSwapper swapper(11, 1E-6);
376 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
377 GeochemicalSystem::ConstraintUserMeaningEnum::KG_SOLVENT_WATER,
378 GeochemicalSystem::ConstraintUserMeaningEnum::FUGACITY,
379 GeochemicalSystem::ConstraintUserMeaningEnum::FUGACITY,
380 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
381 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
382 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
383 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
384 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
385 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
386 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
387 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION};
388 const std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu = {
389 GeochemistryUnitConverter::GeochemistryUnit::KG,
390 GeochemistryUnitConverter::GeochemistryUnit::DIMENSIONLESS,
391 GeochemistryUnitConverter::GeochemistryUnit::DIMENSIONLESS,
392 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
393 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
394 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
395 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
396 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
397 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
398 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
399 GeochemistryUnitConverter::GeochemistryUnit::MOLES};
400 GeochemicalSystem egs(
401 mgd,
402 ac_solver,
403 is_solver,
404 swapper,
405 {"H+", "O2(aq)"},
406 {"CO2(g)", "O2(g)"},
407 "Cl-",
408 {"H2O", "CO2(g)", "O2(g)", "Cl-", "Na+", "SO4--", "Mg++", "Ca++", "K+", "HCO3-", "SiO2(aq)"},
409 {1.0,
410 0.0003162278,
411 0.2,
412 0.5656,
413 0.4850,
414 0.02924,
415 0.05501,
416 0.01063,
417 0.010576055,
418 0.002412,
419 0.00010349},
420 cu,
421 cm,
422 25,
423 0,
424 1E-20,
425 {},
426 {},
427 {});
428
429 // build solver
430 GeochemicalSolver solver(mgd.basis_species_name.size(),
431 mgd.kin_species_name.size(),
432 is_solver,
433 1.0E-15,
434 1.0E-200,
435 100,
436 1E3,
437 0.1,
438 1,
439 {"Antigorite",
440 "Tremolite",
441 "Talc",
442 "Chrysotile",
443 "Sepiolite",
444 "Anthophyllite",
445 "Dolomite",
446 "Dolomite-ord",
447 "Huntite",
448 "Dolomite-dis",
449 "Magnesite",
450 "Calcite",
451 "Aragonite",
452 "Quartz"},
453 3.0,
454 10,
455 false);
456
457 // solve
458 std::stringstream ss;
459 unsigned tot_iter;
460 Real abs_residual;
461 DenseVector<Real> mole_additions(egs.getNumInBasis() + egs.getNumKinetic());
462 DenseMatrix<Real> dmole_additions(egs.getNumInBasis() + egs.getNumKinetic(),
463 egs.getNumInBasis() + egs.getNumKinetic());
464 solver.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
465
466 // check converged
467 EXPECT_LE(tot_iter, (unsigned)100);
468 EXPECT_LE(abs_residual, 1.0E-15);
469
470 // check number in basis, number in redox disequilibrium and number of surface potentials
471 EXPECT_EQ(egs.getNumInBasis(), (unsigned)11);
472 EXPECT_EQ(egs.getNumRedox(), (unsigned)1);
473 EXPECT_EQ(egs.getNumSurfacePotentials(), (unsigned)0);
474 EXPECT_EQ(egs.getNumInEquilibrium(), mgd.eqm_species_name.size());
475 EXPECT_EQ(egs.getNumKinetic(), (unsigned)0);
476
477 // check that the constraints are satisfied
478 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
479 if (mgd.basis_species_name[i] == "H2O")
480 {
481 EXPECT_FALSE(mgd.basis_species_gas[i]);
482 EXPECT_FALSE(mgd.basis_species_mineral[i]);
483 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
484 EXPECT_NEAR(egs.getSolventWaterMass(), 1.0, 1.0E-15);
485 EXPECT_NEAR(egs.getSolventMassAndFreeMolalityAndMineralMoles()[0], 1.0, 1.0E-15);
486 }
487 else if (mgd.basis_species_name[i] == "CO2(g)")
488 {
489 EXPECT_TRUE(mgd.basis_species_gas[i]);
490 EXPECT_FALSE(mgd.basis_species_mineral[i]);
491 EXPECT_TRUE(egs.getBasisActivityKnown()[i]);
492 EXPECT_NEAR(egs.getBasisActivity(i), 0.0003162278, 1.0E-15);
493 }
494 else if (mgd.basis_species_name[i] == "O2(g)")
495 {
496 EXPECT_TRUE(mgd.basis_species_gas[i]);
497 EXPECT_FALSE(mgd.basis_species_mineral[i]);
498 EXPECT_TRUE(egs.getBasisActivityKnown()[i]);
499 EXPECT_NEAR(egs.getBasisActivity(i), 0.2, 1.0E-15);
500 }
501 else if (mgd.basis_species_name[i] == "Cl-")
502 {
503 EXPECT_FALSE(mgd.basis_species_gas[i]);
504 EXPECT_FALSE(mgd.basis_species_mineral[i]);
505 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
506 // do not know the bulk composition as it is dictated by charge neutrality
507 EXPECT_EQ(egs.getChargeBalanceBasisIndex(), i);
508 }
509 else if (mgd.basis_species_name[i] == "Na+")
510 {
511 EXPECT_FALSE(mgd.basis_species_gas[i]);
512 EXPECT_FALSE(mgd.basis_species_mineral[i]);
513 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
514 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.4850, 1.0E-15);
515 }
516 else if (mgd.basis_species_name[i] == "SO4--")
517 {
518 EXPECT_FALSE(mgd.basis_species_gas[i]);
519 EXPECT_FALSE(mgd.basis_species_mineral[i]);
520 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
521 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.02924, 1.0E-15);
522 }
523 else if (mgd.basis_species_name[i] == "Mg++")
524 {
525 EXPECT_FALSE(mgd.basis_species_gas[i]);
526 EXPECT_FALSE(mgd.basis_species_mineral[i]);
527 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
528 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.05501, 1.0E-15);
529 }
530 else if (mgd.basis_species_name[i] == "Ca++")
531 {
532 EXPECT_FALSE(mgd.basis_species_gas[i]);
533 EXPECT_FALSE(mgd.basis_species_mineral[i]);
534 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
535 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.01063, 1.0E-15);
536 }
537 else if (mgd.basis_species_name[i] == "K+")
538 {
539 EXPECT_FALSE(mgd.basis_species_gas[i]);
540 EXPECT_FALSE(mgd.basis_species_mineral[i]);
541 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
542 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.010576055, 1.0E-15);
543 }
544 else if (mgd.basis_species_name[i] == "HCO3-")
545 {
546 EXPECT_FALSE(mgd.basis_species_gas[i]);
547 EXPECT_FALSE(mgd.basis_species_mineral[i]);
548 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
549 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.002412, 1.0E-15);
550 }
551 else if (mgd.basis_species_name[i] == "SiO2(aq)")
552 {
553 EXPECT_FALSE(mgd.basis_species_gas[i]);
554 EXPECT_FALSE(mgd.basis_species_mineral[i]);
555 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
556 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.00010349, 1.0E-15);
557 }
558 else
559 FAIL() << "Incorrect basis species";
560
561 // check total charge
562 EXPECT_NEAR(egs.getTotalChargeOld(), 0.0, 1E-15);
563 // check total charge by summing up basis charges
564 Real tot_charge = 0.0;
565 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
566 tot_charge += egs.getBulkMolesOld()[i] * mgd.basis_species_charge[i];
567 EXPECT_NEAR(tot_charge, 0.0, 1E-15);
568
569 // surface charges
570 for (unsigned sp = 0; sp < egs.getNumSurfacePotentials(); ++sp)
571 {
572 Real charge = 0.0;
573 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
574 if (mgd.surface_sorption_related[j] && mgd.surface_sorption_number[j] == sp)
575 charge += mgd.eqm_species_charge[j] * egs.getEquilibriumMolality(j);
576 EXPECT_NEAR(charge,
577 egs.getSorbingSurfaceArea()[sp] * egs.getSurfaceCharge(sp) /
578 GeochemistryConstants::FARADAY,
579 1E-15);
580 }
581
582 // check basis activity = activity_coefficient * molality
583 for (unsigned i = 1; i < egs.getNumInBasis(); ++i) // don't loop over water
584 if (mgd.basis_species_gas[i] || mgd.basis_species_mineral[i] || egs.getBasisActivityKnown()[i])
585 continue;
586 else
587 EXPECT_NEAR(egs.getBasisActivity(i),
588 egs.getBasisActivityCoefficient(i) *
589 egs.getSolventMassAndFreeMolalityAndMineralMoles()[i],
590 1E-15);
591
592 // check residuals are zero
593 for (unsigned a = 0; a < egs.getNumInAlgebraicSystem(); ++a)
594 EXPECT_LE(std::abs(egs.getResidualComponent(a, mole_additions)), 1E-15);
595 // check residuals are zero by summing molalities, bulk compositions, etc
596 Real nw = egs.getSolventWaterMass();
597 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
598 {
599 Real res = -egs.getBulkMolesOld()[i];
600 if (i == 0)
601 res += nw * GeochemistryConstants::MOLES_PER_KG_WATER;
602 else if (mgd.basis_species_mineral[i])
603 res += egs.getSolventMassAndFreeMolalityAndMineralMoles()[i];
604 else if (mgd.basis_species_gas[i])
605 res += 0.0;
606 else
607 res += nw * egs.getSolventMassAndFreeMolalityAndMineralMoles()[i];
608 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
609 res += nw * mgd.eqm_stoichiometry(j, i) * egs.getEquilibriumMolality(j);
610 EXPECT_LE(std::abs(res), 1E-13);
611 }
612 // surface potentials
613 const Real prefactor = std::sqrt(GeochemistryConstants::GAS_CONSTANT *
614 (25.0 + GeochemistryConstants::CELSIUS_TO_KELVIN) *
615 GeochemistryConstants::PERMITTIVITY_FREE_SPACE *
616 GeochemistryConstants::DIELECTRIC_CONSTANT_WATER *
617 GeochemistryConstants::DENSITY_WATER * egs.getIonicStrength()) /
618 nw;
619 for (unsigned sp = 0; sp < egs.getNumSurfacePotentials(); ++sp)
620 {
621 EXPECT_NEAR(prefactor * std::sinh(egs.getSurfacePotential(sp) * GeochemistryConstants::FARADAY /
622 2.0 / GeochemistryConstants::GAS_CONSTANT /
623 (25.0 + GeochemistryConstants::CELSIUS_TO_KELVIN)),
624 egs.getSurfaceCharge(sp),
625 1.0E-13);
626 }
627
628 // check equilibrium mass balance
629 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
630 {
631 if (mgd.eqm_species_mineral[j] || mgd.eqm_species_gas[j])
632 continue;
633 Real log10ap = 0.0;
634 for (unsigned i = 0; i < mgd.basis_species_name.size(); ++i)
635 log10ap += mgd.eqm_stoichiometry(j, i) * std::log10(egs.getBasisActivity(i));
636 log10ap -= egs.getLog10K(j);
637 if (mgd.surface_sorption_related[j])
638 log10ap -= std::log10(egs.getSurfacePotential(mgd.surface_sorption_number[j]));
639 else
640 log10ap -= std::log10(egs.getEquilibriumActivityCoefficient(j));
641 if (log10ap < -300.0)
642 EXPECT_LE(std::abs(egs.getEquilibriumMolality(j)), 1E-25);
643 else
644 EXPECT_LE(std::abs(egs.getEquilibriumMolality(j) - std::pow(10.0, log10ap)), 1E-15);
645 }
646}
GeochemistryConstants::DIELECTRIC_CONSTANT_WATER
constexpr Real DIELECTRIC_CONSTANT_WATER
Definition GeochemistryConstants.h:24
GeochemistryConstants::PERMITTIVITY_FREE_SPACE
constexpr Real PERMITTIVITY_FREE_SPACE
Definition GeochemistryConstants.h:23
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_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_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::surface_sorption_number
std::vector< unsigned > surface_sorption_number
surface_sorption_number[j] = the index of the surface potential that should be used to modify the equ...
Definition PertinentGeochemicalSystem.h:255
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::surface_sorption_related
std::vector< bool > surface_sorption_related
surface_sorption_related[j] = true iff the j^th equilibrium species is involved in surface sorption
Definition PertinentGeochemicalSystem.h:246
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::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

◆ TEST() [7/16]

TEST ( GeochemicalSolverTest  ,
solve3   
)

Solve realistic case with minerals that precipitate, no redox, no sorption.

Definition at line 649 of file GeochemicalSolverTest.C.

650{
651 // build the model
652 const PertinentGeochemicalSystem model(
653 db_full,
654 {"H2O", "H+", "Cl-", "Na+", "SO4--", "Mg++", "Ca++", "K+", "HCO3-", "SiO2(aq)", "O2(aq)"},
655 {"Antigorite",
656 "Tremolite",
657 "Talc",
658 "Chrysotile",
659 "Sepiolite",
660 "Anthophyllite",
661 "Dolomite",
662 "Dolomite-ord",
663 "Huntite",
664 "Dolomite-dis",
665 "Magnesite",
666 "Calcite",
667 "Aragonite",
668 "Quartz"},
669 {},
670 {},
671 {},
672 {},
673 "O2(aq)",
674 "e-");
675 ModelGeochemicalDatabase mgd = model.modelGeochemicalDatabase();
676
677 // build the equilibrium system
678 GeochemistrySpeciesSwapper swapper(11, 1E-6);
679 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
680 GeochemicalSystem::ConstraintUserMeaningEnum::KG_SOLVENT_WATER,
681 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
682 GeochemicalSystem::ConstraintUserMeaningEnum::FREE_CONCENTRATION,
683 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
684 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
685 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
686 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
687 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
688 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
689 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
690 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION};
691 const std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu = {
692 GeochemistryUnitConverter::GeochemistryUnit::KG,
693 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
694 GeochemistryUnitConverter::GeochemistryUnit::MOLAL,
695 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
696 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
697 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
698 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
699 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
700 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
701 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
702 GeochemistryUnitConverter::GeochemistryUnit::MOLES};
703 GeochemicalSystem egs(
704 mgd,
705 ac_solver,
706 is_solver,
707 swapper,
708 {"H+"},
709 {"MgCO3"},
710 "Cl-",
711 {"H2O", "MgCO3", "O2(aq)", "Cl-", "Na+", "SO4--", "Mg++", "Ca++", "K+", "HCO3-", "SiO2(aq)"},
712 {1.0,
713 0.196E-3,
714 0.2151E-3,
715 0.5656,
716 0.4850,
717 0.02924,
718 0.05501,
719 0.01063,
720 0.010576055,
721 0.002412,
722 0.00010349},
723 cu,
724 cm,
725 25,
726 0,
727 1E-20,
728 {},
729 {},
730 {});
731
732 // build solver (4 swaps are needed to solve this system)
733 GeochemicalSolver solver(mgd.basis_species_name.size(),
734 mgd.kin_species_name.size(),
735 is_solver,
736 1.0E-15,
737 1.0E-200,
738 100,
739 1E3,
740 0.1,
741 4,
742 {"Dolomite-ord", "Dolomite-dis"},
743 3.0,
744 10,
745 false);
746
747 // solve
748 std::stringstream ss;
749 unsigned tot_iter;
750 Real abs_residual;
751 DenseVector<Real> mole_additions(egs.getNumInBasis() + egs.getNumKinetic());
752 DenseMatrix<Real> dmole_additions(egs.getNumInBasis() + egs.getNumKinetic(),
753 egs.getNumInBasis() + egs.getNumKinetic());
754 solver.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
755
756 // check converged
757 EXPECT_LE(tot_iter, (unsigned)100);
758 EXPECT_LE(abs_residual, 1.0E-15);
759
760 // check number in basis, number in redox disequilibrium and number of surface potentials
761 EXPECT_EQ(egs.getNumInBasis(), (unsigned)11);
762 EXPECT_EQ(egs.getNumRedox(), (unsigned)1);
763 EXPECT_EQ(egs.getNumSurfacePotentials(), (unsigned)0);
764 EXPECT_EQ(egs.getNumInEquilibrium(), mgd.eqm_species_name.size());
765
766 // check that the constraints are satisfied
767 Real bulk_dolo = 0.0;
768 Real bulk_ca = 0.0;
769 Real bulk_mg = 0.0;
770 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
771 if (mgd.basis_species_name[i] == "H2O")
772 {
773 EXPECT_FALSE(mgd.basis_species_gas[i]);
774 EXPECT_FALSE(mgd.basis_species_mineral[i]);
775 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
776 EXPECT_NEAR(egs.getSolventWaterMass(), 1.0, 1.0E-15);
777 EXPECT_NEAR(egs.getSolventMassAndFreeMolalityAndMineralMoles()[0], 1.0, 1.0E-15);
778 }
779 else if (mgd.basis_species_name[i] == "O2(aq)")
780 {
781 EXPECT_FALSE(mgd.basis_species_gas[i]);
782 EXPECT_FALSE(mgd.basis_species_mineral[i]);
783 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
784 EXPECT_NEAR(egs.getSolventMassAndFreeMolalityAndMineralMoles()[i], 0.2151E-3, 1.0E-15);
785 }
786 else if (mgd.basis_species_name[i] == "Cl-")
787 {
788 EXPECT_FALSE(mgd.basis_species_gas[i]);
789 EXPECT_FALSE(mgd.basis_species_mineral[i]);
790 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
791 // do not know the bulk composition as it is dictated by charge neutrality
792 EXPECT_EQ(egs.getChargeBalanceBasisIndex(), i);
793 }
794 else if (mgd.basis_species_name[i] == "Na+")
795 {
796 EXPECT_FALSE(mgd.basis_species_gas[i]);
797 EXPECT_FALSE(mgd.basis_species_mineral[i]);
798 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
799 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.4850, 1.0E-15);
800 }
801 else if (mgd.basis_species_name[i] == "SO4--")
802 {
803 EXPECT_FALSE(mgd.basis_species_gas[i]);
804 EXPECT_FALSE(mgd.basis_species_mineral[i]);
805 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
806 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.02924, 1.0E-15);
807 }
808 else if (mgd.basis_species_name[i] == "Mg++")
809 {
810 EXPECT_FALSE(mgd.basis_species_gas[i]);
811 EXPECT_FALSE(mgd.basis_species_mineral[i]);
812 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
813 bulk_mg = egs.getBulkMolesOld()[i];
814 }
815 else if (mgd.basis_species_name[i] == "Ca++")
816 {
817 EXPECT_FALSE(mgd.basis_species_gas[i]);
818 EXPECT_FALSE(mgd.basis_species_mineral[i]);
819 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
820 bulk_ca = egs.getBulkMolesOld()[i];
821 }
822 else if (mgd.basis_species_name[i] == "K+")
823 {
824 EXPECT_FALSE(mgd.basis_species_gas[i]);
825 EXPECT_FALSE(mgd.basis_species_mineral[i]);
826 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
827 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.010576055, 1.0E-15);
828 }
829 else if (mgd.basis_species_name[i] == "Quartz")
830 {
831 EXPECT_FALSE(mgd.basis_species_gas[i]);
832 EXPECT_TRUE(mgd.basis_species_mineral[i]);
833 EXPECT_TRUE(egs.getBasisActivityKnown()[i]);
834 EXPECT_GE(egs.getSolventMassAndFreeMolalityAndMineralMoles()[i], 0.0);
835 // Quartz = SiO2(aq)
836 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.00010349, 1.0);
837 }
838 else if (mgd.basis_species_name[i] == "Dolomite")
839 {
840 EXPECT_FALSE(mgd.basis_species_gas[i]);
841 EXPECT_TRUE(mgd.basis_species_mineral[i]);
842 EXPECT_TRUE(egs.getBasisActivityKnown()[i]);
843 EXPECT_GE(egs.getSolventMassAndFreeMolalityAndMineralMoles()[i], 0.0);
844 bulk_dolo = egs.getBulkMolesOld()[i];
845 }
846 else if (mgd.basis_species_name[i] == "HCO3-")
847 {
848 EXPECT_FALSE(mgd.basis_species_gas[i]);
849 EXPECT_FALSE(mgd.basis_species_mineral[i]);
850 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
851 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.002412, 1.0E-13);
852 }
853 else
854 FAIL() << "Incorrect basis species";
855 // dolomite = 2MgCO3 - Mg + Ca
856 EXPECT_NEAR(bulk_dolo + bulk_ca, 0.01063, 1E-13);
857 EXPECT_NEAR(-bulk_dolo + bulk_mg, 0.05501, 1E-13);
858 EXPECT_NEAR(2.0 * bulk_dolo, 0.196E-3, 1E-13);
859
860 // check total charge
861 EXPECT_NEAR(egs.getTotalChargeOld(), 0.0, 1E-15);
862 // check total charge by summing up basis charges
863 Real tot_charge = 0.0;
864 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
865 tot_charge += egs.getBulkMolesOld()[i] * mgd.basis_species_charge[i];
866 EXPECT_NEAR(tot_charge, 0.0, 1E-15);
867
868 // surface charges
869 for (unsigned sp = 0; sp < egs.getNumSurfacePotentials(); ++sp)
870 {
871 Real charge = 0.0;
872 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
873 if (mgd.surface_sorption_related[j] && mgd.surface_sorption_number[j] == sp)
874 charge += mgd.eqm_species_charge[j] * egs.getEquilibriumMolality(j);
875 EXPECT_NEAR(charge,
876 egs.getSorbingSurfaceArea()[sp] * egs.getSurfaceCharge(sp) /
877 GeochemistryConstants::FARADAY,
878 1E-15);
879 }
880
881 // check basis activity = activity_coefficient * molality
882 for (unsigned i = 1; i < egs.getNumInBasis(); ++i) // don't loop over water
883 if (mgd.basis_species_gas[i] || mgd.basis_species_mineral[i] || egs.getBasisActivityKnown()[i])
884 continue;
885 else
886 EXPECT_NEAR(egs.getBasisActivity(i),
887 egs.getBasisActivityCoefficient(i) *
888 egs.getSolventMassAndFreeMolalityAndMineralMoles()[i],
889 1E-15);
890
891 // check residuals are zero
892 for (unsigned a = 0; a < egs.getNumInAlgebraicSystem(); ++a)
893 EXPECT_LE(std::abs(egs.getResidualComponent(a, mole_additions)), 1E-15);
894 // check residuals are zero by summing molalities, bulk compositions, etc
895 Real nw = egs.getSolventWaterMass();
896 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
897 {
898 Real res = -egs.getBulkMolesOld()[i];
899 if (i == 0)
900 res += nw * GeochemistryConstants::MOLES_PER_KG_WATER;
901 else if (mgd.basis_species_mineral[i])
902 res += egs.getSolventMassAndFreeMolalityAndMineralMoles()[i];
903 else if (mgd.basis_species_gas[i])
904 res += 0.0;
905 else
906 res += nw * egs.getSolventMassAndFreeMolalityAndMineralMoles()[i];
907 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
908 res += nw * mgd.eqm_stoichiometry(j, i) * egs.getEquilibriumMolality(j);
909 EXPECT_LE(std::abs(res), 1E-13);
910 }
911 // surface potentials
912 const Real prefactor = std::sqrt(GeochemistryConstants::GAS_CONSTANT *
913 (25.0 + GeochemistryConstants::CELSIUS_TO_KELVIN) *
914 GeochemistryConstants::PERMITTIVITY_FREE_SPACE *
915 GeochemistryConstants::DIELECTRIC_CONSTANT_WATER *
916 GeochemistryConstants::DENSITY_WATER * egs.getIonicStrength()) /
917 nw;
918 for (unsigned sp = 0; sp < egs.getNumSurfacePotentials(); ++sp)
919 {
920 EXPECT_NEAR(prefactor * std::sinh(egs.getSurfacePotential(sp) * GeochemistryConstants::FARADAY /
921 2.0 / GeochemistryConstants::GAS_CONSTANT /
922 (25.0 + GeochemistryConstants::CELSIUS_TO_KELVIN)),
923 egs.getSurfaceCharge(sp),
924 1.0E-13);
925 }
926
927 // check equilibrium mass balance
928 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
929 {
930 if (mgd.eqm_species_mineral[j] || mgd.eqm_species_gas[j])
931 continue;
932 Real log10ap = 0.0;
933 for (unsigned i = 0; i < mgd.basis_species_name.size(); ++i)
934 log10ap += mgd.eqm_stoichiometry(j, i) * std::log10(egs.getBasisActivity(i));
935 log10ap -= egs.getLog10K(j);
936 if (mgd.surface_sorption_related[j])
937 log10ap -= std::log10(egs.getSurfacePotential(mgd.surface_sorption_number[j]));
938 else
939 log10ap -= std::log10(egs.getEquilibriumActivityCoefficient(j));
940 if (log10ap < -300.0)
941 EXPECT_LE(std::abs(egs.getEquilibriumMolality(j)), 1E-25);
942 else
943 EXPECT_LE(std::abs(egs.getEquilibriumMolality(j) - std::pow(10.0, log10ap)), 1E-15);
944 }
945
946 // check equilibrium species have negative saturation indices, except for prevent_precipitation
947 // minerals
948 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
949 {
950 if (!mgd.eqm_species_mineral[j])
951 continue;
952 if (mgd.eqm_species_name[j] == "Dolomite-ord" || mgd.eqm_species_name[j] == "Dolomite-dis")
953 continue;
954 Real log10ap = 0.0;
955 for (unsigned i = 0; i < mgd.basis_species_name.size(); ++i)
956 log10ap += mgd.eqm_stoichiometry(j, i) * std::log10(egs.getBasisActivity(i));
957 EXPECT_LE(log10ap, egs.getLog10K(j));
958 }
959}

◆ TEST() [8/16]

TEST ( GeochemicalSolverTest  ,
solve3_restore   
)

Solve realistic case with minerals that precipitate, no redox, no sorption, and then "restore" to check it works correctly, and then check copy-assignment of GeochemicalSystem works properly.

Definition at line 962 of file GeochemicalSolverTest.C.

963{
964 // build the model
965 const PertinentGeochemicalSystem model(
966 db_full,
967 {"H2O", "H+", "Cl-", "Na+", "SO4--", "Mg++", "Ca++", "K+", "HCO3-", "SiO2(aq)", "O2(aq)"},
968 {"Antigorite",
969 "Tremolite",
970 "Talc",
971 "Chrysotile",
972 "Sepiolite",
973 "Anthophyllite",
974 "Dolomite",
975 "Dolomite-ord",
976 "Huntite",
977 "Dolomite-dis",
978 "Magnesite",
979 "Calcite",
980 "Aragonite",
981 "Quartz"},
982 {},
983 {},
984 {},
985 {},
986 "O2(aq)",
987 "e-");
988 ModelGeochemicalDatabase mgd = model.modelGeochemicalDatabase();
989
990 // build the equilibrium system
991 GeochemistrySpeciesSwapper swapper(11, 1E-6);
992 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
993 GeochemicalSystem::ConstraintUserMeaningEnum::KG_SOLVENT_WATER,
994 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
995 GeochemicalSystem::ConstraintUserMeaningEnum::FREE_CONCENTRATION,
996 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
997 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
998 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
999 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1000 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1001 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1002 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1003 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION};
1004 const std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu = {
1005 GeochemistryUnitConverter::GeochemistryUnit::KG,
1006 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1007 GeochemistryUnitConverter::GeochemistryUnit::MOLAL,
1008 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1009 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1010 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1011 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1012 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1013 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1014 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1015 GeochemistryUnitConverter::GeochemistryUnit::MOLES};
1016 GeochemicalSystem egs(
1017 mgd,
1018 ac_solver,
1019 is_solver,
1020 swapper,
1021 {"H+"},
1022 {"MgCO3"},
1023 "Cl-",
1024 {"H2O", "MgCO3", "O2(aq)", "Cl-", "Na+", "SO4--", "Mg++", "Ca++", "K+", "HCO3-", "SiO2(aq)"},
1025 {1.0,
1026 0.196E-3,
1027 0.2151E-3,
1028 0.5656,
1029 0.4850,
1030 0.02924,
1031 0.05501,
1032 0.01063,
1033 0.010576055,
1034 0.002412,
1035 0.00010349},
1036 cu,
1037 cm,
1038 25,
1039 0,
1040 1E-20,
1041 {},
1042 {},
1043 {});
1044
1045 // build solver
1046 GeochemicalSolver solver(mgd.basis_species_name.size(),
1047 mgd.kin_species_name.size(),
1048 is_solver,
1049 1.0E-15,
1050 1.0E-200,
1051 100,
1052 1E3,
1053 0.1,
1054 5,
1055 {"Dolomite-ord", "Dolomite-dis"},
1056 3.0,
1057 10,
1058 false);
1059
1060 // solve
1061 std::stringstream ss;
1062 unsigned tot_iter;
1063 Real abs_residual;
1064 DenseVector<Real> mole_additions(egs.getNumInBasis() + egs.getNumKinetic());
1065 DenseMatrix<Real> dmole_additions(egs.getNumInBasis() + egs.getNumKinetic(),
1066 egs.getNumInBasis() + egs.getNumKinetic());
1067 solver.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
1068 const Real old_residual = abs_residual;
1069
1070 // check converged
1071 EXPECT_LE(tot_iter, (unsigned)100);
1072 EXPECT_LE(abs_residual, 1.0E-15);
1073
1074 // retrieve the molalities, and set them into egs: this should result in no change if the
1075 // "restore" is working correctly
1076 std::vector<std::string> names = mgd.basis_species_name;
1077 names.insert(names.end(), mgd.eqm_species_name.begin(), mgd.eqm_species_name.end());
1078 std::vector<Real> molal = egs.getSolventMassAndFreeMolalityAndMineralMoles();
1079 for (unsigned j = 0; j < egs.getNumInEquilibrium(); ++j)
1080 molal.push_back(egs.getEquilibriumMolality(j));
1081 // form constraints_from_molalities, which are all false
1082 const std::vector<bool> com_false(egs.getNumInBasis(), false);
1083 egs.setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(
1084 names, molal, com_false);
1085
1086 // build solver with no ramping of the maximum ionic strength
1087 GeochemicalSolver solver0(mgd.basis_species_name.size(),
1088 mgd.kin_species_name.size(),
1089 is_solver,
1090 1.0E-15,
1091 1.0E-200,
1092 100,
1093 1E3,
1094 0.1,
1095 0,
1096 {"Dolomite-ord", "Dolomite-dis"},
1097 3.0,
1098 0,
1099 false);
1100
1101 // solve this: no swaps should be necessary
1102 solver0.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
1103
1104 // check that the soler thinks this is truly a solution and the residual has not changed
1105 EXPECT_EQ(tot_iter, (unsigned)0);
1106 EXPECT_EQ(abs_residual, old_residual);
1107
1108 // Now use constraints_from_molalities = true
1109 const std::vector<bool> com_true(egs.getNumInBasis(), true);
1110 egs.setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(
1111 names, molal, com_true);
1112
1113 solver0.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
1114
1115 // check that the soler thinks this is truly a solution and the residual has not increased
1116 EXPECT_EQ(tot_iter, (unsigned)0);
1117 EXPECT_LE(abs_residual, old_residual);
1118
1119 // now check the copy-assignment of GeochemicalSystem
1120 const PertinentGeochemicalSystem model_dest(
1121 db_full,
1122 {"H2O", "H+", "Cl-", "Na+", "SO4--", "Mg++", "Ca++", "K+", "HCO3-", "SiO2(aq)", "O2(aq)"},
1123 {"Antigorite",
1124 "Tremolite",
1125 "Talc",
1126 "Chrysotile",
1127 "Sepiolite",
1128 "Anthophyllite",
1129 "Dolomite",
1130 "Dolomite-ord",
1131 "Huntite",
1132 "Dolomite-dis",
1133 "Magnesite",
1134 "Calcite",
1135 "Aragonite",
1136 "Quartz"},
1137 {},
1138 {},
1139 {},
1140 {},
1141 "O2(aq)",
1142 "e-");
1143 ModelGeochemicalDatabase mgd_dest = model_dest.modelGeochemicalDatabase();
1144 GeochemicalSystem egs_dest(
1145 mgd_dest,
1146 ac_solver,
1147 is_solver,
1148 swapper,
1149 {"H+"},
1150 {"MgCO3"},
1151 "Cl-",
1152 {"H2O", "MgCO3", "O2(aq)", "Cl-", "Na+", "SO4--", "Mg++", "Ca++", "K+", "HCO3-", "SiO2(aq)"},
1153 {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0},
1154 cu,
1155 cm,
1156 25,
1157 0,
1158 1E-20,
1159 {},
1160 {},
1161 {});
1162
1163 // change constraint meanings in egs to provide a more thorough check of the copy assignment
1164 egs.closeSystem();
1165 // copy assignment
1166 egs_dest = egs;
1167 EXPECT_EQ(egs_dest, egs);
1168}
GeochemicalSystem::closeSystem
void closeSystem()
Changes a KG_SOLVENT_WATER constraint to MOLES_BULK_WATER (if there is such a constraint) and all FRE...
Definition GeochemicalSystem.C:2100

◆ TEST() [9/16]

TEST ( GeochemicalSolverTest  ,
solve4   
)

Solve realistic case with redox disequilibrium, no minerals, no sorption.

Definition at line 1171 of file GeochemicalSolverTest.C.

1172{
1173 // build the model
1174 const PertinentGeochemicalSystem model(db_full,
1175 {"H2O",
1176 "H+",
1177 "Cl-",
1178 "O2(aq)",
1179 "HCO3-",
1180 "Ca++",
1181 "Mg++",
1182 "Na+",
1183 "K+",
1184 "Fe++",
1185 "Fe+++",
1186 "Mn++",
1187 "Zn++",
1188 "SO4--"},
1189 {},
1190 {},
1191 {},
1192 {},
1193 {},
1194 "O2(aq)",
1195 "e-");
1196 ModelGeochemicalDatabase mgd = model.modelGeochemicalDatabase();
1197
1198 // build the equilibrium system
1199 GeochemistrySpeciesSwapper swapper(14, 1E-6);
1200 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
1201 GeochemicalSystem::ConstraintUserMeaningEnum::KG_SOLVENT_WATER,
1202 GeochemicalSystem::ConstraintUserMeaningEnum::ACTIVITY,
1203 GeochemicalSystem::ConstraintUserMeaningEnum::FREE_CONCENTRATION,
1204 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1205 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1206 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1207 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1208 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1209 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1210 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1211 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1212 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1213 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1214 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION};
1215 const std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu = {
1216 GeochemistryUnitConverter::GeochemistryUnit::KG,
1217 GeochemistryUnitConverter::GeochemistryUnit::DIMENSIONLESS,
1218 GeochemistryUnitConverter::GeochemistryUnit::MOLAL,
1219 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1220 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1221 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1222 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1223 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1224 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1225 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1226 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1227 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1228 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1229 GeochemistryUnitConverter::GeochemistryUnit::MOLES};
1230 GeochemicalSystem egs(mgd,
1231 ac_solver,
1232 is_solver,
1233 swapper,
1234 {},
1235 {},
1236 "Cl-",
1237 {"H2O",
1238 "H+",
1239 "O2(aq)",
1240 "Cl-",
1241 "HCO3-",
1242 "Ca++",
1243 "Mg++",
1244 "Na+",
1245 "K+",
1246 "Fe++",
1247 "Fe+++",
1248 "Mn++",
1249 "Zn++",
1250 "SO4--"},
1251 {1.0,
1252 0.8913E-6,
1253 0.13438E-3,
1254 3.041E-5,
1255 0.0295E-3,
1256 0.005938E-3,
1257 0.01448E-3,
1258 0.0018704E-3,
1259 0.005115E-3,
1260 0.012534E-3,
1261 0.0005372E-3,
1262 0.005042E-3,
1263 0.001897E-3,
1264 0.01562E-3},
1265 cu,
1266 cm,
1267 25,
1268 0,
1269 1E-20,
1270 {},
1271 {},
1272 {});
1273
1274 // build solver
1275 GeochemicalSolver solver(mgd.basis_species_name.size(),
1276 mgd.kin_species_name.size(),
1277 is_solver,
1278 1.0E-15,
1279 1.0E-200,
1280 100,
1281 1E-2,
1282 0.1,
1283 1,
1284 {},
1285 3.0,
1286 10,
1287 false);
1288
1289 // solve
1290 std::stringstream ss;
1291 unsigned tot_iter;
1292 Real abs_residual;
1293 DenseVector<Real> mole_additions(egs.getNumInBasis() + egs.getNumKinetic());
1294 DenseMatrix<Real> dmole_additions(egs.getNumInBasis() + egs.getNumKinetic(),
1295 egs.getNumInBasis() + egs.getNumKinetic());
1296 solver.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
1297
1298 // check converged
1299 EXPECT_LE(tot_iter, (unsigned)100);
1300 EXPECT_LE(abs_residual, 1.0E-15);
1301
1302 // check number in basis, number in redox disequilibrium and number of surface potentials
1303 EXPECT_EQ(egs.getNumInBasis(), (unsigned)14);
1304 EXPECT_EQ(egs.getNumRedox(), (unsigned)2);
1305 EXPECT_EQ(egs.getNumSurfacePotentials(), (unsigned)0);
1306 EXPECT_EQ(egs.getNumInEquilibrium(), mgd.eqm_species_name.size());
1307
1308 // check that the constraints are satisfied
1309 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
1310 if (mgd.basis_species_name[i] == "H2O")
1311 {
1312 EXPECT_FALSE(mgd.basis_species_gas[i]);
1313 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1314 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1315 EXPECT_NEAR(egs.getSolventWaterMass(), 1.0, 1.0E-15);
1316 EXPECT_NEAR(egs.getSolventMassAndFreeMolalityAndMineralMoles()[0], 1.0, 1.0E-15);
1317 }
1318 else if (mgd.basis_species_name[i] == "H+")
1319 {
1320 EXPECT_FALSE(mgd.basis_species_gas[i]);
1321 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1322 EXPECT_TRUE(egs.getBasisActivityKnown()[i]);
1323 EXPECT_NEAR(egs.getBasisActivity(i), 0.8913E-6, 1.0E-15);
1324 }
1325 else if (mgd.basis_species_name[i] == "O2(aq)")
1326 {
1327 EXPECT_FALSE(mgd.basis_species_gas[i]);
1328 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1329 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1330 EXPECT_NEAR(egs.getSolventMassAndFreeMolalityAndMineralMoles()[i], 0.13438E-3, 1.0E-15);
1331 }
1332 else if (mgd.basis_species_name[i] == "Cl-")
1333 {
1334 EXPECT_FALSE(mgd.basis_species_gas[i]);
1335 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1336 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1337 // do not know the bulk composition as it is dictated by charge neutrality
1338 EXPECT_EQ(egs.getChargeBalanceBasisIndex(), i);
1339 }
1340 else if (mgd.basis_species_name[i] == "HCO3-")
1341 {
1342 EXPECT_FALSE(mgd.basis_species_gas[i]);
1343 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1344 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1345 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.0295E-3, 1.0E-15);
1346 }
1347 else if (mgd.basis_species_name[i] == "Ca++")
1348 {
1349 EXPECT_FALSE(mgd.basis_species_gas[i]);
1350 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1351 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1352 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.005938E-3, 1.0E-15);
1353 }
1354 else if (mgd.basis_species_name[i] == "Mg++")
1355 {
1356 EXPECT_FALSE(mgd.basis_species_gas[i]);
1357 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1358 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1359 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.01448E-3, 1.0E-15);
1360 }
1361 else if (mgd.basis_species_name[i] == "Na+")
1362 {
1363 EXPECT_FALSE(mgd.basis_species_gas[i]);
1364 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1365 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1366 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.0018704E-3, 1.0E-15);
1367 }
1368 else if (mgd.basis_species_name[i] == "K+")
1369 {
1370 EXPECT_FALSE(mgd.basis_species_gas[i]);
1371 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1372 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1373 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.005115E-3, 1.0E-15);
1374 }
1375 else if (mgd.basis_species_name[i] == "Fe++")
1376 {
1377 EXPECT_FALSE(mgd.basis_species_gas[i]);
1378 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1379 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1380 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.012534E-3, 1.0E-15);
1381 }
1382 else if (mgd.basis_species_name[i] == "Fe+++")
1383 {
1384 EXPECT_FALSE(mgd.basis_species_gas[i]);
1385 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1386 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1387 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.0005372E-3, 1.0E-15);
1388 }
1389 else if (mgd.basis_species_name[i] == "Mn++")
1390 {
1391 EXPECT_FALSE(mgd.basis_species_gas[i]);
1392 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1393 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1394 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.005042E-3, 1.0E-15);
1395 }
1396 else if (mgd.basis_species_name[i] == "Zn++")
1397 {
1398 EXPECT_FALSE(mgd.basis_species_gas[i]);
1399 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1400 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1401 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.001897E-3, 1.0E-15);
1402 }
1403 else if (mgd.basis_species_name[i] == "SO4--")
1404 {
1405 EXPECT_FALSE(mgd.basis_species_gas[i]);
1406 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1407 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1408 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.01562E-3, 1.0E-15);
1409 }
1410 else
1411 FAIL() << "Incorrect basis species";
1412
1413 // check total charge
1414 EXPECT_NEAR(egs.getTotalChargeOld(), 0.0, 1E-15);
1415 // check total charge by summing up basis charges
1416 Real tot_charge = 0.0;
1417 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
1418 tot_charge += egs.getBulkMolesOld()[i] * mgd.basis_species_charge[i];
1419 EXPECT_NEAR(tot_charge, 0.0, 1E-15);
1420
1421 // surface charges
1422 for (unsigned sp = 0; sp < egs.getNumSurfacePotentials(); ++sp)
1423 {
1424 Real charge = 0.0;
1425 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
1426 if (mgd.surface_sorption_related[j] && mgd.surface_sorption_number[j] == sp)
1427 charge += mgd.eqm_species_charge[j] * egs.getEquilibriumMolality(j);
1428 EXPECT_NEAR(charge,
1429 egs.getSorbingSurfaceArea()[sp] * egs.getSurfaceCharge(sp) /
1430 GeochemistryConstants::FARADAY,
1431 1E-15);
1432 }
1433
1434 // check basis activity = activity_coefficient * molality
1435 for (unsigned i = 1; i < egs.getNumInBasis(); ++i) // don't loop over water
1436 if (mgd.basis_species_gas[i] || mgd.basis_species_mineral[i] || egs.getBasisActivityKnown()[i])
1437 continue;
1438 else
1439 EXPECT_NEAR(egs.getBasisActivity(i),
1440 egs.getBasisActivityCoefficient(i) *
1441 egs.getSolventMassAndFreeMolalityAndMineralMoles()[i],
1442 1E-15);
1443
1444 // check residuals are zero
1445 for (unsigned a = 0; a < egs.getNumInAlgebraicSystem(); ++a)
1446 EXPECT_LE(std::abs(egs.getResidualComponent(a, mole_additions)), 1E-15);
1447 // check residuals are zero by summing molalities, bulk compositions, etc
1448 Real nw = egs.getSolventWaterMass();
1449 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
1450 {
1451 Real res = -egs.getBulkMolesOld()[i];
1452 if (i == 0)
1453 res += nw * GeochemistryConstants::MOLES_PER_KG_WATER;
1454 else if (mgd.basis_species_mineral[i])
1455 res += egs.getSolventMassAndFreeMolalityAndMineralMoles()[i];
1456 else if (mgd.basis_species_gas[i])
1457 res += 0.0;
1458 else
1459 res += nw * egs.getSolventMassAndFreeMolalityAndMineralMoles()[i];
1460 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
1461 res += nw * mgd.eqm_stoichiometry(j, i) * egs.getEquilibriumMolality(j);
1462 EXPECT_LE(std::abs(res), 1E-14);
1463 }
1464 // surface potentials
1465 const Real prefactor = std::sqrt(GeochemistryConstants::GAS_CONSTANT *
1466 (25.0 + GeochemistryConstants::CELSIUS_TO_KELVIN) *
1467 GeochemistryConstants::PERMITTIVITY_FREE_SPACE *
1468 GeochemistryConstants::DIELECTRIC_CONSTANT_WATER *
1469 GeochemistryConstants::DENSITY_WATER * egs.getIonicStrength()) /
1470 nw;
1471 for (unsigned sp = 0; sp < egs.getNumSurfacePotentials(); ++sp)
1472 {
1473 EXPECT_NEAR(prefactor * std::sinh(egs.getSurfacePotential(sp) * GeochemistryConstants::FARADAY /
1474 2.0 / GeochemistryConstants::GAS_CONSTANT /
1475 (25.0 + GeochemistryConstants::CELSIUS_TO_KELVIN)),
1476 egs.getSurfaceCharge(sp),
1477 1.0E-13);
1478 }
1479
1480 // check equilibrium mass balance
1481 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
1482 {
1483 if (mgd.eqm_species_mineral[j] || mgd.eqm_species_gas[j])
1484 continue;
1485 Real log10ap = 0.0;
1486 for (unsigned i = 0; i < mgd.basis_species_name.size(); ++i)
1487 log10ap += mgd.eqm_stoichiometry(j, i) * std::log10(egs.getBasisActivity(i));
1488 log10ap -= egs.getLog10K(j);
1489 if (mgd.surface_sorption_related[j])
1490 log10ap -= std::log10(egs.getSurfacePotential(mgd.surface_sorption_number[j]));
1491 else
1492 log10ap -= std::log10(egs.getEquilibriumActivityCoefficient(j));
1493 if (log10ap < -300.0)
1494 EXPECT_LE(std::abs(egs.getEquilibriumMolality(j)), 1E-25);
1495 else
1496 EXPECT_LE(std::abs(egs.getEquilibriumMolality(j) - std::pow(10.0, log10ap)), 1E-15);
1497 }
1498
1499 // check equilibrium species have negative saturation indices, except for prevent_precipitation
1500 // minerals
1501 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
1502 {
1503 if (!mgd.eqm_species_mineral[j])
1504 continue;
1505 Real log10ap = 0.0;
1506 for (unsigned i = 0; i < mgd.basis_species_name.size(); ++i)
1507 log10ap += mgd.eqm_stoichiometry(j, i) * std::log10(egs.getBasisActivity(i));
1508 EXPECT_LE(log10ap, egs.getLog10K(j));
1509 }
1510}

◆ TEST() [10/16]

TEST ( GeochemicalSolverTest  ,
solve4_restore   
)

Solve realistic case with redox disequilibrium, no minerals, no sorption, then "restore" the solution and check.

Definition at line 1513 of file GeochemicalSolverTest.C.

1514{
1515 // build the model
1516 const PertinentGeochemicalSystem model(db_full,
1517 {"H2O",
1518 "H+",
1519 "Cl-",
1520 "O2(aq)",
1521 "HCO3-",
1522 "Ca++",
1523 "Mg++",
1524 "Na+",
1525 "K+",
1526 "Fe++",
1527 "Fe+++",
1528 "Mn++",
1529 "Zn++",
1530 "SO4--"},
1531 {},
1532 {},
1533 {},
1534 {},
1535 {},
1536 "O2(aq)",
1537 "e-");
1538 ModelGeochemicalDatabase mgd = model.modelGeochemicalDatabase();
1539
1540 // build the equilibrium system
1541 GeochemistrySpeciesSwapper swapper(14, 1E-6);
1542 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
1543 GeochemicalSystem::ConstraintUserMeaningEnum::KG_SOLVENT_WATER,
1544 GeochemicalSystem::ConstraintUserMeaningEnum::ACTIVITY,
1545 GeochemicalSystem::ConstraintUserMeaningEnum::FREE_CONCENTRATION,
1546 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1547 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1548 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1549 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1550 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1551 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1552 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1553 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1554 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1555 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1556 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION};
1557 const std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu = {
1558 GeochemistryUnitConverter::GeochemistryUnit::KG,
1559 GeochemistryUnitConverter::GeochemistryUnit::DIMENSIONLESS,
1560 GeochemistryUnitConverter::GeochemistryUnit::MOLAL,
1561 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1562 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1563 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1564 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1565 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1566 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1567 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1568 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1569 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1570 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1571 GeochemistryUnitConverter::GeochemistryUnit::MOLES};
1572 GeochemicalSystem egs(mgd,
1573 ac_solver,
1574 is_solver,
1575 swapper,
1576 {},
1577 {},
1578 "Cl-",
1579 {"H2O",
1580 "H+",
1581 "O2(aq)",
1582 "Cl-",
1583 "HCO3-",
1584 "Ca++",
1585 "Mg++",
1586 "Na+",
1587 "K+",
1588 "Fe++",
1589 "Fe+++",
1590 "Mn++",
1591 "Zn++",
1592 "SO4--"},
1593 {1.0,
1594 0.8913E-6,
1595 0.13438E-3,
1596 3.041E-5,
1597 0.0295E-3,
1598 0.005938E-3,
1599 0.01448E-3,
1600 0.0018704E-3,
1601 0.005115E-3,
1602 0.012534E-3,
1603 0.0005372E-3,
1604 0.005042E-3,
1605 0.001897E-3,
1606 0.01562E-3},
1607 cu,
1608 cm,
1609 25,
1610 0,
1611 1E-20,
1612 {},
1613 {},
1614 {});
1615
1616 // build solver
1617 GeochemicalSolver solver(mgd.basis_species_name.size(),
1618 mgd.kin_species_name.size(),
1619 is_solver,
1620 1.0E-15,
1621 1.0E-200,
1622 100,
1623 1E-2,
1624 0.1,
1625 1,
1626 {},
1627 3.0,
1628 10,
1629 false);
1630
1631 // solve
1632 std::stringstream ss;
1633 unsigned tot_iter;
1634 Real abs_residual;
1635 DenseVector<Real> mole_additions(egs.getNumInBasis() + egs.getNumKinetic());
1636 DenseMatrix<Real> dmole_additions(egs.getNumInBasis() + egs.getNumKinetic(),
1637 egs.getNumInBasis() + egs.getNumKinetic());
1638 solver.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
1639 const Real old_residual = abs_residual;
1640
1641 // check converged
1642 EXPECT_LE(tot_iter, (unsigned)100);
1643 EXPECT_LE(abs_residual, 1.0E-15);
1644
1645 // retrieve the molalities, and set them into egs: this should result in no change if the
1646 // "restore" is working correctly
1647 std::vector<std::string> names = mgd.basis_species_name;
1648 names.insert(names.end(), mgd.eqm_species_name.begin(), mgd.eqm_species_name.end());
1649 std::vector<Real> molal = egs.getSolventMassAndFreeMolalityAndMineralMoles();
1650 for (unsigned j = 0; j < egs.getNumInEquilibrium(); ++j)
1651 molal.push_back(egs.getEquilibriumMolality(j));
1652 // form constraints_from_molalities, which are all false
1653 const std::vector<bool> com_false(egs.getNumInBasis(), false);
1654 egs.setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(
1655 names, molal, com_false);
1656
1657 // build solver with no ramping of the maximum ionic strength
1658 GeochemicalSolver solver0(mgd.basis_species_name.size(),
1659 mgd.kin_species_name.size(),
1660 is_solver,
1661 1.0E-15,
1662 1.0E-200,
1663 100,
1664 1E-2,
1665 0.1,
1666 1,
1667 {},
1668 3.0,
1669 0,
1670 false);
1671
1672 solver0.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
1673
1674 // check that the soler thinks this is truly a solution and the residual has not changed (up to
1675 // precision-loss)
1676 EXPECT_EQ(tot_iter, (unsigned)0);
1677 EXPECT_LE(abs_residual, old_residual);
1678
1679 // Now use constraints_from_molalities = true
1680 const std::vector<bool> com_true(egs.getNumInBasis(), true);
1681 egs.setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(
1682 names, molal, com_true);
1683
1684 solver0.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
1685
1686 // check that the soler thinks this is truly a solution and the residual has not increased
1687 EXPECT_EQ(tot_iter, (unsigned)0);
1688 EXPECT_LE(abs_residual, old_residual);
1689}

◆ TEST() [11/16]

TEST ( GeochemicalSolverTest  ,
solve5   
)

Solve realistic case with sorption and minerals, no redox.

Definition at line 1692 of file GeochemicalSolverTest.C.

1693{
1694 // build the model
1695 const PertinentGeochemicalSystem model(
1696 db_ferric,
1697 {"H2O", "H+", "Na+", "Cl-", "Hg++", "Pb++", "SO4--", "Fe+++", ">(s)FeOH", ">(w)FeOH"},
1698 {"Fe(OH)3(ppd)"},
1699 {},
1700 {},
1701 {},
1702 {},
1703 "O2(aq)",
1704 "e-");
1705 ModelGeochemicalDatabase mgd = model.modelGeochemicalDatabase();
1706
1707 // build the equilibrium system
1708 GeochemistrySpeciesSwapper swapper(10, 1E-6);
1709 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
1710 GeochemicalSystem::ConstraintUserMeaningEnum::KG_SOLVENT_WATER,
1711 GeochemicalSystem::ConstraintUserMeaningEnum::ACTIVITY,
1712 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1713 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1714 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1715 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1716 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1717 GeochemicalSystem::ConstraintUserMeaningEnum::FREE_MINERAL,
1718 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1719 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION};
1720 const std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu = {
1721 GeochemistryUnitConverter::GeochemistryUnit::KG,
1722 GeochemistryUnitConverter::GeochemistryUnit::DIMENSIONLESS,
1723 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1724 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1725 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1726 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1727 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1728 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1729 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1730 GeochemistryUnitConverter::GeochemistryUnit::MOLES};
1731 GeochemicalSystem egs(
1732 mgd,
1733 ac_solver,
1734 is_solver,
1735 swapper,
1736 {"Fe+++"},
1737 {"Fe(OH)3(ppd)"},
1738 "Cl-",
1739 {"H2O", "H+", "Na+", "Cl-", "Hg++", "Pb++", "SO4--", "Fe(OH)3(ppd)", ">(s)FeOH", ">(w)FeOH"},
1740 {1.0, 1.0E-4, 10E-3, 10E-3, 0.1E-3, 0.1E-3, 0.2E-3, 9.3573E-3, 4.6786E-5, 1.87145E-3},
1741 cu,
1742 cm,
1743 25,
1744 0,
1745 1E-20,
1746 {},
1747 {},
1748 {});
1749
1750 // build solver
1751 GeochemicalSolver solver(mgd.basis_species_name.size(),
1752 mgd.kin_species_name.size(),
1753 is_solver,
1754 1.0E-15,
1755 1.0E-200,
1756 100,
1757 1.0,
1758 0.1,
1759 1,
1760 {},
1761 3.0,
1762 10,
1763 false);
1764
1765 // solve
1766 std::stringstream ss;
1767 unsigned tot_iter;
1768 Real abs_residual;
1769 DenseVector<Real> mole_additions(egs.getNumInBasis() + egs.getNumKinetic());
1770 DenseMatrix<Real> dmole_additions(egs.getNumInBasis() + egs.getNumKinetic(),
1771 egs.getNumInBasis() + egs.getNumKinetic());
1772 solver.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
1773
1774 // check converged
1775 EXPECT_LE(tot_iter, (unsigned)100);
1776 EXPECT_LE(abs_residual, 1.0E-15);
1777
1778 // check number in basis, number in redox disequilibrium and number of surface potentials
1779 EXPECT_EQ(egs.getNumInBasis(), (unsigned)10);
1780 EXPECT_EQ(egs.getNumRedox(), (unsigned)0);
1781 EXPECT_EQ(egs.getNumSurfacePotentials(), (unsigned)1);
1782 EXPECT_EQ(egs.getNumInEquilibrium(), mgd.eqm_species_name.size());
1783
1784 // check that the constraints are satisfied
1785 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
1786 if (mgd.basis_species_name[i] == "H2O")
1787 {
1788 EXPECT_FALSE(mgd.basis_species_gas[i]);
1789 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1790 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1791 EXPECT_NEAR(egs.getSolventWaterMass(), 1.0, 1.0E-15);
1792 EXPECT_NEAR(egs.getSolventMassAndFreeMolalityAndMineralMoles()[0], 1.0, 1.0E-15);
1793 }
1794 else if (mgd.basis_species_name[i] == "H+")
1795 {
1796 EXPECT_FALSE(mgd.basis_species_gas[i]);
1797 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1798 EXPECT_TRUE(egs.getBasisActivityKnown()[i]);
1799 EXPECT_NEAR(egs.getBasisActivity(i), 1.0E-4, 1.0E-15);
1800 }
1801 else if (mgd.basis_species_name[i] == "Na+")
1802 {
1803 EXPECT_FALSE(mgd.basis_species_gas[i]);
1804 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1805 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1806 EXPECT_NEAR(egs.getBulkMolesOld()[i], 10E-3, 1.0E-15);
1807 }
1808 else if (mgd.basis_species_name[i] == "Cl-")
1809 {
1810 EXPECT_FALSE(mgd.basis_species_gas[i]);
1811 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1812 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1813 // do not know the bulk composition as it is dictated by charge neutrality
1814 EXPECT_EQ(egs.getChargeBalanceBasisIndex(), i);
1815 }
1816 else if (mgd.basis_species_name[i] == "Hg++")
1817 {
1818 EXPECT_FALSE(mgd.basis_species_gas[i]);
1819 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1820 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1821 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.1E-3, 1.0E-15);
1822 }
1823 else if (mgd.basis_species_name[i] == "Pb++")
1824 {
1825 EXPECT_FALSE(mgd.basis_species_gas[i]);
1826 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1827 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1828 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.1E-3, 1.0E-15);
1829 }
1830 else if (mgd.basis_species_name[i] == "SO4--")
1831 {
1832 EXPECT_FALSE(mgd.basis_species_gas[i]);
1833 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1834 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1835 EXPECT_NEAR(egs.getBulkMolesOld()[i], 0.2E-3, 1.0E-15);
1836 }
1837 else if (mgd.basis_species_name[i] == "Fe(OH)3(ppd)")
1838 {
1839 EXPECT_FALSE(mgd.basis_species_gas[i]);
1840 EXPECT_TRUE(mgd.basis_species_mineral[i]);
1841 EXPECT_TRUE(egs.getBasisActivityKnown()[i]);
1842 EXPECT_NEAR(egs.getSolventMassAndFreeMolalityAndMineralMoles()[i], 9.3573E-3, 1.0E-15);
1843 }
1844 else if (mgd.basis_species_name[i] == ">(s)FeOH")
1845 {
1846 EXPECT_FALSE(mgd.basis_species_gas[i]);
1847 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1848 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1849 EXPECT_NEAR(egs.getBulkMolesOld()[i], 4.6786E-5, 1.0E-15);
1850 }
1851 else if (mgd.basis_species_name[i] == ">(w)FeOH")
1852 {
1853 EXPECT_FALSE(mgd.basis_species_gas[i]);
1854 EXPECT_FALSE(mgd.basis_species_mineral[i]);
1855 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
1856 EXPECT_NEAR(egs.getBulkMolesOld()[i], 1.87145E-3, 1.0E-15);
1857 }
1858 else
1859 FAIL() << "Incorrect basis species";
1860
1861 // check total charge
1862 EXPECT_NEAR(egs.getTotalChargeOld(), 0.0, 1E-15);
1863 // check total charge by summing up basis charges
1864 Real tot_charge = 0.0;
1865 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
1866 tot_charge += egs.getBulkMolesOld()[i] * mgd.basis_species_charge[i];
1867 EXPECT_NEAR(tot_charge, 0.0, 1E-15);
1868
1869 // surface charges
1870 for (unsigned sp = 0; sp < egs.getNumSurfacePotentials(); ++sp)
1871 {
1872 Real charge = 0.0;
1873 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
1874 if (mgd.surface_sorption_related[j] && mgd.surface_sorption_number[j] == sp)
1875 charge += mgd.eqm_species_charge[j] * egs.getEquilibriumMolality(j);
1876 EXPECT_NEAR(charge,
1877 egs.getSorbingSurfaceArea()[sp] * egs.getSurfaceCharge(sp) /
1878 GeochemistryConstants::FARADAY,
1879 1E-15);
1880 }
1881
1882 // check basis activity = activity_coefficient * molality
1883 for (unsigned i = 1; i < egs.getNumInBasis(); ++i) // don't loop over water
1884 if (mgd.basis_species_gas[i] || mgd.basis_species_mineral[i] || egs.getBasisActivityKnown()[i])
1885 continue;
1886 else
1887 EXPECT_NEAR(egs.getBasisActivity(i),
1888 egs.getBasisActivityCoefficient(i) *
1889 egs.getSolventMassAndFreeMolalityAndMineralMoles()[i],
1890 1E-15);
1891
1892 // check residuals are zero
1893 for (unsigned a = 0; a < egs.getNumInAlgebraicSystem(); ++a)
1894 EXPECT_LE(std::abs(egs.getResidualComponent(a, mole_additions)), 1E-15);
1895 // check residuals are zero by summing molalities, bulk compositions, etc
1896 Real nw = egs.getSolventWaterMass();
1897 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
1898 {
1899 Real res = -egs.getBulkMolesOld()[i];
1900 if (i == 0)
1901 res += nw * GeochemistryConstants::MOLES_PER_KG_WATER;
1902 else if (mgd.basis_species_mineral[i])
1903 res += egs.getSolventMassAndFreeMolalityAndMineralMoles()[i];
1904 else if (mgd.basis_species_gas[i])
1905 res += 0.0;
1906 else
1907 res += nw * egs.getSolventMassAndFreeMolalityAndMineralMoles()[i];
1908 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
1909 res += nw * mgd.eqm_stoichiometry(j, i) * egs.getEquilibriumMolality(j);
1910 EXPECT_LE(std::abs(res), 1E-14);
1911 }
1912 // surface potentials
1913 const Real prefactor = std::sqrt(8.0 * GeochemistryConstants::GAS_CONSTANT *
1914 (25.0 + GeochemistryConstants::CELSIUS_TO_KELVIN) *
1915 GeochemistryConstants::PERMITTIVITY_FREE_SPACE *
1916 GeochemistryConstants::DIELECTRIC_CONSTANT_WATER *
1917 GeochemistryConstants::DENSITY_WATER * egs.getIonicStrength()) /
1918 nw;
1919 for (unsigned sp = 0; sp < egs.getNumSurfacePotentials(); ++sp)
1920 {
1921 EXPECT_NEAR(prefactor * std::sinh(egs.getSurfacePotential(sp) * GeochemistryConstants::FARADAY /
1922 2.0 / GeochemistryConstants::GAS_CONSTANT /
1923 (25.0 + GeochemistryConstants::CELSIUS_TO_KELVIN)),
1924 egs.getSurfaceCharge(sp),
1925 1.0E-13);
1926 }
1927
1928 // check equilibrium mass balance
1929 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
1930 {
1931 if (mgd.eqm_species_mineral[j] || mgd.eqm_species_gas[j])
1932 continue;
1933 Real log10ap = 0.0;
1934 for (unsigned i = 0; i < mgd.basis_species_name.size(); ++i)
1935 log10ap += mgd.eqm_stoichiometry(j, i) * std::log10(egs.getBasisActivity(i));
1936 log10ap -= egs.getLog10K(j);
1937 if (mgd.surface_sorption_related[j])
1938 log10ap -= std::log10(std::exp(mgd.eqm_species_charge[j] * GeochemistryConstants::FARADAY *
1939 egs.getSurfacePotential(mgd.surface_sorption_number[j]) /
1940 GeochemistryConstants::GAS_CONSTANT /
1941 (25.0 + GeochemistryConstants::CELSIUS_TO_KELVIN)));
1942 else
1943 log10ap -= std::log10(egs.getEquilibriumActivityCoefficient(j));
1944 if (log10ap < -300.0)
1945 EXPECT_LE(std::abs(egs.getEquilibriumMolality(j)), 1E-25);
1946 else
1947 EXPECT_LE(std::abs(egs.getEquilibriumMolality(j) - std::pow(10.0, log10ap)), 1E-15);
1948 }
1949
1950 // check equilibrium species have negative saturation indices, except for prevent_precipitation
1951 // minerals
1952 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
1953 {
1954 if (!mgd.eqm_species_mineral[j])
1955 continue;
1956 Real log10ap = 0.0;
1957 for (unsigned i = 0; i < mgd.basis_species_name.size(); ++i)
1958 log10ap += mgd.eqm_stoichiometry(j, i) * std::log10(egs.getBasisActivity(i));
1959 EXPECT_LE(log10ap, egs.getLog10K(j));
1960 }
1961}
db_ferric
const GeochemicalDatabaseReader db_ferric("../test/database/ferric_hydroxide_sorption.json", true, true)

◆ TEST() [12/16]

TEST ( GeochemicalSolverTest  ,
solve5_restore   
)

Repeat the realistic case with sorption and minerals, no redox, this time with a "restore", and then use the copy assignment operator of GeochemicalSystem to create a new geochemical system.

Definition at line 1964 of file GeochemicalSolverTest.C.

1965{
1966 // build the model
1967 const PertinentGeochemicalSystem model(
1968 db_ferric,
1969 {"H2O", "H+", "Na+", "Cl-", "Hg++", "Pb++", "SO4--", "Fe+++", ">(s)FeOH", ">(w)FeOH"},
1970 {"Fe(OH)3(ppd)"},
1971 {},
1972 {},
1973 {},
1974 {},
1975 "O2(aq)",
1976 "e-");
1977 ModelGeochemicalDatabase mgd = model.modelGeochemicalDatabase();
1978
1979 // build the equilibrium system
1980 GeochemistrySpeciesSwapper swapper(10, 1E-6);
1981 const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm = {
1982 GeochemicalSystem::ConstraintUserMeaningEnum::KG_SOLVENT_WATER,
1983 GeochemicalSystem::ConstraintUserMeaningEnum::ACTIVITY,
1984 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1985 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1986 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1987 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1988 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1989 GeochemicalSystem::ConstraintUserMeaningEnum::FREE_MINERAL,
1990 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
1991 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION};
1992 const std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu = {
1993 GeochemistryUnitConverter::GeochemistryUnit::KG,
1994 GeochemistryUnitConverter::GeochemistryUnit::DIMENSIONLESS,
1995 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1996 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1997 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1998 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
1999 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
2000 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
2001 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
2002 GeochemistryUnitConverter::GeochemistryUnit::MOLES};
2003 GeochemicalSystem egs(
2004 mgd,
2005 ac_solver,
2006 is_solver,
2007 swapper,
2008 {"Fe+++"},
2009 {"Fe(OH)3(ppd)"},
2010 "Cl-",
2011 {"H2O", "H+", "Na+", "Cl-", "Hg++", "Pb++", "SO4--", "Fe(OH)3(ppd)", ">(s)FeOH", ">(w)FeOH"},
2012 {1.0, 1.0E-4, 10E-3, 10E-3, 0.1E-3, 0.1E-3, 0.2E-3, 9.3573E-3, 4.6786E-5, 1.87145E-3},
2013 cu,
2014 cm,
2015 25,
2016 0,
2017 1E-20,
2018 {},
2019 {},
2020 {});
2021
2022 // build solver
2023 GeochemicalSolver solver(mgd.basis_species_name.size(),
2024 mgd.kin_species_name.size(),
2025 is_solver,
2026 1.0E-15,
2027 1.0E-200,
2028 100,
2029 1.0,
2030 0.1,
2031 1,
2032 {},
2033 3.0,
2034 0,
2035 false);
2036
2037 // solve
2038 std::stringstream ss;
2039 unsigned tot_iter;
2040 Real abs_residual;
2041 DenseVector<Real> mole_additions(egs.getNumInBasis() + egs.getNumKinetic());
2042 DenseMatrix<Real> dmole_additions(egs.getNumInBasis() + egs.getNumKinetic(),
2043 egs.getNumInBasis() + egs.getNumKinetic());
2044 solver.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
2045 const Real old_residual = abs_residual;
2046
2047 // check converged
2048 EXPECT_LE(tot_iter, (unsigned)100);
2049 EXPECT_LE(abs_residual, 1.0E-15);
2050
2051 // now setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles to the solution
2052 // obtained This should actually do nothing: we are just checking the "restore" doesn't muck
2053 // something up! Of course, activity coefficients will slightly change due to making the
2054 // configuration slightly more consistent (activities and molalities match) so there will be very
2055 // minor changes
2056 std::vector<std::string> names = mgd.basis_species_name;
2057 names.insert(names.end(), mgd.eqm_species_name.begin(), mgd.eqm_species_name.end());
2058 names.push_back("Fe(OH)3(ppd)_surface_potential_expr");
2059 std::vector<Real> molal = egs.getSolventMassAndFreeMolalityAndMineralMoles();
2060 for (unsigned j = 0; j < egs.getNumInEquilibrium(); ++j)
2061 molal.push_back(egs.getEquilibriumMolality(j));
2062 molal.push_back(std::exp(egs.getSurfacePotential(0) * GeochemistryConstants::FARADAY /
2063 (-2.0 * GeochemistryConstants::GAS_CONSTANT *
2064 (25.0 + GeochemistryConstants::CELSIUS_TO_KELVIN))));
2065 // form constraints_from_molalities, which are all false
2066 const std::vector<bool> com_false(egs.getNumInBasis(), false);
2067 egs.setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(
2068 names, molal, com_false);
2069
2070 solver.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
2071 EXPECT_EQ(tot_iter, (unsigned)0);
2072 EXPECT_LE(abs_residual, 2.0 * old_residual);
2073
2074 // Now use constraints_from_molalities = true
2075 const std::vector<bool> com_true(egs.getNumInBasis(), true);
2076 egs.setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(
2077 names, molal, com_true);
2078
2079 solver.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
2080
2081 EXPECT_EQ(tot_iter, (unsigned)0);
2082 EXPECT_LE(abs_residual, 2.0 * old_residual);
2083
2084 // now check the copy-assignment operator of GeochemicalSystem
2085
2086 // build the destination model
2087 const PertinentGeochemicalSystem dest_model(
2088 db_ferric,
2089 {"H2O", "H+", "Na+", "Cl-", "Hg++", "Pb++", "SO4--", "Fe+++", ">(s)FeOH", ">(w)FeOH"},
2090 {"Fe(OH)3(ppd)"},
2091 {},
2092 {},
2093 {},
2094 {},
2095 "O2(aq)",
2096 "e-");
2097 ModelGeochemicalDatabase dest_mgd = dest_model.modelGeochemicalDatabase();
2098
2099 // build the destination equilibrium system
2100 GeochemicalSystem dest_egs(
2101 dest_mgd,
2102 ac_solver,
2103 is_solver,
2104 swapper,
2105 {"Fe+++"},
2106 {"Fe(OH)3(ppd)"},
2107 "Cl-",
2108 {"H2O", "H+", "Na+", "Cl-", "Hg++", "Pb++", "SO4--", "Fe(OH)3(ppd)", ">(s)FeOH", ">(w)FeOH"},
2109 {2.0, 2.0E-4, 20E-3, 20E-3, 0.2E-3, 0.2E-3, 0.4E-3, 19.3573E-3, 14.6786E-5, 2.87145E-3},
2110 cu,
2111 cm,
2112 250,
2113 0,
2114 1E-20,
2115 {},
2116 {},
2117 {});
2118
2119 // change constraint meanings in egs to provide a more thorough check of the copy assignment
2120 egs.closeSystem();
2121 // copy assignment
2122 dest_egs = egs;
2123 EXPECT_EQ(dest_egs, egs);
2124}

◆ TEST() [13/16]

TEST ( GeochemicalSolverTest  ,
solve_addH   
)

Test progressively adding H+ to a system.

Definition at line 2127 of file GeochemicalSolverTest.C.

2128{
2129 ModelGeochemicalDatabase mgd = model_simplest.modelGeochemicalDatabase();
2130 GeochemicalSystem egs(mgd,
2131 ac_solver,
2132 is_solver,
2133 swapper2,
2134 {},
2135 {},
2136 "H+",
2137 {"H2O", "H+"},
2138 {1.75, 1},
2139 cu2,
2140 cm2,
2141 25,
2142 0,
2143 1E-20,
2144 {},
2145 {},
2146 {});
2147 GeochemicalSolver solver(mgd.basis_species_name.size(),
2148 mgd.kin_species_name.size(),
2149 is_solver,
2150 1.0E-12,
2151 1.0E-200,
2152 100,
2153 10.0,
2154 0.1,
2155 1,
2156 {},
2157 3.0,
2158 0,
2159 false);
2160
2161 std::stringstream ss;
2162 unsigned tot_iter;
2163 Real abs_residual;
2164 DenseVector<Real> mole_additions(egs.getNumInBasis() + egs.getNumKinetic());
2165 DenseMatrix<Real> dmole_additions(egs.getNumInBasis() + egs.getNumKinetic(),
2166 egs.getNumInBasis() + egs.getNumKinetic());
2167 solver.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
2168
2169 // check Newton has converged
2170 EXPECT_LE(tot_iter, (unsigned)100);
2171 EXPECT_LE(abs_residual, 1.0E-12);
2172
2173 // check Bulk is as expected
2174 const std::vector<Real> bulk = egs.getBulkMolesOld();
2175 EXPECT_NEAR(bulk[0],
2176 1.75 * (GeochemistryConstants::MOLES_PER_KG_WATER + egs.getEquilibriumMolality(0)),
2177 1.0E-9);
2178 EXPECT_NEAR(bulk[1], 0.0, 1.0E-9);
2179
2180 // close the system by setting to fixed number of moles
2181 egs.changeConstraintToBulk(0);
2182
2183 // check Bulk is as expected
2184 std::vector<Real> bulk_new = egs.getBulkMolesOld();
2185 for (unsigned i = 0; i < 2; ++i)
2186 EXPECT_EQ(bulk_new[i], bulk[i]);
2187
2188 // add 1 mol of water
2189 mole_additions(0) = 1.0;
2190 solver.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
2191
2192 // check Newton has converged
2193 EXPECT_LE(tot_iter, (unsigned)100);
2194 EXPECT_LE(abs_residual, 1.0E-12);
2195
2196 // check Bulk is as expected
2197 bulk_new = egs.getBulkMolesOld();
2198 EXPECT_EQ(bulk_new[0], bulk[0] + 1.0);
2199 EXPECT_EQ(bulk_new[1], bulk[1]);
2200
2201 // add 1 mol of H+ : this shouldn't make any difference because charge-balance instantly removes
2202 // it!
2203 mole_additions(0) = 0.0;
2204 mole_additions(1) = 1.0;
2205 solver.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
2206
2207 // check Newton has converged
2208 EXPECT_LE(tot_iter, (unsigned)100);
2209 EXPECT_LE(abs_residual, 1.0E-12);
2210
2211 // check Bulk is as expected
2212 bulk_new = egs.getBulkMolesOld();
2213 EXPECT_EQ(bulk_new[0], bulk[0] + 1.0);
2214 EXPECT_EQ(bulk_new[1], bulk[1]);
2215}

◆ TEST() [14/16]

TEST ( GeochemicalSolverTest  ,
solve_kinetic1   
)

Solve case that involves kinetic species with zero rates (so kinetic species should have no impact except to modify the bulk composition)

Definition at line 2389 of file GeochemicalSolverTest.C.

2390{
2391 const PertinentGeochemicalSystem model(db_solver,
2392 {"H2O", "H+", "Fe+++", "HCO3-"},
2393 {},
2394 {},
2395 {"Something", "Fe(OH)3(ppd)fake"},
2396 {},
2397 {},
2398 "O2(aq)",
2399 "e-");
2400 ModelGeochemicalDatabase mgd = model.modelGeochemicalDatabase();
2401 std::vector<GeochemistryUnitConverter::GeochemistryUnit> ku;
2402 ku.push_back(GeochemistryUnitConverter::GeochemistryUnit::MOLES);
2403 ku.push_back(GeochemistryUnitConverter::GeochemistryUnit::MOLES);
2404 GeochemicalSystem egs(mgd,
2405 ac_solver,
2406 is_solver,
2407 swapper_kin,
2408 {},
2409 {},
2410 "HCO3-",
2411 {"H2O", "H+", "Fe+++", "HCO3-"},
2412 {1.75, 1E-5, 1E-5, 4E-5},
2413 cu4,
2414 cm4,
2415 25,
2416 0,
2417 1E-20,
2418 {"Something", "Fe(OH)3(ppd)fake"},
2419 {1.0E-6, 2.0E-6},
2420 ku);
2421 GeochemicalSolver solver(mgd.basis_species_name.size(),
2422 mgd.kin_species_name.size(),
2423 is_solver,
2424 1.0E-15,
2425 1.0E-200,
2426 100,
2427 1E-5,
2428 1E-16,
2429 1,
2430 {},
2431 3.0,
2432 0,
2433 true);
2434
2435 std::stringstream ss;
2436 unsigned tot_iter;
2437 Real abs_residual;
2438 DenseVector<Real> mole_additions(egs.getNumInBasis() + egs.getNumKinetic());
2439 DenseMatrix<Real> dmole_additions(egs.getNumInBasis() + egs.getNumKinetic(),
2440 egs.getNumInBasis() + egs.getNumKinetic());
2441 solver.solveSystem(egs, ss, tot_iter, abs_residual, 0.0, mole_additions, dmole_additions);
2442
2443 // check Newton has converged
2444 EXPECT_LE(tot_iter, (unsigned)100);
2445 EXPECT_LE(abs_residual, 1.0E-15);
2446
2447 // check numbers are correct
2448 EXPECT_EQ(egs.getNumInBasis(), (unsigned)4);
2449 EXPECT_EQ(egs.getNumRedox(), (unsigned)0);
2450 EXPECT_EQ(egs.getNumSurfacePotentials(), (unsigned)0);
2451 EXPECT_EQ(egs.getNumInEquilibrium(), mgd.eqm_species_name.size());
2452 EXPECT_EQ(egs.getNumKinetic(), (unsigned)2);
2453
2454 // check that kinetic moles have not changed (kinetic rates are zero)
2455 for (unsigned k = 0; k < egs.getNumKinetic(); ++k)
2456 if (mgd.kin_species_name[k] == "Something")
2457 EXPECT_EQ(egs.getKineticMoles(k), 1.0E-6);
2458 else if (mgd.kin_species_name[k] == "Fe(OH)3(ppd)fake")
2459 EXPECT_EQ(egs.getKineticMoles(k), 2.0E-6);
2460 else
2461 FAIL() << "Incorrect kinetic species";
2462
2463 // check that the constraints are satisfied
2464 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
2465 if (mgd.basis_species_name[i] == "H2O")
2466 {
2467 EXPECT_FALSE(mgd.basis_species_gas[i]);
2468 EXPECT_FALSE(mgd.basis_species_mineral[i]);
2469 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
2470 EXPECT_NEAR(egs.getSolventWaterMass(), 1.75, 1.0E-15);
2471 EXPECT_NEAR(egs.getSolventMassAndFreeMolalityAndMineralMoles()[0], 1.75, 1.0E-15);
2472 }
2473 else if (mgd.basis_species_name[i] == "H+")
2474 {
2475 EXPECT_FALSE(mgd.basis_species_gas[i]);
2476 EXPECT_FALSE(mgd.basis_species_mineral[i]);
2477 EXPECT_TRUE(egs.getBasisActivityKnown()[i]);
2478 EXPECT_NEAR(egs.getBasisActivity(i), 1E-5, 1.0E-15);
2479 }
2480 else if (mgd.basis_species_name[i] == "Fe+++")
2481 {
2482 EXPECT_FALSE(mgd.basis_species_gas[i]);
2483 EXPECT_FALSE(mgd.basis_species_mineral[i]);
2484 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
2485 EXPECT_NEAR(egs.getBulkMolesOld()[i], 1E-5 + 1E-6 + 2 * 2E-6, 1.0E-15);
2486 }
2487 else if (mgd.basis_species_name[i] == "HCO3-")
2488 {
2489 EXPECT_FALSE(mgd.basis_species_gas[i]);
2490 EXPECT_FALSE(mgd.basis_species_mineral[i]);
2491 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
2492 // do not know the bulk composition as it is dictated by charge neutrality
2493 EXPECT_EQ(egs.getChargeBalanceBasisIndex(), i);
2494 }
2495 else
2496 FAIL() << "Incorrect basis species";
2497
2498 // check total charge
2499 EXPECT_NEAR(egs.getTotalChargeOld(), 0.0, 1E-15);
2500 // check total charge by summing up basis charges
2501 Real tot_charge = 0.0;
2502 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
2503 tot_charge += egs.getBulkMolesOld()[i] * mgd.basis_species_charge[i];
2504 EXPECT_NEAR(tot_charge, 0.0, 1E-15);
2505
2506 // check basis activity = activity_coefficient * molality
2507 for (unsigned i = 1; i < egs.getNumInBasis(); ++i) // don't loop over water
2508 if (mgd.basis_species_gas[i] || mgd.basis_species_mineral[i] || egs.getBasisActivityKnown()[i])
2509 continue;
2510 else
2511 EXPECT_NEAR(egs.getBasisActivity(i),
2512 egs.getBasisActivityCoefficient(i) *
2513 egs.getSolventMassAndFreeMolalityAndMineralMoles()[i],
2514 1E-15);
2515
2516 // check residuals are zero
2517 for (unsigned a = 0; a < egs.getNumInAlgebraicSystem(); ++a)
2518 EXPECT_LE(std::abs(egs.getResidualComponent(a, mole_additions)), 1E-15);
2519 // check residuals are zero by summing molalities, bulk compositions, etc
2520 Real nw = egs.getSolventWaterMass();
2521 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
2522 {
2523 Real res = -egs.getBulkMolesOld()[i];
2524 if (i == 0)
2525 res += nw * GeochemistryConstants::MOLES_PER_KG_WATER;
2526 else if (mgd.basis_species_mineral[i])
2527 res += egs.getSolventMassAndFreeMolalityAndMineralMoles()[i];
2528 else if (mgd.basis_species_gas[i])
2529 res += 0.0;
2530 else
2531 res += nw * egs.getSolventMassAndFreeMolalityAndMineralMoles()[i];
2532 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
2533 res += nw * mgd.eqm_stoichiometry(j, i) * egs.getEquilibriumMolality(j);
2534 for (unsigned k = 0; k < mgd.kin_species_name.size(); ++k)
2535 res += mgd.kin_stoichiometry(k, i) *
2536 egs.getKineticMoles(
2537 k); // this should be the only important kinetic contribution to this test!
2538 EXPECT_LE(std::abs(res), 1E-13);
2539 }
2540
2541 // check equilibrium mass balance
2542 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
2543 {
2544 if (mgd.eqm_species_mineral[j] || mgd.eqm_species_gas[j])
2545 continue;
2546 Real log10ap = 0.0;
2547 for (unsigned i = 0; i < mgd.basis_species_name.size(); ++i)
2548 log10ap += mgd.eqm_stoichiometry(j, i) * std::log10(egs.getBasisActivity(i));
2549 log10ap -= egs.getLog10K(j);
2550 if (mgd.surface_sorption_related[j])
2551 log10ap -= std::log10(egs.getSurfacePotential(mgd.surface_sorption_number[j]));
2552 else
2553 log10ap -= std::log10(egs.getEquilibriumActivityCoefficient(j));
2554 if (log10ap < -300.0)
2555 EXPECT_LE(std::abs(egs.getEquilibriumMolality(j)), 1E-25);
2556 else
2557 EXPECT_LE(std::abs(egs.getEquilibriumMolality(j) - std::pow(10.0, log10ap)), 1E-15);
2558 }
2559}
cu4
const std::vector< GeochemistryUnitConverter::GeochemistryUnit > cu4
Definition GeochemicalSolverTest.C:35
cm4
const std::vector< GeochemicalSystem::ConstraintUserMeaningEnum > cm4
Definition GeochemicalSolverTest.C:30
db_solver
const GeochemicalDatabaseReader db_solver("database/moose_testdb.json", true, true, false)
swapper_kin
GeochemistrySpeciesSwapper swapper_kin(4, 1E-6)
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

◆ TEST() [15/16]

TEST ( GeochemicalSolverTest  ,
solve_kinetic2   
)

Solve case that involves kinetic species with constant rates.

Definition at line 2562 of file GeochemicalSolverTest.C.

2563{
2564 PertinentGeochemicalSystem model(db_solver,
2565 {"H2O", "H+", "Fe+++", "HCO3-"},
2566 {},
2567 {},
2568 {"Something", "Fe(OH)3(ppd)fake"},
2569 {},
2570 {},
2571 "O2(aq)",
2572 "e-");
2573 // Define constant rates
2574 // Something = -3H+ + Fe+++ + 2H2O + 1.5HCO3-, so Q=1E15*1E-5*1E-7=1E3 (approx) >> K, so
2575 // rate_Something produces a negative rate, so mole_additions > 0, so Something will increase in
2576 // mole number
2577 KineticRateUserDescription rate_Something("Something",
2578 1.0E-7,
2579 1.0,
2580 false,
2581 0.0,
2582 0.0,
2583 0.0,
2584 {},
2585 {},
2586 {},
2587 {},
2588 1.0,
2589 0.0,
2590 0.0,
2591 0.0,
2592 DirectionChoiceEnum::BOTH,
2593 "H2O",
2594 0.0,
2595 -1.0,
2596 0.0);
2597 // Fe(OH)3(ppd)fake = -3H+ + 2Fe+++ + 3H2O, so Q=1E15*1E-10=1E5 (approx) < K, so rate_Fe produces
2598 // a positive rate, so mole_additions < 0, so Something will decrease in mole number
2599 KineticRateUserDescription rate_Fe("Fe(OH)3(ppd)fake",
2600 1.0E-7,
2601 1.0,
2602 false,
2603 0.0,
2604 0.0,
2605 0.0,
2606 {},
2607 {},
2608 {},
2609 {},
2610 1.0,
2611 0.0,
2612 0.0,
2613 0.0,
2614 DirectionChoiceEnum::BOTH,
2615 "H2O",
2616 0.0,
2617 -1.0,
2618 0.0);
2619 model.addKineticRate(rate_Something);
2620 model.addKineticRate(rate_Fe);
2621
2622 ModelGeochemicalDatabase mgd = model.modelGeochemicalDatabase();
2623 std::vector<GeochemistryUnitConverter::GeochemistryUnit> ku;
2624 ku.push_back(GeochemistryUnitConverter::GeochemistryUnit::MOLES);
2625 ku.push_back(GeochemistryUnitConverter::GeochemistryUnit::MOLES);
2626 GeochemicalSystem egs(mgd,
2627 ac_solver,
2628 is_solver,
2629 swapper_kin,
2630 {},
2631 {},
2632 "HCO3-",
2633 {"H2O", "H+", "Fe+++", "HCO3-"},
2634 {1.75, 1E-5, 1E-5, 4E-5},
2635 cu4,
2636 cm4,
2637 25,
2638 0,
2639 1E-20,
2640 {"Something", "Fe(OH)3(ppd)fake"},
2641 {1.0E-6, 2.0E-6},
2642 ku);
2643 GeochemicalSolver solver(mgd.basis_species_name.size(),
2644 mgd.kin_species_name.size(),
2645 is_solver,
2646 1.0E-15,
2647 1.0E-200,
2648 100,
2649 1E-5,
2650 1E-16,
2651 1,
2652 {},
2653 3.0,
2654 0,
2655 true);
2656
2657 const unsigned index_something = mgd.kin_species_index.at("Something");
2658 const unsigned index_fe = mgd.kin_species_index.at("Fe(OH)3(ppd)fake");
2659
2660 std::stringstream ss;
2661 unsigned tot_iter;
2662 Real abs_residual;
2663 DenseVector<Real> mole_additions(egs.getNumInBasis() + egs.getNumKinetic());
2664 mole_additions(egs.getNumInBasis() + index_something) = 1.5E-6; // external input to "Something"
2665 DenseMatrix<Real> dmole_additions(egs.getNumInBasis() + egs.getNumKinetic(),
2666 egs.getNumInBasis() + egs.getNumKinetic());
2667 // solve with a time-step size of 10.0
2668 solver.solveSystem(egs, ss, tot_iter, abs_residual, 10.0, mole_additions, dmole_additions);
2669
2670 // check Newton has converged
2671 EXPECT_LE(tot_iter, (unsigned)100);
2672 EXPECT_LE(abs_residual, 1.0E-15);
2673
2674 // check numbers are correct
2675 EXPECT_EQ(egs.getNumInBasis(), (unsigned)4);
2676 EXPECT_EQ(egs.getNumRedox(), (unsigned)0);
2677 EXPECT_EQ(egs.getNumSurfacePotentials(), (unsigned)0);
2678 EXPECT_EQ(egs.getNumInEquilibrium(), mgd.eqm_species_name.size());
2679 EXPECT_EQ(egs.getNumKinetic(), (unsigned)2);
2680
2681 // check activity products and log10K are ordered as calculated above
2682 EXPECT_GT(egs.log10KineticActivityProduct(index_something),
2683 egs.getKineticLog10K(index_something));
2684 EXPECT_LT(egs.log10KineticActivityProduct(index_fe), egs.getKineticLog10K(index_fe));
2685
2686 // check that mole additions is as calculated above
2687 EXPECT_EQ(mole_additions(egs.getNumInBasis() + index_something), 1.5E-6 + 1.0E-6);
2688 EXPECT_EQ(mole_additions(egs.getNumInBasis() + index_fe), -1.0E-6);
2689
2690 // check that kinetic moles have changed according to the constant rates
2691 for (unsigned k = 0; k < egs.getNumKinetic(); ++k)
2692 if (mgd.kin_species_name[k] == "Something")
2693 EXPECT_NEAR(egs.getKineticMoles(k), 2.0E-6 + 1.5E-6, 1.0E-12);
2694 else if (mgd.kin_species_name[k] == "Fe(OH)3(ppd)fake")
2695 EXPECT_EQ(egs.getKineticMoles(k), 1.0E-6);
2696 else
2697 FAIL() << "Incorrect kinetic species";
2698
2699 // check that the constraints are satisfied
2700 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
2701 if (mgd.basis_species_name[i] == "H2O")
2702 {
2703 EXPECT_FALSE(mgd.basis_species_gas[i]);
2704 EXPECT_FALSE(mgd.basis_species_mineral[i]);
2705 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
2706 EXPECT_NEAR(egs.getSolventWaterMass(), 1.75, 1.0E-15);
2707 EXPECT_NEAR(egs.getSolventMassAndFreeMolalityAndMineralMoles()[0], 1.75, 1.0E-15);
2708 }
2709 else if (mgd.basis_species_name[i] == "H+")
2710 {
2711 EXPECT_FALSE(mgd.basis_species_gas[i]);
2712 EXPECT_FALSE(mgd.basis_species_mineral[i]);
2713 EXPECT_TRUE(egs.getBasisActivityKnown()[i]);
2714 EXPECT_NEAR(egs.getBasisActivity(i), 1E-5, 1.0E-15);
2715 }
2716 else if (mgd.basis_species_name[i] == "Fe+++")
2717 {
2718 EXPECT_FALSE(mgd.basis_species_gas[i]);
2719 EXPECT_FALSE(mgd.basis_species_mineral[i]);
2720 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
2721 EXPECT_NEAR(egs.getBulkMolesOld()[i], 1E-5 + 1E-6 + 2 * 2E-6, 1.0E-15);
2722 }
2723 else if (mgd.basis_species_name[i] == "HCO3-")
2724 {
2725 EXPECT_FALSE(mgd.basis_species_gas[i]);
2726 EXPECT_FALSE(mgd.basis_species_mineral[i]);
2727 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
2728 // do not know the bulk composition as it is dictated by charge neutrality
2729 EXPECT_EQ(egs.getChargeBalanceBasisIndex(), i);
2730 }
2731 else
2732 FAIL() << "Incorrect basis species";
2733
2734 // check total charge
2735 EXPECT_NEAR(egs.getTotalChargeOld(), 0.0, 1E-15);
2736 // check total charge by summing up basis charges
2737 Real tot_charge = 0.0;
2738 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
2739 tot_charge += egs.getBulkMolesOld()[i] * mgd.basis_species_charge[i];
2740 EXPECT_NEAR(tot_charge, 0.0, 1E-15);
2741
2742 // check basis activity = activity_coefficient * molality
2743 for (unsigned i = 1; i < egs.getNumInBasis(); ++i) // don't loop over water
2744 if (mgd.basis_species_gas[i] || mgd.basis_species_mineral[i] || egs.getBasisActivityKnown()[i])
2745 continue;
2746 else
2747 EXPECT_NEAR(egs.getBasisActivity(i),
2748 egs.getBasisActivityCoefficient(i) *
2749 egs.getSolventMassAndFreeMolalityAndMineralMoles()[i],
2750 1E-15);
2751
2752 // check residuals are zero
2753 mole_additions.zero(); // to remove the kinetic rate contributions
2754 for (unsigned a = 0; a < egs.getNumInAlgebraicSystem(); ++a)
2755 EXPECT_LE(std::abs(egs.getResidualComponent(a, mole_additions)), 1E-15);
2756 // check residuals are zero by summing molalities, bulk compositions, etc
2757 Real nw = egs.getSolventWaterMass();
2758 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
2759 {
2760 Real res = -egs.getBulkMolesOld()[i];
2761 if (i == 0)
2762 res += nw * GeochemistryConstants::MOLES_PER_KG_WATER;
2763 else if (mgd.basis_species_mineral[i])
2764 res += egs.getSolventMassAndFreeMolalityAndMineralMoles()[i];
2765 else if (mgd.basis_species_gas[i])
2766 res += 0.0;
2767 else
2768 res += nw * egs.getSolventMassAndFreeMolalityAndMineralMoles()[i];
2769 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
2770 res += nw * mgd.eqm_stoichiometry(j, i) * egs.getEquilibriumMolality(j);
2771 for (unsigned k = 0; k < mgd.kin_species_name.size(); ++k)
2772 res += mgd.kin_stoichiometry(k, i) * egs.getKineticMoles(k);
2773 EXPECT_LE(std::abs(res), 1E-13);
2774 }
2775
2776 // check equilibrium mass balance
2777 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
2778 {
2779 if (mgd.eqm_species_mineral[j] || mgd.eqm_species_gas[j])
2780 continue;
2781 Real log10ap = 0.0;
2782 for (unsigned i = 0; i < mgd.basis_species_name.size(); ++i)
2783 log10ap += mgd.eqm_stoichiometry(j, i) * std::log10(egs.getBasisActivity(i));
2784 log10ap -= egs.getLog10K(j);
2785 if (mgd.surface_sorption_related[j])
2786 log10ap -= std::log10(egs.getSurfacePotential(mgd.surface_sorption_number[j]));
2787 else
2788 log10ap -= std::log10(egs.getEquilibriumActivityCoefficient(j));
2789 if (log10ap < -300.0)
2790 EXPECT_LE(std::abs(egs.getEquilibriumMolality(j)), 1E-25);
2791 else
2792 EXPECT_LE(std::abs(egs.getEquilibriumMolality(j) - std::pow(10.0, log10ap)), 1E-15);
2793 }
2794}
PertinentGeochemicalSystem::addKineticRate
void addKineticRate(const KineticRateUserDescription &description)
Adds a rate description for kinetic_species.
Definition PertinentGeochemicalSystem.C:856
KineticRateUserDescription
Holds a user-specified description of a kinetic rate.
Definition GeochemistryKineticRateCalculator.h:76
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

◆ TEST() [16/16]

TEST ( GeochemicalSolverTest  ,
solve_kinetic3   
)

Solve case that involves kinetic species with promoting indices and implicit solve.

Definition at line 2797 of file GeochemicalSolverTest.C.

2798{
2799 PertinentGeochemicalSystem model(
2800 db_solver, {"H2O", "H+", "Fe+++", "HCO3-"}, {}, {}, {"Something"}, {}, {}, "O2(aq)", "e-");
2801 // Define rate. The following produces
2802 // rate = 1.0E-2 * moles_something * 88.8537 * a_{H+} * exp(1E5 / R * (1/303.15 - 1/T))
2803 // = 4.56796E-06 * moles_something
2804 // Something = -3H+ + Fe+++ + 2H2O + 1.5HCO3-, so Q=1E15*1E-5*1E-7=1E3 (approx) >> K, so
2805 // rate_Something produces a negative rate, so mole_additions > 0, so Something will increase in
2806 // mole number
2807 KineticRateUserDescription rate_Something("Something",
2808 1.0E-2,
2809 1.0,
2810 true,
2811 0.0,
2812 0.0,
2813 0.0,
2814 {"H+"},
2815 {1.0},
2816 {0.0},
2817 {0.0},
2818 1.0,
2819 0.0,
2820 1E5,
2821 1.0 / 303.15,
2822 DirectionChoiceEnum::BOTH,
2823 "H2O",
2824 0.0,
2825 -1.0,
2826 0.0);
2827 model.addKineticRate(rate_Something);
2828
2829 ModelGeochemicalDatabase mgd = model.modelGeochemicalDatabase();
2830 std::vector<GeochemistryUnitConverter::GeochemistryUnit> ku;
2831 ku.push_back(GeochemistryUnitConverter::GeochemistryUnit::MOLES);
2832 GeochemicalSystem egs(mgd,
2833 ac_solver,
2834 is_solver,
2835 swapper_kin,
2836 {},
2837 {},
2838 "HCO3-",
2839 {"H2O", "H+", "Fe+++", "HCO3-"},
2840 {1.75, 1E-5, 1E-5, 4E-5},
2841 cu4,
2842 cm4,
2843 25,
2844 0,
2845 1E-20,
2846 {"Something"},
2847 {1.0E-6},
2848 ku);
2849 GeochemicalSolver solver(mgd.basis_species_name.size(),
2850 mgd.kin_species_name.size(),
2851 is_solver,
2852 1.0E-15,
2853 1.0E-200,
2854 100,
2855 1E-5,
2856 1E-16,
2857 1,
2858 {},
2859 3.0,
2860 0,
2861 true);
2862
2863 std::stringstream ss;
2864 unsigned tot_iter;
2865 Real abs_residual;
2866 DenseVector<Real> mole_additions(egs.getNumInBasis() + egs.getNumKinetic());
2867 DenseMatrix<Real> dmole_additions(egs.getNumInBasis() + egs.getNumKinetic(),
2868 egs.getNumInBasis() + egs.getNumKinetic());
2869 // solve with a time-step size of 1E4
2870 solver.solveSystem(egs, ss, tot_iter, abs_residual, 1.0E4, mole_additions, dmole_additions);
2871
2872 // check Newton has converged
2873 EXPECT_LE(tot_iter, (unsigned)100);
2874 EXPECT_LE(abs_residual, 1.0E-15);
2875
2876 // check numbers are correct
2877 EXPECT_EQ(egs.getNumInBasis(), (unsigned)4);
2878 EXPECT_EQ(egs.getNumRedox(), (unsigned)0);
2879 EXPECT_EQ(egs.getNumSurfacePotentials(), (unsigned)0);
2880 EXPECT_EQ(egs.getNumInEquilibrium(), mgd.eqm_species_name.size());
2881 EXPECT_EQ(egs.getNumKinetic(), (unsigned)1);
2882
2883 // check activity products and log10K are ordered as calculated above
2884 const unsigned index_something = 0;
2885 EXPECT_GT(egs.log10KineticActivityProduct(index_something),
2886 egs.getKineticLog10K(index_something));
2887
2888 // check that mole additions is as calculated above
2889 EXPECT_NEAR(mole_additions(egs.getNumInBasis() + index_something) /
2890 (1.0E4 * 4.56796204026978E-06 * egs.getKineticMoles(0)),
2891 1.0,
2892 1.0E-8);
2893
2894 // check that kinetic moles have changed according to the constant rates
2895 for (unsigned k = 0; k < egs.getNumKinetic(); ++k)
2896 if (mgd.kin_species_name[k] == "Something")
2897 EXPECT_NEAR(egs.getKineticMoles(k),
2898 1.0E-6 + mole_additions(egs.getNumInBasis() + index_something),
2899 1.0E-12);
2900 else
2901 FAIL() << "Incorrect kinetic species";
2902
2903 // check that the constraints are satisfied
2904 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
2905 if (mgd.basis_species_name[i] == "H2O")
2906 {
2907 EXPECT_FALSE(mgd.basis_species_gas[i]);
2908 EXPECT_FALSE(mgd.basis_species_mineral[i]);
2909 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
2910 EXPECT_NEAR(egs.getSolventWaterMass(), 1.75, 1.0E-15);
2911 EXPECT_NEAR(egs.getSolventMassAndFreeMolalityAndMineralMoles()[0], 1.75, 1.0E-15);
2912 }
2913 else if (mgd.basis_species_name[i] == "H+")
2914 {
2915 EXPECT_FALSE(mgd.basis_species_gas[i]);
2916 EXPECT_FALSE(mgd.basis_species_mineral[i]);
2917 EXPECT_TRUE(egs.getBasisActivityKnown()[i]);
2918 EXPECT_NEAR(egs.getBasisActivity(i), 1E-5, 1.0E-15);
2919 }
2920 else if (mgd.basis_species_name[i] == "Fe+++")
2921 {
2922 EXPECT_FALSE(mgd.basis_species_gas[i]);
2923 EXPECT_FALSE(mgd.basis_species_mineral[i]);
2924 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
2925 EXPECT_NEAR(egs.getBulkMolesOld()[i], 1E-5 + 1E-6, 1E-15); // 1E-6 from Something
2926 }
2927 else if (mgd.basis_species_name[i] == "HCO3-")
2928 {
2929 EXPECT_FALSE(mgd.basis_species_gas[i]);
2930 EXPECT_FALSE(mgd.basis_species_mineral[i]);
2931 EXPECT_FALSE(egs.getBasisActivityKnown()[i]);
2932 // do not know the bulk composition as it is dictated by charge neutrality
2933 EXPECT_EQ(egs.getChargeBalanceBasisIndex(), i);
2934 }
2935 else
2936 FAIL() << "Incorrect basis species";
2937
2938 // check total charge
2939 EXPECT_NEAR(egs.getTotalChargeOld(), 0.0, 1E-15);
2940 // check total charge by summing up basis charges
2941 Real tot_charge = 0.0;
2942 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
2943 tot_charge += egs.getBulkMolesOld()[i] * mgd.basis_species_charge[i];
2944 EXPECT_NEAR(tot_charge, 0.0, 1E-15);
2945
2946 // check basis activity = activity_coefficient * molality
2947 for (unsigned i = 1; i < egs.getNumInBasis(); ++i) // don't loop over water
2948 if (mgd.basis_species_gas[i] || mgd.basis_species_mineral[i] || egs.getBasisActivityKnown()[i])
2949 continue;
2950 else
2951 EXPECT_NEAR(egs.getBasisActivity(i),
2952 egs.getBasisActivityCoefficient(i) *
2953 egs.getSolventMassAndFreeMolalityAndMineralMoles()[i],
2954 1E-15);
2955
2956 // check residuals are zero
2957 mole_additions.zero(); // to remove the kinetic rate contributions
2958 for (unsigned a = 0; a < egs.getNumInAlgebraicSystem(); ++a)
2959 EXPECT_LE(std::abs(egs.getResidualComponent(a, mole_additions)), 1E-15);
2960 // check residuals are zero by summing molalities, bulk compositions, etc
2961 Real nw = egs.getSolventWaterMass();
2962 for (unsigned i = 0; i < egs.getNumInBasis(); ++i)
2963 {
2964 Real res = -egs.getBulkMolesOld()[i];
2965 if (i == 0)
2966 res += nw * GeochemistryConstants::MOLES_PER_KG_WATER;
2967 else if (mgd.basis_species_mineral[i])
2968 res += egs.getSolventMassAndFreeMolalityAndMineralMoles()[i];
2969 else if (mgd.basis_species_gas[i])
2970 res += 0.0;
2971 else
2972 res += nw * egs.getSolventMassAndFreeMolalityAndMineralMoles()[i];
2973 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
2974 res += nw * mgd.eqm_stoichiometry(j, i) * egs.getEquilibriumMolality(j);
2975 for (unsigned k = 0; k < mgd.kin_species_name.size(); ++k)
2976 res += mgd.kin_stoichiometry(k, i) * egs.getKineticMoles(k);
2977 EXPECT_LE(std::abs(res), 1E-13);
2978 }
2979
2980 // check equilibrium mass balance
2981 for (unsigned j = 0; j < mgd.eqm_species_name.size(); ++j)
2982 {
2983 if (mgd.eqm_species_mineral[j] || mgd.eqm_species_gas[j])
2984 continue;
2985 Real log10ap = 0.0;
2986 for (unsigned i = 0; i < mgd.basis_species_name.size(); ++i)
2987 log10ap += mgd.eqm_stoichiometry(j, i) * std::log10(egs.getBasisActivity(i));
2988 log10ap -= egs.getLog10K(j);
2989 if (mgd.surface_sorption_related[j])
2990 log10ap -= std::log10(egs.getSurfacePotential(mgd.surface_sorption_number[j]));
2991 else
2992 log10ap -= std::log10(egs.getEquilibriumActivityCoefficient(j));
2993 if (log10ap < -300.0)
2994 EXPECT_LE(std::abs(egs.getEquilibriumMolality(j)), 1E-25);
2995 else
2996 EXPECT_LE(std::abs(egs.getEquilibriumMolality(j) - std::pow(10.0, log10ap)), 1E-15);
2997 }
2998}

Variable Documentation

◆ ac_solver

GeochemistryActivityCoefficientsDebyeHuckel ac_solver(is_solver, db_solver) ( is_solver  ,
db_solver   
)

Referenced by TEST(), TEST(), TEST(), TEST(), TEST(), TEST(), TEST(), TEST(), TEST(), TEST(), TEST(), TEST(), TEST(), TEST(), TEST(), and TEST().

◆ cm2

const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm2
Initial value:
= {
GeochemicalSystem::ConstraintUserMeaningEnum::KG_SOLVENT_WATER,
GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION}

Definition at line 24 of file GeochemicalSolverTest.C.

24 {
25 GeochemicalSystem::ConstraintUserMeaningEnum::KG_SOLVENT_WATER,
26 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION};

Referenced by TEST(), TEST(), TEST(), TEST(), TEST(), TEST_F(), TEST_F(), TEST_F(), TEST_F(), and TEST_F().

◆ cm4

const std::vector<GeochemicalSystem::ConstraintUserMeaningEnum> cm4
Initial value:
= {
GeochemicalSystem::ConstraintUserMeaningEnum::KG_SOLVENT_WATER,
GeochemicalSystem::ConstraintUserMeaningEnum::ACTIVITY,
GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION_WITH_KINETIC}

Definition at line 30 of file GeochemicalSolverTest.C.

30 {
31 GeochemicalSystem::ConstraintUserMeaningEnum::KG_SOLVENT_WATER,
32 GeochemicalSystem::ConstraintUserMeaningEnum::ACTIVITY,
33 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION,
34 GeochemicalSystem::ConstraintUserMeaningEnum::BULK_COMPOSITION_WITH_KINETIC};

Referenced by TEST(), TEST(), and TEST().

◆ cu2

const std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu2
Initial value:
= {
GeochemistryUnitConverter::GeochemistryUnit::KG,
GeochemistryUnitConverter::GeochemistryUnit::MOLES}

Definition at line 27 of file GeochemicalSolverTest.C.

27 {
28 GeochemistryUnitConverter::GeochemistryUnit::KG,
29 GeochemistryUnitConverter::GeochemistryUnit::MOLES};

Referenced by TEST(), TEST(), TEST(), TEST(), TEST(), TEST_F(), TEST_F(), TEST_F(), TEST_F(), and TEST_F().

◆ cu4

const std::vector<GeochemistryUnitConverter::GeochemistryUnit> cu4
Initial value:
= {
GeochemistryUnitConverter::GeochemistryUnit::KG,
GeochemistryUnitConverter::GeochemistryUnit::DIMENSIONLESS,
GeochemistryUnitConverter::GeochemistryUnit::MOLES,
GeochemistryUnitConverter::GeochemistryUnit::MOLES}

Definition at line 35 of file GeochemicalSolverTest.C.

35 {
36 GeochemistryUnitConverter::GeochemistryUnit::KG,
37 GeochemistryUnitConverter::GeochemistryUnit::DIMENSIONLESS,
38 GeochemistryUnitConverter::GeochemistryUnit::MOLES,
39 GeochemistryUnitConverter::GeochemistryUnit::MOLES};

Referenced by TEST(), TEST(), and TEST().

◆ db_ferric

const GeochemicalDatabaseReader db_ferric("../test/database/ferric_hydroxide_sorption.json", true, true) ( "../test/database/ferric_hydroxide_sorption.json"  ,
true  ,
true   
)

Referenced by TEST(), and TEST().

◆ db_full

const GeochemicalDatabaseReader db_full("../database/moose_geochemdb.json", true, true, false) ( "../database/moose_geochemdb.json"  ,
true  ,
true  ,
false   
)

Referenced by TEST(), TEST(), TEST(), TEST(), TEST(), and TEST().

◆ db_solver

const GeochemicalDatabaseReader db_solver("database/moose_testdb.json", true, true, false) ( "database/moose_testdb.json"  ,
true  ,
true  ,
false   
)

Referenced by TEST(), TEST(), and TEST().

◆ is_solver

GeochemistryIonicStrength is_solver(3.0, 3.0, false, false) ( 3.  0,
3.  0,
false  ,
false   
)

Referenced by TEST(), TEST(), TEST(), TEST(), TEST(), TEST(), TEST(), TEST(), TEST(), TEST(), TEST(), TEST(), TEST(), TEST(), TEST(), and TEST().

◆ model_simplest

const PertinentGeochemicalSystem model_simplest(db_solver, {"H2O", "H+"}, {}, {}, {}, {}, {}, "O2(aq)", "e-") ( db_solver  ,
{"H2O", "H+"}  ,
{}  ,
{}  ,
{}  ,
{}  ,
{}  ,
"O2(aq)"  ,
"e-"   
)

Referenced by TEST(), TEST(), TEST(), TEST(), and TEST().

◆ swapper2

GeochemistrySpeciesSwapper swapper2(2, 1E-6) ( 2  ,
1E-  6 
)

Referenced by TEST(), TEST(), TEST(), TEST(), and TEST().

◆ swapper_kin

GeochemistrySpeciesSwapper swapper_kin(4, 1E-6) ( 4  ,
1E-  6 
)

Referenced by TEST(), TEST(), and TEST().

Generated on Fri Sep 11 2026 20:37:17 for https://mooseframework.inl.gov by doxygen 1.9.8