[ck2yaml] Add support for format with single intermediate temperature

This format, used in the original NASA thermo database by McBride et al stores the molecular weight in the last section of the first line of each thermo entry. Also, correct handling of cases where the inferred intermediate temperature is outside the range (Tmin, Tmax) and the species only has a single set of non-zero coefficients.
2025-02-25 18:55:29 -06:00 · 2022-04-20 11:38:11 -04:00
parent a41522c84e
commit 05ed4b3d64
3 changed files with 111 additions and 30 deletions
--- a/interfaces/cython/cantera/ck2yaml.py
+++ b/interfaces/cython/cantera/ck2yaml.py
@@ -15,6 +15,7 @@ Usage:
            [--name=<name>]
            [--extra=<filename>]
            [--output=<filename>]
+            [--single-intermediate-temperature]
            [--permissive]
            [--quiet]
            [--no-validate]
@@ -34,6 +35,10 @@ For the case of a surface mechanism, the gas phase input file should be
 specified as 'input' and the surface phase input file should be specified as
 'surface'.

+The '--single-intermediate-temperature' option should be used with thermo data where
+only a single break temperature is used and the last value in the first line of each
+species thermo entry is the molecular weight instead.
+
 The '--permissive' option allows certain recoverable parsing errors (e.g.
 duplicate transport data) to be ignored. The '--name=<name>' option
 is used to override default phase names (i.e. 'gas').
@@ -217,8 +222,12 @@ class Nasa7:
    @classmethod
    def to_yaml(cls, representer, node):
        out = BlockMap([('model', 'NASA7')])
-        out['temperature-ranges'] = FlowList([node.Tmin, node.Tmid, node.Tmax])
-        out['data'] = [FlowList(node.low_coeffs), FlowList(node.high_coeffs)]
+        if node.Tmid is not None:
+            out['temperature-ranges'] = FlowList([node.Tmin, node.Tmid, node.Tmax])
+            out['data'] = [FlowList(node.low_coeffs), FlowList(node.high_coeffs)]
+        else:
+            out['temperature-ranges'] = FlowList([node.Tmin, node.Tmax])
+            out['data'] = [FlowList(node.low_coeffs)]
        if node.note:
            note = textwrap.dedent(node.note.rstrip())
            if '\n' in note:
@@ -761,6 +770,7 @@ class Parser:
        self.quantity_units = 'mol'  # for the current REACTIONS section
        self.output_quantity_units = 'mol'  # for the output file
        self.motz_wise = None
+        self.single_intermediate_temperature = False
        self.warning_as_error = True

        self.elements = []
@@ -901,10 +911,16 @@ class Parser:
        # Remember that the high-T polynomial comes first!
        Tmin = fortFloat(lines[0][45:55])
        Tmax = fortFloat(lines[0][55:65])
-        try:
-            Tint = fortFloat(lines[0][65:75])
-        except ValueError:
-            Tint = TintDefault
+        if self.single_intermediate_temperature:
+            # Intermediate temperature is shared across all species, except if the
+            # species only has one temperature range
+            Tint = TintDefault if Tmin < TintDefault < Tmax else None
+        else:
+            # Non-default intermediate temperature can be provided
+            try:
+                Tint = fortFloat(lines[0][65:75])
+            except ValueError:
+                Tint = TintDefault if Tmin < TintDefault < Tmax else None

        high_coeffs = [fortFloat(lines[i][j:k])
                       for i,j,k in [(1,0,15), (1,15,30), (1,30,45), (1,45,60),
@@ -913,16 +929,31 @@ class Parser:
                       for i,j,k in [(2,30,45), (2,45,60), (2,60,75), (3,0,15),
                                     (3,15,30), (3,30,45), (3,45,60)]]

-        # Duplicate the valid set of coefficients if only one range is provided
-        if all(c == 0 for c in low_coeffs) and Tmin == Tint:
-            low_coeffs = high_coeffs
-        elif all(c == 0 for c in high_coeffs) and Tmax == Tint:
-            high_coeffs = low_coeffs
+        # Cases where only one temperature range is needed
+        if Tint == Tmin or Tint == Tmax or high_coeffs == low_coeffs:
+            Tint = None

-        # Construct and return the thermodynamics model
-        thermo = Nasa7(Tmin=Tmin, Tmax=Tmax, Tmid=Tint,
-                           low_coeffs=low_coeffs, high_coeffs=high_coeffs,
-                           note=note)
+        # Duplicate the valid set of coefficients if only one range is provided
+        if Tint is None:
+            if all(c == 0 for c in low_coeffs):
+                # Use the first set of coefficients if the second is all zeros
+                coeffs = high_coeffs
+            elif all(c == 0 for c in high_coeffs):
+                # Use the second set of coefficients if the first is all zeros
+                coeffs = low_coeffs
+            elif high_coeffs == low_coeffs:
+                # If the coefficients are duplicated, that's fine too
+                coeffs = low_coeffs
+            else:
+                raise InputError(
+                    "Only one temperature range defined but two distinct sets of "
+                    "coefficients given for species thermo entry:\n{}\n",
+                    "".join(lines))
+            thermo = Nasa7(Tmin=Tmin, Tmax=Tmax, Tmid=None, low_coeffs=coeffs,
+                           high_coeffs=None, note=note)
+        else:
+            thermo = Nasa7(Tmin=Tmin, Tmax=Tmax, Tmid=Tint, low_coeffs=low_coeffs,
+                           high_coeffs=high_coeffs, note=note)

        return species, thermo, composition

@@ -2007,9 +2038,11 @@ class Parser:
    @staticmethod
    def convert_mech(input_file, thermo_file=None, transport_file=None,
                     surface_file=None, phase_name='gas', extra_file=None,
-                     out_name=None, quiet=False, permissive=None):
+                     out_name=None, single_intermediate_temperature=False, quiet=False,
+                     permissive=None):

        parser = Parser()
