https://mooseframework.inl.gov
Functions | Variables
GeochemistryKineticRateCalculatorTest.C File Reference

Go to the source code of this file.

Functions

 TEST (KineticRateUserDescriptionTest, exceptions)
 Test exceptions in KineticRateUserDescription. More...
 
 TEST (GeochemistryKineticRateCalculatorTest, exceptions)
 Test exceptions. More...
 
 TEST (GeochemistryKineticRateCalculatorTest, rate1)
 Test rate calculations in various non-monod scenarios. More...
 
 TEST (GeochemistryKineticRateCalculatorTest, rate_monod)
 Test rate calculations in various monod scenarios. More...
 
 TEST (GeochemistryKineticRateCalculatorTest, rate_energy_captured)
 Test rate calculations for nonzero energy captured. More...
 
 TEST (GeochemistryKineticRateCalculatorTest, rate_eta_theta)
 Test rate calculations for eta and theta edge cases. More...
 

Variables

const GeochemicalDatabaseReader db_kin ("database/moose_testdb.json")
 
PertinentGeochemicalSystem model_kin (db_kin, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++"}, {}, {}, {"Calcite"}, {"CH4(aq)"}, {}, "O2(aq)", "e-")
 
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)
 
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)
 
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)
 
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)
 
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)
 
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)
 
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)
 
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)
 
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)
 
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)
 
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)
 
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)
 
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)
 
const ModelGeochemicalDatabase & mgd_kin = model_kin.modelGeochemicalDatabase()
 

Function Documentation

◆ TEST() [1/6]

TEST ( KineticRateUserDescriptionTest  ,
exceptions   
)

Test exceptions in KineticRateUserDescription.

Definition at line 288 of file GeochemistryKineticRateCalculatorTest.C.

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 }
KineticRateUserDescription
Holds a user-specified description of a kinetic rate.
Definition: GeochemistryKineticRateCalculator.h:75

◆ TEST() [2/6]

TEST ( GeochemistryKineticRateCalculatorTest  ,
exceptions   
)

Test exceptions.

Definition at line 425 of file GeochemistryKineticRateCalculatorTest.C.

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 }
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
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
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::addKineticRate
void addKineticRate(const KineticRateUserDescription &description)
Adds a rate description for kinetic_species.
Definition: PertinentGeochemicalSystem.C:856
Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
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::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

◆ TEST() [3/6]

TEST ( GeochemistryKineticRateCalculatorTest  ,
rate1   
)

Test rate calculations in various non-monod scenarios.

Definition at line 992 of file GeochemistryKineticRateCalculatorTest.C.

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 }
GeochemistryConstants::CELSIUS_TO_KELVIN
constexpr Real CELSIUS_TO_KELVIN
Definition: GeochemistryConstants.h:22
eps
const Real eps
Definition: GeochemistrySpeciesSwapperTest.C:14
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)
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
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)
model_kin
PertinentGeochemicalSystem model_kin(db_kin, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++"}, {}, {}, {"Calcite"}, {"CH4(aq)"}, {}, "O2(aq)", "e-")
PertinentGeochemicalSystem::addKineticRate
void addKineticRate(const KineticRateUserDescription &description)
Adds a rate description for kinetic_species.
Definition: PertinentGeochemicalSystem.C:856
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
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)
std::pow
MooseUnits pow(const MooseUnits &, int)

◆ TEST() [4/6]

TEST ( GeochemistryKineticRateCalculatorTest  ,
rate_monod   
)

Test rate calculations in various monod scenarios.

Definition at line 1568 of file GeochemistryKineticRateCalculatorTest.C.

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 }
GeochemistryConstants::CELSIUS_TO_KELVIN
constexpr Real CELSIUS_TO_KELVIN
Definition: GeochemistryConstants.h:22
eps
const Real eps
Definition: GeochemistrySpeciesSwapperTest.C:14
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)
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)
model_kin
PertinentGeochemicalSystem model_kin(db_kin, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++"}, {}, {}, {"Calcite"}, {"CH4(aq)"}, {}, "O2(aq)", "e-")
PertinentGeochemicalSystem::addKineticRate
void addKineticRate(const KineticRateUserDescription &description)
Adds a rate description for kinetic_species.
Definition: PertinentGeochemicalSystem.C:856
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
EM::j
static const std::complex< double > j(0, 1)
Complex number "j" (also known as "i")
std::pow
MooseUnits pow(const MooseUnits &, int)

◆ TEST() [5/6]

TEST ( GeochemistryKineticRateCalculatorTest  ,
rate_energy_captured   
)

Test rate calculations for nonzero energy captured.

Definition at line 2031 of file GeochemistryKineticRateCalculatorTest.C.

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 }
GeochemistryConstants::CELSIUS_TO_KELVIN
constexpr Real CELSIUS_TO_KELVIN
Definition: GeochemistryConstants.h:22
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
model_kin
PertinentGeochemicalSystem model_kin(db_kin, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++"}, {}, {}, {"Calcite"}, {"CH4(aq)"}, {}, "O2(aq)", "e-")
PertinentGeochemicalSystem::addKineticRate
void addKineticRate(const KineticRateUserDescription &description)
Adds a rate description for kinetic_species.
Definition: PertinentGeochemicalSystem.C:856
Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
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)

◆ TEST() [6/6]

TEST ( GeochemistryKineticRateCalculatorTest  ,
rate_eta_theta   
)

Test rate calculations for eta and theta edge cases.

Definition at line 2110 of file GeochemistryKineticRateCalculatorTest.C.

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
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
model_kin
PertinentGeochemicalSystem model_kin(db_kin, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++"}, {}, {}, {"Calcite"}, {"CH4(aq)"}, {}, "O2(aq)", "e-")
PertinentGeochemicalSystem::addKineticRate
void addKineticRate(const KineticRateUserDescription &description)
Adds a rate description for kinetic_species.
Definition: PertinentGeochemicalSystem.C:856
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)
Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
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_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)

Variable Documentation

◆ db_kin

const GeochemicalDatabaseReader db_kin("database/moose_testdb.json")

◆ mgd_kin

const ModelGeochemicalDatabase& mgd_kin = model_kin.modelGeochemicalDatabase()

Definition at line 285 of file GeochemistryKineticRateCalculatorTest.C.

Referenced by TEST(), and TEST_F().

◆ model_kin

PertinentGeochemicalSystem model_kin(db_kin, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++"}, {}, {}, {"Calcite"}, {"CH4(aq)"}, {}, "O2(aq)", "e-")

Referenced by TEST().

◆ 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)

Referenced by TEST(), and TEST_F().

◆ 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)

Referenced by TEST().

◆ 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)

Referenced by TEST().

◆ 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)

Referenced by TEST().

◆ 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)

Referenced by TEST().

◆ 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)

Referenced by TEST().

◆ 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)

Referenced by TEST(), and TEST_F().

◆ 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)

Referenced by TEST(), and TEST_F().

◆ 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)

Referenced by TEST().

◆ 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)

Referenced by TEST().

◆ 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)

Referenced by TEST().

◆ 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)

Referenced by TEST().

◆ 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)

Referenced by TEST().

Generated on Fri Jul 18 2025 11:57:53 for https://mooseframework.inl.gov by   doxygen 1.8.14