# Calling Cantera using PyCall from within DifferentialEquations function

**URL:** <https://discourse.julialang.org/t/calling-cantera-using-pycall-from-within-differentialequations-function/104942>\
**Category:** Chemistry\
**Created:** [October 13, 2023, 11:35pm UTC](https://discourse.julialang.org/t/calling-cantera-using-pycall-from-within-differentialequations-function/104942 "2023-10-13T23:35:45Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![NikoBiele](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nikobiele/32/208184_2.png) [@NikoBiele](https://discourse.julialang.org/u/NikoBiele)\
**Post date:** [October 13, 2023, 11:35pm UTC](https://discourse.julialang.org/t/calling-cantera-using-pycall-from-within-differentialequations-function/104942/1 "2023-10-13T23:35:45Z")

</div>

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](https://cantera.org/examples/python/reactors/custom.py.html)) Python script to Julia. In order to install Cantera to run the code below, please follow:

1. [Installing Cantera | Cantera](https://cantera.org/install/index.html)
2. [Calling chemical kinetics solvers like CHEMKIN or Cantera from Julia - #2 by R\_Surya\_Narayan](https://discourse.julialang.org/t/calling-chemical-kinetics-solvers-like-chemkin-or-cantera-from-julia/50568/2)

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

```julia
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](https://global.discourse-cdn.com/julialang/original/3X/2/e/2ed4e338ac5a98510701a89083bfac32f5c93c9c.png)

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

```julia
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. 🙂

---

<div class="post-metadata">

**Author:** ![contradict](https://avatars.discourse-cdn.com/v4/letter/c/ac91a4/32.png) [@contradict](https://discourse.julialang.org/u/contradict)\
**Post date:** [October 14, 2023, 12:04am UTC](https://discourse.julialang.org/t/calling-cantera-using-pycall-from-within-differentialequations-function/104942/2 "2023-10-14T00:04:39Z")

</div>

> [@NikoBiele](#):
>
> `du = [dTdt; dYdt]`

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

```julia
du .= [dTdt; dYdt]

```

---

<div class="post-metadata">

**Author:** ![NikoBiele](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nikobiele/32/208184_2.png) [@NikoBiele](https://discourse.julialang.org/u/NikoBiele)\
**Post date:** [October 14, 2023, 6:14am UTC](https://discourse.julialang.org/t/calling-cantera-using-pycall-from-within-differentialequations-function/104942/3 "2023-10-14T06:14:27Z")

</div>

It works, thanks a lot for your help 🙂