<commit_before># Python 2 compatibility
from __future__ import division
from __future__ import print_function

# Standard libraries
import os
import sys
import subprocess
from argparse import ArgumentParser

# More Python 2 compatibility
if sys.version_info.major == 3:
    from itertools import zip
elif sys.version_info.major == 2:
    from itertools import izip as zip

# Related modules
import numpy as np

try:
    import cantera as ct
    from cantera import ck2cti
except ImportError:
    print('Error: Cantera must be installed.')
    raise

# Local imports
import utils
from create_jacobian import create_jacobian

# Compiler based on language
cmd_compile = dict(c='gcc',
                   cuda='nvcc',
                   fortran='gfortran'
                   )

# Flags based on language
flags = dict(c=['-std=c99'],
             cuda=['-arch=sm_20',
                   '-I/usr/local/cuda/include/',
                   '-I/usr/local/cuda/samples/common/inc/',
                   '-dc'],
             fortran='')

libs = dict(c='-lm',
            cuda='-arch=sm_20',
            fortran='')


class ReactorConstPres(object):

    def __init__(self, gas):
        # Parameters of the ODE system and auxiliary data are stored in the
        # ReactorConstPres object.
        self.gas = gas
        self.P = gas.P

    def __call__(self, t=None, y=None):
        """the ODE function, y' = f(t,y) """

        # State vector is [T, Y_1, Y_2, ... Y_K]
        # Only set state if something is input
        if y is not None:
            self.gas.set_unnormalized_mass_fractions(y[1:])
            self.gas.TP = y[0], self.P
        rho = self.gas.density

        wdot = self.gas.net_production_rates
        dTdt = - (np.dot(self.gas.partial_molar_enthalpies, wdot) /
                  (rho * self.gas.cp)
                  )
        dYdt = wdot * self.gas.molecular_weights / rho

        return np.hstack((dTdt, dYdt))


class ReactorConstVol(object):

    def __init__(self, gas):
        # Parameters of the ODE system and auxiliary data are stored in the
        # ReactorConstVol object.
        self.gas = gas
        self.density = gas.density

    def __call__(self, t=None, y=None):
        """the ODE function, y' = f(t,y) """

        # State vector is [T, Y_1, Y_2, ... Y_K]
        # Only set state if something is input
        if y is not None:
            self.gas.set_unnormalized_mass_fractions(y[1:])
            self.gas.TD = y[0], self.density

        wdot = self.gas.net_production_rates
        dTdt = - (np.dot(self.gas.partial_molar_int_energies, wdot) /
                  (self.density * self.gas.cv)
                  )
        dYdt = wdot * self.gas.molecular_weights / self.density

        return np.hstack((dTdt, dYdt))


def convert_mech(mech_filename, therm_filename=None):
    """Convert a mechanism and return a string with the filename.

    Convert a CHEMKIN format mechanism to the Cantera CTI format using
    the Cantera built-in script `ck2cti`.

    :param mech_filename:
        Filename of the input CHEMKIN format mechanism. The converted
        CTI file will have the same name, but with `.cti` extension.
    :param thermo_filename:
        Filename of the thermodynamic database. Optional if the
        thermodynamic database is present in the mechanism input.
    """

    arg = ['--input=' + mech_filename]
    if therm_filename is not None:
        arg.append('--thermo=' + therm_filename)
    arg.append('--permissive')

    # Convert the mechanism
    ck2cti.main(arg)
    mech_filename = mech_filename[:-4] + '.cti'
    print('Mechanism conversion successful, written to '
          '{}'.format(mech_filename)
          )
    return mech_filename

