# Optimizing an ODEProblem with Threads?

I have a sizable (402) ODE system which comes from the discretization of a parabolic PDE (Beeler-Reuter model in 1D); in the past I used ModelingToolkit to optimize my poorly written in-place function by evaluating it over Julia’s threads and to generate a likewise threaded Jacobian (the system is pretty stiff).

However, some updates to the package seem to have removed the functions I was using and replaced them with a simpler method which does not work on my system.

``````using DifferentialEquations, ModelingToolkit

const N = 50
const D = 0.001/(0.02*0.02)

function ab(C,V)
# eq (13) from original paper
return (C*exp(C*(V+C))+C*(V+C))/(exp(C*(V+C))+C)
end

function stimulationcurrent(oscillatorcomponent,p)
return p.*oscillatorcomponent.^p
end

function Casmoother(Ca; ep=1.5e-10)
return ep*0.5*(1.0 + tanh(1.0-(Ca/ep))) + Ca*0.5*(1.0 + tanh((Ca/ep)-1.0))
end

function BRN!(dx,x,p,t)

V = view(  x,(1+0N):(1N))
C = view(  x,(1+1N):(2N))
X = view(  x,(1+2N):(3N))
m = view(  x,(1+3N):(4N))
h = view(  x,(1+4N):(5N))
j = view(  x,(1+5N):(6N))
d = view(  x,(1+6N):(7N))
f = view(  x,(1+7N):(8N))

dV = view(dx,(1+0N):(1N))
dC = view(dx,(1+1N):(2N))
dX = view(dx,(1+2N):(3N))
dm = view(dx,(1+3N):(4N))
dh = view(dx,(1+4N):(5N))
dj = view(dx,(1+5N):(6N))
dd = view(dx,(1+6N):(7N))
df = view(dx,(1+7N):(8N))

# iterate over space
@inbounds for n in 1:N

# spatially local currents
IK = (exp(0.08*(V[n]+53.0)) + exp(0.04*(V[n]+53.0)))
IK = 4.0*(exp(0.04*(V[n]+85.0)) - 1.0)/IK
IK = IK+0.2*(V[n]+23.0)/(1.0-exp(-0.04*(V[n]+23.0)))
IK = 0.35*IK
Ix = X[n]*0.8*(exp(0.04*(V[n]+77.0))-1.0)/exp(0.04*(V[n]+35.0))
INa= (4.0*m[n]*m[n]*m[n]*h[n]*j[n] + 0.003)*(V[n]-50.0)
Is = 0.09*d[n]*f[n]*(V[n]+82.3+13.0287*log(Casmoother(C[n])))

# these from Beeler & Reuter table:
ax = ab([ 0.0005, 0.083, 50.0, 0.0, 0.0, 0.057, 1.0],V[n])
bx = ab([ 0.0013,-0.06 , 20.0, 0.0, 0.0,-0.04 , 1.0],V[n])
am = ab([ 0.0   , 0.0  , 47.0,-1.0,47.0,-0.1  ,-1.0],V[n])
bm = ab([40.0   ,-0.056, 72.0, 0.0, 0.0, 0.0  , 0.0],V[n])
ah = ab([ 0.126 ,-0.25 , 77.0, 0.0, 0.0, 0.0  , 0.0],V[n])
bh = ab([ 1.7   , 0.0  , 22.5, 0.0, 0.0,-0.082, 1.0],V[n])
aj = ab([ 0.055 ,-0.25 , 78.0, 0.0, 0.0,-0.2  , 1.0],V[n])
bj = ab([ 0.3   , 0.0  , 32.0, 0.0, 0.0,-0.1  , 1.0],V[n])
ad = ab([ 0.095 ,-0.01 , -5.0, 0.0, 0.0,-0.072, 1.0],V[n])
bd = ab([ 0.07  ,-0.017, 44.0, 0.0, 0.0, 0.05 , 1.0],V[n])
af = ab([ 0.012 ,-0.008, 28.0, 0.0, 0.0, 0.15 , 1.0],V[n])
bf = ab([ 0.0065,-0.02 , 30.0, 0.0, 0.0,-0.2  , 1.0],V[n])

# BR dynamics
dV[n] = -(IK + Ix + INa + Is)/p
dC[n] = -10^-7 * Is + 0.07*(10^-7 - C[n])
dX[n] = ax*(1.0 - X[n]) - bx*X[n]
dm[n] = am*(1.0 - m[n]) - bm*m[n]
dh[n] = ah*(1.0 - h[n]) - bh*h[n]
dj[n] = aj*(1.0 - j[n]) - bj*j[n]
dd[n] = ad*(1.0 - d[n]) - bd*d[n]
df[n] = af*(1.0 - f[n]) - bf*f[n]

# diffusion
if N > 1
if n==1
dV[n] = dV[n] + D*(2V[n+1] - 2V[n])
elseif n==N
dV[n] = dV[n] + D*(2V[n-1] - 2V[n])
else
dV[n] = dV[n] + D*(V[n+1] + V[n-1] - 2V[n])
end
end

end

# this is the forcing current, left boundary point only
dx = dx + stimulationcurrent(x[8N+1],p)/p

# nonlinear oscillator dynamics
ω = 2*pi*p/1000.0
dx[8*N+1] = ω*x[8*N+2] + x[8*N+1]*(1.0 - x[8*N+1]^2 - x[8*N+2]^2)
dx[8*N+2] =-ω*x[8*N+1] + x[8*N+2]*(1.0 - x[8*N+1]^2 - x[8*N+2]^2)

return nothing

u0 = zeros(Float64,8N+2)
u0[(0N+1):(1N)] .= -84.0
u0[(1N+1):(2N)] .= 10^-7
u0[(2N+1):(3N)] .= 0.01
u0[(3N+1):(4N)] .= 0.01
u0[(4N+1):(5N)] .= 0.99
u0[(5N+1):(6N)] .= 0.99
u0[(6N+1):(7N)] .= 0.01
u0[(7N+1):(8N)] .= 0.99
u0[8N+1]         = 0.0
u0[8N+2]         = 1.0

# parameters
#	uF/cm2    	uA/cm2       	Hz		power
p = [	1.0,		2.3,		1.25, 		1.0]

# tspan
tspan = (0.0,300.0)

``````

The current method of optimizing the model with threading appears to be making an `ODEProblem`, toolkitizing it, and then generating appropriately simplified and threaded functions for the ODE and Jacobian:

``````prob = ODEProblem(BRN!, u0, tspan, p)
sys = modelingtoolkitize(prob)
However, the call for `J!` fails due to
`ERROR: MethodError: no method matching +(::AbstractAlgebra.Generic.MPoly{BigInt})`, and when calling the successful `f!`, it is still single-threaded (despite `Threads.nthreads() = 8`). Oddly enough, setting `const N=1` and building a serial version of `J!` fails with a similar error, so I think this must be an issue with my right-hand-side.
Since this functionality is not expressly demonstrated in the documentation, I wonder if there is a better way to do this or if there are some limitations to account for the failure of `J!`? Obviously adding `Threads.@threads` to the spatial loop in `BRN!` would parallelize it, but the old simplify/parallel method was considerably faster than even that. Many thanks.