https://mooseframework.inl.gov
GeochemistryKineticRateCalculatorTest.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://mooseframework.inl.gov
3 //*
4 //* All rights reserved, see COPYRIGHT for full restrictions
5 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 //*
7 //* Licensed under LGPL 2.1, please see LICENSE for details
8 //* https://www.gnu.org/licenses/lgpl-2.1.html
9 
10 #include "gtest/gtest.h"
11 
12 #include "GeochemicalSystem.h"
13 #include "GeochemistryKineticRateCalculator.h"
14 
15 const GeochemicalDatabaseReader db_kin("database/moose_testdb.json");
16 PertinentGeochemicalSystem model_kin(db_kin,
17  {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++"},
18  {},
19  {},
20  {"Calcite"},
21  {"CH4(aq)"},
22  {},
23  "O2(aq)",
24  "e-");
25 KineticRateUserDescription rate_ch4("CH4(aq)",
26  1.5,
27  2.0,
28  true,
29  0.0,
30  0.0,
31  0.0,
32  {"H2O", "OH-", "O2(aq)", "CO2(aq)", "CaCO3"},
33  {3.0, 3.1, 3.2, 3.3, 3.4},
34  {0.0, 0.0, 0.0, 0.0, 0.0},
35  {0.0, 0.0, 0.0, 0.0, 0.0},
36  0.8,
37  2.5,
38  66.0,
39  0.003,
40  DirectionChoiceEnum::BOTH,
41  "H2O",
42  0.0,
43  -1.0,
44  0.0);
45 KineticRateUserDescription rate_ch4_dissolution("CH4(aq)",
46  1.5,
47  2.0,
48  true,
49  0.0,
50  0.0,
51  0.0,
52  {"H2O", "OH-", "O2(aq)", "CO2(aq)", "CaCO3"},
53  {3.0, 3.1, 3.2, 3.3, 3.4},
54  {0.0, 0.0, 0.0, 0.0, 0.0},
55  {0.0, 0.0, 0.0, 0.0, 0.0},
56  0.8,
57  2.5,
58  66.0,
59  0.003,
60  DirectionChoiceEnum::DISSOLUTION,
61  "H2O",
62  0.0,
63  -1.0,
64  0.0);
65 KineticRateUserDescription rate_ch4_precipitation("CH4(aq)",
66  1.5,
67  2.0,
68  true,
69  0.0,
70  0.0,
71  0.0,
72  {"H2O", "OH-", "O2(aq)", "CO2(aq)", "CaCO3"},
73  {3.0, 3.1, 3.2, 3.3, 3.4},
74  {0.0, 0.0, 0.0, 0.0, 0.0},
75  {0.0, 0.0, 0.0, 0.0, 0.0},
76  0.8,
77  2.5,
78  66.0,
79  0.003,
80  DirectionChoiceEnum::PRECIPITATION,
81  "H2O",
82  0.0,
83  -1.0,
84  0.0);
85 KineticRateUserDescription rate_ch4_raw("CH4(aq)",
86  1.5,
87  2.0,
88  true,
89  0.0,
90  0.0,
91  0.0,
92  {"H2O", "OH-", "O2(aq)", "CO2(aq)", "CaCO3"},
93  {3.0, 3.1, 3.2, 3.3, 3.4},
94  {0.0, 0.0, 0.0, 0.0, 0.0},
95  {0.0, 0.0, 0.0, 0.0, 0.0},
96  0.8,
97  2.5,
98  66.0,
99  0.003,
100  DirectionChoiceEnum::RAW,
101  "H2O",
102  0.0,
103  -1.0,
104  0.0);
105 KineticRateUserDescription rate_ch4_death("CH4(aq)",
106  1.5,
107  2.0,
108  true,
109  0.0,
110  0.0,
111  0.0,
112  {"H2O", "OH-", "O2(aq)", "CO2(aq)", "CaCO3"},
113  {3.0, 3.1, 3.2, 3.3, 3.4},
114  {0.0, 0.0, 0.0, 0.0, 0.0},
115  {0.0, 0.0, 0.0, 0.0, 0.0},
116  0.8,
117  2.5,
118  66.0,
119  0.003,
120  DirectionChoiceEnum::DEATH,
121  "H2O",
122  0.0,
123  -1.0,
124  0.0);
125 KineticRateUserDescription rate_ch4_kin_mon("CH4(aq)",
126  1.5,
127  2.0,
128  true,
129  1.75,
130  2.125,
131  0.875,
132  {"H2O", "OH-", "O2(aq)", "CO2(aq)", "CaCO3"},
133  {3.0, 3.1, 3.2, 3.3, 3.4},
134  {0.0, 0.0, 0.0, 0.0, 0.0},
135  {0.0, 0.0, 0.0, 0.0, 0.0},
136  0.8,
137  2.5,
138  66.0,
139  0.003,
140  DirectionChoiceEnum::BOTH,
141  "H2O",
142  0.0,
143  -1.0,
144  0.0);
145 KineticRateUserDescription rate_ch4_energy_captured("CH4(aq)",
146  1.5,
147  2.0,
148  true,
149  0.0,
150  0.0,
151  0.0,
152  {"H2O", "OH-", "O2(aq)", "CO2(aq)", "CaCO3"},
153  {3.0, 3.1, 3.2, 3.3, 3.4},
154  {0.0, 0.0, 0.0, 0.0, 0.0},
155  {0.0, 0.0, 0.0, 0.0, 0.0},
156  0.8,
157  2.5,
158  66.0,
159  0.003,
160  DirectionChoiceEnum::BOTH,
161  "H2O",
162  0.0,
163  -1.0,
164  1.0);
165 KineticRateUserDescription rate_cal("Calcite",
166  7.0,
167  6.0,
168  false,
169  0.0,
170  0.0,
171  0.0,
172  {"H+"},
173  {-3.0},
174  {0.0},
175  {0.0},
176  2.5,
177  0.8,
178  55.0,
179  0.00315,
180  DirectionChoiceEnum::BOTH,
181  "H2O",
182  0.0,
183  -1.0,
184  0.0);
185 KineticRateUserDescription rate_cal_theta0_eta2("Calcite",
186  7.0,
187  6.0,
188  false,
189  0.0,
190  0.0,
191  0.0,
192  {"H+"},
193  {-3.0},
194  {0.0},
195  {0.0},
196  0.0,
197  2.0,
198  55.0,
199  0.00315,
200  DirectionChoiceEnum::BOTH,
201  "H2O",
202  0.0,
203  -1.0,
204  0.0);
205 KineticRateUserDescription rate_cal_theta0_eta1("Calcite",
206  7.0,
207  6.0,
208  false,
209  0.0,
210  0.0,
211  0.0,
212  {"H+"},
213  {-3.0},
214  {0.0},
215  {0.0},
216  0.0,
217  1.0,
218  55.0,
219  0.00315,
220  DirectionChoiceEnum::BOTH,
221  "H2O",
222  0.0,
223  -1.0,
224  0.0);
225 KineticRateUserDescription rate_cal_theta2_eta1("Calcite",
226  7.0,
227  6.0,
228  false,
229  0.0,
230  0.0,
231  0.0,
232  {"H+"},
233  {-3.0},
234  {0.0},
235  {0.0},
236  2.0,
237  1.0,
238  55.0,
239  0.00315,
240  DirectionChoiceEnum::BOTH,
241  "H2O",
242  0.0,
243  -1.0,
244  0.0);
245 KineticRateUserDescription rate_cal_theta0_eta0("Calcite",
246  7.0,
247  6.0,
248  false,
249  0.0,
250  0.0,
251  0.0,
252  {"H+"},
253  {-3.0},
254  {0.0},
255  {0.0},
256  0.0,
257  0.0,
258  55.0,
259  0.00315,
260  DirectionChoiceEnum::BOTH,
261  "H2O",
262  0.0,
263  -1.0,
264  0.0);
265 KineticRateUserDescription rate_cal_theta0_eta05("Calcite",
266  7.0,
267  6.0,
268  false,
269  0.0,
270  0.0,
271  0.0,
272  {"H+"},
273  {-3.0},
274  {0.0},
275  {0.0},
276  0.0,
277  0.5,
278  55.0,
279  0.00315,
280  DirectionChoiceEnum::BOTH,
281  "H2O",
282  0.0,
283  -1.0,
284  0.0);
285 const ModelGeochemicalDatabase & mgd_kin = model_kin.modelGeochemicalDatabase();
286 
288 TEST(KineticRateUserDescriptionTest, exceptions)
289 {
290  try
291  {
292  KineticRateUserDescription rate("CH4(aq)",
293  1.0,
294  2.0,
295  true,
296  0.0,
297  0.0,
298  0.0,
299  {"H2O", "H+"},
300  {3.0},
301  {0.0, 0.0},
302  {0.0, 0.0},
303  4.0,
304  5.0,
305  6.0,
306  7.0,
307  DirectionChoiceEnum::BOTH,
308  "H2O",
309  0.0,
310  -1.0,
311  0.0);
312  FAIL() << "Missing expected exception.";
313  }
314  catch (const std::exception & e)
315  {
316  std::string msg(e.what());
317  ASSERT_TRUE(
318  msg.find("The promoting_species and promoting_indices vectors must be the same size") !=
319  std::string::npos)
320  << "Failed with unexpected error message: " << msg;
321  }
322 
323  try
324  {
325  KineticRateUserDescription rate("CH4(aq)",
326  1.0,
327  2.0,
328  true,
329  0.0,
330  0.0,
331  0.0,
332  {"H2O", "H+"},
333  {3.0, 0.0},
334  {3.0},
335  {0.0, 0.0},
336  4.0,
337  5.0,
338  6.0,
339  7.0,
340  DirectionChoiceEnum::BOTH,
341  "H2O",
342  0.0,
343  -1.0,
344  0.0);
345  FAIL() << "Missing expected exception.";
346  }
347  catch (const std::exception & e)
348  {
349  std::string msg(e.what());
350  ASSERT_TRUE(
351  msg.find(
352  "The promoting_species and promoting_monod_indices vectors must be the same size") !=
353  std::string::npos)
354  << "Failed with unexpected error message: " << msg;
355  }
356 
357  try
358  {
359  KineticRateUserDescription rate("CH4(aq)",
360  1.0,
361  2.0,
362  true,
363  0.0,
364  0.0,
365  0.0,
366  {"H2O", "H+"},
367  {3.0, 0.0},
368  {3.0, 0.0},
369  {1.0},
370  4.0,
371  5.0,
372  6.0,
373  7.0,
374  DirectionChoiceEnum::BOTH,
375  "H2O",
376  0.0,
377  -1.0,
378  0.0);
379  FAIL() << "Missing expected exception.";
380  }
381  catch (const std::exception & e)
382  {
383  std::string msg(e.what());
384  ASSERT_TRUE(
385  msg.find(
386  "The promoting_species and promoting_half_saturation vectors must be the same size") !=
387  std::string::npos)
388  << "Failed with unexpected error message: " << msg;
389  }
390 
391  try
392  {
393  KineticRateUserDescription rate("CH4(aq)",
394  1.0,
395  2.0,
396  true,
397  0.0,
398  0.0,
399  0.0,
400  {"H2O", "OH-", "H2O"},
401  {3.0, 1.0, -1.0},
402  {0.0, 0.0, 0.0},
403  {0.0, 0.0, 0.0},
404  4.0,
405  5.0,
406  6.0,
407  7.0,
408  DirectionChoiceEnum::BOTH,
409  "H2O",
410  0.0,
411  -1.0,
412  0.0);
413  FAIL() << "Missing expected exception.";
414  }
415  catch (const std::exception & e)
416  {
417  std::string msg(e.what());
418  ASSERT_TRUE(msg.find("Promoting species H2O has already been provided with an exponent") !=
419  std::string::npos)
420  << "Failed with unexpected error message: " << msg;
421  }
422 }
423 
425 TEST(GeochemistryKineticRateCalculatorTest, exceptions)
426 {
427  model_kin.addKineticRate(rate_ch4);
428  const std::vector<Real> promoting_indices(5 + 6, 0.0);
429  const std::vector<std::string> basis_name = {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++"};
430  const std::vector<bool> & basis_species_gas = mgd_kin.basis_species_gas;
431  const std::vector<Real> basis_molality = {1.0, 2.0, 3.0, 4.0, 5.0};
432  const std::vector<Real> basis_activity = {0.1, 0.2, 0.3, 0.4, 0.5};
433  const std::vector<bool> basis_activity_known(5, false);
434  const std::vector<std::string> & eqm_name = mgd_kin.eqm_species_name;
435  const std::vector<bool> & eqm_species_gas = mgd_kin.eqm_species_gas;
436  const std::vector<Real> eqm_molality = {1.5, 2.5, 3.5, 4.5, 5.5, 7.5};
437  const std::vector<Real> eqm_activity = {1.1, 1.2, 1.3, 1.4, 1.5, 1.7};
438  const DenseMatrix<Real> & eqm_stoichiometry = mgd_kin.eqm_stoichiometry;
439  const DenseMatrix<Real> & kin_stoichiometry = mgd_kin.kin_stoichiometry;
440 
441  Real rate;
442  Real drate_dkin;
443  std::vector<Real> drate_dmol(5);
444 
445  try
446  {
447  GeochemistryKineticRateCalculator::calculateRate({},
448  std::vector<Real>(5 + 6, 0),
449  std::vector<Real>(5 + 6, 0),
450  rate_ch4,
451  basis_name,
452  basis_species_gas,
453  basis_molality,
454  basis_activity,
455  basis_activity_known,
456  eqm_name,
457  eqm_species_gas,
458  eqm_molality,
459  eqm_activity,
460  eqm_stoichiometry,
461  1.0,
462  2.0,
463  3.0,
464  4.0,
465  kin_stoichiometry,
466  0,
467  25.0,
468  rate,
469  drate_dkin,
470  drate_dmol);
471  FAIL() << "Missing expected exception.";
472  }
473  catch (const std::exception & e)
474  {
475  std::string msg(e.what());
476  ASSERT_TRUE(msg.find("kinetic_rate: promoting_indices incorrectly sized 0") !=
477  std::string::npos)
478  << "Failed with unexpected error message: " << msg;
479  }
480 
481  try
482  {
483  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
484  {},
485  std::vector<Real>(5 + 6, 0),
486  rate_ch4,
487  basis_name,
488  basis_species_gas,
489  basis_molality,
490  basis_activity,
491  basis_activity_known,
492  eqm_name,
493  eqm_species_gas,
494  eqm_molality,
495  eqm_activity,
496  eqm_stoichiometry,
497  1.0,
498  2.0,
499  3.0,
500  4.0,
501  kin_stoichiometry,
502  0,
503  25.0,
504  rate,
505  drate_dkin,
506  drate_dmol);
507  FAIL() << "Missing expected exception.";
508  }
509  catch (const std::exception & e)
510  {
511  std::string msg(e.what());
512  ASSERT_TRUE(msg.find("kinetic_rate: promoting_monod_indices incorrectly sized 0") !=
513  std::string::npos)
514  << "Failed with unexpected error message: " << msg;
515  }
516 
517  try
518  {
519  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
520  std::vector<Real>(5 + 6, 0),
521  {},
522  rate_ch4,
523  basis_name,
524  basis_species_gas,
525  basis_molality,
526  basis_activity,
527  basis_activity_known,
528  eqm_name,
529  eqm_species_gas,
530  eqm_molality,
531  eqm_activity,
532  eqm_stoichiometry,
533  1.0,
534  2.0,
535  3.0,
536  4.0,
537  kin_stoichiometry,
538  0,
539  25.0,
540  rate,
541  drate_dkin,
542  drate_dmol);
543  FAIL() << "Missing expected exception.";
544  }
545  catch (const std::exception & e)
546  {
547  std::string msg(e.what());
548  ASSERT_TRUE(msg.find("kinetic_rate: promoting_half_saturation incorrectly sized 0") !=
549  std::string::npos)
550  << "Failed with unexpected error message: " << msg;
551  }
552 
553  try
554  {
555  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
556  std::vector<Real>(5 + 6, 0),
557  std::vector<Real>(5 + 6, 0),
558  rate_ch4,
559  basis_name,
560  {},
561  basis_molality,
562  basis_activity,
563  basis_activity_known,
564  eqm_name,
565  eqm_species_gas,
566  eqm_molality,
567  eqm_activity,
568  eqm_stoichiometry,
569  1.0,
570  2.0,
571  3.0,
572  4.0,
573  kin_stoichiometry,
574  0,
575  25.0,
576  rate,
577  drate_dkin,
578  drate_dmol);
579  FAIL() << "Missing expected exception.";
580  }
581  catch (const std::exception & e)
582  {
583  std::string msg(e.what());
584  ASSERT_TRUE(msg.find("kinetic_rate: incorrectly sized basis-species vectors 5 0 5 5 5 5") !=
585  std::string::npos)
586  << "Failed with unexpected error message: " << msg;
587  }
588 
589  try
590  {
591  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
592  std::vector<Real>(5 + 6, 0),
593  std::vector<Real>(5 + 6, 0),
594  rate_ch4,
595  basis_name,
596  basis_species_gas,
597  {},
598  basis_activity,
599  basis_activity_known,
600  eqm_name,
601  eqm_species_gas,
602  eqm_molality,
603  eqm_activity,
604  eqm_stoichiometry,
605  1.0,
606  2.0,
607  3.0,
608  4.0,
609  kin_stoichiometry,
610  0,
611  25.0,
612  rate,
613  drate_dkin,
614  drate_dmol);
615  FAIL() << "Missing expected exception.";
616  }
617  catch (const std::exception & e)
618  {
619  std::string msg(e.what());
620  ASSERT_TRUE(msg.find("kinetic_rate: incorrectly sized basis-species vectors 5 5 0 5 5 5") !=
621  std::string::npos)
622  << "Failed with unexpected error message: " << msg;
623  }
624 
625  try
626  {
627  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
628  std::vector<Real>(5 + 6, 0),
629  std::vector<Real>(5 + 6, 0),
630  rate_ch4,
631  basis_name,
632  basis_species_gas,
633  basis_molality,
634  {},
635  basis_activity_known,
636  eqm_name,
637  eqm_species_gas,
638  eqm_molality,
639  eqm_activity,
640  eqm_stoichiometry,
641  1.0,
642  2.0,
643  3.0,
644  4.0,
645  kin_stoichiometry,
646  0,
647  25.0,
648  rate,
649  drate_dkin,
650  drate_dmol);
651  FAIL() << "Missing expected exception.";
652  }
653  catch (const std::exception & e)
654  {
655  std::string msg(e.what());
656  ASSERT_TRUE(msg.find("kinetic_rate: incorrectly sized basis-species vectors 5 5 5 0 5 5") !=
657  std::string::npos)
658  << "Failed with unexpected error message: " << msg;
659  }
660 
661  try
662  {
663  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
664  std::vector<Real>(5 + 6, 0),
665  std::vector<Real>(5 + 6, 0),
666  rate_ch4,
667  basis_name,
668  basis_species_gas,
669  basis_molality,
670  basis_activity,
671  {},
672  eqm_name,
673  eqm_species_gas,
674  eqm_molality,
675  eqm_activity,
676  eqm_stoichiometry,
677  1.0,
678  2.0,
679  3.0,
680  4.0,
681  kin_stoichiometry,
682  0,
683  25.0,
684  rate,
685  drate_dkin,
686  drate_dmol);
687  FAIL() << "Missing expected exception.";
688  }
689  catch (const std::exception & e)
690  {
691  std::string msg(e.what());
692  ASSERT_TRUE(msg.find("kinetic_rate: incorrectly sized basis-species vectors 5 5 5 5 0 5") !=
693  std::string::npos)
694  << "Failed with unexpected error message: " << msg;
695  }
696 
697  try
698  {
699  std::vector<Real> drdm(4);
700  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
701  std::vector<Real>(5 + 6, 0),
702  std::vector<Real>(5 + 6, 0),
703  rate_ch4,
704  basis_name,
705  basis_species_gas,
706  basis_molality,
707  basis_activity,
708  basis_activity_known,
709  eqm_name,
710  eqm_species_gas,
711  eqm_molality,
712  eqm_activity,
713  eqm_stoichiometry,
714  1.0,
715  2.0,
716  3.0,
717  4.0,
718  kin_stoichiometry,
719  0,
720  25.0,
721  rate,
722  drate_dkin,
723  drdm);
724  FAIL() << "Missing expected exception.";
725  }
726  catch (const std::exception & e)
727  {
728  std::string msg(e.what());
729  ASSERT_TRUE(msg.find("kinetic_rate: incorrectly sized basis-species vectors 5 5 5 5 5 4") !=
730  std::string::npos)
731  << "Failed with unexpected error message: " << msg;
732  }
733 
734  try
735  {
736  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
737  std::vector<Real>(5 + 6, 0),
738  std::vector<Real>(5 + 6, 0),
739  rate_ch4,
740  basis_name,
741  basis_species_gas,
742  basis_molality,
743  basis_activity,
744  basis_activity_known,
745  eqm_name,
746  {},
747  eqm_molality,
748  eqm_activity,
749  eqm_stoichiometry,
750  1.0,
751  2.0,
752  3.0,
753  4.0,
754  kin_stoichiometry,
755  0,
756  25.0,
757  rate,
758  drate_dkin,
759  drate_dmol);
760  FAIL() << "Missing expected exception.";
761  }
762  catch (const std::exception & e)
763  {
764  std::string msg(e.what());
765  ASSERT_TRUE(msg.find("kinetic_rate: incorrectly sized equilibrium-species vectors 6 0 6 6") !=
766  std::string::npos)
767  << "Failed with unexpected error message: " << msg;
768  }
769 
770  try
771  {
772  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
773  std::vector<Real>(5 + 6, 0),
774  std::vector<Real>(5 + 6, 0),
775  rate_ch4,
776  basis_name,
777  basis_species_gas,
778  basis_molality,
779  basis_activity,
780  basis_activity_known,
781  eqm_name,
782  eqm_species_gas,
783  {},
784  eqm_activity,
785  eqm_stoichiometry,
786  1.0,
787  2.0,
788  3.0,
789  4.0,
790  kin_stoichiometry,
791  0,
792  25.0,
793  rate,
794  drate_dkin,
795  drate_dmol);
796  FAIL() << "Missing expected exception.";
797  }
798  catch (const std::exception & e)
799  {
800  std::string msg(e.what());
801  ASSERT_TRUE(msg.find("kinetic_rate: incorrectly sized equilibrium-species vectors 6 6 0 6") !=
802  std::string::npos)
803  << "Failed with unexpected error message: " << msg;
804  }
805 
806  try
807  {
808  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
809  std::vector<Real>(5 + 6, 0),
810  std::vector<Real>(5 + 6, 0),
811  rate_ch4,
812  basis_name,
813  basis_species_gas,
814  basis_molality,
815  basis_activity,
816  basis_activity_known,
817  eqm_name,
818  eqm_species_gas,
819  eqm_molality,
820  {},
821  eqm_stoichiometry,
822  1.0,
823  2.0,
824  3.0,
825  4.0,
826  kin_stoichiometry,
827  0,
828  25.0,
829  rate,
830  drate_dkin,
831  drate_dmol);
832  FAIL() << "Missing expected exception.";
833  }
834  catch (const std::exception & e)
835  {
836  std::string msg(e.what());
837  ASSERT_TRUE(msg.find("kinetic_rate: incorrectly sized equilibrium-species vectors 6 6 6 0") !=
838  std::string::npos)
839  << "Failed with unexpected error message: " << msg;
840  }
841 
842  try
843  {
844  DenseMatrix<Real> stoi(7, 4);
845  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
846  std::vector<Real>(5 + 6, 0),
847  std::vector<Real>(5 + 6, 0),
848  rate_ch4,
849  basis_name,
850  basis_species_gas,
851  basis_molality,
852  basis_activity,
853  basis_activity_known,
854  eqm_name,
855  eqm_species_gas,
856  eqm_molality,
857  eqm_activity,
858  stoi,
859  1.0,
860  2.0,
861  3.0,
862  4.0,
863  kin_stoichiometry,
864  0,
865  25.0,
866  rate,
867  drate_dkin,
868  drate_dmol);
869  FAIL() << "Missing expected exception.";
870  }
871  catch (const std::exception & e)
872  {
873  std::string msg(e.what());
874  ASSERT_TRUE(msg.find("kinetic_rate: incorrectly sized eqm stoichiometry matrix 7 4") !=
875  std::string::npos)
876  << "Failed with unexpected error message: " << msg;
877  }
878 
879  try
880  {
881  DenseMatrix<Real> stoi(4, 5);
882  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
883  std::vector<Real>(5 + 6, 0),
884  std::vector<Real>(5 + 6, 0),
885  rate_ch4,
886  basis_name,
887  basis_species_gas,
888  basis_molality,
889  basis_activity,
890  basis_activity_known,
891  eqm_name,
892  eqm_species_gas,
893  eqm_molality,
894  eqm_activity,
895  stoi,
896  1.0,
897  2.0,
898  3.0,
899  4.0,
900  kin_stoichiometry,
901  0,
902  25.0,
903  rate,
904  drate_dkin,
905  drate_dmol);
906  FAIL() << "Missing expected exception.";
907  }
908  catch (const std::exception & e)
909  {
910  std::string msg(e.what());
911  ASSERT_TRUE(msg.find("kinetic_rate: incorrectly sized eqm stoichiometry matrix 4 5") !=
912  std::string::npos)
913  << "Failed with unexpected error message: " << msg;
914  }
915 
916  try
917  {
918  DenseMatrix<Real> stoi(1, 5);
919  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
920  std::vector<Real>(5 + 6, 0),
921  std::vector<Real>(5 + 6, 0),
922  rate_ch4,
923  basis_name,
924  basis_species_gas,
925  basis_molality,
926  basis_activity,
927  basis_activity_known,
928  eqm_name,
929  eqm_species_gas,
930  eqm_molality,
931  eqm_activity,
932  eqm_stoichiometry,
933  1.0,
934  2.0,
935  3.0,
936  4.0,
937  stoi,
938  1,
939  25.0,
940  rate,
941  drate_dkin,
942  drate_dmol);
943  FAIL() << "Missing expected exception.";
944  }
945  catch (const std::exception & e)
946  {
947  std::string msg(e.what());
948  ASSERT_TRUE(msg.find("kinetic_rate: incorrectly sized kinetic stoichiometry matrix 1 5") !=
949  std::string::npos)
950  << "Failed with unexpected error message: " << msg;
951  }
952 
953  try
954  {
955  DenseMatrix<Real> stoi(1, 4);
956  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
957  std::vector<Real>(5 + 6, 0),
958  std::vector<Real>(5 + 6, 0),
959  rate_ch4,
960  basis_name,
961  basis_species_gas,
962  basis_molality,
963  basis_activity,
964  basis_activity_known,
965  eqm_name,
966  eqm_species_gas,
967  eqm_molality,
968  eqm_activity,
969  eqm_stoichiometry,
970  1.0,
971  2.0,
972  3.0,
973  4.0,
974  stoi,
975  0,
976  25.0,
977  rate,
978  drate_dkin,
979  drate_dmol);
980  FAIL() << "Missing expected exception.";
981  }
982  catch (const std::exception & e)
983  {
984  std::string msg(e.what());
985  ASSERT_TRUE(msg.find("kinetic_rate: incorrectly sized kinetic stoichiometry matrix 1 4") !=
986  std::string::npos)
987  << "Failed with unexpected error message: " << msg;
988  }
989 }
990 
992 TEST(GeochemistryKineticRateCalculatorTest, rate1)
993 {
994  model_kin.addKineticRate(rate_ch4);
995  model_kin.addKineticRate(rate_ch4_dissolution);
996  model_kin.addKineticRate(rate_ch4_precipitation);
997  model_kin.addKineticRate(rate_ch4_raw);
998  model_kin.addKineticRate(rate_ch4_death);
999  model_kin.addKineticRate(rate_cal);
1000  std::vector<Real> promoting_indices(5 + 7, 0.0);
1001  const std::vector<std::string> basis_name = {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++"};
1002  const std::vector<bool> basis_species_gas = {false, false, false, true, false};
1003  std::vector<Real> basis_molality = {1.5, 2.0, 3.0, 4.0, 5.0};
1004  const std::vector<Real> basis_activity = {0.1, 0.2, 0.3, 0.4, 0.5};
1005  const std::vector<bool> basis_activity_known = {false, false, false, true, true};
1006  const std::vector<std::string> eqm_name = {"OH-", "n2", "n3", "n4", "n5", "n6", "n7"};
1007  const std::vector<bool> eqm_species_gas = {false, false, false, true, false, false, false};
1008  const std::vector<Real> eqm_molality = {1.5, 2.5, 2.3, 2.1, 1.9, 1.7, 1.3};
1009  const std::vector<Real> eqm_activity = {1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7};
1010  DenseMatrix<Real> eqm_stoichiometry(7, 5);
1011  DenseMatrix<Real> kin_stoichiometry(2, 5);
1012 
1013  Real rate;
1014  Real drate_dkin;
1015  std::vector<Real> drate_dmol(5, 1.0);
1016 
1017  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
1018  std::vector<Real>(12, 0),
1019  std::vector<Real>(12, 0),
1020  rate_ch4,
1021  basis_name,
1022  basis_species_gas,
1023  basis_molality,
1024  basis_activity,
1025  basis_activity_known,
1026  eqm_name,
1027  eqm_species_gas,
1028  eqm_molality,
1029  eqm_activity,
1030  eqm_stoichiometry,
1031  1.25,
1032  2.25,
1033  3.5,
1034  4.0,
1035  kin_stoichiometry,
1036  0,
1037  50.0,
1038  rate,
1039  drate_dkin,
1040  drate_dmol);
1041  EXPECT_NEAR(rate,
1042  -1.5 * 2.0 * 1.25 * 2.25 *
1043  std::pow(std::abs(1 - std::pow(std::pow(10.0, 4.0 - 3.5), 0.8)), 2.5) *
1044  std::exp(66.0 / GeochemistryConstants::GAS_CONSTANT *
1045  (0.003 - 1.0 / (50.0 + GeochemistryConstants::CELSIUS_TO_KELVIN))),
1046  1.0E-6);
1047  EXPECT_NEAR(drate_dkin,
1048  -1.5 * 2.0 * 2.25 *
1049  std::pow(std::abs(1 - std::pow(std::pow(10.0, 4.0 - 3.5), 0.8)), 2.5) *
1050  std::exp(66.0 / GeochemistryConstants::GAS_CONSTANT *
1051  (0.003 - 1.0 / (50.0 + GeochemistryConstants::CELSIUS_TO_KELVIN))),
1052  1.0E-6);
1053  for (unsigned i = 0; i < 5; ++i)
1054  ASSERT_EQ(drate_dmol[i], 0.0);
1055 
1056  // Since rate_ch4_dissolution does dissolution only, so when activity_product > eqm_constant, rate
1057  // = 0:
1058  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
1059  std::vector<Real>(12, 0),
1060  std::vector<Real>(12, 0),
1061  rate_ch4_dissolution,
1062  basis_name,
1063  basis_species_gas,
1064  basis_molality,
1065  basis_activity,
1066  basis_activity_known,
1067  eqm_name,
1068  eqm_species_gas,
1069  eqm_molality,
1070  eqm_activity,
1071  eqm_stoichiometry,
1072  1.25,
1073  2.25,
1074  3.5,
1075  4.0,
1076  kin_stoichiometry,
1077  0,
1078  50.0,
1079  rate,
1080  drate_dkin,
1081  drate_dmol);
1082  EXPECT_EQ(rate, 0.0);
1083  EXPECT_EQ(drate_dkin, 0.0);
1084  for (unsigned i = 0; i < 5; ++i)
1085  ASSERT_EQ(drate_dmol[i], 0.0);
1086 
1087  // Since rate_ch4_dissolution has some dissolution when activity_product < eqm_constant
1088  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
1089  std::vector<Real>(12, 0),
1090  std::vector<Real>(12, 0),
1091  rate_ch4_dissolution,
1092  basis_name,
1093  basis_species_gas,
1094  basis_molality,
1095  basis_activity,
1096  basis_activity_known,
1097  eqm_name,
1098  eqm_species_gas,
1099  eqm_molality,
1100  eqm_activity,
1101  eqm_stoichiometry,
1102  1.25,
1103  2.25,
1104  4.0,
1105  3.5,
1106  kin_stoichiometry,
1107  0,
1108  50.0,
1109  rate,
1110  drate_dkin,
1111  drate_dmol);
1112  EXPECT_NEAR(rate,
1113  1.5 * 2.0 * 1.25 * 2.25 *
1114  std::pow(std::abs(1 - std::pow(std::pow(10.0, 3.5 - 4.0), 0.8)), 2.5) *
1115  std::exp(66.0 / GeochemistryConstants::GAS_CONSTANT *
1116  (0.003 - 1.0 / (50.0 + GeochemistryConstants::CELSIUS_TO_KELVIN))),
1117  1.0E-6);
1118  EXPECT_NEAR(drate_dkin,
1119  1.5 * 2.0 * 2.25 *
1120  std::pow(std::abs(1 - std::pow(std::pow(10.0, 3.5 - 4.0), 0.8)), 2.5) *
1121  std::exp(66.0 / GeochemistryConstants::GAS_CONSTANT *
1122  (0.003 - 1.0 / (50.0 + GeochemistryConstants::CELSIUS_TO_KELVIN))),
1123  1.0E-6);
1124  for (unsigned i = 0; i < 5; ++i)
1125  ASSERT_EQ(drate_dmol[i], 0.0);
1126 
1127  // Since rate_ch4_precipitation does precipitation only, and since activity_product >
1128  // eqm_constant, rate != 0:
1129  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
1130  std::vector<Real>(12, 0),
1131  std::vector<Real>(12, 0),
1132  rate_ch4_precipitation,
1133  basis_name,
1134  basis_species_gas,
1135  basis_molality,
1136  basis_activity,
1137  basis_activity_known,
1138  eqm_name,
1139  eqm_species_gas,
1140  eqm_molality,
1141  eqm_activity,
1142  eqm_stoichiometry,
1143  1.25,
1144  2.25,
1145  3.5,
1146  4.0,
1147  kin_stoichiometry,
1148  0,
1149  50.0,
1150  rate,
1151  drate_dkin,
1152  drate_dmol);
1153  EXPECT_NEAR(rate,
1154  -1.5 * 2.0 * 1.25 * 2.25 *
1155  std::pow(std::abs(1 - std::pow(std::pow(10.0, 4.0 - 3.5), 0.8)), 2.5) *
1156  std::exp(66.0 / GeochemistryConstants::GAS_CONSTANT *
1157  (0.003 - 1.0 / (50.0 + GeochemistryConstants::CELSIUS_TO_KELVIN))),
1158  1.0E-6);
1159  EXPECT_NEAR(drate_dkin,
1160  -1.5 * 2.0 * 2.25 *
1161  std::pow(std::abs(1 - std::pow(std::pow(10.0, 4.0 - 3.5), 0.8)), 2.5) *
1162  std::exp(66.0 / GeochemistryConstants::GAS_CONSTANT *
1163  (0.003 - 1.0 / (50.0 + GeochemistryConstants::CELSIUS_TO_KELVIN))),
1164  1.0E-6);
1165  for (unsigned i = 0; i < 5; ++i)
1166  ASSERT_EQ(drate_dmol[i], 0.0);
1167 
1168  // rate_ch4_precipitation does precipitation only, but because activity_product < eqm_constant,
1169  // rate = 0
1170  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
1171  std::vector<Real>(12, 0),
1172  std::vector<Real>(12, 0),
1173  rate_ch4_precipitation,
1174  basis_name,
1175  basis_species_gas,
1176  basis_molality,
1177  basis_activity,
1178  basis_activity_known,
1179  eqm_name,
1180  eqm_species_gas,
1181  eqm_molality,
1182  eqm_activity,
1183  eqm_stoichiometry,
1184  1.25,
1185  2.25,
1186  4.0,
1187  3.5,
1188  kin_stoichiometry,
1189  0,
1190  50.0,
1191  rate,
1192  drate_dkin,
1193  drate_dmol);
1194  EXPECT_EQ(rate, 0.0);
1195  EXPECT_EQ(drate_dkin, 0.0);
1196  for (unsigned i = 0; i < 5; ++i)
1197  ASSERT_EQ(drate_dmol[i], 0.0);
1198 
1199  // rate_ch4_raw does not depend on activity_product / eqm_constant
1200  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
1201  std::vector<Real>(12, 0),
1202  std::vector<Real>(12, 0),
1203  rate_ch4_raw,
1204  basis_name,
1205  basis_species_gas,
1206  basis_molality,
1207  basis_activity,
1208  basis_activity_known,
1209  eqm_name,
1210  eqm_species_gas,
1211  eqm_molality,
1212  eqm_activity,
1213  eqm_stoichiometry,
1214  1.25,
1215  2.25,
1216  3.5,
1217  4.0,
1218  kin_stoichiometry,
1219  0,
1220  50.0,
1221  rate,
1222  drate_dkin,
1223  drate_dmol);
1224  EXPECT_NEAR(rate,
1225  1.5 * 2.0 * 1.25 * 2.25 *
1226  std::pow(std::abs(1 - std::pow(std::pow(10.0, 4.0 - 3.5), 0.8)), 2.5) *
1227  std::exp(66.0 / GeochemistryConstants::GAS_CONSTANT *
1228  (0.003 - 1.0 / (50.0 + GeochemistryConstants::CELSIUS_TO_KELVIN))),
1229  1.0E-6);
1230  EXPECT_NEAR(drate_dkin,
1231  1.5 * 2.0 * 2.25 *
1232  std::pow(std::abs(1 - std::pow(std::pow(10.0, 4.0 - 3.5), 0.8)), 2.5) *
1233  std::exp(66.0 / GeochemistryConstants::GAS_CONSTANT *
1234  (0.003 - 1.0 / (50.0 + GeochemistryConstants::CELSIUS_TO_KELVIN))),
1235  1.0E-6);
1236  for (unsigned i = 0; i < 5; ++i)
1237  ASSERT_EQ(drate_dmol[i], 0.0);
1238 
1239  // rate_ch4_death does not depend on activity_product / eqm_constant
1240  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
1241  std::vector<Real>(12, 0),
1242  std::vector<Real>(12, 0),
1243  rate_ch4_death,
1244  basis_name,
1245  basis_species_gas,
1246  basis_molality,
1247  basis_activity,
1248  basis_activity_known,
1249  eqm_name,
1250  eqm_species_gas,
1251  eqm_molality,
1252  eqm_activity,
1253  eqm_stoichiometry,
1254  1.25,
1255  2.25,
1256  3.5,
1257  4.0,
1258  kin_stoichiometry,
1259  0,
1260  50.0,
1261  rate,
1262  drate_dkin,
1263  drate_dmol);
1264  EXPECT_NEAR(rate,
1265  1.5 * 2.0 * 1.25 * 2.25 *
1266  std::pow(std::abs(1 - std::pow(std::pow(10.0, 4.0 - 3.5), 0.8)), 2.5) *
1267  std::exp(66.0 / GeochemistryConstants::GAS_CONSTANT *
1268  (0.003 - 1.0 / (50.0 + GeochemistryConstants::CELSIUS_TO_KELVIN))),
1269  1.0E-6);
1270  EXPECT_NEAR(drate_dkin,
1271  1.5 * 2.0 * 2.25 *
1272  std::pow(std::abs(1 - std::pow(std::pow(10.0, 4.0 - 3.5), 0.8)), 2.5) *
1273  std::exp(66.0 / GeochemistryConstants::GAS_CONSTANT *
1274  (0.003 - 1.0 / (50.0 + GeochemistryConstants::CELSIUS_TO_KELVIN))),
1275  1.0E-6);
1276  for (unsigned i = 0; i < 5; ++i)
1277  ASSERT_EQ(drate_dmol[i], 0.0);
1278 
1279  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
1280  std::vector<Real>(12, 0),
1281  std::vector<Real>(12, 0),
1282  rate_ch4,
1283  basis_name,
1284  basis_species_gas,
1285  basis_molality,
1286  basis_activity,
1287  basis_activity_known,
1288  eqm_name,
1289  eqm_species_gas,
1290  eqm_molality,
1291  eqm_activity,
1292  eqm_stoichiometry,
1293  1.25,
1294  2.25,
1295  4.5,
1296  4.0,
1297  kin_stoichiometry,
1298  0,
1299  50.0,
1300  rate,
1301  drate_dkin,
1302  drate_dmol);
1303  EXPECT_NEAR(rate,
1304  1.5 * 2.0 * 1.25 * 2.25 *
1305  std::pow(std::abs(1 - std::pow(std::pow(10.0, 4.0 - 4.5), 0.8)), 2.5) *
1306  std::exp(66.0 / GeochemistryConstants::GAS_CONSTANT *
1307  (0.003 - 1.0 / (50.0 + GeochemistryConstants::CELSIUS_TO_KELVIN))),
1308  1.0E-6);
1309  EXPECT_NEAR(drate_dkin,
1310  1.5 * 2.0 * 2.25 *
1311  std::pow(std::abs(1 - std::pow(std::pow(10.0, 4.0 - 4.5), 0.8)), 2.5) *
1312  std::exp(66.0 / GeochemistryConstants::GAS_CONSTANT *
1313  (0.003 - 1.0 / (50.0 + GeochemistryConstants::CELSIUS_TO_KELVIN))),
1314  1.0E-6);
1315  for (unsigned i = 0; i < 5; ++i)
1316  ASSERT_EQ(drate_dmol[i], 0.0);
1317 
1318  for (unsigned i = 0; i < 5; ++i)
1319  promoting_indices[i] = 0.5 * (i + 1);
1320  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
1321  std::vector<Real>(12, 0),
1322  std::vector<Real>(12, 0),
1323  rate_cal,
1324  basis_name,
1325  basis_species_gas,
1326  basis_molality,
1327  basis_activity,
1328  basis_activity_known,
1329  eqm_name,
1330  eqm_species_gas,
1331  eqm_molality,
1332  eqm_activity,
1333  eqm_stoichiometry,
1334  1.25,
1335  2.25,
1336  3.5,
1337  4.0,
1338  kin_stoichiometry,
1339  0,
1340  40.0,
1341  rate,
1342  drate_dkin,
1343  drate_dmol);
1344  EXPECT_NEAR(rate,
1345  -7.0 * 6.0 * std::pow(1.5, 0.5) * std::pow(0.2, 1.0) * std::pow(3.0, 1.5) *
1346  std::pow(0.4, 2.0) * std::pow(5.0, 2.5) *
1347  std::pow(std::abs(1 - std::pow(std::pow(10.0, 4.0 - 3.5), 2.5)), 0.8) *
1348  std::exp(55.0 / GeochemistryConstants::GAS_CONSTANT *
1349  (0.00315 - 1.0 / (40.0 + GeochemistryConstants::CELSIUS_TO_KELVIN))),
1350  1.0E-6);
1351  EXPECT_EQ(drate_dkin, 0.0);
1352  EXPECT_NEAR(drate_dmol[0], 0.5 * rate / 1.5, 1.0E-6);
1353  EXPECT_NEAR(drate_dmol[1], 1.0 * rate / 2.0, 1.0E-6);
1354  EXPECT_NEAR(drate_dmol[2], 1.5 * rate / 3.0, 1.0E-6);
1355  EXPECT_EQ(drate_dmol[3], 0.0);
1356  EXPECT_NEAR(drate_dmol[4], 2.5 * rate / 5.0, 1.0E-6);
1357 
1358  for (unsigned i = 5; i < 12; ++i)
1359  promoting_indices[i] = -0.1 * (i + 1);
1360  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
1361  std::vector<Real>(12, 0),
1362  std::vector<Real>(12, 0),
1363  rate_cal,
1364  basis_name,
1365  basis_species_gas,
1366  basis_molality,
1367  basis_activity,
1368  basis_activity_known,
1369  eqm_name,
1370  eqm_species_gas,
1371  eqm_molality,
1372  eqm_activity,
1373  eqm_stoichiometry,
1374  1.25,
1375  2.25,
1376  3.5,
1377  4.0,
1378  kin_stoichiometry,
1379  0,
1380  40.0,
1381  rate,
1382  drate_dkin,
1383  drate_dmol);
1384  EXPECT_NEAR(rate,
1385  -7.0 * 6.0 * std::pow(1.5, 0.5) * std::pow(0.2, 1.0) * std::pow(3.0, 1.5) *
1386  std::pow(0.4, 2.0) * std::pow(5.0, 2.5) * std::pow(1.1, -0.6) *
1387  std::pow(2.5, -0.7) * std::pow(2.3, -0.8) * std::pow(1.4, -0.9) *
1388  std::pow(1.9, -1.0) * std::pow(1.7, -1.1) * std::pow(1.3, -1.2) *
1389  std::pow(std::abs(1 - std::pow(std::pow(10.0, 4.0 - 3.5), 2.5)), 0.8) *
1390  std::exp(55.0 / GeochemistryConstants::GAS_CONSTANT *
1391  (0.00315 - 1.0 / (40.0 + GeochemistryConstants::CELSIUS_TO_KELVIN))),
1392  1.0E-6);
1393  EXPECT_EQ(drate_dkin, 0.0);
1394  EXPECT_NEAR(drate_dmol[0], 0.5 * rate / 1.5, 1.0E-6);
1395  EXPECT_NEAR(drate_dmol[1], 1.0 * rate / 2.0, 1.0E-6);
1396  EXPECT_NEAR(drate_dmol[2], 1.5 * rate / 3.0, 1.0E-6);
1397  EXPECT_EQ(drate_dmol[3], 0.0);
1398  EXPECT_NEAR(drate_dmol[4], 2.5 * rate / 5.0, 1.0E-6);
1399 
1400  // derivatives of activity coefficients are zero
1401  eqm_stoichiometry(0, 0) = 1.0; // OH- depends on H2O
1402  eqm_stoichiometry(0, 1) = -1.25; // OH- depends on H+
1403  eqm_stoichiometry(0, 2) = 1.5; // OH- depends on another thing
1404  eqm_stoichiometry(0, 3) = 3.5; // OH- depends on the gas
1405  eqm_stoichiometry(3, 0) = 2.0; // gas depends on H2O
1406  eqm_stoichiometry(3, 1) = -2.25; // gas depends on H+
1407  eqm_stoichiometry(3, 2) = -2.5; // gas depends on another thing
1408  eqm_stoichiometry(3, 3) = -33.5; // gas depends on the gas
1409  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
1410  std::vector<Real>(12, 0),
1411  std::vector<Real>(12, 0),
1412  rate_cal,
1413  basis_name,
1414  basis_species_gas,
1415  basis_molality,
1416  basis_activity,
1417  basis_activity_known,
1418  eqm_name,
1419  eqm_species_gas,
1420  eqm_molality,
1421  eqm_activity,
1422  eqm_stoichiometry,
1423  1.25,
1424  2.25,
1425  3.5,
1426  4.0,
1427  kin_stoichiometry,
1428  0,
1429  40.0,
1430  rate,
1431  drate_dkin,
1432  drate_dmol);
1433  EXPECT_NEAR(rate,
1434  -7.0 * 6.0 * std::pow(1.5, 0.5) * std::pow(0.2, 1.0) * std::pow(3.0, 1.5) *
1435  std::pow(0.4, 2.0) * std::pow(5.0, 2.5) * std::pow(1.1, -0.6) *
1436  std::pow(2.5, -0.7) * std::pow(2.3, -0.8) * std::pow(1.4, -0.9) *
1437  std::pow(1.9, -1.0) * std::pow(1.7, -1.1) * std::pow(1.3, -1.2) *
1438  std::pow(std::abs(1 - std::pow(std::pow(10.0, 4.0 - 3.5), 2.5)), 0.8) *
1439  std::exp(55.0 / GeochemistryConstants::GAS_CONSTANT *
1440  (0.00315 - 1.0 / (40.0 + GeochemistryConstants::CELSIUS_TO_KELVIN))),
1441  1.0E-6);
1442  EXPECT_EQ(drate_dkin, 0.0);
1443  EXPECT_NEAR(drate_dmol[0], 0.5 * rate / 1.5, 1.0E-6);
1444  EXPECT_NEAR(drate_dmol[1],
1445  1.0 * rate / 2.0 - 0.6 * (-1.25) * rate / 2.0 - 0.9 * (-2.25) * rate / 2.0,
1446  1.0E-6);
1447  EXPECT_NEAR(drate_dmol[2],
1448  1.5 * rate / 3.0 - 0.6 * (1.5) * rate / 3.0 - 0.9 * (-2.5) * rate / 3.0,
1449  1.0E-6);
1450  EXPECT_EQ(drate_dmol[3], 0.0);
1451  EXPECT_NEAR(drate_dmol[4], 2.5 * rate / 5.0, 1.0E-6);
1452 
1453  // nonzero contributions from kinetic stoichiometry via the activity product
1454  eqm_stoichiometry.zero();
1455  kin_stoichiometry(0, 0) = -1.0;
1456  kin_stoichiometry(0, 1) = 1.25;
1457  kin_stoichiometry(0, 2) = -1.5;
1458  kin_stoichiometry(0, 3) = 1.75;
1459  kin_stoichiometry(1, 1) = -11.0;
1460  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
1461  std::vector<Real>(12, 0),
1462  std::vector<Real>(12, 0),
1463  rate_cal,
1464  basis_name,
1465  basis_species_gas,
1466  basis_molality,
1467  basis_activity,
1468  basis_activity_known,
1469  eqm_name,
1470  eqm_species_gas,
1471  eqm_molality,
1472  eqm_activity,
1473  eqm_stoichiometry,
1474  1.25,
1475  2.25,
1476  3.5,
1477  4.0,
1478  kin_stoichiometry,
1479  0,
1480  40.0,
1481  rate,
1482  drate_dkin,
1483  drate_dmol);
1484  const Real theta_term = std::pow(std::pow(10.0, 4.0 - 3.5), 2.5);
1485  EXPECT_NEAR(rate,
1486  -7.0 * 6.0 * std::pow(1.5, 0.5) * std::pow(0.2, 1.0) * std::pow(3.0, 1.5) *
1487  std::pow(0.4, 2.0) * std::pow(5.0, 2.5) * std::pow(1.1, -0.6) *
1488  std::pow(2.5, -0.7) * std::pow(2.3, -0.8) * std::pow(1.4, -0.9) *
1489  std::pow(1.9, -1.0) * std::pow(1.7, -1.1) * std::pow(1.3, -1.2) *
1490  std::pow(std::abs(1 - theta_term), 0.8) *
1491  std::exp(55.0 / GeochemistryConstants::GAS_CONSTANT *
1492  (0.00315 - 1.0 / (40.0 + GeochemistryConstants::CELSIUS_TO_KELVIN))),
1493  1.0E-6);
1494  EXPECT_EQ(drate_dkin, 0.0);
1495  EXPECT_NEAR(drate_dmol[0], 0.5 * rate / 1.5, 1.0E-6);
1496  EXPECT_NEAR(drate_dmol[1],
1497  1.0 * rate / 2.0 +
1498  0.8 * rate / std::abs(1 - theta_term) * (-1) * (-2.5) * theta_term * 1.25 / 2.0,
1499  1.0E-6);
1500  EXPECT_NEAR(drate_dmol[2],
1501  1.5 * rate / 3.0 +
1502  0.8 * rate / std::abs(1 - theta_term) * (-1) * (-2.5) * theta_term * (-1.5 / 3.0),
1503  1.0E-6);
1504  EXPECT_EQ(drate_dmol[3], 0.0);
1505  EXPECT_NEAR(drate_dmol[4], 2.5 * rate / 5.0, 1.0E-6);
1506 
1507  // do a finite-difference of the previous test.
1508  // Assume activity coefficients = 1
1509  Real ap = std::pow(0.1, -1.0) * std::pow(basis_activity[1], 1.25) *
1510  std::pow(basis_molality[2], -1.5) * std::pow(0.4, 1.75);
1511  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
1512  std::vector<Real>(12, 0),
1513  std::vector<Real>(12, 0),
1514  rate_cal,
1515  basis_name,
1516  basis_species_gas,
1517  basis_molality,
1518  basis_activity,
1519  basis_activity_known,
1520  eqm_name,
1521  eqm_species_gas,
1522  eqm_molality,
1523  eqm_activity,
1524  eqm_stoichiometry,
1525  1.25,
1526  2.25,
1527  -1.5,
1528  std::log10(ap),
1529  kin_stoichiometry,
1530  0,
1531  40.0,
1532  rate,
1533  drate_dkin,
1534  drate_dmol);
1535  const Real eps = 1.0E-8;
1536  basis_molality[2] += eps;
1537  ap = std::pow(0.1, -1.0) * std::pow(basis_activity[1], 1.25) * std::pow(basis_molality[2], -1.5) *
1538  std::pow(0.4, 1.75);
1539  Real rate_new;
1540  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
1541  std::vector<Real>(12, 0),
1542  std::vector<Real>(12, 0),
1543  rate_cal,
1544  basis_name,
1545  basis_species_gas,
1546  basis_molality,
1547  basis_activity,
1548  basis_activity_known,
1549  eqm_name,
1550  eqm_species_gas,
1551  eqm_molality,
1552  eqm_activity,
1553  eqm_stoichiometry,
1554  1.25,
1555  2.25,
1556  -1.5,
1557  std::log10(ap),
1558  kin_stoichiometry,
1559  0,
1560  40.0,
1561  rate_new,
1562  drate_dkin,
1563  drate_dmol);
1564  EXPECT_NEAR(drate_dmol[2], (rate_new - rate) / eps, 1E-5);
1565 }
1566 
1568 TEST(GeochemistryKineticRateCalculatorTest, rate_monod)
1569 {
1570  model_kin.addKineticRate(rate_ch4);
1571  model_kin.addKineticRate(rate_ch4_kin_mon);
1572  model_kin.addKineticRate(rate_cal);
1573  std::vector<Real> promoting_indices(5 + 7, 0.0);
1574  std::vector<Real> promoting_monod(5 + 7, 0.0);
1575  std::vector<Real> half_saturation(5 + 7, 0.0);
1576  const std::vector<std::string> basis_name = {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++"};
1577  const std::vector<bool> basis_species_gas = {false, false, false, true, false};
1578  const std::vector<Real> basis_molality = {1.5, 2.0, 3.0, 4.0, 5.0};
1579  const std::vector<Real> basis_activity = {0.1, 0.2, 0.3, 0.4, 0.5};
1580  const std::vector<bool> basis_activity_known = {false, false, false, true, true};
1581  const std::vector<std::string> eqm_name = {"OH-", "n2", "n3", "n4", "n5", "n6", "n7"};
1582  const std::vector<bool> eqm_species_gas = {false, false, false, true, false, false, false};
1583  std::vector<Real> eqm_molality = {1.5, 2.5, 2.3, 2.1, 1.9, 1.7, 1.3};
1584  std::vector<Real> eqm_activity = {1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7};
1585  DenseMatrix<Real> eqm_stoichiometry(7, 5);
1586  DenseMatrix<Real> kin_stoichiometry(2, 5);
1587 
1588  Real rate;
1589  Real drate_dkin;
1590  std::vector<Real> drate_dmol(5, 1.0);
1591 
1592  // rate with kinetic molal and monod terms active
1593  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
1594  std::vector<Real>(12, 0),
1595  std::vector<Real>(12, 0),
1596  rate_ch4_kin_mon,
1597  basis_name,
1598  basis_species_gas,
1599  basis_molality,
1600  basis_activity,
1601  basis_activity_known,
1602  eqm_name,
1603  eqm_species_gas,
1604  eqm_molality,
1605  eqm_activity,
1606  eqm_stoichiometry,
1607  1.25,
1608  2.25,
1609  3.5,
1610  4.0,
1611  kin_stoichiometry,
1612  0,
1613  50.0,
1614  rate,
1615  drate_dkin,
1616  drate_dmol);
1617  EXPECT_NEAR(rate,
1618  -1.5 * 2.0 * 1.25 * 2.25 * std::pow(1.25 / 1.5, 1.75) /
1619  std::pow(std::pow(1.25 / 1.5, 1.75) + std::pow(0.875, 1.75), 2.125) *
1620  std::pow(std::abs(1 - std::pow(std::pow(10.0, 4.0 - 3.5), 0.8)), 2.5) *
1621  std::exp(66.0 / GeochemistryConstants::GAS_CONSTANT *
1622  (0.003 - 1.0 / (50.0 + GeochemistryConstants::CELSIUS_TO_KELVIN))),
1623  1.0E-6);
1624  EXPECT_NEAR(drate_dkin,
1625  -1.5 * 2.0 * 2.25 *
1626  (std::pow(1.25 / 1.5, 1.75) /
1627  std::pow(std::pow(1.25 / 1.5, 1.75) + std::pow(0.875, 1.75), 2.125) +
1628  1.25 * 1.75 * std::pow(1.25 / 1.5, 1.75 - 1.0) / 1.5 /
1629  std::pow(std::pow(1.25 / 1.5, 1.75) + std::pow(0.875, 1.75), 2.125) +
1630  1.25 * std::pow(1.25 / 1.5, 1.75) * (-2.125) /
1631  std::pow(std::pow(1.25 / 1.5, 1.75) + std::pow(0.875, 1.75), 2.125 + 1.0) *
1632  1.75 * std::pow(1.25 / 1.5, 1.75 - 1.0) / 1.5) *
1633  std::pow(std::abs(1 - std::pow(std::pow(10.0, 4.0 - 3.5), 0.8)), 2.5) *
1634  std::exp(66.0 / GeochemistryConstants::GAS_CONSTANT *
1635  (0.003 - 1.0 / (50.0 + GeochemistryConstants::CELSIUS_TO_KELVIN))),
1636  1.0E-6);
1637  EXPECT_NEAR(drate_dmol[0],
1638  -1.5 * 2.0 * 1.25 * 2.25 *
1639  (1.75 * std::pow(1.25 / 1.5, 1.75 - 1.0) * (-1.25 / 1.5 / 1.5) /
1640  std::pow(std::pow(1.25 / 1.5, 1.75) + std::pow(0.875, 1.75), 2.125) +
1641  std::pow(1.25 / 1.5, 1.75) * (-2.125) /
1642  std::pow(std::pow(1.25 / 1.5, 1.75) + std::pow(0.875, 1.75), 2.125 + 1.0) *
1643  1.75 * std::pow(1.25 / 1.5, 1.75 - 1.0) * (-1.25 / 1.5 / 1.5)) *
1644  std::pow(std::abs(1 - std::pow(std::pow(10.0, 4.0 - 3.5), 0.8)), 2.5) *
1645  std::exp(66.0 / GeochemistryConstants::GAS_CONSTANT *
1646  (0.003 - 1.0 / (50.0 + GeochemistryConstants::CELSIUS_TO_KELVIN))),
1647  1.0E-6);
1648  for (unsigned i = 1; i < 5; ++i)
1649  ASSERT_EQ(drate_dmol[i], 0.0);
1650 
1651  // with basis promoting indices and monod terms
1652  for (unsigned i = 0; i < 5; ++i)
1653  {
1654  promoting_indices[i] = 0.5 * (i + 1);
1655  promoting_monod[i] = 0.25 * i;
1656  half_saturation[i] = 0.125 * i;
1657  }
1658  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
1659  promoting_monod,
1660  half_saturation,
1661  rate_cal,
1662  basis_name,
1663  basis_species_gas,
1664  basis_molality,
1665  basis_activity,
1666  basis_activity_known,
1667  eqm_name,
1668  eqm_species_gas,
1669  eqm_molality,
1670  eqm_activity,
1671  eqm_stoichiometry,
1672  1.25,
1673  2.25,
1674  3.5,
1675  4.0,
1676  kin_stoichiometry,
1677  0,
1678  40.0,
1679  rate,
1680  drate_dkin,
1681  drate_dmol);
1682  EXPECT_NEAR(rate,
1683  -7.0 * 6.0 *
1684  (std::pow(1.5, 0.5) / std::pow(std::pow(1.5, 0.5) + std::pow(0.0, 0.5), 0.0)) *
1685  (std::pow(0.2, 1.0) / std::pow(std::pow(0.2, 1.0) + std::pow(0.125, 1.0), 0.25)) *
1686  (std::pow(3.0, 1.5) / std::pow(std::pow(3.0, 1.5) + std::pow(0.25, 1.5), 0.5)) *
1687  (std::pow(0.4, 2.0) / std::pow(std::pow(0.4, 2.0) + std::pow(0.375, 2.0), 0.75)) *
1688  (std::pow(5.0, 2.5) / std::pow(std::pow(5.0, 2.5) + std::pow(0.5, 2.5), 1.0)) *
1689  std::pow(std::abs(1 - std::pow(std::pow(10.0, 4.0 - 3.5), 2.5)), 0.8) *
1690  std::exp(55.0 / GeochemistryConstants::GAS_CONSTANT *
1691  (0.00315 - 1.0 / (40.0 + GeochemistryConstants::CELSIUS_TO_KELVIN))),
1692  1.0E-6);
1693  EXPECT_EQ(drate_dkin, 0.0);
1694  EXPECT_NEAR(drate_dmol[0], 0.5 * rate / 1.5, 1.0E-6);
1695  EXPECT_NEAR(drate_dmol[1],
1696  1.0 * rate / 2.0 + (-0.25) / (std::pow(0.2, 1.0) + std::pow(0.125, 1.0)) * rate *
1697  1.0 * std::pow(0.2, 1.0 - 1.0) * 0.2 / 2.0,
1698  1.0E-6);
1699  EXPECT_NEAR(drate_dmol[2],
1700  1.5 * rate / 3.0 + (-0.5) / (std::pow(3.0, 1.5) + std::pow(0.25, 1.5)) * rate * 1.5 *
1701  std::pow(3.0, 1.5 - 1.0),
1702  1.0E-6);
1703  EXPECT_EQ(drate_dmol[3], 0.0);
1704  EXPECT_NEAR(drate_dmol[4],
1705  2.5 * rate / 5.0 + (-1.0) / (std::pow(5.0, 2.5) + std::pow(0.5, 2.5)) * rate * 2.5 *
1706  std::pow(5.0, 2.5 - 1.0),
1707  1.0E-6);
1708 
1709  // with basis and equilibrium promoting indices and monod terms, but with eqm_stoichiometry = 0
1710  for (unsigned i = 5; i < 12; ++i)
1711  {
1712  promoting_indices[i] = -0.1 * (i + 1);
1713  promoting_monod[i] = 0.25 * i;
1714  half_saturation[i] = 0.125 * i;
1715  }
1716  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
1717  promoting_monod,
1718  half_saturation,
1719  rate_cal,
1720  basis_name,
1721  basis_species_gas,
1722  basis_molality,
1723  basis_activity,
1724  basis_activity_known,
1725  eqm_name,
1726  eqm_species_gas,
1727  eqm_molality,
1728  eqm_activity,
1729  eqm_stoichiometry,
1730  1.25,
1731  2.25,
1732  3.5,
1733  4.0,
1734  kin_stoichiometry,
1735  0,
1736  40.0,
1737  rate,
1738  drate_dkin,
1739  drate_dmol);
1740  EXPECT_NEAR(
1741  rate,
1742  -7.0 * 6.0 * (std::pow(1.5, 0.5) / std::pow(std::pow(1.5, 0.5) + std::pow(0.0, 0.5), 0.0)) *
1743  (std::pow(0.2, 1.0) / std::pow(std::pow(0.2, 1.0) + std::pow(0.125, 1.0), 0.25)) *
1744  (std::pow(3.0, 1.5) / std::pow(std::pow(3.0, 1.5) + std::pow(0.25, 1.5), 0.5)) *
1745  (std::pow(0.4, 2.0) / std::pow(std::pow(0.4, 2.0) + std::pow(0.375, 2.0), 0.75)) *
1746  (std::pow(5.0, 2.5) / std::pow(std::pow(5.0, 2.5) + std::pow(0.5, 2.5), 1.0)) *
1747  (std::pow(1.1, -0.6) / std::pow(std::pow(1.1, -0.6) + std::pow(0.625, -0.6), 1.25)) *
1748  (std::pow(2.5, -0.7) / std::pow(std::pow(2.5, -0.7) + std::pow(0.75, -0.7), 1.5)) *
1749  (std::pow(2.3, -0.8) / std::pow(std::pow(2.3, -0.8) + std::pow(0.875, -0.8), 1.75)) *
1750  (std::pow(1.4, -0.9) / std::pow(std::pow(1.4, -0.9) + std::pow(1.0, -0.9), 2.0)) *
1751  (std::pow(1.9, -1.0) / std::pow(std::pow(1.9, -1.0) + std::pow(1.125, -1.0), 2.25)) *
1752  (std::pow(1.7, -1.1) / std::pow(std::pow(1.7, -1.1) + std::pow(1.25, -1.1), 2.5)) *
1753  (std::pow(1.3, -1.2) / std::pow(std::pow(1.3, -1.2) + std::pow(1.375, -1.2), 2.75)) *
1754  std::pow(std::abs(1 - std::pow(std::pow(10.0, 4.0 - 3.5), 2.5)), 0.8) *
1755  std::exp(55.0 / GeochemistryConstants::GAS_CONSTANT *
1756  (0.00315 - 1.0 / (40.0 + GeochemistryConstants::CELSIUS_TO_KELVIN))),
1757  1.0E-6);
1758  EXPECT_EQ(drate_dkin, 0.0);
1759  EXPECT_NEAR(drate_dmol[0], 0.5 * rate / 1.5, 1.0E-6);
1760  EXPECT_NEAR(drate_dmol[1],
1761  1.0 * rate / 2.0 + (-0.25) / (std::pow(0.2, 1.0) + std::pow(0.125, 1.0)) * rate *
1762  1.0 * std::pow(0.2, 1.0 - 1.0) * 0.2 / 2.0,
1763  1.0E-6);
1764  EXPECT_NEAR(drate_dmol[2],
1765  1.5 * rate / 3.0 + (-0.5) / (std::pow(3.0, 1.5) + std::pow(0.25, 1.5)) * rate * 1.5 *
1766  std::pow(3.0, 1.5 - 1.0),
1767  1.0E-6);
1768  EXPECT_EQ(drate_dmol[3], 0.0);
1769  EXPECT_NEAR(drate_dmol[4],
1770  2.5 * rate / 5.0 + (-1.0) / (std::pow(5.0, 2.5) + std::pow(0.5, 2.5)) * rate * 2.5 *
1771  std::pow(5.0, 2.5 - 1.0),
1772  1.0E-6);
1773 
1774  // do finite-difference version of the previous test
1775  // note that since the derivatives of activity coefficients are zero, all activity coefficients
1776  // are set to 1 in the following
1777  for (unsigned i = 0; i < 5; ++i)
1778  {
1779  Real rate_new;
1780  const Real eps = 1.0E-8;
1781  if (basis_species_gas[i])
1782  continue; // molality is undefined
1783 
1784  // original
1785  std::vector<Real> fd_basis_molality = basis_molality;
1786  std::vector<Real> fd_basis_activity =
1787  fd_basis_molality; // because activity product is assumed to be constant
1788  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
1789  promoting_monod,
1790  half_saturation,
1791  rate_cal,
1792  basis_name,
1793  basis_species_gas,
1794  fd_basis_molality,
1795  fd_basis_activity,
1796  basis_activity_known,
1797  eqm_name,
1798  eqm_species_gas,
1799  eqm_molality,
1800  eqm_activity,
1801  eqm_stoichiometry,
1802  1.25,
1803  2.25,
1804  3.5,
1805  4.0,
1806  kin_stoichiometry,
1807  0,
1808  40.0,
1809  rate,
1810  drate_dkin,
1811  drate_dmol);
1812  // original + eps
1813  fd_basis_molality[i] += eps;
1814  fd_basis_activity = fd_basis_molality; // because activity product is assumed to be constant
1815  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
1816  promoting_monod,
1817  half_saturation,
1818  rate_cal,
1819  basis_name,
1820  basis_species_gas,
1821  fd_basis_molality,
1822  fd_basis_activity,
1823  basis_activity_known,
1824  eqm_name,
1825  eqm_species_gas,
1826  eqm_molality,
1827  eqm_activity,
1828  eqm_stoichiometry,
1829  1.25,
1830  2.25,
1831  3.5,
1832  4.0,
1833  kin_stoichiometry,
1834  0,
1835  40.0,
1836  rate_new,
1837  drate_dkin,
1838  drate_dmol);
1839  const Real fd = (rate_new - rate) / eps;
1840  EXPECT_TRUE(std::abs(drate_dmol[i] / fd - 1.0) < 1E-4);
1841  }
1842 
1843  // with basis and equilibrium promoting indices and monod terms, and with eqm_stoichiometry != 0
1844  eqm_stoichiometry(0, 0) = 1.0; // OH- depends on H2O
1845  eqm_stoichiometry(0, 1) = -1.25; // OH- depends on H+
1846  eqm_stoichiometry(0, 2) = 1.5; // OH- depends on HCO3-
1847  eqm_stoichiometry(0, 3) = 3.5; // OH- depends on the gas
1848  eqm_stoichiometry(1, 0) = 1.25; // n2 depends on H2O
1849  eqm_stoichiometry(1, 1) = 2.75; // n2 depends on H+
1850  eqm_stoichiometry(1, 2) = 0.75; // n2 depends on HCO3-
1851  eqm_stoichiometry(1, 3) = -1.5; // n2 depends on the gas
1852  eqm_stoichiometry(3, 0) = 2.0; // gas depends on H2O
1853  eqm_stoichiometry(3, 1) = -2.25; // gas depends on H+
1854  eqm_stoichiometry(3, 2) = -2.5; // gas depends on HCO3-
1855  eqm_stoichiometry(3, 3) = -33.5; // gas depends on the gas
1856  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
1857  promoting_monod,
1858  half_saturation,
1859  rate_cal,
1860  basis_name,
1861  basis_species_gas,
1862  basis_molality,
1863  basis_activity,
1864  basis_activity_known,
1865  eqm_name,
1866  eqm_species_gas,
1867  eqm_molality,
1868  eqm_activity,
1869  eqm_stoichiometry,
1870  1.25,
1871  2.25,
1872  3.5,
1873  4.0,
1874  kin_stoichiometry,
1875  0,
1876  40.0,
1877  rate,
1878  drate_dkin,
1879  drate_dmol);
1880  EXPECT_NEAR(
1881  rate,
1882  -7.0 * 6.0 * (std::pow(1.5, 0.5) / std::pow(std::pow(1.5, 0.5) + std::pow(0.0, 0.5), 0.0)) *
1883  (std::pow(0.2, 1.0) / std::pow(std::pow(0.2, 1.0) + std::pow(0.125, 1.0), 0.25)) *
1884  (std::pow(3.0, 1.5) / std::pow(std::pow(3.0, 1.5) + std::pow(0.25, 1.5), 0.5)) *
1885  (std::pow(0.4, 2.0) / std::pow(std::pow(0.4, 2.0) + std::pow(0.375, 2.0), 0.75)) *
1886  (std::pow(5.0, 2.5) / std::pow(std::pow(5.0, 2.5) + std::pow(0.5, 2.5), 1.0)) *
1887  (std::pow(1.1, -0.6) / std::pow(std::pow(1.1, -0.6) + std::pow(0.625, -0.6), 1.25)) *
1888  (std::pow(2.5, -0.7) / std::pow(std::pow(2.5, -0.7) + std::pow(0.75, -0.7), 1.5)) *
1889  (std::pow(2.3, -0.8) / std::pow(std::pow(2.3, -0.8) + std::pow(0.875, -0.8), 1.75)) *
1890  (std::pow(1.4, -0.9) / std::pow(std::pow(1.4, -0.9) + std::pow(1.0, -0.9), 2.0)) *
1891  (std::pow(1.9, -1.0) / std::pow(std::pow(1.9, -1.0) + std::pow(1.125, -1.0), 2.25)) *
1892  (std::pow(1.7, -1.1) / std::pow(std::pow(1.7, -1.1) + std::pow(1.25, -1.1), 2.5)) *
1893  (std::pow(1.3, -1.2) / std::pow(std::pow(1.3, -1.2) + std::pow(1.375, -1.2), 2.75)) *
1894  std::pow(std::abs(1 - std::pow(std::pow(10.0, 4.0 - 3.5), 2.5)), 0.8) *
1895  std::exp(55.0 / GeochemistryConstants::GAS_CONSTANT *
1896  (0.00315 - 1.0 / (40.0 + GeochemistryConstants::CELSIUS_TO_KELVIN))),
1897  1.0E-6);
1898  EXPECT_EQ(drate_dkin, 0.0);
1899  EXPECT_NEAR(drate_dmol[0], 0.5 * rate / 1.5, 1.0E-6); // derivative of water activity is ignored
1900  EXPECT_NEAR(
1901  drate_dmol[1],
1902  (1.0 * rate / 2.0 + (-0.25) / (std::pow(0.2, 1.0) + std::pow(0.125, 1.0)) * rate * 1.0 *
1903  std::pow(0.2, 1.0 - 1.0) * 0.2 / 2.0) +
1904  (-1.25 *
1905  (-0.6 * rate / 2.0 + (-1.25) / (std::pow(1.1, -0.6) + std::pow(0.625, -0.6)) * rate *
1906  (-0.6) * std::pow(1.1, -0.6 - 1.0) * 1.1 / 2.0)) +
1907  (2.75 * (-0.7 * rate / 2.0 + (-1.5) / (std::pow(2.5, -0.7) + std::pow(0.75, -0.7)) *
1908  rate * (-0.7) * std::pow(2.5, -0.7 - 1.0) * 2.5 / 2.0)) +
1909  (-2.25 * (-0.9 * rate / 2.0 + (-2.0) / (std::pow(1.4, -0.9) + std::pow(1.0, -0.9)) *
1910  rate * (-0.9) * std::pow(1.4, -0.9 - 1.0) * 1.4 / 2.0)),
1911  1.0E-6);
1912  EXPECT_NEAR(
1913  drate_dmol[2],
1914  (1.5 * rate / 3.0 + (-0.5) / (std::pow(3.0, 1.5) + std::pow(0.25, 1.5)) * rate * 1.5 *
1915  std::pow(3.0, 1.5 - 1.0)) +
1916  (1.5 *
1917  (-0.6 * rate / 1.1 * 1.1 / 1.5 +
1918  (-1.25) / (std::pow(1.1, -0.6) + std::pow(0.625, -0.6)) * rate * (-0.6) *
1919  std::pow(1.1, -0.6 - 1.0) * 1.1 / 1.5) *
1920  1.5 / 3.0) +
1921  (0.75 * (-0.7 * rate / 3.0 + (-1.5) / (std::pow(2.5, -0.7) + std::pow(0.75, -0.7)) *
1922  rate * (-0.7) * std::pow(2.5, -0.7 - 1.0) * 2.5 / 3.0)) +
1923  (-2.5 *
1924  (-0.9 * rate / 1.4 + (-2.0) / (std::pow(1.4, -0.9) + std::pow(1.0, -0.9)) * rate *
1925  (-0.9) * std::pow(1.4, -0.9 - 1.0)) *
1926  1.4 / 3.0),
1927  1.0E-6);
1928 
1929  EXPECT_EQ(drate_dmol[3], 0.0);
1930  EXPECT_NEAR(drate_dmol[4],
1931  2.5 * rate / 5.0 + (-1.0) / (std::pow(5.0, 2.5) + std::pow(0.5, 2.5)) * rate * 2.5 *
1932  std::pow(5.0, 2.5 - 1.0),
1933  1.0E-6);
1934 
1935  // do finite-difference version of the previous test
1936  // note that since the derivatives of activity coefficients are zero, all activity coefficients
1937  // are set to 1 in the following
1938  for (unsigned i = 0; i < 5; ++i)
1939  {
1940  Real rate_new;
1941  const Real eps = 1.0E-8;
1942  if (basis_species_gas[i])
1943  continue; // molality is undefined
1944 
1945  // original
1946  std::vector<Real> fd_basis_molality = basis_molality;
1947  std::vector<Real> fd_basis_activity =
1948  fd_basis_molality; // because activity product is assumed to be constant
1949  for (unsigned j = 0; j < 7; ++j)
1950  {
1951  eqm_molality[j] = std::pow(fd_basis_activity[0], eqm_stoichiometry(j, 0));
1952  for (unsigned basis = 1; basis < 5; ++basis)
1953  {
1954  if (basis_species_gas[basis])
1955  continue;
1956  eqm_molality[j] *= std::pow(fd_basis_molality[basis], eqm_stoichiometry(j, basis));
1957  }
1958  }
1959  eqm_activity = eqm_molality;
1960  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
1961  promoting_monod,
1962  half_saturation,
1963  rate_cal,
1964  basis_name,
1965  basis_species_gas,
1966  fd_basis_molality,
1967  fd_basis_activity,
1968  basis_activity_known,
1969  eqm_name,
1970  eqm_species_gas,
1971  eqm_molality,
1972  eqm_activity,
1973  eqm_stoichiometry,
1974  1.25,
1975  2.25,
1976  3.5,
1977  4.0,
1978  kin_stoichiometry,
1979  0,
1980  40.0,
1981  rate,
1982  drate_dkin,
1983  drate_dmol);
1984  // original + eps
1985  fd_basis_molality[i] += eps;
1986  for (unsigned basis = 1; basis < 5; ++basis)
1987  fd_basis_activity[basis] =
1988  fd_basis_molality[basis]; // because activity product is assumed to be constant, but
1989  // derivs wrt water activity are not included
1990  for (unsigned j = 0; j < 7; ++j)
1991  {
1992  eqm_molality[j] = std::pow(fd_basis_activity[0], eqm_stoichiometry(j, 0));
1993  for (unsigned basis = 1; basis < 5; ++basis)
1994  {
1995  if (basis_species_gas[basis])
1996  continue;
1997  eqm_molality[j] *= std::pow(fd_basis_molality[basis], eqm_stoichiometry(j, basis));
1998  }
1999  }
2000  eqm_activity = eqm_molality;
2001  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
2002  promoting_monod,
2003  half_saturation,
2004  rate_cal,
2005  basis_name,
2006  basis_species_gas,
2007  fd_basis_molality,
2008  fd_basis_activity,
2009  basis_activity_known,
2010  eqm_name,
2011  eqm_species_gas,
2012  eqm_molality,
2013  eqm_activity,
2014  eqm_stoichiometry,
2015  1.25,
2016  2.25,
2017  3.5,
2018  4.0,
2019  kin_stoichiometry,
2020  0,
2021  40.0,
2022  rate_new,
2023  drate_dkin,
2024  drate_dmol);
2025  const Real fd = (rate_new - rate) / eps;
2026  EXPECT_TRUE(std::abs(drate_dmol[i] / fd - 1.0) < 1E-4);
2027  }
2028 }
2029 
2031 TEST(GeochemistryKineticRateCalculatorTest, rate_energy_captured)
2032 {
2033  model_kin.addKineticRate(rate_ch4_energy_captured);
2034  std::vector<Real> promoting_indices(5 + 7, 0.0);
2035  const std::vector<std::string> basis_name = {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++"};
2036  const std::vector<bool> basis_species_gas = {false, false, false, true, false};
2037  const std::vector<Real> basis_molality = {1.5, 2.0, 3.0, 4.0, 5.0};
2038  const std::vector<Real> basis_activity = {0.1, 0.2, 0.3, 0.4, 0.5};
2039  const std::vector<bool> basis_activity_known = {false, false, false, true, true};
2040  const std::vector<std::string> eqm_name = {"OH-", "n2", "n3", "n4", "n5", "n6", "n7"};
2041  const std::vector<bool> eqm_species_gas = {false, false, false, true, false, false, false};
2042  const std::vector<Real> eqm_molality = {1.5, 2.5, 2.3, 2.1, 1.9, 1.7, 1.3};
2043  const std::vector<Real> eqm_activity = {1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7};
2044  DenseMatrix<Real> eqm_stoichiometry(7, 5);
2045  DenseMatrix<Real> kin_stoichiometry(2, 5);
2046 
2047  Real rate;
2048  Real drate_dkin;
2049  std::vector<Real> drate_dmol(5, 1.0);
2050 
2051  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
2052  std::vector<Real>(12, 0),
2053  std::vector<Real>(12, 0),
2054  rate_ch4_energy_captured,
2055  basis_name,
2056  basis_species_gas,
2057  basis_molality,
2058  basis_activity,
2059  basis_activity_known,
2060  eqm_name,
2061  eqm_species_gas,
2062  eqm_molality,
2063  eqm_activity,
2064  eqm_stoichiometry,
2065  1.25,
2066  2.25,
2067  3.5,
2068  4.0,
2069  kin_stoichiometry,
2070  0,
2071  50.0,
2072  rate,
2073  drate_dkin,
2074  drate_dmol);
2075  EXPECT_NEAR(
2076  rate,
2077  -1.5 * 2.0 * 1.25 * 2.25 *
2078  std::pow(
2079  std::abs(1 - std::pow(std::pow(10.0,
2080  4.0 - (3.5 -
2081  1.0 / GeochemistryConstants::GAS_CONSTANT /
2082  (50.0 +
2083  GeochemistryConstants::CELSIUS_TO_KELVIN) /
2084  GeochemistryConstants::LOGTEN)),
2085  0.8)),
2086  2.5) *
2087  std::exp(66.0 / GeochemistryConstants::GAS_CONSTANT *
2088  (0.003 - 1.0 / (50.0 + GeochemistryConstants::CELSIUS_TO_KELVIN))),
2089  1.0E-6);
2090  EXPECT_NEAR(
2091  drate_dkin,
2092  -1.5 * 2.0 * 2.25 *
2093  std::pow(
2094  std::abs(1 - std::pow(std::pow(10.0,
2095  4.0 - (3.5 -
2096  1.0 / GeochemistryConstants::GAS_CONSTANT /
2097  (50.0 +
2098  GeochemistryConstants::CELSIUS_TO_KELVIN) /
2099  GeochemistryConstants::LOGTEN)),
2100  0.8)),
2101  2.5) *
2102  std::exp(66.0 / GeochemistryConstants::GAS_CONSTANT *
2103  (0.003 - 1.0 / (50.0 + GeochemistryConstants::CELSIUS_TO_KELVIN))),
2104  1.0E-6);
2105  for (unsigned i = 0; i < 5; ++i)
2106  ASSERT_EQ(drate_dmol[i], 0.0);
2107 }
2108 
2110 TEST(GeochemistryKineticRateCalculatorTest, rate_eta_theta)
2111 {
2112  model_kin.addKineticRate(rate_cal_theta0_eta2);
2113  model_kin.addKineticRate(rate_cal_theta0_eta1);
2114  model_kin.addKineticRate(rate_cal_theta2_eta1);
2115  model_kin.addKineticRate(rate_cal_theta0_eta0);
2116  model_kin.addKineticRate(rate_cal_theta0_eta05);
2117  std::vector<Real> promoting_indices(5 + 7, 0.0);
2118  const std::vector<std::string> basis_name = {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++"};
2119  const std::vector<bool> basis_species_gas = {false, false, false, true, false};
2120  const std::vector<Real> basis_molality = {1.5, 2.0, 3.0, 4.0, 5.0};
2121  const std::vector<Real> basis_activity = {0.1, 0.2, 0.3, 0.4, 0.5};
2122  const std::vector<bool> basis_activity_known = {false, false, false, true, true};
2123  const std::vector<std::string> eqm_name = {"OH-", "n2", "n3", "n4", "n5", "n6", "n7"};
2124  const std::vector<bool> eqm_species_gas = {false, false, false, true, false, false, false};
2125  const std::vector<Real> eqm_molality = {1.5, 2.5, 2.3, 2.1, 1.9, 1.7, 1.3};
2126  const std::vector<Real> eqm_activity = {1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7};
2127  DenseMatrix<Real> eqm_stoichiometry(7, 5);
2128  DenseMatrix<Real> kin_stoichiometry(2, 5);
2129 
2130  Real rate;
2131  Real drate_dkin;
2132  std::vector<Real> drate_dmol(5, 1.0);
2133 
2134  for (unsigned i = 0; i < 5; ++i)
2135  promoting_indices[i] = 0.5 * (i + 1);
2136 
2137  kin_stoichiometry(0, 0) = -1.0;
2138  kin_stoichiometry(0, 1) = 1.25;
2139  kin_stoichiometry(0, 2) = -1.5;
2140  kin_stoichiometry(0, 3) = 1.75;
2141  kin_stoichiometry(1, 1) = -11.0;
2142  kin_stoichiometry(1, 2) = 0.5;
2143  kin_stoichiometry(1, 3) = 4.0;
2144 
2145  // theta = 0 and eta = 2
2146  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
2147  std::vector<Real>(12, 0),
2148  std::vector<Real>(12, 0),
2149  rate_cal_theta0_eta2,
2150  basis_name,
2151  basis_species_gas,
2152  basis_molality,
2153  basis_activity,
2154  basis_activity_known,
2155  eqm_name,
2156  eqm_species_gas,
2157  eqm_molality,
2158  eqm_activity,
2159  eqm_stoichiometry,
2160  1.25,
2161  2.25,
2162  3.5,
2163  4.0,
2164  kin_stoichiometry,
2165  1,
2166  40.0,
2167  rate,
2168  drate_dkin,
2169  drate_dmol);
2170  EXPECT_EQ(rate, 0.0);
2171  EXPECT_EQ(drate_dkin, 0.0);
2172  EXPECT_NEAR(drate_dmol[0], 0.5 * rate / 1.5, 1.0E-6);
2173  EXPECT_NEAR(drate_dmol[1], 1.0 * rate / 2.0 + 0.0, 1E-6);
2174  EXPECT_NEAR(drate_dmol[2], 1.5 * rate / 3.0 + 0.0, 1E-6);
2175  EXPECT_EQ(drate_dmol[3], 0.0);
2176  EXPECT_NEAR(drate_dmol[4], 2.5 * rate / 5.0, 1.0E-6);
2177 
2178  // theta = 0 and eta = 1
2179  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
2180  std::vector<Real>(12, 0),
2181  std::vector<Real>(12, 0),
2182  rate_cal_theta0_eta1,
2183  basis_name,
2184  basis_species_gas,
2185  basis_molality,
2186  basis_activity,
2187  basis_activity_known,
2188  eqm_name,
2189  eqm_species_gas,
2190  eqm_molality,
2191  eqm_activity,
2192  eqm_stoichiometry,
2193  1.25,
2194  2.25,
2195  3.5,
2196  4.0,
2197  kin_stoichiometry,
2198  1,
2199  40.0,
2200  rate,
2201  drate_dkin,
2202  drate_dmol);
2203  EXPECT_EQ(rate, 0.0);
2204  EXPECT_EQ(drate_dkin, 0.0);
2205  for (unsigned i = 0; i < 5; ++i)
2206  EXPECT_EQ(drate_dmol[i], 0.0);
2207 
2208  // theta = 2 and eta = 1 and Q = K
2209  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
2210  std::vector<Real>(12, 0),
2211  std::vector<Real>(12, 0),
2212  rate_cal_theta2_eta1,
2213  basis_name,
2214  basis_species_gas,
2215  basis_molality,
2216  basis_activity,
2217  basis_activity_known,
2218  eqm_name,
2219  eqm_species_gas,
2220  eqm_molality,
2221  eqm_activity,
2222  eqm_stoichiometry,
2223  1.25,
2224  2.25,
2225  4.0,
2226  4.0,
2227  kin_stoichiometry,
2228  1,
2229  40.0,
2230  rate,
2231  drate_dkin,
2232  drate_dmol);
2233  EXPECT_EQ(rate, 0.0);
2234  EXPECT_EQ(drate_dkin, 0.0);
2235  /*
2236  const Real rate_no_theta =
2237  7.0 * 6.0 * std::pow(1.5, 0.5) * std::pow(0.2, 1.0) * std::pow(3.0, 1.5) *
2238  std::pow(0.4, 2.0) * std::pow(5.0, 2.5) *
2239  std::exp(55.0 / GeochemistryConstants::GAS_CONSTANT *
2240  (0.00315 - 1.0 / (40.0 + GeochemistryConstants::CELSIUS_TO_KELVIN)));
2241  const Real deriv_ap_term = rate_no_theta * 2.0 * 1.0;
2242  */
2243  EXPECT_EQ(drate_dmol[0], 0.0);
2244  EXPECT_NEAR(drate_dmol[1], 0.0, 1E-6); // deriv_ap_term * (-11.0) / 2.0, 1E-6);
2245  EXPECT_NEAR(drate_dmol[2], 0.0, 1E-6); // deriv_ap_term * (0.5) / 3.0, 1E-6);
2246  EXPECT_EQ(drate_dmol[3], 0.0);
2247  EXPECT_EQ(drate_dmol[4], 0.0);
2248 
2249  // theta = 0 and eta = 0
2250  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
2251  std::vector<Real>(12, 0),
2252  std::vector<Real>(12, 0),
2253  rate_cal_theta0_eta0,
2254  basis_name,
2255  basis_species_gas,
2256  basis_molality,
2257  basis_activity,
2258  basis_activity_known,
2259  eqm_name,
2260  eqm_species_gas,
2261  eqm_molality,
2262  eqm_activity,
2263  eqm_stoichiometry,
2264  1.25,
2265  2.25,
2266  3.5,
2267  4.0,
2268  kin_stoichiometry,
2269  1,
2270  40.0,
2271  rate,
2272  drate_dkin,
2273  drate_dmol);
2274  EXPECT_NEAR(rate,
2275  -7.0 * 6.0 * std::pow(1.5, 0.5) * std::pow(0.2, 1.0) * std::pow(3.0, 1.5) *
2276  std::pow(0.4, 2.0) * std::pow(5.0, 2.5) *
2277  std::exp(55.0 / GeochemistryConstants::GAS_CONSTANT *
2278  (0.00315 - 1.0 / (40.0 + GeochemistryConstants::CELSIUS_TO_KELVIN))),
2279  1.0E-6);
2280  EXPECT_EQ(drate_dkin, 0.0);
2281  EXPECT_NEAR(drate_dmol[0], 0.5 * rate / 1.5, 1.0E-6);
2282  EXPECT_NEAR(drate_dmol[1], 1.0 * rate / 2.0 + 0.0, 1E-6);
2283  EXPECT_NEAR(drate_dmol[2], 1.5 * rate / 3.0 + 0.0, 1E-6);
2284  EXPECT_EQ(drate_dmol[3], 0.0);
2285  EXPECT_NEAR(drate_dmol[4], 2.5 * rate / 5.0, 1.0E-6);
2286 
2287  // theta = 0 and eta = 0.5
2288  GeochemistryKineticRateCalculator::calculateRate(promoting_indices,
2289  std::vector<Real>(12, 0),
2290  std::vector<Real>(12, 0),
2291  rate_cal_theta0_eta05,
2292  basis_name,
2293  basis_species_gas,
2294  basis_molality,
2295  basis_activity,
2296  basis_activity_known,
2297  eqm_name,
2298  eqm_species_gas,
2299  eqm_molality,
2300  eqm_activity,
2301  eqm_stoichiometry,
2302  1.25,
2303  2.25,
2304  3.5,
2305  4.0,
2306  kin_stoichiometry,
2307  1,
2308  40.0,
2309  rate,
2310  drate_dkin,
2311  drate_dmol);
2312  EXPECT_EQ(rate, 0.0);
2313  EXPECT_EQ(drate_dkin, 0.0);
2314  EXPECT_NEAR(drate_dmol[0], 0.5 * rate / 1.5, 1.0E-6);
2315  EXPECT_LE(drate_dmol[1], -std::numeric_limits<Real>::max());
2316  EXPECT_GE(drate_dmol[2], 0.1 * std::numeric_limits<Real>::max());
2317  EXPECT_EQ(drate_dmol[3], 0.0);
2318  EXPECT_NEAR(drate_dmol[4], 2.5 * rate / 5.0, 1.0E-6);
2319 }
GeochemistryConstants::CELSIUS_TO_KELVIN
constexpr Real CELSIUS_TO_KELVIN
Definition: GeochemistryConstants.h:22
eps
const Real eps
Definition: GeochemistrySpeciesSwapperTest.C:14
DenseMatrix< Real >::zero
virtual void zero() override final
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:199
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:330
rate_ch4_raw
KineticRateUserDescription rate_ch4_raw("CH4(aq)", 1.5, 2.0, true, 0.0, 0.0, 0.0, {"H2O", "OH-", "O2(aq)", "CO2(aq)", "CaCO3"}, {3.0, 3.1, 3.2, 3.3, 3.4}, {0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0}, 0.8, 2.5, 66.0, 0.003, DirectionChoiceEnum::RAW, "H2O", 0.0, -1.0, 0.0)
rate_cal_theta0_eta0
KineticRateUserDescription rate_cal_theta0_eta0("Calcite", 7.0, 6.0, false, 0.0, 0.0, 0.0, {"H+"}, {-3.0}, {0.0}, {0.0}, 0.0, 0.0, 55.0, 0.00315, DirectionChoiceEnum::BOTH, "H2O", 0.0, -1.0, 0.0)
GeochemistryKineticRateCalculator::calculateRate
void calculateRate(const std::vector< Real > &promoting_indices, const std::vector< Real > &promoting_monod_indices, const std::vector< Real > &promoting_half_saturation, const KineticRateUserDescription &description, const std::vector< std::string > &basis_species_name, const std::vector< bool > &basis_species_gas, const std::vector< Real > &basis_molality, const std::vector< Real > &basis_activity, const std::vector< bool > &basis_activity_known, const std::vector< std::string > &eqm_species_name, const std::vector< bool > &eqm_species_gas, const std::vector< Real > &eqm_molality, const std::vector< Real > &eqm_activity, const DenseMatrix< Real > &eqm_stoichiometry, Real kin_moles, Real kin_species_molecular_weight, Real log10K, Real log10_activity_product, const DenseMatrix< Real > &kin_stoichiometry, unsigned kin, Real temp_degC, Real &rate, Real &drate_dkin, std::vector< Real > &drate_dmol)
Calclates a kinetic rate and its derivative.
Definition: GeochemistryKineticRateCalculator.C:15
rate_ch4_kin_mon
KineticRateUserDescription rate_ch4_kin_mon("CH4(aq)", 1.5, 2.0, true, 1.75, 2.125, 0.875, {"H2O", "OH-", "O2(aq)", "CO2(aq)", "CaCO3"}, {3.0, 3.1, 3.2, 3.3, 3.4}, {0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0}, 0.8, 2.5, 66.0, 0.003, DirectionChoiceEnum::BOTH, "H2O", 0.0, -1.0, 0.0)
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:220
rate_ch4
KineticRateUserDescription rate_ch4("CH4(aq)", 1.5, 2.0, true, 0.0, 0.0, 0.0, {"H2O", "OH-", "O2(aq)", "CO2(aq)", "CaCO3"}, {3.0, 3.1, 3.2, 3.3, 3.4}, {0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0}, 0.8, 2.5, 66.0, 0.003, DirectionChoiceEnum::BOTH, "H2O", 0.0, -1.0, 0.0)
mgd_kin
const ModelGeochemicalDatabase & mgd_kin
Definition: GeochemistryKineticRateCalculatorTest.C:285
model_kin
PertinentGeochemicalSystem model_kin(db_kin, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++"}, {}, {}, {"Calcite"}, {"CH4(aq)"}, {}, "O2(aq)", "e-")
PertinentGeochemicalSystem::modelGeochemicalDatabase
const ModelGeochemicalDatabase & modelGeochemicalDatabase() const
Return a reference to the ModelGeochemicalDatabase structure.
Definition: PertinentGeochemicalSystem.C:850
PertinentGeochemicalSystem::addKineticRate
void addKineticRate(const KineticRateUserDescription &description)
Adds a rate description for kinetic_species.
Definition: PertinentGeochemicalSystem.C:856
PertinentGeochemicalSystem
Constructs and stores a minimal amount of information that is pertinent to the user-defined geochemic...
Definition: PertinentGeochemicalSystem.h:382
rate_cal_theta2_eta1
KineticRateUserDescription rate_cal_theta2_eta1("Calcite", 7.0, 6.0, false, 0.0, 0.0, 0.0, {"H+"}, {-3.0}, {0.0}, {0.0}, 2.0, 1.0, 55.0, 0.00315, DirectionChoiceEnum::BOTH, "H2O", 0.0, -1.0, 0.0)
rate_cal_theta0_eta05
KineticRateUserDescription rate_cal_theta0_eta05("Calcite", 7.0, 6.0, false, 0.0, 0.0, 0.0, {"H+"}, {-3.0}, {0.0}, {0.0}, 0.0, 0.5, 55.0, 0.00315, DirectionChoiceEnum::BOTH, "H2O", 0.0, -1.0, 0.0)
TEST
TEST(KineticRateUserDescriptionTest, exceptions)
Test exceptions in KineticRateUserDescription.
Definition: GeochemistryKineticRateCalculatorTest.C:288
rate_cal
KineticRateUserDescription rate_cal("Calcite", 7.0, 6.0, false, 0.0, 0.0, 0.0, {"H+"}, {-3.0}, {0.0}, {0.0}, 2.5, 0.8, 55.0, 0.00315, DirectionChoiceEnum::BOTH, "H2O", 0.0, -1.0, 0.0)
Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
db_kin
const GeochemicalDatabaseReader db_kin("database/moose_testdb.json")
KineticRateUserDescription
Holds a user-specified description of a kinetic rate.
Definition: GeochemistryKineticRateCalculator.h:75
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:167
ModelGeochemicalDatabase
Data structure to hold all relevant information from the database file.
Definition: PertinentGeochemicalSystem.h:99
EM::j
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
rate_ch4_death
KineticRateUserDescription rate_ch4_death("CH4(aq)", 1.5, 2.0, true, 0.0, 0.0, 0.0, {"H2O", "OH-", "O2(aq)", "CO2(aq)", "CaCO3"}, {3.0, 3.1, 3.2, 3.3, 3.4}, {0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0}, 0.8, 2.5, 66.0, 0.003, DirectionChoiceEnum::DEATH, "H2O", 0.0, -1.0, 0.0)
rate_ch4_precipitation
KineticRateUserDescription rate_ch4_precipitation("CH4(aq)", 1.5, 2.0, true, 0.0, 0.0, 0.0, {"H2O", "OH-", "O2(aq)", "CO2(aq)", "CaCO3"}, {3.0, 3.1, 3.2, 3.3, 3.4}, {0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0}, 0.8, 2.5, 66.0, 0.003, DirectionChoiceEnum::PRECIPITATION, "H2O", 0.0, -1.0, 0.0)
rate_ch4_dissolution
KineticRateUserDescription rate_ch4_dissolution("CH4(aq)", 1.5, 2.0, true, 0.0, 0.0, 0.0, {"H2O", "OH-", "O2(aq)", "CO2(aq)", "CaCO3"}, {3.0, 3.1, 3.2, 3.3, 3.4}, {0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0}, 0.8, 2.5, 66.0, 0.003, DirectionChoiceEnum::DISSOLUTION, "H2O", 0.0, -1.0, 0.0)
rate_cal_theta0_eta1
KineticRateUserDescription rate_cal_theta0_eta1("Calcite", 7.0, 6.0, false, 0.0, 0.0, 0.0, {"H+"}, {-3.0}, {0.0}, {0.0}, 0.0, 1.0, 55.0, 0.00315, DirectionChoiceEnum::BOTH, "H2O", 0.0, -1.0, 0.0)
std::pow
MooseUnits pow(const MooseUnits &, int)
rate_ch4_energy_captured
KineticRateUserDescription rate_ch4_energy_captured("CH4(aq)", 1.5, 2.0, true, 0.0, 0.0, 0.0, {"H2O", "OH-", "O2(aq)", "CO2(aq)", "CaCO3"}, {3.0, 3.1, 3.2, 3.3, 3.4}, {0.0, 0.0, 0.0, 0.0, 0.0}, {0.0, 0.0, 0.0, 0.0, 0.0}, 0.8, 2.5, 66.0, 0.003, DirectionChoiceEnum::BOTH, "H2O", 0.0, -1.0, 1.0)
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:193
GeochemicalDatabaseReader
Class for reading geochemical reactions from a MOOSE geochemical database.
Definition: GeochemicalDatabaseReader.h:213
rate_cal_theta0_eta2
KineticRateUserDescription rate_cal_theta0_eta2("Calcite", 7.0, 6.0, false, 0.0, 0.0, 0.0, {"H+"}, {-3.0}, {0.0}, {0.0}, 0.0, 2.0, 55.0, 0.00315, DirectionChoiceEnum::BOTH, "H2O", 0.0, -1.0, 0.0)
Generated on Thu Jul 17 2025 01:43:52 for https://mooseframework.inl.gov by   doxygen 1.8.14