+        parser.single_intermediate_temperature = single_intermediate_temperature
        if quiet:
            logger.setLevel(level=logging.ERROR)
        else:
@@ -2126,18 +2159,19 @@ class Parser:

 def convert_mech(input_file, thermo_file=None, transport_file=None,
                 surface_file=None, phase_name='gas', extra_file=None,
-                 out_name=None, quiet=False, permissive=None):
+                 out_name=None, single_intermediate_temperature=False, quiet=False,
+                 permissive=None):
    _, surface_names = Parser.convert_mech(
        input_file, thermo_file, transport_file, surface_file, phase_name,
-        extra_file, out_name, quiet, permissive)
+        extra_file, out_name, single_intermediate_temperature, quiet, permissive)
    return surface_names


 def main(argv):

    longOptions = ['input=', 'thermo=', 'transport=', 'surface=', 'name=',
-                   'extra=', 'output=', 'permissive', 'help', 'debug', 'quiet',
-                   'no-validate', 'id=']
+                   'extra=', 'output=', 'permissive', 'help', 'debug',
+                   'single-intermediate-temperature', 'quiet', 'no-validate', 'id=']

    try:
        optlist, args = getopt.getopt(argv, 'dh', longOptions)
@@ -2161,6 +2195,7 @@ def main(argv):

    input_file = options.get('--input')
    thermo_file = options.get('--thermo')
+    single_intermediate_temperature = '--single-intermediate-temperature' in options
    permissive = '--permissive' in options
    quiet = '--quiet' in options
    transport_file = options.get('--transport')
@@ -2191,7 +2226,7 @@ def main(argv):

    parser, surfaces = Parser.convert_mech(input_file, thermo_file,
            transport_file, surface_file, phase_name, extra_file, out_name,
-            quiet, permissive)
+            single_intermediate_temperature, quiet, permissive)

    # Do full validation by importing the resulting mechanism
    if not input_file:
--- a/interfaces/cython/cantera/test/test_convert.py
+++ b/interfaces/cython/cantera/test/test_convert.py
@@ -12,8 +12,7 @@ from cantera import ck2cti, ck2yaml, cti2yaml, ctml2yaml

 class converterTestCommon:
    def convert(self, inputFile, thermo=None, transport=None,
-                surface=None, output=None, extra=None,
-                quiet=True, permissive=None):
+                surface=None, output=None, extra=None, **kwargs):
        if output is None:
            output = Path(inputFile).stem  # strip '.inp'
        if inputFile is not None:
@@ -31,8 +30,7 @@ class converterTestCommon:
        if output.is_file():
            output.unlink()
        self._convert(inputFile, thermo=thermo, transport=transport,
-            surface=surface, output=output, extra=extra,
-            quiet=quiet, permissive=permissive)
+            surface=surface, output=output, extra=extra, **kwargs)
        return output

    def checkConversion(self, refFile, testFile):
@@ -482,10 +480,9 @@ class ck2ctiTest(converterTestCommon, utilities.CanteraTest):
        ct.make_deprecation_warnings_fatal()

    def _convert(self, inputFile, *, thermo, transport, surface, output, extra,
-                 quiet, permissive):
+                 **kwargs):
        ck2cti.convertMech(inputFile, thermoFile=thermo,
-            transportFile=transport, surfaceFile=surface, outName=output,
-            quiet=quiet, permissive=permissive)
+            transportFile=transport, surfaceFile=surface, outName=output, **kwargs)

    def test_missingElement(self):
        with self.assertRaisesRegex(self.InputError, 'Undefined elements'):
