Modeling HyChem

[jiweiqi]

Note down my questions about modeling HyChem

• during fitting the pyrolysis steps, do we need to consider the core mechanism? how much does the core mechanism participate in the breakdown of large hydrocarbon molecular?
• how do we hybrid the neural networks and the core mechanism? we probably have to bring in the package of RMS.jl. I am not sure how much the workload will be. Maybe we should collaborate with the RMS.jl team.

[jiweiqi]

A quick working around solution is to disable the core mechanism when testing CRNN.

[jiweiqi]

Useful tips

• start from 2 experiments for good initialization
• initial small lr=0.001 and then 0.005 (warm up)
  

[jiweiqi]

example reaction pathflux analysis

It seems essential to couple with the core mechanism.

[jiweiqi]

Code for path flux analysis

import sys, os
import numpy as np
import matplotlib.pyplot as plt
from skopt.space import Space
from skopt.sampler import Lhs

import cantera as ct

np.random.seed(1)

gas = ct.Solution('../../mech/JP10skeletal.cti')
i_var = [gas.species_index(s) for s in [
    "C10H16", "H", "H2", "CH3", "CH4", "aC3H5", "C2H4", "C3H6", "C5H6", "C6H6", "C6H5CH3", "N2",
    "OH", "H2O", "O2", "HO2", "O", "H2O2"]]

# for i in range(3,232):
#     gas.set_multiplier(i, 0)

lhs = Lhs(lhs_type="classic", criterion=None)
space = Space([(1100., 1400.), (1., 5.), ('0', '1', '2', '3', '4')])

nsamples = 10
x = lhs.generate(space.dimensions, nsamples)

comp = ['C10H16:0.01,n2:0.99',
        'C10H16:0.02,n2:0.98',
        'C10H16:0.01,o2:0.014,h2o:0.1,n2:0.9',
        'C10H16:0.01,o2:0.028,h2o:0.1,n2:0.9',
        'C10H16:0.01,o2:0.007,h2o:0.1,n2:0.9']

for i in range(1):

    gas.TPX = x[i][0], x[i][1] * ct.one_atm, comp[np.int16(x[i][2])]
    
    Y_fuel_0 = gas.X[gas.species_index('C10H16')]

    r = ct.IdealGasConstPressureReactor(gas)
    sim = ct.ReactorNet([r])
    time = 0.0
    states = ct.SolutionArray(gas, extra=['t'])
    
    print(x[i][0], x[i][1], comp)

    print('%10s %10s %10s %14s' % ('t [s]', 'T [K]', 'P [Pa]', 'u [J/kg]'))
    for n in range(10000):
        states.append(r.thermo.state, t=sim.time)
        # time += 1.e-4
        # sim.advance(time)
        sim.step()

        if r.thermo.X[gas.species_index('C10H16')] < Y_fuel_0 * 0.75:
            break

        # print('%10.3e %10.3f %10.3f %14.6e' % (sim.time, r.T,
        #                                        r.thermo.P, r.thermo.u))
        
    element = 'H'

    diagram = ct.ReactionPathDiagram(gas, element)
    diagram.title = 'Reaction path diagram following {} @ {:.1f} K {:.1f} atm {}'.format(element, x[i][0], x[i][1], comp[np.int16(x[i][2])])
    diagram.label_threshold = 0.05
    diagram.show_details = True
    
    dot_file = 'rxnpath.dot'
    img_file = 'rxnpath.png'
    img_path = os.path.join(os.getcwd(), img_file)
    
    diagram.write_dot(dot_file)
    # print(diagram.get_data())
    
    print("Wrote graphviz input file to '{0}'.".format(os.path.join(os.getcwd(), dot_file)))
    # dot .\rxnpath.dot -Tpng -o rxnpath.png -Gdpi=300
    
    os.system('dot {0} -Tpng -o{1} -Gdpi=200'.format(dot_file, img_file))
    print("Wrote graphviz output file to '{0}'.".format(img_path))

[jiweiqi]

