# Calling Cantera using PyCall from within DifferentialEquations function

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:

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:

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.

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

It works, thanks a lot for your help