https://mooseframework.inl.gov
Loading...
Searching...
No Matches
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
15const GeochemicalDatabaseReader db_kin("database/moose_testdb.json");
16PertinentGeochemicalSystem model_kin(db_kin,
17 {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++"},
18 {},
19 {},
20 {"Calcite"},
21 {"CH4(aq)"},
22 {},
23 "O2(aq)",
24 "e-");
25KineticRateUserDescription 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);
45KineticRateUserDescription 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);
65KineticRateUserDescription 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);
85KineticRateUserDescription 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);
105KineticRateUserDescription 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);
125KineticRateUserDescription 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);
145KineticRateUserDescription 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);
165KineticRateUserDescription 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);
185KineticRateUserDescription 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);
205KineticRateUserDescription 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);
225KineticRateUserDescription 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);
245KineticRateUserDescription 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);
265KineticRateUserDescription 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);
285const ModelGeochemicalDatabase & mgd_kin = model_kin.modelGeochemicalDatabase();
286
288TEST(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
425TEST(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
992TEST(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
1568TEST(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
2031TEST(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
2110TEST(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}
model_kin
PertinentGeochemicalSystem model_kin(db_kin, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++"}, {}, {}, {"Calcite"}, {"CH4(aq)"}, {}, "O2(aq)", "e-")
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)
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)
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)
TEST
TEST(KineticRateUserDescriptionTest, exceptions)
Test exceptions in KineticRateUserDescription.
Definition GeochemistryKineticRateCalculatorTest.C:288
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)
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)
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)
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_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)
mgd_kin
const ModelGeochemicalDatabase & mgd_kin
Definition GeochemistryKineticRateCalculatorTest.C:285
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)
db_kin
const GeochemicalDatabaseReader db_kin("database/moose_testdb.json")
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_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_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)
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)
GeochemicalDatabaseReader
Class for reading geochemical reactions from a MOOSE geochemical database.
Definition GeochemicalDatabaseReader.h:216
PertinentGeochemicalSystem
Constructs and stores a minimal amount of information that is pertinent to the user-defined geochemic...
Definition PertinentGeochemicalSystem.h:387
PertinentGeochemicalSystem::modelGeochemicalDatabase
const ModelGeochemicalDatabase & modelGeochemicalDatabase() const
Return a reference to the ModelGeochemicalDatabase structure.
Definition PertinentGeochemicalSystem.C:850
PertinentGeochemicalSystem::addKineticRate
void addKineticRate(const KineticRateUserDescription &description)
Adds a rate description for kinetic_species.
Definition PertinentGeochemicalSystem.C:856
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
KineticRateUserDescription
Holds a user-specified description of a kinetic rate.
Definition GeochemistryKineticRateCalculator.h:76
ModelGeochemicalDatabase
Data structure to hold all relevant information from the database file.
Definition PertinentGeochemicalSystem.h:102
ModelGeochemicalDatabase::eqm_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::kin_stoichiometry
DenseMatrix< Real > kin_stoichiometry
kin_stoichiometry(i, j) = stoichiometric coefficient for kinetic species "i" in terms of the basis sp...
Definition PertinentGeochemicalSystem.h:334
ModelGeochemicalDatabase::basis_species_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_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
Generated on Fri Aug 21 2026 13:41:03 for https://mooseframework.inl.gov by doxygen 1.9.8