Good hyper-parameters

include("header.jl")

is_restart = true;
n_epoch = 100000;
ntotal = 20;
batch_size = 4;
n_plot = 5;
grad_max = 10.0^(1);
maxiters = 1000;
n_exp = 2;
nr = 3;
opt = ADAMW(0.01, (0.9, 0.999), 1.e-6);
atol = 1.e-8;
rtol = 1.e-3;
lb = atol;
ub = 1.0;

include("dataset.jl")
include("network.jl")
include("callback.jl")

# opt = ADAMW(5.e-3, (0.9, 0.999), 1.f-6)

epochs = ProgressBar(iter:n_epoch);
loss_epoch = zeros(Float64, n_exp);
grad_norm = zeros(Float64, n_exp);
for epoch in epochs
    global p
    for i_exp in randperm(n_exp)
        sample = rand(batch_size:ntotal)
        # _, ind = loss_n_ode(p, i_exp, sample; get_ind = true)
        grad = ForwardDiff.gradient(x -> loss_n_ode(x, i_exp, sample), p)
        grad_norm[i_exp] = norm(grad, 2)
        if grad_norm[i_exp] > grad_max
            grad = grad ./ grad_norm[i_exp] .* grad_max
        end
        update!(opt, p, grad)
    end
    for i_exp = 1:n_exp
        loss_epoch[i_exp] = loss_n_ode(p, i_exp)
    end
    _loss = mean(loss_epoch)
    _gnorm = mean(grad_norm)
    set_description(
        epochs,
        string(
            @sprintf("Loss %.4e grad %.2e lr %.2e", _loss, _gnorm, opt[1].eta)
        ),
    )
    cb(p, _loss, _gnorm)
end

[jiweiqi]

varnames = ["C10H16", "H", "H2", "CH3", "CH4", "aC3H5", 
            "C2H4", "C3H6", "C5H6", "C6H6", "C6H5CH3"];

ns = size(varnames)[1];

gas = CreateSolution("./mechanism/JP10skeletal.yaml")
const MW = gas.MW

ind_crnn = []
for var in varnames
    push!(ind_crnn, species_index(gas, var))
end

ind_obs = []
ind_obs_data = []
varnames_obs = ["C10H16", "H2", "CH4", "C2H4", "C3H6", "C5H6", "C6H6", "C6H5CH3"]
# varnames = varnames
for var in varnames_obs
    push!(ind_obs, species_index(gas, var))
    push!(ind_obs_data, findfirst(varnames .== var))
end

