Merge pull request #191 from NCAR/kpp_to_micm

Kpp to micm

etc/configs/kpp/arr.eqn

@@ -0,0 +1,15 @@
+#EQUATIONS
+
+<1> NO2 + hv = NO + O3P : 				6.69e-1*(SUN/60.0e0);
+<2> O3P + O2 + AIR = O3 : 				ARR_ac(5.68e-34,  -2.80e0);
+<3> O3P + O3 = 2O2 : 					ARR_ab(8.00e-12, 2060.0e0);
+<4> O3P + NO + AIR = NO2 : 				ARR_ac(1.00e-31,  -1.60e0);
+<5> O3P + NO2 = NO : 					ARR_ab(6.50e-12,- 120.0e0);
+<7> O3 + NO = NO2 : 					ARR_ab(1.80e-12, 1370.0e0);
+<8> O3 + NO2 = NO3 : 					ARR_ab(1.40e-13, 2470.0e0);
+<9> NO + NO3 = 2NO2 : 					ARR_ab(1.80e-11, -110.0e0);
+<10> NO + NO + O2 = 2NO2 : 				ARR_ab(3.30e-39, -530.0e0);
+<14> NO2 + NO3 = NO + NO2 : 				ARR_ab(4.50e-14, 1260.0e0);
+<15> NO3 + hv = NO : 					1.59e0*(SUN/60.0e0);
+<16> NO3 + hv = NO2 + O3P : 				1.50e+1*(SUN/60.0e0);
+<17> O3 + hv = O3P : 					3.76e-2*(SUN/60.0e0);

etc/configs/kpp/arr.spc

@@ -0,0 +1,19 @@
+#INCLUDE atoms.kpp
+
+#DEFVAR   
+
+        O3        	= 3O;
+	H2O2		= 2H + 2O;        
+	NO		= N + O;
+	NO2		= N + 2O;
+	NO3		= N + 3O;
+	O3P		= O;
+
+
+
+#DEFFIX	
+	AIR		= IGNORE;
+	O2		= 2O;
+	H2O		= 2H + O;
+	H2		= 2H;	
+	CH4		= C + 4H;	

etc/scripts/kpp_to_micm.py

@@ -13,24 +13,24 @@
 Description:
     kpp_to_micm.py translates KPP config files to MICM JSON config files
     (desginated by the suffixes .kpp, .spc, .eqn, .def)
-    from a single directory specified by the --kpp_dir argument.
+    from a single directory specified by the
+    --kpp_dir and --kpp_name arguments.
 
     In the initial implementation,
     the KPP sections #ATOMS (not yet used), #DEFFIX,
     #DEFVAR, and #EQUATIONS are read and parsed.
     Equations with the hv reactant are MICM PHOTOLYSIS reactions,
-    all others are assumed to be ARRHENIUS reactions.
+    equations with a single coefficient are assumed to be ARRHENIUS reactions.
 
 TODO:
-    (1) Parse both A and C Arrhenius coefficients from KPP equations
-    Currently a single rate coeffient is assigned to A and C is set to 0.
-    (2) Translate stoichiometric coefficients in the equation string
+    (1) Translate stoichiometric coefficients in the equation string
     with more than one digit.
-    (3) Add method unit tests with pytest.
-    (4) Add support for many more reaction types ...
+    (2) Add pytest unit test for method micm_equation_json.
+    (3+) Add support for many more reaction types ...
 
 Revision History:
-    2023/08/03 Initial implementation
+    v1.00 2023/08/03 Initial implementation
+    v1.01 2023/08/16 Added method parse_kpp_arrhenius
 """
 
 import os
@@ -40,10 +40,10 @@
 import json
 from glob import glob
 
-__version__ = 1.0
+__version__ = 'v1.01'
 
 
-def read_kpp_config(kpp_dir):
+def read_kpp_config(kpp_dir, kpp_name):
     """
     Read all KPP config files in a directory
 
@@ -59,7 +59,7 @@ def read_kpp_config(kpp_dir):
     lines = list()
 
     for suffix in suffixes:
-        files = glob(os.path.join(kpp_dir, '*' + suffix))
+        files = glob(os.path.join(kpp_dir, kpp_name + '*' + suffix))
         logging.debug(files)
         for filename in files:
             with open(filename, 'r') as f:
@@ -129,6 +129,69 @@ def micm_species_json(lines, fixed=False, tolerance=1.0e-12):
     return species_json
 
 
