From a77b83d15606489b549ac869151987002aad924a Mon Sep 17 00:00:00 2001 From: ZhiyuShi <66129026+ZhiyuShi@users.noreply.github.com> Date: Fri, 10 Jul 2020 15:28:07 +0800 Subject: [PATCH 1/2] Update __init__.py --- reactorch/__init__.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/reactorch/__init__.py b/reactorch/__init__.py index 7a49306..9359263 100644 --- a/reactorch/__init__.py +++ b/reactorch/__init__.py @@ -1 +1,3 @@ from .solution import * +from .flame_1d import * +from .fit_transport_data import * From 25a89622f2cc7a5d78b92bb6f69678c449c142a8 Mon Sep 17 00:00:00 2001 From: ZhiyuShi <66129026+ZhiyuShi@users.noreply.github.com> Date: Fri, 10 Jul 2020 15:29:24 +0800 Subject: [PATCH 2/2] Add files via upload --- reactorch/fit_transport_data.py | 100 ++++++++++++++++++++++++++++++++ reactorch/flame_1d.py | 91 +++++++++++++++++++++++++++++ 2 files changed, 191 insertions(+) create mode 100644 reactorch/fit_transport_data.py create mode 100644 reactorch/flame_1d.py diff --git a/reactorch/fit_transport_data.py b/reactorch/fit_transport_data.py new file mode 100644 index 0000000..fed0a8d --- /dev/null +++ b/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) diff --git a/reactorch/flame_1d.py b/reactorch/flame_1d.py new file mode 100644 index 0000000..02c52cb --- /dev/null +++ b/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__ = "weiqiji@mit.edu" +__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) +