E_C = [10, 0, 0, 1, 1, 3, 2, 3, 5, 6, 7];
E_H = [16, 1, 2, 3, 4, 5, 4, 6, 6, 6, 8];
E_O = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0];
E_ = hcat(E_C, E_H)[1:ns, :];
E_null = nullspace(E_')';
nse = size(E_null)[1];

yall = zeros(n_exp, ns + 3, ntotal);
for i_exp in 1:n_exp
    rawdata = readdlm("data/data_$i_exp")'
    ts = rawdata[1, :]
    tend = ts[end]

    _ts = 10 .^ range(log10(tend/1e4), log10(tend/1.01), length=ntotal);
    _ts[1] = 0.0;

    _ylabel = zeros(ns+2, ntotal);
    for i in 1:ns + 2
        itp = LinearInterpolation(ts, rawdata[i+1, :])
        _ylabel[i, :] .= itp.(_ts)
    end

    yall[i_exp, 1, :] .= _ts
    yall[i_exp, 2:end, :] .= _ylabel
end
yalls = yall[:, 4:end, :];
# yscale = (maximum(yalls, dims=[1, 3]) - minimum(yalls, dims=[1, 3]) .+ lb)[1, :, 1];

println("include ", @__FILE__)

np = nr * (nse + 3) + 1;
p = randn(Float64, np) .* 0.01 .- 0.005;
p[end] = 0.1;  # slope
function p2vec(p)
    slope = p[end] .* 10.0
    w_b = p[1:nr] .* slope
    w_in_b = p[nr + 1:nr * 2]
    w_in_Ea = p[nr * 2 + 1:nr * 3] .* slope
    w_out = E_null' * reshape(p[nr * 3 + 1:nr * (nse + 3)], nse, nr)
    w_in = vcat(clamp.(-w_out, 0.0, 4.0), w_in_Ea', w_in_b')
    return w_in, w_b, w_out
end

function display_p(p)
    w_in, w_b, w_out = p2vec(p)
    println("\n species (column) reaction (row)")
    println("w_in | Ea | b | logA")
    show(stdout, "text/plain", round.(hcat(w_in', w_b), digits=2))
    println("\n w_out")
    show(stdout, "text/plain", round.(w_out', digits=2))
    println("\n")
end

function crnn!(du, u, p, t)
    T::typeof(u[1]) = itpT(t)
    P = itpP(t)
    Y = clamp.(u, -lb, 1.0)
    mean_MW = 1. / dot(Y, 1 ./ MW)
    ρ_mass = P / R / T * mean_MW
    X = Y2X(gas, Y, mean_MW)
    C = Y2C(gas, Y, ρ_mass)
    cp_mole, cp_mass = get_cp(gas, T, X, mean_MW)
    h_mole = get_H(gas, T, Y, X)
    S0 = get_S(gas, T, P, X)
    wdot = wdot_func(gas.reaction, T, C, S0, h_mole)
    crnn_in = vcat(log.(clamp.(@view(C[ind_crnn]), 1.e-12, 10.0)),
                   - 1.0 / T, log(T))
    wdot[ind_crnn] .+= w_out * exp.(w_in' * crnn_in + w_b)
    @. du = wdot * MW / ρ_mass
end

tspan = [0, 0.01];
u0 = zeros(gas.n_species);
prob = ODEProblem(crnn!, u0, tspan)

ode_solver = TRBDF2();
sense = ForwardDiffSensitivity()
function predict_n_ode(p, i_exp, sample=ntotal; dense=false)
    global w_in, w_b, w_out = p2vec(p)
    ylabel = @view(yall[i_exp, :, :])
    ts = @view(ylabel[1, :])
    global itpT = LinearInterpolation(ts, @view(ylabel[2, :]); extrapolation_bc=Flat())
    global itpP = LinearInterpolation(ts, @view(ylabel[3, :]); extrapolation_bc=Flat())
    u0 = zeros(gas.n_species)
    u0[ind_crnn] .= @view(ylabel[4:end, 1])
    _prob = remake(prob, u0=u0, p=p, tspan=[0, ts[sample]])
    if dense
        _ts = []
    else
        _ts = @view(ts[1:sample])
    end
    sol = solve(_prob, ode_solver, saveat=_ts,
                atol=atol, rtol=rtol, sensalg=sense, maxiters=maxiters)
    pred = Array(sol)
end
@time pred = predict_n_ode(p, 1; dense=false);

function loss_n_ode(p, i_exp, sample=ntotal; get_ind=false)
    pred = clamp.(predict_n_ode(p, i_exp, sample), -ub, ub)[ind_obs, :]
    ind = size(pred)[2]
    ylabel = clamp.(@view(yall[i_exp, ind_obs_data .+ 3, 1:ind]), -ub, ub)
    yscale_ = maximum(ylabel, dims=2) - minimum(ylabel, dims=2) .+ 1.e-6
    loss = mae(pred ./ yscale_, ylabel ./ yscale_)

    if get_ind
        return loss, ind
    else
        return loss
    end
end
@time loss_n_ode(p, 1; get_ind=true)

# @time grad = ForwardDiff.gradient(x -> loss_n_ode(x, 1), p)

[jiweiqi]

Sampling conditions
No / T / P / Composition

0 1226.1335674919455 4.557842665401539 C10H16:0.01,o2:0.014,h2o:0.1,n2:0.9
1 1165.5878063413302 3.479240595633952 C10H16:0.02,n2:0.98
2 1112.5106601410773 3.9253696712636974 C10H16:0.01,o2:0.014,h2o:0.1,n2:0.9
3 1366.2916745688813 1.9382242908172191 C10H16:0.01,o2:0.028,h2o:0.1,n2:0.9
4 1371.1716434969865 4.667932167825827 C10H16:0.01,n2:0.99
5 1139.0699771789552 2.367677805761318 C10H16:0.01,n2:0.99
6 1270.114025305352 3.166921920946851 C10H16:0.01,o2:0.028,h2o:0.1,n2:0.9
7 1284.211608157857 1.2881297973768633 C10H16:0.01,o2:0.028,h2o:0.1,n2:0.9
8 1206.1645020201008 2.951246974556378 C10H16:0.02,n2:0.98
9 1339.047847271582 1.4587023563268453 C10H16:0.01,o2:0.007,h2o:0.1,n2:0.9

[jiweiqi]

Working version

@inbounds function dudt!(du, u, p, t)
    Y = @view(u[1:ns])
    T = u[end]
    mean_MW = 1.0 / dot(Y, 1 ./ gas.MW)
    ρ_mass = P / R / T * mean_MW
    X = Y2X(gas, Y, mean_MW)
    C = Y2C(gas, Y, ρ_mass)
    cp_mole, cp_mass = get_cp(gas, T, X, mean_MW)
    h_mole = get_H(gas, T, Y, X)
    S0 = get_S(gas, T, P, X)
    # _p = reshape(p, nr, 3)
    # kp = @. @views(exp(_p[:, 1] + _p[:, 2] * log(T) - _p[:, 3] * 4184.0 / R / T))
    kp = exp.(p)
    qdot = wdot_func(gas.reaction, T, C, S0, h_mole; get_qdot = true) .* kp
    qdot[[2,3,4,6,7]] .*= 0.0
    wdot = gas.reaction.vk * qdot
    Ydot = wdot / ρ_mass .* gas.MW
    Tdot = -dot(h_mole, wdot) / ρ_mass / cp_mass
    du .= vcat(Ydot, Tdot)
    return du
end


@inbounds function dudt!(du, u, p, t)
    w_in, w_b, w_out = p2vec(p)
    Y = clamp.(@view(u[1:ns]), -1.e-3, 1.0)
    T = clamp(u[end], 600.0, 3500.0) * ones(eltype(p), 1)[1]
    mean_MW = 1.0 / dot(Y, 1 ./ gas.MW)
    ρ_mass = P / R / T * mean_MW
    X = Y2X(gas, Y, mean_MW)
    C = Y2C(gas, Y, ρ_mass)
    cp_mole, cp_mass = get_cp(gas, T, X, mean_MW)
    h_mole = get_H(gas, T, Y, X)
    S0 = get_S(gas, T, P, X)
    qdot = wdot_func(gas.reaction, T, C, S0, h_mole; get_qdot=true)
    qdot[2:7] .*= 0.0
    wdot = gas.reaction.vk * qdot
    crnn_in = vcat(log.(clamp.(@view(C[ind_crnn]), lb, 1.0)), -1.0 / T, log(T))
    wdot[ind_crnn] .+= w_out * exp.(w_in' * crnn_in + w_b)
    Ydot = wdot / ρ_mass .* gas.MW
    Tdot = -dot(h_mole, wdot) / ρ_mass / cp_mass
    du .= vcat(Ydot, Tdot)
    return du
end


np = crnn_nr * (nse + 3) + 1;
p = randn(Float64, np) .* 0.1 .- 0.05;
p[1:crnn_nr] .+= 1.0;  #lnA
p[crnn_nr * 2 + 1:crnn_nr * 3] .+= 0.0;  # Ea
p[end] = 0.1;  # slope
@inbounds function p2vec(p)
    slope = p[end] .* 10.0
    w_b = @view(p[1:crnn_nr]) .* slope
    w_in_b = @view(p[crnn_nr + 1:crnn_nr * 2])
    w_in_Ea = @view(p[crnn_nr * 2 + 1:crnn_nr * 3]) .* slope
    w_out = E_null' *
            reshape(@view(p[crnn_nr * 3 + 1:crnn_nr * (nse + 3)]),
                    nse, crnn_nr)
    # w_out = reshape(@view(p[crnn_nr * 3 + 1:crnn_nr * (crnn_ns + 3)]),
    #         crnn_ns, crnn_nr)
    # w_in = reshape(@view(p[crnn_nr * (crnn_ns + 3)+1:crnn_nr * (crnn_ns*2 + 3)]),
    #         crnn_ns, crnn_nr)
    # w_out = @. -w_in * (10 ^ w_out)
    w_in = vcat(clamp.(-w_out, 0.0, 2.5), w_in_Ea', w_in_b')
    return w_in, w_b, w_out
end

[jiweiqi]

[jiweiqi]

- equation: C10H16 => 1.79 H + 0.01 CH3 + 0.2 aC3H5 + 0.866256 C2H4 + 0.0693
    C3H6 + 0.952882 C5H6 + 0.407143 H2 + 0.351373 C6H6 + 0.08242 C6H5CH3  # Reaction 1
  rate-constant: {A: 3.2e+15, b: 0.0, Ea: 7.25e+04}
- equation: C10H16 + H => H2 + 0.43 H + 0.26 CH3 + 0.31 aC3H5 + 0.781607
    C2H4 + 0.062529 C3H6 + 0.859765 C5H6 + 0.367353 H2 + 0.361213 C6H6 +
    0.084728 C6H5CH3  # Reaction 2
  rate-constant: {A: 1.83e+04, b: 2.765, Ea: 2551.0}
- equation: C10H16 + CH3 => CH4 + 0.43 H + 0.26 CH3 + 0.31 aC3H5 + 0.781607
    C2H4 + 0.062529 C3H6 + 0.859765 C5H6 + 0.367353 H2 + 0.361213 C6H6 +
    0.084728 C6H5CH3  # Reaction 3
  rate-constant: {A: 1.84e-05, b: 4.891, Ea: 3621.0}
- equation: C10H16 + OH => H2O + 0.43 H + 0.26 CH3 + 0.31 aC3H5 + 0.781607
    C2H4 + 0.062529 C3H6 + 0.859765 C5H6 + 0.367353 H2 + 0.361213 C6H6 +
    0.084728 C6H5CH3  # Reaction 4
  rate-constant: {A: 1.42e+04, b: 2.878, Ea: -1052.0}
- equation: C10H16 + O2 => HO2 + 0.43 H + 0.26 CH3 + 0.31 aC3H5 + 0.781607
    C2H4 + 0.062529 C3H6 + 0.859765 C5H6 + 0.367353 H2 + 0.361213 C6H6 +
    0.084728 C6H5CH3  # Reaction 5
  rate-constant: {A: 4.37e+12, b: 0.949, Ea: 4.3495e+04}
- equation: C10H16 + O => OH + 0.43 H + 0.26 CH3 + 0.31 aC3H5 + 0.781607
    C2H4 + 0.062529 C3H6 + 0.859765 C5H6 + 0.367353 H2 + 0.361213 C6H6 +
    0.084728 C6H5CH3  # Reaction 6
  rate-constant: {A: 1.03e+06, b: 2.731, Ea: 1158.0}
- equation: C10H16 + HO2 => H2O2 + 0.43 H + 0.26 CH3 + 0.31 aC3H5 + 0.781607
    C2H4 + 0.062529 C3H6 + 0.859765 C5H6 + 0.367353 H2 + 0.361213 C6H6 +
    0.084728 C6H5CH3  # Reaction 7

[jiweiqi]

Reference profiles

[jiweiqi]

With the optimized HyChem model, the system does show non-monotonic species profiles.

[jiweiqi]

This is an interesting phenomenon that is similar to NTC. Might merit some efforts.

[jiweiqi]

[jiweiqi]