@@ -497,10 +494,10 @@ class ck2yamlTest(converterTestCommon, utilities.CanteraTest):
    InputError = ck2yaml.InputError

    def _convert(self, inputFile, *, thermo, transport, surface, output, extra,
-                 quiet, permissive):
+                 **kwargs):
        ck2yaml.convert_mech(inputFile, thermo_file=thermo,
            transport_file=transport, surface_file=surface, out_name=output,
-            extra_file=extra, quiet=quiet, permissive=permissive)
+            extra_file=extra, **kwargs)

    @utilities.slow_test
    def test_extra(self):
@@ -555,6 +552,36 @@ class ck2yamlTest(converterTestCommon, utilities.CanteraTest):
        for token in ('FAILED', 'lines 12 and 14', 'R1A', 'R1B'):
            self.assertIn(token, message)

+    def test_single_Tint(self):
+        output = self.convert(None, thermo="thermo_single_Tint.dat",
+                              output="thermo_single_Tint",
+                              single_intermediate_temperature=True)
+        mech = utilities.load_yaml(output)
+
+        # Al(cr)
+        thermo = mech["species"][0]["thermo"]
+        assert thermo["temperature-ranges"] == [200.0, 933.61]
+        assert len(thermo["data"]) == 1
+        assert thermo["data"][0][0] == 1.01040191
+
+        # AlBr3(L)
+        thermo = mech["species"][1]["thermo"]
+        assert thermo["temperature-ranges"] == [370.6, 5000.0]
+        assert len(thermo["data"]) == 1
+        assert thermo["data"][0][0] == 15.02975
+
+        # AlF3(b)
+        thermo = mech["species"][2]["thermo"]
+        assert thermo["temperature-ranges"] == [728.0, 1000.0, 2523.0]
+        assert len(thermo["data"]) == 2
+        assert thermo["data"][1][0] == 10.41947
+
+        # AlF3(L)
+        thermo = mech["species"][3]["thermo"]
+        assert thermo["temperature-ranges"] == [2523.0, 5000.0]
+        assert len(thermo["data"]) == 1
+        assert thermo["data"][0][0] == 15.096679
+
    def test_error_for_big_element_number(self):
        with self.assertRaisesRegex(self.InputError, 'Element amounts can have no more than 3 digits.'):
            self.convert('big_element_num_err.inp')
--- a/test/data/thermo_single_Tint.dat
+++ b/test/data/thermo_single_Tint.dat
@@ -0,0 +1,19 @@
+THERMO
+ 300.0   1000.0   5000.0
+AL(cr)            CODA89AL 1.   0.   0.   0.C   200.000   933.610     26.98154 1
+ 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00    2
+ 0.00000000E+00 0.00000000E+00 1.01040191E+00 1.20769743E-02-2.62083556E-05    3
+ 2.64282413E-08-9.01916513E-12-6.54454196E+02-5.00471254E+00 0.00000000E+00    4
+ALBr3(L)          J 9/79AL 1.BR 3.   0.   0.C   370.600  5000.000    266.69354 1
+ 1.50297500E+01 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00    2
+-6.47837290E+04-6.07991010E+01 1.50297500E+01 0.00000000E+00 0.00000000E+00    3
+ 0.00000000E+00 0.00000000E+00-6.47837290E+04-6.07991010E+01 0.00000000E+00    4
+ALF3(b)           J 9/79AL 1.F  3.   0.   0.C   728.000  2523.000     83.97675 1
+ 1.04194700E+01 2.33765010E-03-8.80830770E-07 2.85578830E-10-3.46072630E-14    2
+-1.84922050E+05-5.23714020E+01 9.50345050E+00 5.13025090E-03-3.71167640E-06    3
+ 1.20523570E-09 0.00000000E+00-1.84695550E+05-4.77361470E+01 0.00000000E+00    4
+ALF3(L)           J 9/79AL 1.F  3.   0.   0.C  2523.000  5000.000     83.97675 1
+ 1.50966790E+01 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00    2
+-1.79986860E+05-8.00491030E+01 1.50966790E+01 0.00000000E+00 0.00000000E+00    3
+ 0.00000000E+00 0.00000000E+00-1.79986860E+05-8.00491030E+01 0.00000000E+00    4
+END