Calling Cantera using PyCall from within DifferentialEquations function

Specific Domains Chemistry
1

Hi,

I’m trying to simulate the combustion of hydrogen using Cantera from Julia via PyCall. But I’m having trouble using PyCall from within my DifferentialEquations function. I have translated this (custom.py | Cantera) Python script to Julia. In order to install Cantera to run the code below, please follow:

  1. Installing Cantera | Cantera
  2. Calling chemical kinetics solvers like CHEMKIN or Cantera from Julia - #2 by R_Surya_Narayan

First a version of the script which uses a for-loop with simple Euler integration:

using PyCall
using LinearAlgebra
using Plots

# import Cantera
ct = pyimport("cantera")

# create a gas
gas = ct.Solution("h2o2.yaml")

# Initial condition
P = ct.one_atm
gas.TPX = 1250.0, P, "H2:2,O2:1,N2:4"
u0 = [gas.T; gas.Y]
u = zeros(length(u0),500_000)
u[:,1] = u0

# timestep
dt = 1e-9

for i = 1:500_000-1
    # set state
    gas.TPY = u[1,i], P, u[2:end,i]

    # State vector is [T, Y_1, Y_2, ... Y_K]
    gas.set_unnormalized_mass_fractions(u[2:end,i])
    gas.TP = u[1], P
    rho = gas.density

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

    du = [dTdt; dYdt]

    u[:,i+1] = u[:,i] + du[:] * dt
end
plot(legend=false, xlabel="time / s", ylabel="molefraction")
for i = 2:length(u0)-1
    display(plot!([1:100:length(u[1,:])]*dt,u[i,1:100:end]))
end

Which produces a plot that makes sense:

Screenshot 2023-10-14 011550
Screenshot 2023-10-14 0115501370×911 28.2 KB

Next, my attempt at implementing the same code with DifferentialEquations.jl:

using PyCall
using DifferentialEquations
using LinearAlgebra

# import Cantera
ct = pyimport("cantera")

# create a gas
gas = ct.Solution("h2o2.yaml")

# Initial condition
P = ct.one_atm
gas.TPX = 1250.0, P, "H2:2,O2:1,N2:4"
u0 = [gas.T; gas.Y]

function react!(du, u, params, t)
    # parameters
    P = params

    # set state
    gas.TPY = u[1], P, u[2:end]

    # State vector is [T, Y_1, Y_2, ... Y_K]
    gas.set_unnormalized_mass_fractions(u[2:end])
    gas.TP = u[1], P
    rho = gas.density

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

    du = [dTdt; dYdt]
end

tspan = (0.0, 1e-3)
params = P
prob = ODEProblem(react!, u0, tspan, params)
sol = solve(prob, Rodas5(autodiff=false), abstol = 1e-12, reltol = 1e-12)

This code will run, but produces zero change in concentrations of reactants. Any help is highly appreciated. :slight_smile:

2

This line creates a new variable du with the provided value. In order to modify the output paramter du, you should do

du .= [dTdt; dYdt]
2 Likes
3

It works, thanks a lot for your help :slight_smile: