Fix bug in parser (#723)

* fix tutorial discrepancy

* update rate constant tutorials

* update tests to cover original underflow error conditions

README.md

@@ -167,6 +167,7 @@ int main(const int argc, const char *argv[])
 
   state.conditions_[0].temperature_ = 287.45;  // K
   state.conditions_[0].pressure_ = 101319.9;   // Pa
+  state.conditions_[0].CalculateIdealAirDensity();
   state.SetConcentration(foo, 20.0);           // mol m-3
 
   state.PrintHeader();

docs/source/user_guide/rate_constant_tutorial.rst

@@ -69,13 +69,13 @@ rosenbrock solver at the top of the file.
 
   .. literalinclude:: ../../../test/tutorial/test_rate_constants_no_user_defined_by_hand.cpp
     :language: cpp
-    :lines: 1-13
+    :lines: 1-14
 
 After that, we'll use the ``micm`` namespace so that we don't have to repeat it everywhere we need it.
 
   .. literalinclude:: ../../../test/tutorial/test_rate_constants_no_user_defined_by_hand.cpp
     :language: cpp
-    :lines: 14-15
+    :lines: 16-17
 
 To create a :cpp:class:`micm::RosenbrockSolver`, we have to define a chemical system (:cpp:class:`micm::System`)
 and our reactions, which will be a vector of :cpp:class:`micm::Process` We will use the species to define these.
@@ -90,22 +90,24 @@ and our reactions, which will be a vector of :cpp:class:`micm::Process` We will
 
         .. literalinclude:: ../../../test/tutorial/test_rate_constants_no_user_defined_by_hand.cpp
           :language: cpp
-          :lines: 19-30
+          :lines: 25-36
 
         Now that we have a gas phase and our species, we can start building the reactions. Two things to note are that
         stoichiemtric coefficients for reactants are represented by repeating that product as many times as you need.
         To specify the yield of a product, we've created a typedef :cpp:type:`micm::Yield` 
-        and a function :cpp:func:`micm::Yields` that produces these.
+        and a function :cpp:func:`micm::Yields` that produces these. Note that we add a conversion for
+        some rate constant parameters to be consistent with the configuration file that expects rate constants
+        to be in cm^3/molecule/s. (All units will be mks in the next version of the configuration file format.)
 
         .. literalinclude:: ../../../test/tutorial/test_rate_constants_no_user_defined_by_hand.cpp
           :language: cpp
-          :lines: 32-96
+          :lines: 38-102
 
         And finally we define our chemical system and reactions
 
         .. literalinclude:: ../../../test/tutorial/test_rate_constants_no_user_defined_by_hand.cpp
           :language: cpp
-          :lines: 98-99
+          :lines: 104-105
 
     .. tab-item:: OpenAtmos Configuration reading
 
@@ -114,15 +116,15 @@ and our reactions, which will be a vector of :cpp:class:`micm::Process` We will
 
         .. literalinclude:: ../../../test/tutorial/test_rate_constants_no_user_defined_with_config.cpp
           :language: cpp
-          :lines: 20-32
+          :lines: 21-32
 
 Now that we have a chemical system and a list of reactions, we can create the RosenbrockSolver.
 There are several ways to configure the solver. Here we are using a three stage solver. More options
 can be found in the :cpp:class:`micm::RosenbrockSolverParameters` and in the :ref:`Solver Configurations` tutorial.
 
   .. literalinclude:: ../../../test/tutorial/test_rate_constants_no_user_defined_by_hand.cpp
     :language: cpp
-    :lines: 101
+    :lines: 107-110
 
 The rosenbrock solver will provide us a state, which we can use to set the concentrations,
 custom rate parameters, and temperature and pressure. Note that setting the custom rate paramters is different depending
@@ -140,20 +142,20 @@ Initializing the state
 
         .. literalinclude:: ../../../test/tutorial/test_rate_constants_no_user_defined_by_hand.cpp
           :language: cpp
-          :lines: 102-116
+          :lines: 111-126
 
     .. tab-item:: OpenAtmos Configuration reading
 
         .. literalinclude:: ../../../test/tutorial/test_rate_constants_no_user_defined_with_config.cpp
           :language: cpp
-          :lines: 36-54
+          :lines: 44-63
 
 
 Finally, we are ready to pick a timestep and solve the system.
 
   .. literalinclude:: ../../../test/tutorial/test_rate_constants_no_user_defined_by_hand.cpp
     :language: cpp
-    :lines: 118-142
+    :lines: 129-151
 
 
 This is the output:
@@ -163,14 +165,14 @@ This is the output:
    :header: "time", "A", "B", "C", "D", "E", "F", "G"
    :widths: 10, 15, 15, 15, 15, 15, 15, 15
 
-      0,   1.00e+00,   0.00e+00,   0.00e+00,   0.00e+00,   0.00e+00,   0.00e+00,   0.00e+00
-    500,   3.18e-09,   3.66e-09,   9.83e-01,   3.88e-14,   1.41e-03,   2.02e-13,   7.92e-03
-    1000,   1.14e-14,   1.31e-14,   9.66e-01,   1.39e-19,   1.40e-03,   7.24e-19,   1.64e-02
-    1500,   4.09e-20,   4.71e-20,   9.49e-01,   4.98e-25,   1.39e-03,   2.59e-24,   2.48e-02
-    2000,   1.47e-25,   1.69e-25,   9.33e-01,   1.79e-30,   1.38e-03,   9.30e-30,   3.30e-02
-    2500,   5.26e-31,   6.05e-31,   9.17e-01,   6.40e-36,   1.37e-03,   3.33e-35,   4.11e-02
-    3000,   1.89e-36,   2.17e-36,   9.01e-01,   2.30e-41,   1.36e-03,   1.20e-40,   4.90e-02
-    3500,   6.77e-42,   7.78e-42,   8.85e-01,   8.23e-47,   1.34e-03,   4.29e-46,   5.68e-02
-    4000,   2.43e-47,   2.79e-47,   8.70e-01,   2.95e-52,   1.33e-03,   1.54e-51,   6.44e-02
-    4500,   8.70e-53,   1.00e-52,   8.55e-01,   1.06e-57,   1.32e-03,   5.51e-57,   7.20e-02
-    5000,   3.12e-58,   3.59e-58,   8.40e-01,   3.80e-63,   1.31e-03,   1.98e-62,   7.94e-02
+     0,   1.00e+00,   0.00e+00,   0.00e+00,   0.00e+00,   0.00e+00,   0.00e+00,   0.00e+00
+   500,   8.54e-01,   4.57e-04,   1.44e-01,   1.55e-04,   6.47e-14,   1.23e-22,   6.44e-04
+  1000,   7.30e-01,   3.90e-04,   2.65e-01,   2.89e-04,   2.53e-13,   2.28e-22,   2.44e-03
+  1500,   6.23e-01,   3.33e-04,   3.66e-01,   4.02e-04,   2.98e-13,   3.18e-22,   5.20e-03
+  2000,   5.32e-01,   2.85e-04,   4.49e-01,   5.00e-04,   3.30e-13,   3.95e-22,   8.77e-03
+  2500,   4.55e-01,   2.43e-04,   5.18e-01,   5.83e-04,   3.55e-13,   4.61e-22,   1.30e-02
+  3000,   3.88e-01,   2.08e-04,   5.75e-01,   6.54e-04,   3.74e-13,   5.17e-22,   1.78e-02
+  3500,   3.32e-01,   1.77e-04,   6.21e-01,   7.14e-04,   3.88e-13,   5.65e-22,   2.30e-02
+  4000,   2.83e-01,   1.52e-04,   6.59e-01,   7.66e-04,   4.00e-13,   6.06e-22,   2.86e-02
+  4500,   2.42e-01,   1.29e-04,   6.88e-01,   8.10e-04,   4.09e-13,   6.41e-22,   3.45e-02
+  5000,   2.07e-01,   1.11e-04,   7.11e-01,   8.48e-04,   4.15e-13,   6.71e-22,   4.06e-02

docs/source/user_guide/user_defined_rate_constant_tutorial.rst

@@ -153,12 +153,12 @@ Finally, set and upate the rate constants as needed:
           // solving until we finish
           double elapsed_solve_time = 0;
       +   state.SetCustomRateParameter("my photolysis rate", photo_rate);
+          solver.CalculateRateConstants(state);
 
           while (elapsed_solve_time < time_step)
           {
             auto result = solver.Solve(time_step - elapsed_solve_time, state);
             elapsed_solve_time = result.final_time_;
-            state.variables_[0] = result.result_.AsVector();
           }
 
           print_state(time_step * (i + 1), state);
@@ -179,6 +179,7 @@ Finally, set and upate the rate constants as needed:
       + // these rates are constant through the simulation
       + state.SetCustomRateParameter("EMIS.my emission rate", emission_rate);
       + state.SetCustomRateParameter("LOSS.my loss rate", loss);
+
         // solve for ten iterations
         for (int i = 0; i < 10; ++i)
         {
@@ -188,12 +189,12 @@ Finally, set and upate the rate constants as needed:
           // solving until we finish
           double elapsed_solve_time = 0;
       +   state.SetCustomRateParameter("PHOTO.my photolysis rate", photo_rate);
+          solver.CalculateRateConstants(state);
 
           while (elapsed_solve_time < time_step)
           {
             auto result = solver.Solve(time_step - elapsed_solve_time, state);
             elapsed_solve_time = result.final_time_;
-            state.variables_[0] = result.result_.AsVector();
           }
 
           print_state(time_step * (i + 1), state);
@@ -207,14 +208,14 @@ the :ref:`Rate constants` tutorial's result.
    :header: "time", "A", "B", "C", "D", "E", "F", "G"
    :widths: 10, 15, 15, 15, 15, 15, 15, 15
 
-   "0", "1.00e+00", "0.00e+00", "0.00e+00", "0.00e+00", "0.00e+00", "0.00e+00", "0.00e+00"
-   "500", "3.18e-09", "3.66e-09", "9.83e-01", "3.88e-14", "1.41e-03", "2.02e-13", "7.92e-03"
-   "1000", "1.14e-14", "1.31e-14", "9.66e-01", "1.39e-19", "1.40e-03", "7.24e-19", "1.64e-02"
-   "1500", "7.27e-20", "6.40e-20", "9.49e-01", "6.53e-25", "1.39e-03", "3.19e-24", "2.48e-02"
-   "2000", "3.17e-20", "1.70e-20", "9.33e-01", "1.55e-25", "1.38e-03", "5.92e-25", "3.30e-02"
-   "2500", "3.17e-20", "1.70e-20", "9.17e-01", "1.55e-25", "1.37e-03", "5.92e-25", "4.11e-02"
-   "3000", "3.17e-20", "1.70e-20", "9.01e-01", "1.55e-25", "1.36e-03", "5.92e-25", "4.90e-02"
-   "3500", "3.17e-20", "1.70e-20", "8.85e-01", "1.55e-25", "1.34e-03", "5.92e-25", "5.68e-02"
-   "4000", "3.17e-20", "1.70e-20", "8.70e-01", "1.55e-25", "1.33e-03", "5.92e-25", "6.44e-02"
-   "4500", "3.17e-20", "1.70e-20", "8.55e-01", "1.55e-25", "1.32e-03", "5.92e-25", "7.20e-02"
-   "5000", "3.17e-20", "1.70e-20", "8.40e-01", "1.55e-25", "1.31e-03", "5.92e-25", "7.94e-02"
+     0,   1.00e+00,   0.00e+00,   0.00e+00,   0.00e+00,   0.00e+00,   0.00e+00,   0.00e+00
+   500,   8.54e-01,   4.57e-04,   1.44e-01,   1.55e-04,   6.47e-14,   1.23e-22,   6.44e-04
+  1000,   7.30e-01,   3.90e-04,   2.65e-01,   2.89e-04,   2.53e-13,   2.28e-22,   2.44e-03
+  1500,   6.23e-01,   3.33e-04,   3.66e-01,   4.02e-04,   2.98e-13,   3.18e-22,   5.20e-03
+  2000,   5.32e-01,   2.85e-04,   4.49e-01,   5.00e-04,   3.30e-13,   3.95e-22,   8.77e-03
+  2500,   4.55e-01,   2.43e-04,   5.18e-01,   5.83e-04,   3.55e-13,   4.61e-22,   1.30e-02
+  3000,   3.88e-01,   2.08e-04,   5.75e-01,   6.54e-04,   3.74e-13,   5.17e-22,   1.78e-02
+  3500,   3.32e-01,   1.77e-04,   6.21e-01,   7.14e-04,   3.88e-13,   5.65e-22,   2.30e-02
+  4000,   2.83e-01,   1.52e-04,   6.59e-01,   7.66e-04,   4.00e-13,   6.06e-22,   2.86e-02
+  4500,   2.42e-01,   1.29e-04,   6.88e-01,   8.10e-04,   4.09e-13,   6.41e-22,   3.45e-02
+  5000,   2.07e-01,   1.11e-04,   7.11e-01,   8.48e-04,   4.15e-13,   6.71e-22,   4.06e-02

include/micm/configure/solver_config.hpp

@@ -94,7 +94,7 @@ inline std::error_code make_error_code(MicmConfigErrc e)
 namespace micm
 {
 
-  constexpr double MOLES_M3_TO_MOLECULES_CM3 = 1.0e-6 * 6.02214076e23;
+  constexpr double MOLES_M3_TO_MOLECULES_CM3 = 1.0e-6 * constants::AVOGADRO_CONSTANT;
 
   // Solver parameters
   struct SolverParameters
@@ -538,7 +538,8 @@ namespace micm
       {
         parameters.A_ = object["A"].as<double>();
       }
-      parameters.A_ *= std::pow(MOLES_M3_TO_MOLECULES_CM3, reactants.size() - 1);
+      int n_reactants = reactants.size();
+      parameters.A_ *= std::pow(MOLES_M3_TO_MOLECULES_CM3, n_reactants - 1);
       if (object["B"])
       {
         parameters.B_ = object["B"].as<double>();
@@ -587,7 +588,8 @@ namespace micm
         parameters.k0_A_ = object["k0_A"].as<double>();
       }
       // Account for the conversion of reactant concentrations (including M) to molecules cm-3
-      parameters.k0_A_ *= std::pow(MOLES_M3_TO_MOLECULES_CM3, reactants.size());
+      int n_reactants = reactants.size();
+      parameters.k0_A_ *= std::pow(MOLES_M3_TO_MOLECULES_CM3, n_reactants);
       if (object["k0_B"])
       {
         parameters.k0_B_ = object["k0_B"].as<double>();
@@ -601,7 +603,7 @@ namespace micm
         parameters.kinf_A_ = object["kinf_A"].as<double>();
       }
       // Account for terms in denominator and exponent that include [M] but not other reactants
-      parameters.kinf_A_ *= std::pow(MOLES_M3_TO_MOLECULES_CM3, reactants.size() - 1);
+      parameters.kinf_A_ *= std::pow(MOLES_M3_TO_MOLECULES_CM3, n_reactants - 1);
       if (object["kinf_B"])
       {
         parameters.kinf_B_ = object["kinf_B"].as<double>();
@@ -640,7 +642,8 @@ namespace micm
         parameters.k0_A_ = object["k0_A"].as<double>();
       }
       // Account for the conversion of reactant concentrations (including M) to molecules cm-3
-      parameters.k0_A_ *= std::pow(MOLES_M3_TO_MOLECULES_CM3, reactants.size() - 1);
+      int n_reactants = reactants.size();
+      parameters.k0_A_ *= std::pow(MOLES_M3_TO_MOLECULES_CM3, n_reactants - 1);
       if (object["k0_B"])
       {
         parameters.k0_B_ = object["k0_B"].as<double>();
@@ -654,7 +657,7 @@ namespace micm
         parameters.kinf_A_ = object["kinf_A"].as<double>();
       }
       // Account for terms in denominator and exponent that include [M] but not other reactants
-      parameters.kinf_A_ *= std::pow(MOLES_M3_TO_MOLECULES_CM3, reactants.size() - 2);
+      parameters.kinf_A_ *= std::pow(MOLES_M3_TO_MOLECULES_CM3, n_reactants - 2);
       if (object["kinf_B"])
       {
         parameters.kinf_B_ = object["kinf_B"].as<double>();
@@ -697,7 +700,8 @@ namespace micm
       BranchedRateConstantParameters parameters;
       parameters.X_ = object[X].as<double>();
       // Account for the conversion of reactant concentrations to molecules cm-3
-      parameters.X_ *= std::pow(MOLES_M3_TO_MOLECULES_CM3, reactants.size() - 1);
+      int n_reactants = reactants.size();
+      parameters.X_ *= std::pow(MOLES_M3_TO_MOLECULES_CM3, n_reactants - 1);
       parameters.Y_ = object[Y].as<double>();
       parameters.a0_ = object[A0].as<double>();
       parameters.n_ = object[N].as<int>();
@@ -731,7 +735,8 @@ namespace micm
         parameters.A_ = object["A"].as<double>();
       }
       // Account for the conversion of reactant concentrations to molecules cm-3
-      parameters.A_ *= std::pow(MOLES_M3_TO_MOLECULES_CM3, reactants.size() - 1);
+      int n_reactants = reactants.size();
+      parameters.A_ *= std::pow(MOLES_M3_TO_MOLECULES_CM3, n_reactants - 1);
       if (object["B"])
       {
         parameters.B_ = object["B"].as<double>();

include/micm/system/conditions.hpp

@@ -2,13 +2,21 @@
 // SPDX-License-Identifier: Apache-2.0
 #pragma once
 
+#include <micm/util/constants.hpp>
+
 namespace micm
 {
   /// @brief Environemental conditions
   struct Conditions
   {
     double temperature_{ 0.0 };  // K
     double pressure_{ 0.0 };     // Pa
-    double air_density_{ 1.0 };  // mol m-3
+    double air_density_{ 0.0 };  // mol m-3
+    void CalculateIdealAirDensity();
   };
+
+  inline void Conditions::CalculateIdealAirDensity()
+  {
+    air_density_ = pressure_ / (constants::GAS_CONSTANT * temperature_);
+  }
 }  // namespace micm

test/tutorial/configs/rate_constants_no_user_defined/reactions.yaml

@@ -7,7 +7,7 @@ camp-data:
           A: {}
         products:
           B: {}
-        A: 2.15e-1
+        A: 2.15e-4
         B: 0
         C: 110.0
       - type: "BRANCHED"

test/tutorial/configs/rate_constants_user_defined/reactions.json

@@ -12,7 +12,7 @@
           "products": {
             "B": {}
           },
-          "A": 2.15e-1,
+          "A": 2.15e-4,
           "B": 0,
           "C": 110.0
         },

test/tutorial/test_README_example.cpp

@@ -40,6 +40,7 @@ int main(const int argc, const char *argv[])
 
   state.conditions_[0].temperature_ = 287.45;  // K
   state.conditions_[0].pressure_ = 101319.9;   // Pa
+  state.conditions_[0].CalculateIdealAirDensity();
   state.SetConcentration(foo, 20.0);           // mol m-3
 
   state.PrintHeader();

test/tutorial/test_rate_constants_no_user_defined_by_hand.cpp

@@ -7,6 +7,7 @@
 #include <micm/process/tunneling_rate_constant.hpp>
 #include <micm/solver/rosenbrock.hpp>
 #include <micm/solver/solver_builder.hpp>
+#include <micm/util/constants.hpp>
 
 #include <iomanip>
 #include <iostream>
@@ -15,6 +16,10 @@
 // Use our namespace so that this example is easier to read
 using namespace micm;
 
+// Conversion factor from moles m-3 to molecules cm-3 for consistency
+// with the configuraion file
+constexpr double MOLES_M3_TO_MOLECULES_CM3 = 1.0e-6 * constants::AVOGADRO_CONSTANT;
+
 int main(const int argc, const char* argv[])
 {
   auto a = Species("A");
@@ -33,7 +38,7 @@ int main(const int argc, const char* argv[])
   Process r1 = Process::Create()
                    .SetReactants({ a })
                    .SetProducts({ Yields(b, 1) })
-                   .SetRateConstant(ArrheniusRateConstant({ .A_ = 2.15e-1, .B_ = 0, .C_ = 110 }))
+                   .SetRateConstant(ArrheniusRateConstant({ .A_ = 2.15e-4, .B_ = 0, .C_ = 110 }))
                    .SetPhase(gas_phase);
 
   // a branched reaction has two output pathways
@@ -68,7 +73,7 @@ int main(const int argc, const char* argv[])
                    .SetRateConstant(TernaryChemicalActivationRateConstant({ .k0_A_ = 1.2,
                                                                           .k0_B_ = 2.3,
                                                                           .k0_C_ = 302.3,
-                                                                          .kinf_A_ = 2.6,
+                                                                          .kinf_A_ = 2.6 / MOLES_M3_TO_MOLECULES_CM3,
                                                                           .kinf_B_ = -3.1,
                                                                           .kinf_C_ = 402.1,
                                                                           .Fc_ = 0.9,
@@ -80,10 +85,10 @@ int main(const int argc, const char* argv[])
   Process r6 = Process::Create()
                    .SetReactants({ e, e })
                    .SetProducts({ Yields(g, 1) })
-                   .SetRateConstant(TroeRateConstant({ .k0_A_ = 1.2e4,
+                   .SetRateConstant(TroeRateConstant({ .k0_A_ = 1.2e4 * MOLES_M3_TO_MOLECULES_CM3 * MOLES_M3_TO_MOLECULES_CM3,
                                                      .k0_B_ = 167.0,
                                                      .k0_C_ = 3.0,
-                                                     .kinf_A_ = 136.0,
+                                                     .kinf_A_ = 136.0 * MOLES_M3_TO_MOLECULES_CM3,
                                                      .kinf_B_ = 5.0,
                                                      .kinf_C_ = 24.0,
                                                      .Fc_ = 0.9,
@@ -107,6 +112,7 @@ int main(const int argc, const char* argv[])
 
   state.conditions_[0].temperature_ = 287.45;  // K
   state.conditions_[0].pressure_ = 101319.9;   // Pa
+  state.conditions_[0].CalculateIdealAirDensity();
 
   state.SetConcentration(a, 1.0);  // mol m-3
   state.SetConcentration(b, 0.0);  // mol m-3

test/tutorial/test_rate_constants_no_user_defined_with_config.cpp

@@ -45,6 +45,7 @@ int main(const int argc, const char* argv[])
 
   state.conditions_[0].temperature_ = 287.45;  // K
   state.conditions_[0].pressure_ = 101319.9;   // Pa
+  state.conditions_[0].CalculateIdealAirDensity();
 
   std::unordered_map<std::string, std::vector<double>> intial_concentration = {
     { "A", { 1.0 } },  // mol m-3