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

Go to the source code of this file.

Functions

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

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 ( 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}
model_kin
PertinentGeochemicalSystem model_kin(db_kin, {"H2O", "H+", "HCO3-", "O2(aq)", "Ca++"}, {}, {}, {"Calcite"}, {"CH4(aq)"}, {}, "O2(aq)", "e-")
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
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
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
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

◆ TEST() [2/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}
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_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)
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)

◆ TEST() [3/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}
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() [4/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}
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_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_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_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)

◆ TEST() [5/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}
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)

◆ TEST() [6/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:76

Variable Documentation

◆ db_kin

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

◆ mgd_kin

const ModelGeochemicalDatabase& mgd_kin = model_kin.modelGeochemicalDatabase()

Definition at line 285 of file GeochemistryKineticRateCalculatorTest.C.

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

◆ model_kin

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

Referenced by TEST(), TEST(), TEST(), TEST(), and 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) ( "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(), 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) ( "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) ( "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) ( "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) ( "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) ( "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) ( "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(), TEST(), TEST(), TEST_F(), TEST_F(), 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) ( "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) ( "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) ( "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) ( "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) ( "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) ( "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 Sep 11 2026 20:37:17 for https://mooseframework.inl.gov by doxygen 1.9.8