def __write_fd_jacob(file, lang):
    file.write(
    """
    #define FD_ORD 6

    void eval_fd_jacob (const double t, const double pres, double * y, double * jac) {

      double ydot[NN] = {0};
      dydt (t, pres, y, ydot);

      // Finite difference coefficients
      double x_coeffs[FD_ORD];
      double y_coeffs[FD_ORD];

      if (FD_ORD == 2) {
        // 2nd order central difference
        x_coeffs[0] = -1.0;
        x_coeffs[1] = 1.0;
        y_coeffs[0] = -0.5;
        y_coeffs[1] = 0.5;
      } else if (FD_ORD == 4) {
        // 4th order central difference
        x_coeffs[0] = -2.0;
        x_coeffs[1] = -1.0;
        x_coeffs[2] = 1.0;
        x_coeffs[3] = 2.0;
        y_coeffs[0] = 1.0 / 12.0;
        y_coeffs[1] = -2.0 / 3.0;
        y_coeffs[2] = 2.0 / 3.0;
        y_coeffs[3] = -1.0 / 12.0;
      } else {
        // 6th order central difference
        x_coeffs[0] = -3.0;
        x_coeffs[1] = -2.0;
        x_coeffs[2] = -1.0;
        x_coeffs[3] = 1.0;
        x_coeffs[4] = 2.0;
        x_coeffs[5] = 3.0;

        y_coeffs[0] = -1.0 / 60.0;
        y_coeffs[1] = 3.0 / 20.0;
        y_coeffs[2] = -3.0 / 4.0;
        y_coeffs[3] = 3.0 / 4.0;
        y_coeffs[4] = -3.0 / 20.0;
        y_coeffs[5] = 1.0 / 60.0;
      }

      double ewt[NN];
      #pragma unroll
      for (int i = 0; i < NN; ++i) {
        ewt[i] = ATOL + (RTOL * fabs(y[i]));
      }

      // unit roundoff of machine
      double srur = sqrt(DBL_EPSILON);

      double sum = 0.0;
      #pragma unroll
      for (int i = 0; i < NN; ++i) {
        sum += (ewt[i] * ydot[i]) * (ewt[i] * ydot[i]);
      }
      double fac = sqrt(sum / ((double)(NN)));
      double r0 = 1000.0 * 1e-8 * DBL_EPSILON * ((double)(NN)) * fac;

      double ftemp[NN] = {0};

      #pragma unroll
      for (int j = 0; j < NN; ++j) {
        double yj_orig = y[j];
        double r = fmax(srur * fabs(yj_orig), r0 / ewt[j]);

        #pragma unroll
        for (int i = 0; i < NN; ++i) {
          jac[i + NN*j] = 0.0;
        }

        #pragma unroll
        for (int k = 0; k < FD_ORD; ++k) {
          y[j] = yj_orig + x_coeffs[k] * r;
          dydt (t, pres, y, ftemp);    
          #pragma unroll
          for (int i = 0; i < NN; ++i) {
            jac[i + NN*j] += y_coeffs[k] * ftemp[i];
          }
        }

        #pragma unroll
        for (int i = 0; i < NN; ++i) {
          jac[i + NN*j] /= r;
        }

        y[j] = yj_orig;<commit_msg>add atol/rtol to FD jacob<commit_after># Python 2 compatibility
from __future__ import division
from __future__ import print_function

# Standard libraries
import os
import sys
import subprocess
from argparse import ArgumentParser

# More Python 2 compatibility
if sys.version_info.major == 3:
    from itertools import zip
elif sys.version_info.major == 2:
    from itertools import izip as zip

# Related modules
import numpy as np

try:
    import cantera as ct
    from cantera import ck2cti
except ImportError:
    print('Error: Cantera must be installed.')
    raise

# Local imports
import utils
from create_jacobian import create_jacobian

# Compiler based on language
cmd_compile = dict(c='gcc',
                   cuda='nvcc',
                   fortran='gfortran'
                   )

# Flags based on language
flags = dict(c=['-std=c99'],
             cuda=['-arch=sm_20',
                   '-I/usr/local/cuda/include/',
                   '-I/usr/local/cuda/samples/common/inc/',
                   '-dc'],
             fortran='')

libs = dict(c='-lm',
            cuda='-arch=sm_20',
            fortran='')


class ReactorConstPres(object):

    def __init__(self, gas):
        # Parameters of the ODE system and auxiliary data are stored in the
        # ReactorConstPres object.
        self.gas = gas
        self.P = gas.P

    def __call__(self, t=None, y=None):
        """the ODE function, y' = f(t,y) """

        # State vector is [T, Y_1, Y_2, ... Y_K]
        # Only set state if something is input
        if y is not None:
            self.gas.set_unnormalized_mass_fractions(y[1:])
            self.gas.TP = y[0], self.P
        rho = self.gas.density

        wdot = self.gas.net_production_rates
        dTdt = - (np.dot(self.gas.partial_molar_enthalpies, wdot) /
                  (rho * self.gas.cp)
                  )
        dYdt = wdot * self.gas.molecular_weights / rho

        return np.hstack((dTdt, dYdt))


class ReactorConstVol(object):

    def __init__(self, gas):
        # Parameters of the ODE system and auxiliary data are stored in the
        # ReactorConstVol object.
        self.gas = gas
        self.density = gas.density

    def __call__(self, t=None, y=None):
        """the ODE function, y' = f(t,y) """

        # State vector is [T, Y_1, Y_2, ... Y_K]
        # Only set state if something is input
        if y is not None:
            self.gas.set_unnormalized_mass_fractions(y[1:])
            self.gas.TD = y[0], self.density

        wdot = self.gas.net_production_rates
        dTdt = - (np.dot(self.gas.partial_molar_int_energies, wdot) /
                  (self.density * self.gas.cv)
                  )
        dYdt = wdot * self.gas.molecular_weights / self.density

        return np.hstack((dTdt, dYdt))


def convert_mech(mech_filename, therm_filename=None):
    """Convert a mechanism and return a string with the filename.

    Convert a CHEMKIN format mechanism to the Cantera CTI format using
    the Cantera built-in script `ck2cti`.

    :param mech_filename:
        Filename of the input CHEMKIN format mechanism. The converted
        CTI file will have the same name, but with `.cti` extension.
    :param thermo_filename:
        Filename of the thermodynamic database. Optional if the
        thermodynamic database is present in the mechanism input.
    """

    arg = ['--input=' + mech_filename]
    if therm_filename is not None:
        arg.append('--thermo=' + therm_filename)
    arg.append('--permissive')

    # Convert the mechanism
    ck2cti.main(arg)
    mech_filename = mech_filename[:-4] + '.cti'
    print('Mechanism conversion successful, written to '
          '{}'.format(mech_filename)
          )
    return mech_filename

def __write_fd_jacob(file, lang, atol=1e-16, rtol=1e-8):
    file.write(
    """
    #define FD_ORD 6

    void eval_fd_jacob (const double t, const double pres, double * y, double * jac) {

      double ydot[NN] = {0};
      dydt (t, pres, y, ydot);

      // Finite difference coefficients
      double x_coeffs[FD_ORD];
      double y_coeffs[FD_ORD];

      if (FD_ORD == 2) {
        // 2nd order central difference
        x_coeffs[0] = -1.0;
        x_coeffs[1] = 1.0;
        y_coeffs[0] = -0.5;
        y_coeffs[1] = 0.5;
      } else if (FD_ORD == 4) {
        // 4th order central difference
        x_coeffs[0] = -2.0;
        x_coeffs[1] = -1.0;
        x_coeffs[2] = 1.0;
        x_coeffs[3] = 2.0;
        y_coeffs[0] = 1.0 / 12.0;
        y_coeffs[1] = -2.0 / 3.0;
        y_coeffs[2] = 2.0 / 3.0;
        y_coeffs[3] = -1.0 / 12.0;
      } else {
        // 6th order central difference
        x_coeffs[0] = -3.0;
        x_coeffs[1] = -2.0;
        x_coeffs[2] = -1.0;
        x_coeffs[3] = 1.0;
        x_coeffs[4] = 2.0;
        x_coeffs[5] = 3.0;

        y_coeffs[0] = -1.0 / 60.0;
        y_coeffs[1] = 3.0 / 20.0;
        y_coeffs[2] = -3.0 / 4.0;
        y_coeffs[3] = 3.0 / 4.0;
        y_coeffs[4] = -3.0 / 20.0;
        y_coeffs[5] = 1.0 / 60.0;
      }

      double ewt[NN];
      #pragma unroll
      """
    )
    file.write(
    """
      for (int i = 0; i < NN; ++i) {{
        ewt[i] = {} + ({} * fabs(y[i]));
      }}

      // unit roundoff of machine
      double srur = sqrt(DBL_EPSILON);

      double sum = 0.0;
      #pragma unroll
      for (int i = 0; i < NN; ++i) {{
        sum += (ewt[i] * ydot[i]) * (ewt[i] * ydot[i]);
      }}
      double fac = sqrt(sum / ((double)(NN)));
      double r0 = 1000.0 * {} * DBL_EPSILON * ((double)(NN)) * fac;
      """.format(atol, rtol, rtol)
    )
    file.write(
    """
      double ftemp[NN] = {0};

      #pragma unroll
      for (int j = 0; j < NN; ++j) {
        double yj_orig = y[j];
        double r = fmax(srur * fabs(yj_orig), r0 / ewt[j]);

        #pragma unroll
        for (int i = 0; i < NN; ++i) {
          jac[i + NN*j] = 0.0;
        }

        #pragma unroll
        for (int k = 0; k < FD_ORD; ++k) {
          y[j] = yj_orig + x_coeffs[k] * r;
          dydt (t, pres, y, ftemp);    
          #pragma unroll
          for (int i = 0; i < NN; ++i) {
            jac[i + NN*j] += y_coeffs[k] * ftemp[i];
          }
        }

        #pragma unroll
        for (int i = 0; i < NN; ++i) {
          jac[i + NN*j] /= r;
        }

        y[j] = yj_orig;
{'fraction_non_alphanumeric': 0.12916432527895963, 'fraction_numerical': 0.024002568836798587, 'mean_word_length': 1.8307202908429903, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2378753', 'n_tokens_mistral': 4732, 'n_tokens_neox': 4442, 'n_words': 1153}