+def parse_kpp_arrhenius(kpp_str):
+    """
+    Parse KPP Arrhenius reaction
+
+    Parameters
+        (str) kpp_str: Arrhenius reaction string
+
+    Returns
+        (dict): MICM Arrhenius reaction coefficients
+
+
+    Arrhenius formula from KPP
+    --------------------------
+
+    KPP_REAL ARR_abc( float A0, float B0, float C0 )
+    {
+      double ARR_RES;
+
+      ARR_RES = (double)A0
+        * exp( -(double)B0/TEMP )
+        * pow( (TEMP/300.0), (double)C0 );
+
+    return (KPP_REAL)ARR_RES;
+    }
+
+    Arrhenius formula from MICM
+    ---------------------------
+
+    inline double ArrheniusRateConstant::calculate(
+      const double& temperature, const double& pressure) const
+    {
+    return parameters_.A_ * std::exp(parameters_.C_ / temperature)
+      * pow(temperature / parameters_.D_, parameters_.B_) *
+      (1.0 + parameters_.E_ * pressure);
+    }
+    """
+    logging.debug(kpp_str)
+    coeffs = [float(coeff.replace(' ', '')) for coeff in
+        kpp_str.split('(')[1].split(')')[0].split(',')]
+    logging.debug(coeffs)
+    arr_dict = dict()
+    arr_dict['type'] = 'ARRHENIUS' 
+    # note the interchange of B and C, and change of sign
+    # in the KPP and MICM conventions
+    if ('_abc(' in kpp_str):
+        arr_dict['A'] = coeffs[0]
+        arr_dict['B'] = coeffs[2]
+        arr_dict['C'] = - coeffs[1]
+        arr_dict['D'] = 300.0
+    elif ('_ab(' in kpp_str):
+        arr_dict['A'] = coeffs[0]
+        arr_dict['C'] = - coeffs[1]
+        arr_dict['D'] = 300.0
+    elif ('_ac(' in kpp_str):
+        arr_dict['A'] = coeffs[0]
+        arr_dict['B'] = coeffs[1]
+        arr_dict['D'] = 300.0
+    else:
+        logging.error('unrecognized KPP Arrhenius syntax')
+    logging.debug(arr_dict)
+    return arr_dict
+
+
 def micm_equation_json(lines):
     """
     Generate MICM equation JSON
@@ -143,18 +206,23 @@ def micm_equation_json(lines):
     equations = list() # list of dict
 
     for line in lines:
+        logging.debug(line)
+
+        # split on equal sign into left hand and right hand sides 
         lhs, rhs = tuple(line.split('='))
+
+        # extract reaction coefficients
+        rhs, coeffs = tuple(rhs.split(':'))
+        coeffs = coeffs.replace(';', '')
+
+        # get reactants and products
         reactants = lhs.split('+')
         products = rhs.split('+')
 
         # extract equation label delimited by < >
         label, reactants[0] = tuple(reactants[0].split('>'))
         label = label.lstrip('<')
 
-        # extract equation coefficients delimited by :
-        products[-1], coeffs = tuple(products[-1].split(':'))
-        coeffs = coeffs.replace('(', '').replace(')', '').replace(';', '')
-
         # remove trailing and leading whitespace
         reactants = [reactant.strip().lstrip() for reactant in reactants]
         products = [product.strip().lstrip() for product in products]
@@ -163,12 +231,13 @@ def micm_equation_json(lines):
 
         if 'SUN' in coeffs:
             equation_dict['type'] = 'PHOTOLYSIS' 
+        elif 'ARR' in coeffs:
+            equation_dict = parse_kpp_arrhenius(coeffs)
         else:
+            # default to Arrhenius with a single coefficient
+            coeffs = coeffs.replace('(', '').replace(')', '')
             equation_dict['type'] = 'ARRHENIUS' 
-            # assuming a single coefficient here
             equation_dict['A'] = float(coeffs)
-            # need to generalize to parse both A and C from KPP
-            equation_dict['C'] = 0.0
 
         equation_dict['reactants'] = dict()
         equation_dict['products'] = dict()
@@ -208,6 +277,9 @@ def micm_equation_json(lines):
     parser.add_argument('--kpp_dir', type=str,
         default=os.path.join('..', 'configs', 'kpp'),
         help='KPP input config directory')
+    parser.add_argument('--kpp_name', type=str,
+        default='small_strato',
+        help='KPP config name')
     parser.add_argument('--micm_dir', type=str,
         default=os.path.join('..', 'configs', 'micm'),
         help='MICM output species config file')
@@ -227,7 +299,7 @@ def micm_equation_json(lines):
     """
     Read KPP config files
     """
-    lines = read_kpp_config(args.kpp_dir)
+    lines = read_kpp_config(args.kpp_dir, args.kpp_name)
 
     """
     Split KPP config by section
@@ -277,6 +349,8 @@ def micm_equation_json(lines):
     Write MICM JSON
     """
     micm_mechanism_dir = os.path.join(args.micm_dir, args.mechanism)
+    if not os.path.exists(args.micm_dir):
+        os.mkdir(args.micm_dir)
     if not os.path.exists(micm_mechanism_dir):
         os.mkdir(micm_mechanism_dir)
     with open(os.path.join(micm_mechanism_dir, 'species.json'), 'w') as f:

etc/scripts/test_kpp_to_micm.py

@@ -0,0 +1,40 @@
+"""
+Copyright (C) 2023
+National Center for Atmospheric Research,
+SPDX-License-Identifier: Apache-2.0
+
+File:
+    test_kpp_to_micm.py
+
+Usage:
+    pytest test_kpp_to_micm.py --log-cli-level=DEBUG
+"""
+
+import kpp_to_micm
+
+def test_parse_kpp_arrhenius():
+    """
+    examples
+    ARR_ab(1.0e-12, 2000.0)
+    ARR_ac(1.0e-12, -3.0)
+    ARR_abc(1.0e-12, 2000.0, -3.0)
+    """
+
+    kpp_A, kpp_B, kpp_C = 1.0e-12, 2000.0, -3.0
+
+    arr_dict = kpp_to_micm.parse_kpp_arrhenius(
+        'ARR_ab(%.2e, %.2f)' % (kpp_A, kpp_B))
+    assert arr_dict['A'] == kpp_A
+    assert arr_dict['C'] == - kpp_B
+
+    arr_dict = kpp_to_micm.parse_kpp_arrhenius(
+        'ARR_ac(%.2e, %.2f)' % (kpp_A, kpp_C))
+    assert arr_dict['A'] == kpp_A
+    assert arr_dict['B'] == kpp_C
+
+    arr_dict = kpp_to_micm.parse_kpp_arrhenius(
+        'ARR_abc(%.2e, %.2f, %.2f)' % (kpp_A, kpp_B, kpp_C))
+    assert arr_dict['A'] == kpp_A
+    assert arr_dict['C'] == - kpp_B
+    assert arr_dict['B'] == kpp_C
+