add class Flame_1D and function fit_transport_data

reactorch/__init__.py

@@ -1 +1,3 @@
 from .solution import *
+from .flame_1d import *
+from .fit_transport_data import *

reactorch/fit_transport_data.py

@@ -0,0 +1,100 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Apr 15 21:04:05 2020
+
+@author: weiqi
+"""
+
+import cantera as ct
+import numpy as np
+import matplotlib.pyplot as plt
+import torch
+
+
+poly_order = 6
+
+
+def species_viscosities(gas, T_list, p, comp):
+
+    species_viscosities_list = np.zeros((T_list.shape[0], gas.n_species))
+
+    for i, T in enumerate(T_list):
+        gas.TPX = T, p, comp
+        species_viscosities_list[i, :] = gas.species_viscosities
+
+    poly = np.polyfit(np.log(T_list), species_viscosities_list, deg=poly_order)
+
+    return species_viscosities_list, poly
+
+
+def binary_diff_coeffs(gas, T_list, p, comp):
+    binary_diff_coeffs_list = np.zeros((T_list.shape[0], gas.n_species*gas.n_species))
+
+    for i, T in enumerate(T_list):
+        gas.TPX = T, p, comp
+        binary_diff_coeffs_list[i, :] = gas.binary_diff_coeffs.flatten()
+
+    poly = np.polyfit(np.log(T_list), binary_diff_coeffs_list, deg=poly_order)
+    return binary_diff_coeffs_list, poly
+
+
+def thermal_conductivity(gas, T_list, p, comp):
+    thermal_conductivity_list = np.zeros((T_list.shape[0], gas.n_species))
+
+    arr_thermal_cond = np.zeros(gas.n_species)
+
+    for i, T in enumerate(T_list):
+        gas.TPX = T, p, comp
+        for j in range(gas.n_species):
+            Y = np.zeros(gas.n_species)
+            Y[j] = 1
+            gas.set_unnormalized_mass_fractions(Y)
+            arr_thermal_cond[j] = gas.thermal_conductivity
+
+        thermal_conductivity_list[i, :] = arr_thermal_cond
+
+    poly = np.polyfit(np.log(T_list), thermal_conductivity_list, deg=poly_order)
+    return thermal_conductivity_list, poly
+
+def fit_transport_data(mech_yaml, TPX):
+
+    gas = ct.Solution(mech_yaml)
+
+    T = TPX[0]
+    p = TPX[1]
+    comp = TPX[2]
+
+    n_points = 2700
+    T_min = 300
+    T_max = 3500
+    T_list = np.linspace(start=T_min, stop=T_max, num=n_points, endpoint=True)
+
+    gas.TPX = 1000, p, comp
+
+    species_viscosities_list, species_viscosities_poly = species_viscosities(gas, T_list, p, comp)
+    binary_diff_coeffs_list, binary_diff_coeffs_poly = binary_diff_coeffs(gas, T_list, p, comp)
+    thermal_conductivity_list, thermal_conductivity_poly = thermal_conductivity(gas, T_list, p, comp)
+
+    np.save('mech/transfit/species_viscosities_poly', species_viscosities_poly)
+    np.save('mech/transfit/binary_diff_coeffs_poly', binary_diff_coeffs_poly)
+    np.save('mech/transfit/thermal_conductivity_poly', thermal_conductivity_poly)
+
+    max_error_relative_list = np.zeros(gas.n_species)
+    for i in range(gas.n_species):
+        poly = np.poly1d(species_viscosities_poly[:, i])
+        max_error_relative_list[i] = np.abs(poly(np.log(T_list)) / species_viscosities_list[:, i] - 1).max()
+
+
+    np.set_printoptions(precision=3)
+    print(gas.species_names)
+    print(max_error_relative_list * 100)
+
+
+    max_error_relative_list = np.zeros(gas.n_species)
+    for i in range(gas.n_species):
+        poly = np.poly1d(thermal_conductivity_poly[:, i])
+        max_error_relative_list[i] = np.abs(poly(np.log(T_list)) / (thermal_conductivity_list[:, i] + 1e-12) - 1).max()
+
+    print(gas.species_names)
+    print(max_error_relative_list * 100)

reactorch/flame_1d.py

@@ -0,0 +1,91 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+__author__ = "Weiqi Ji"
+__copyright__ = "Copyright 2020, DENG"
+
+__version__ = "0.1"
+__email__ = "[email protected]"
+__status__ = "Development"
+
+import torch
+
+
+torch.set_default_tensor_type("torch.DoubleTensor")
+
+
+class Flame_1d:
+    def __init__(self, tct, device=None):
+        super().__init__()
+
+        if device is None:
+            self.device = torch.device('cpu')
+        else:
+            self.device = device
+
+        self.tct = tct
+
+    def set_states(self, TPY):
+
+        self.tct.set_states(TPY)
+
+        self.update_u()
+        self.update_diffusive_fluxes()
+        self.update_dTdx()
+        self.update_dYdx()
+        self.update_rhs()
+
+    def set_x(self, x):
+        self.x = x
+
+    def set_mdot(self, m_dot):
+        self.m_dot = torch.Tensor([m_dot]).to(self.device)
+
+    def update_u(self):
+        self.u = self.m_dot / self.tct.density_mass
+
+    def grad_to_x(self, output):
+
+        if output.shape[1] == 1:
+            dydt = torch.autograd.grad(outputs=output.sum(),
+                                       inputs=self.x,
+                                       retain_graph=True,
+                                       create_graph=True,
+                                       allow_unused=True)[0]
+
+        else:
+            dydt = torch.zeros_like(output)
+            for i in range(output.shape[1]):
+                dydt[:, i] = torch.autograd.grad(outputs=output[:, i].sum(),
+                                                 inputs=self.x,
+                                                 retain_graph=True,
+                                                 create_graph=True,
+                                                 allow_unused=True)[0].view(-1)
+
+        return dydt
+
+    def update_diffusive_fluxes(self):
+        # following https://cantera.org/science/flames.html
+        self.tct.update_transport()
+        self.diffusive_fluxes_star = self.tct.density_mass * \
+            (self.tct.molecular_weights.T / self.tct.mean_molecular_weight) \
+            * self.tct.mix_diff_coeffs * self.grad_to_x(self.tct.X)
+        self.diffusive_fluxes = self.diffusive_fluxes_star - \
+            self.tct.Y.clone() * self.diffusive_fluxes_star.sum(dim=1, keepdim=True)
+
+    def update_dYdx(self):
+        # following https://cantera.org/science/flames.html
+        self.dYdx = - (self.grad_to_x(self.diffusive_fluxes) + self.tct.molecular_weights.T * self.tct.wdot) / self.m_dot
+
+    def update_dTdx(self):
+        # following https://cantera.org/science/flames.html
+        self.dT2dx = self.grad_to_x(self.tct.T)
+
+        self.dTdx0 = self.grad_to_x(self.dT2dx * self.tct.thermal_conductivity)
+        self.dTdx1 = self.dTdx0 - self.dT2dx * (self.diffusive_fluxes * self.tct.cp).sum(dim=1, keepdim=True)
+        self.dTdx2 = self.dTdx1 - (self.tct.h * self.tct.molecular_weights.T * self.tct.wdot).sum(dim=1, keepdim=True)
+        self.dTdx = self.dTdx2/self.m_dot / self.tct.cp_mole.unsqueeze(-1)
+
+    def update_rhs(self):
+        self.rhs = torch.cat((self.dTdx, self.dYdx), dim=1)
+