Test kpp micm (#295)

* Start on test_kpp_to_micm program. * Start on CMakeLists for test_kpp_to_micm. * Added kpp test subdir. * Renamed test. * Set config path relative to build dir. * Added kpp_to_micm.py generated JSON. * Copy kpp_to_micm JSON to builds/configs for testing. * Check ConfigParseStatus. * Include regression/RosenbrockChapman/chapman_ode_solver. * Instantiate micm::RosenbrockSolver. * Get solver state and initialize T and p. * Added initial concentrations for Chapman test from Seinfeld and Pandis. * Added solver.Solve iteration. * Added --to_si_units option to kpp_to_micm.py. * Updated Chapman reactions JSON to SI units. * Use SI units. * Use SI units. * Added print_state to test_kpp_to_micm. * Improved print_state format. * Made a KPP config for Chapman without NOx, for comparison to RosenbrockChapman. * Removed NOx from Chapman KPP test. * Added M to print_state. * Set photolysis rates. * Print custom param labels. * Use state.SetCustomRateParameter to set photolysis rates. * Use state.SetConcentration for each species. * Changed print precision to 3 decimal places. * Changed print precision to 3 decimal places. * Set initial O3 concentration to 0. * Cleaned up comments. * Change to SI conversion factor to (N_Avogradro 1e-6)^(N_reactants-1).
NCAR · Oct 12, 2023 · 88d6b5b · 88d6b5b
1 parent 0996e58
commit 88d6b5b
Show file tree

Hide file tree

Showing 8 changed files with 341 additions and 8 deletions.
diff --git a/etc/configs/kpp/chapman.eqn b/etc/configs/kpp/chapman.eqn
@@ -0,0 +1,10 @@
+#EQUATIONS { Small Stratospheric Mechanism }
+
+<R1>  O2   + hv = 2O        : (2.643E-10) * SUN*SUN*SUN;
+<R2>  O    + O2 = O3        : (8.018E-17);
+<R3>  O3   + hv = O + O2    : (6.120E-04) * SUN;
+<R4>  O    + O3 = 2O2       : (1.576E-15);
+<R5>  O3   + hv = O1D + O2  : (1.070E-03) * SUN*SUN;
+<R6>  O1D  + M  = O   + M   : (7.110E-11);
+<R7>  O1D  + O3 = 2O2       : (1.200E-10);
+
diff --git a/etc/configs/kpp/chapman.spc b/etc/configs/kpp/chapman.spc
@@ -0,0 +1,11 @@
+#INCLUDE atoms.kpp
+
+#DEFVAR
+O   = O;          { Oxygen atomic ground state }
+O1D = O;          { Oxygen atomic excited state }
+O3  = O + O + O;  { Ozone }
+
+#DEFFIX
+M   = O + O + N + N;  { Atmospheric generic molecule }
+O2  = O + O;          { Molecular oxygen }
+
diff --git a/etc/scripts/kpp_to_micm.py b/etc/scripts/kpp_to_micm.py
@@ -31,6 +31,7 @@
 Revision History:
     v1.00 2023/08/03 Initial implementation
     v1.01 2023/08/16 Added method parse_kpp_arrhenius
+    v1.02 2023/10/05 Added unit conversion option
 """
 
 import os
@@ -40,7 +41,12 @@
 import json
 from glob import glob
 
-__version__ = 'v1.01'
+__version__ = 'v1.02'
+
+"""
+Physical Constants
+"""
+N_Avogadro = 6.02214076e23 
 
 
 def read_kpp_config(kpp_dir, kpp_name):
@@ -129,17 +135,21 @@ def micm_species_json(lines, fixed=False, tolerance=1.0e-12):
     return species_json
 
 
-def parse_kpp_arrhenius(kpp_str):
+def parse_kpp_arrhenius(kpp_str, to_si_units=False,
+    N_reactants=2):
     """
     Parse KPP Arrhenius reaction
 
     Parameters
         (str) kpp_str: Arrhenius reaction string
+        (bool) to_si_units: convert A coefficient
+        from (cm^3 molecule-1)^(N-1) s-1 to (m^3 mol-1)^(N-1) s-1
+        where N is the number of reactants
+        (int) N_reactants: number of reactants
 
     Returns
         (dict): MICM Arrhenius reaction coefficients
 
-
     Arrhenius formula from KPP
     --------------------------
 
@@ -188,16 +198,20 @@ def parse_kpp_arrhenius(kpp_str):
         arr_dict['D'] = 300.0
     else:
         logging.error('unrecognized KPP Arrhenius syntax')
+    if to_si_units:
+        arr_dict['A'] *= (N_Avogadro * 1.0e-6)**(N_reactants - 1)
     logging.debug(arr_dict)
     return arr_dict
 
 
-def micm_equation_json(lines):
+def micm_equation_json(lines, to_si_units=False):
     """
     Generate MICM equation JSON
 
     Parameters
         (list of str) lines: lines of equation section
+        (bool) to_si_units: convert A coefficient
+        from cm^3 molecule-1 s-1 to m^3 mol-1 s-1
 
     Returns
         (list of dict): list of MICM equation entries
@@ -227,17 +241,22 @@ def micm_equation_json(lines):
         reactants = [reactant.strip().lstrip() for reactant in reactants]
         products = [product.strip().lstrip() for product in products]
 
+        N_reactants = len(reactants)
+
         equation_dict = dict()
 
         if 'SUN' in coeffs:
-            equation_dict['type'] = 'PHOTOLYSIS' 
+            equation_dict['type'] = 'PHOTOLYSIS'
         elif 'ARR' in coeffs:
-            equation_dict = parse_kpp_arrhenius(coeffs)
+            equation_dict = parse_kpp_arrhenius(coeffs,
+                to_si_units=to_si_units, N_reactants=N_reactants)
         else:
             # default to Arrhenius with a single coefficient
             coeffs = coeffs.replace('(', '').replace(')', '')
-            equation_dict['type'] = 'ARRHENIUS' 
+            equation_dict['type'] = 'ARRHENIUS'
             equation_dict['A'] = float(coeffs)
+            if to_si_units:
+                equation_dict['A'] *= (N_Avogadro * 1.0e-6)**(N_reactants - 1)
 
         equation_dict['reactants'] = dict()
         equation_dict['products'] = dict()
@@ -286,6 +305,9 @@ def micm_equation_json(lines):
     parser.add_argument('--mechanism', type=str,
         default='Chapman',
         help='mechanism name')
+    parser.add_argument('--to_si_units', action='store_true',
+        default=False,
+        help='convert units from molecules cm-3 to mol m-3')
     parser.add_argument('--debug', action='store_true',
         help='set logging level to debug')
     args = parser.parse_args()
@@ -324,7 +346,8 @@ def micm_equation_json(lines):
     """
     Generate MICM equations JSON from KPP #EQUATIONS section
     """
-    equations_json = micm_equation_json(sections['#EQUATIONS'])
+    equations_json = micm_equation_json(sections['#EQUATIONS'],
+        to_si_units=args.to_si_units)
 
     """
     Assemble MICM species JSON

diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
@@ -5,3 +5,4 @@ add_subdirectory(unit)
 add_subdirectory(regression)
 add_subdirectory(integration)
 add_subdirectory(tutorial)
+add_subdirectory(kpp)
diff --git a/test/kpp/CMakeLists.txt b/test/kpp/CMakeLists.txt
@@ -0,0 +1,21 @@
+################################################################################
+# Test utilities
+
+include(test_util)
+
+################################################################################
+# Tests
+
+if(ENABLE_JSON)
+  create_standard_test(NAME kpp_to_micm SOURCES test_kpp_to_micm.cpp)
+endif()
+
+################################################################################
+
+################################################################################
+# Copy test data
+
+add_custom_target(copy_kpp_configs ALL ${CMAKE_COMMAND} -E copy_directory
+  ${CMAKE_CURRENT_SOURCE_DIR}/configs ${CMAKE_BINARY_DIR}/configs)
+
+################################################################################
diff --git a/test/kpp/configs/kpp_chapman/reactions.json b/test/kpp/configs/kpp_chapman/reactions.json
@@ -0,0 +1,97 @@
+{
+    "camp-data": [
+        {
+            "name": "Chapman",
+            "type": "MECHANISM",
+            "reactions": [
+                {
+                    "type": "PHOTOLYSIS",
+                    "reactants": {
+                        "O2": {}
+                    },
+                    "products": {
+                        "O": {
+                            "yield": 2.0
+                        }
+                    },
+                    "MUSICA name": "R1"
+                },
+                {
+                    "type": "ARRHENIUS",
+                    "A": 48.28552461368,
+                    "reactants": {
+                        "O": {},
+                        "O2": {}
+                    },
+                    "products": {
+                        "O3": {}
+                    },
+                    "MUSICA name": "R2"
+                },
+                {
+                    "type": "PHOTOLYSIS",
+                    "reactants": {
+                        "O3": {}
+                    },
+                    "products": {
+                        "O": {},
+                        "O2": {}
+                    },
+                    "MUSICA name": "R3"
+                },
+                {
+                    "type": "ARRHENIUS",
+                    "A": 949.089383776,
+                    "reactants": {
+                        "O": {},
+                        "O3": {}
+                    },
+                    "products": {
+                        "O2": {
+                            "yield": 2.0
+                        }
+                    },
+                    "MUSICA name": "R4"
+                },
+                {
+                    "type": "PHOTOLYSIS",
+                    "reactants": {
+                        "O3": {}
+                    },
+                    "products": {
+                        "O1D": {},
+                        "O2": {}
+                    },
+                    "MUSICA name": "R5"
+                },
+                {
+                    "type": "ARRHENIUS",
+                    "A": 42817420.803600006,
+                    "reactants": {
+                        "O1D": {},
+                        "M": {}
+                    },
+                    "products": {
+                        "O": {},
+                        "M": {}
+                    },
+                    "MUSICA name": "R6"
+                },
+                {
+                    "type": "ARRHENIUS",
+                    "A": 72265689.12,
+                    "reactants": {
+                        "O1D": {},
+                        "O3": {}
+                    },
+                    "products": {
+                        "O2": {
+                            "yield": 2.0
+                        }
+                    },
+                    "MUSICA name": "R7"
+                }
+            ]
+        }
+    ]
+}
diff --git a/test/kpp/configs/kpp_chapman/species.json b/test/kpp/configs/kpp_chapman/species.json
@@ -0,0 +1,29 @@
+{
+    "camp-data": [
+        {
+            "name": "M",
+            "type": "CHEM_SPEC",
+            "tracer type": "CONSTANT"
+        },
+        {
+            "name": "O2",
+            "type": "CHEM_SPEC",
+            "tracer type": "CONSTANT"
+        },
+        {
+            "name": "O",
+            "type": "CHEM_SPEC",
+            "absolute tolerance": 1e-12
+        },
+        {
+            "name": "O1D",
+            "type": "CHEM_SPEC",
+            "absolute tolerance": 1e-12
+        },
+        {
+            "name": "O3",
+            "type": "CHEM_SPEC",
+            "absolute tolerance": 1e-12
+        }
+    ]
+}