Julia: A nice benchmark for compiler performance
I was trying to provide myself with a large collection (>500) of test functions for optimization and came across an interesting benchmark for compiler performance. It stresses the ability to infer long function bodies. To understand what the body looks like, follow the URL in the download statement (and look at the expression immediately after the obj =e= statement).
Here's the script:
if !isfile("chebyqad.gms")
download("http://www.gamsworld.org/performance/princetonlib/raw/cute/chebyqad.gms", "chebyqad.gms")
end
# Some definitions we'll need
sqr(x) = x*x
arctan(x) = atan(x)
const pi2 = π/2
## Add type-annotations to variables that start with 'x'. (Doesn't help.)
function annotate!(ex::Expr)
for i = 1:length(ex.args)
ex.args[i] = annotate!(ex.args[i])
end
return ex
end
function annotate!(sym::Symbol)
if startswith(string(sym), 'x')
return :($sym::Float64)
end
return sym
end
annotate!(arg) = arg
println("parsing time:")
@time open("chebyqad.gms") do io
while true
line = readline(io)
if !startswith(line, " obj =e=")
continue
end
str = readuntil(io, ';')
body = Meta.parse(str)
if isa(body, Expr) && body.head == :toplevel
@assert length(body.args) == 1
body = body.args[1]
end
# body = annotate!(body)
vars = [Symbol("x",i) for i = 1:50]
destruct = Expr(:tuple, vars...)
@eval function f(x)
$destruct = x
$body
end
break
end
end
println("done parsing, about to call f")
@time f(rand(50))
@time f(rand(50))
Results (note: this was with an older version of the script above):
julia> include("compiler_perf.jl")
56.592909 seconds (17.62 M allocations: 729.870 MiB, 0.95% gc time)
0.000420 seconds (5.21 k allocations: 122.188 KiB)
0.4243939871551012
So of course the performance is quite nice once it's compiled, but the compilation takes a long time. I tried profiling on 0.7 and got no stacktraces (didn't figure out why), but on 0.6 most of the time was in effect_free, occurs_outside_getfield and replace_getfield! (all inference.jl). For reference 0.6 is not quite twice as fast as 0.7 to build this.
Adding a @noinline or @nospecialize to f's declaration doesn't appear to help at all.
This is not the worst case I've encountered: that mantle belongs to http://www.gamsworld.org/performance/princetonlib/htm/nnls/nnls.htm, which I have yet to see finish inference.
timholy
❤2
All 11 comments
https://github.com/Keno/NewOptimizer.jl might help. It's getfield elim pass is linear time.
Keno
on 7 Feb 2018
❤3
Would it be possible/feasible/desirable to add this to BaseBenchmarks?
ararslan
on 7 Feb 2018
We should set up a compiler benchmark suite. Not sure it needs to be part of BaseBenchmarks, and we probably need some infrastructure in base to do it (e.g. compile this subtree without using any of the caches).
Keno
on 7 Feb 2018
(e.g. compile this subtree without using any of the caches).
That part we have, since you can invoke jl_get_llvmf_defn with params = Base.CodegenParams(cached=false,...) which should ignore all the caches.
vchuravy
on 7 Feb 2018
Yeah, but that runs into https://github.com/JuliaGPU/CUDAnative.jl/issues/122.
Keno
on 7 Feb 2018
Now that the new optimizer is in, I gave this another whirl (script above slightly updated):
julia 0.6.3:
julia> include("compiler_perf.jl")
parsing time:
2.074501 seconds (336.90 k allocations: 21.286 MiB)
done parsing, about to call f
25.365329 seconds (19.31 M allocations: 704.055 MiB, 1.67% gc time)
0.000391 seconds (4.96 k allocations: 78.031 KiB)
0.46427294808313835
julia 0.7.0-beta.129: still waiting after 20 minutes...
I trimmed the file size down to about 10% of the full size and collected this profile. The biggest offender by far seems to be construct_domtree!, specifically its usage of reduce. Hope that's useful.
timholy
on 3 Jul 2018
❤5
construct_domtree! is O(N^3) in the size of the function. There's significantly better algorithms, but it didn't feature prominently in the benchmarks we looked at, so we didn't fix that yet. We should use one of them.
Keno
on 3 Jul 2018
👍3
Here is another example (https://gist.github.com/KristofferC/2c37b616e0cc6769069690b23fbd3310) taken from the FastGaussQuadrature package which also has the construct_domtree! showing up in the profile.
julia> @time feval_asy1(n, a, b, t, idx);
11.324220 seconds (30.48 M allocations: 1.346 GiB, 5.93% gc time)
julia> @time feval_asy1(n, a, b, t, idx);
0.000816 seconds (1.56 k allocations: 537.313 KiB)
It does create kinda nasty types though...
│ %7958 = new(Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{2},Nothing,typeof(*),Tuple{Int64,Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{2},Nothing,typeof(-),Tuple{Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{0},Nothing,typeof(-),Tuple{Float64,Float64}},Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{2},Nothing,typeof(*),Tuple{Float64,Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{2},Nothing,typeof(cos),Tuple{Adjoint{Float64,Array{Float64,1}}}}}}}}}}, *, %7957, nothing)::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{2},Nothing,typeof(*),Tuple{Int64,Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{2},Nothing,typeof(-),Tuple{Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{0},Nothing,typeof(-),Tuple{Float64,Float64}},Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{2},Nothing,typeof(*),Tuple{Float64,Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{2},Nothing,typeof(cos),Tuple{Adjoint{Float64,Array{Float64,1}}}}}}}}}}
KristofferC
on 6 Jul 2018
https://github.com/JuliaLang/julia/issues/27488 seems to be another example of this.
P2D2c = function(dy,y,pars,t)
Dbulk = pars[1]
Dsn = pars[2]
Dsp = pars[3]
F = 96487.000000
R = 8.314300
Rpn = pars[4]
Rpp = pars[5]
Tr = 298.150000
brugn = pars[6]
brugp = pars[7]
brugs = pars[8]
c0 = pars[9]
ctn = pars[10]
ctp = pars[11]
efn = pars[12]
efp = pars[13]
en = pars[14]
ep = pars[15]
es = pars[16]
kn = pars[17]
kp = pars[18]
ln1 = pars[19]
lp = pars[20]
ls = pars[21]
sigmn = pars[22]
sigmp = pars[23]
socn = pars[24]
socp = pars[25]
t1 = pars[26]
current = -1.0
capacity=17.5
rate=-2.
iapp=capacity*rate
up=[0.,0.]
Kr = pars[27]
ap = 3/Rpp*(1-ep-efp)
an = 3/Rpn*(1-en-efn)
sigmap = sigmp*(1-ep-efp)
sigman = sigmn*(1-en-efn)
D2pos = ep^brugp*Dbulk
D2sep = es^brugs*Dbulk
D2neg = en^brugn*Dbulk
dy[1] = 1/ep*(64.0*D2pos/lp^2*(y[60]-2.0*y[1]+y[2])+2.0*ap*(1.0-1.0*t1)*kp*(y[1]*c0)^.5*(-1.0*y[110]*ctp+ctp)^.5*(y[110]*ctp)^.5*sinh(.5*F/R/Tr*(y[86]-1.0*y[65]-1.0*ocp_cathode(y[110],up)))/c0)
dy[2] = 1/ep*(64.0*D2pos/lp^2*(y[1]-2.0*y[2]+y[3])+2.0*ap*(1.0-1.0*t1)*kp*(y[2]*c0)^.5*(-1.0*y[111]*ctp+ctp)^.5*(y[111]*ctp)^.5*sinh(.5*F/R/Tr*(y[87]-1.0*y[66]-1.0*ocp_cathode(y[111],up)))/c0)
dy[3] = 1/ep*(64.0*D2pos/lp^2*(y[2]-2.0*y[3]+y[4])+2.0*ap*(1.0-1.0*t1)*kp*(y[3]*c0)^.5*(-1.0*y[112]*ctp+ctp)^.5*(y[112]*ctp)^.5*sinh(.5*F/R/Tr*(y[88]-1.0*y[67]-1.0*ocp_cathode(y[112],up)))/c0)
dy[4] = 1/ep*(64.0*D2pos/lp^2*(y[3]-2.0*y[4]+y[5])+2.0*ap*(1.0-1.0*t1)*kp*(y[4]*c0)^.5*(-1.0*y[113]*ctp+ctp)^.5*(y[113]*ctp)^.5*sinh(.5*F/R/Tr*(y[89]-1.0*y[68]-1.0*ocp_cathode(y[113],up)))/c0)
dy[5] = 1/ep*(64.0*D2pos/lp^2*(y[4]-2.0*y[5]+y[6])+2.0*ap*(1.0-1.0*t1)*kp*(y[5]*c0)^.5*(-1.0*y[114]*ctp+ctp)^.5*(y[114]*ctp)^.5*sinh(.5*F/R/Tr*(y[90]-1.0*y[69]-1.0*ocp_cathode(y[114],up)))/c0)
dy[6] = 1/ep*(64.0*D2pos/lp^2*(y[5]-2.0*y[6]+y[7])+2.0*ap*(1.0-1.0*t1)*kp*(y[6]*c0)^.5*(-1.0*y[115]*ctp+ctp)^.5*(y[115]*ctp)^.5*sinh(.5*F/R/Tr*(y[91]-1.0*y[70]-1.0*ocp_cathode(y[115],up)))/c0)
dy[7] = 1/ep*(64.0*D2pos/lp^2*(y[6]-2.0*y[7]+y[61])+2.0*ap*(1.0-1.0*t1)*kp*(y[7]*c0)^.5*(-1.0*y[116]*ctp+ctp)^.5*(y[116]*ctp)^.5*sinh(.5*F/R/Tr*(y[92]-1.0*y[71]-1.0*ocp_cathode(y[116],up)))/c0)
dy[8] = 16.0/es*D2sep/ls^2*(y[61]-2.0*y[8]+y[9])
dy[9] = 16.0/es*D2sep/ls^2*(y[8]-2.0*y[9]+y[10])
dy[10] = 16.0/es*D2sep/ls^2*(y[9]-2.0*y[10]+y[62])
dy[11] = 1/en*(64.0*D2neg/ln1^2*(y[62]-2.0*y[11]+y[12])+2.0*an*(1.0-1.0*t1)*kn*(y[11]*c0)^.5*(-1.0*y[124]*ctn+ctn)^.5*(y[124]*ctn)^.5*sinh(.5*F/R/Tr*(y[95]-1.0*y[77]+.5e-1-2.325*exp(-100.0*y[124]^1.15)+.1721*tanh(20.000000*y[124]-21.000000)+.25e-2*tanh(39.347962*y[124]-37.241379)+.34e-1*tanh(89.411765*y[124]-6.0294118)+.2e-2*tanh(7.0422535*y[124]-1.3661972)+.155e-1*tanh(77.519380*y[124]-8.1395349)))/c0)
dy[12] = 1/en*(64.0*D2neg/ln1^2*(y[11]-2.0*y[12]+y[13])+2.0*an*(1.0-1.0*t1)*kn*(y[12]*c0)^.5*(-1.0*y[125]*ctn+ctn)^.5*(y[125]*ctn)^.5*sinh(.5*F/R/Tr*(y[96]-1.0*y[78]+.5e-1-2.325*exp(-100.0*y[125]^1.15)+.1721*tanh(20.000000*y[125]-21.000000)+.25e-2*tanh(39.347962*y[125]-37.241379)+.34e-1*tanh(89.411765*y[125]-6.0294118)+.2e-2*tanh(7.0422535*y[125]-1.3661972)+.155e-1*tanh(77.519380*y[125]-8.1395349)))/c0)
dy[13] = 1/en*(64.0*D2neg/ln1^2*(y[12]-2.0*y[13]+y[14])+2.0*an*(1.0-1.0*t1)*kn*(y[13]*c0)^.5*(-1.0*y[126]*ctn+ctn)^.5*(y[126]*ctn)^.5*sinh(.5*F/R/Tr*(y[97]-1.0*y[79]+.5e-1-2.325*exp(-100.0*y[126]^1.15)+.1721*tanh(20.000000*y[126]-21.000000)+.25e-2*tanh(39.347962*y[126]-37.241379)+.34e-1*tanh(89.411765*y[126]-6.0294118)+.2e-2*tanh(7.0422535*y[126]-1.3661972)+.155e-1*tanh(77.519380*y[126]-8.1395349)))/c0)
dy[14] = 1/en*(64.0*D2neg/ln1^2*(y[13]-2.0*y[14]+y[15])+2.0*an*(1.0-1.0*t1)*kn*(y[14]*c0)^.5*(-1.0*y[127]*ctn+ctn)^.5*(y[127]*ctn)^.5*sinh(.5*F/R/Tr*(y[98]-1.0*y[80]+.5e-1-2.325*exp(-100.0*y[127]^1.15)+.1721*tanh(20.000000*y[127]-21.000000)+.25e-2*tanh(39.347962*y[127]-37.241379)+.34e-1*tanh(89.411765*y[127]-6.0294118)+.2e-2*tanh(7.0422535*y[127]-1.3661972)+.155e-1*tanh(77.519380*y[127]-8.1395349)))/c0)
dy[15] = 1/en*(64.0*D2neg/ln1^2*(y[14]-2.0*y[15]+y[16])+2.0*an*(1.0-1.0*t1)*kn*(y[15]*c0)^.5*(-1.0*y[128]*ctn+ctn)^.5*(y[128]*ctn)^.5*sinh(.5*F/R/Tr*(y[99]-1.0*y[81]+.5e-1-2.325*exp(-100.0*y[128]^1.15)+.1721*tanh(20.000000*y[128]-21.000000)+.25e-2*tanh(39.347962*y[128]-37.241379)+.34e-1*tanh(89.411765*y[128]-6.0294118)+.2e-2*tanh(7.0422535*y[128]-1.3661972)+.155e-1*tanh(77.519380*y[128]-8.1395349)))/c0)
dy[16] = 1/en*(64.0*D2neg/ln1^2*(y[15]-2.0*y[16]+y[17])+2.0*an*(1.0-1.0*t1)*kn*(y[16]*c0)^.5*(-1.0*y[129]*ctn+ctn)^.5*(y[129]*ctn)^.5*sinh(.5*F/R/Tr*(y[100]-1.0*y[82]+.5e-1-2.325*exp(-100.0*y[129]^1.15)+.1721*tanh(20.000000*y[129]-21.000000)+.25e-2*tanh(39.347962*y[129]-37.241379)+.34e-1*tanh(89.411765*y[129]-6.0294118)+.2e-2*tanh(7.0422535*y[129]-1.3661972)+.155e-1*tanh(77.519380*y[129]-8.1395349)))/c0)
dy[17] = 1/en*(64.0*D2neg/ln1^2*(y[16]-2.0*y[17]+y[63])+2.0*an*(1.0-1.0*t1)*kn*(y[17]*c0)^.5*(-1.0*y[130]*ctn+ctn)^.5*(y[130]*ctn)^.5*sinh(.5*F/R/Tr*(y[101]-1.0*y[83]+.5e-1-2.325*exp(-100.0*y[130]^1.15)+.1721*tanh(20.000000*y[130]-21.000000)+.25e-2*tanh(39.347962*y[130]-37.241379)+.34e-1*tanh(89.411765*y[130]-6.0294118)+.2e-2*tanh(7.0422535*y[130]-1.3661972)+.155e-1*tanh(77.519380*y[130]-8.1395349)))/c0)
dy[18] = 16.0*(Dsp*(y[25]-1.0*y[103])+Dsp*(y[103]-2.0*y[18]+y[25]))/Rpp^2
dy[19] = 16.0*(Dsp*(y[26]-1.0*y[104])+Dsp*(y[104]-2.0*y[19]+y[26]))/Rpp^2
dy[20] = 16.0*(Dsp*(y[27]-1.0*y[105])+Dsp*(y[105]-2.0*y[20]+y[27]))/Rpp^2
dy[21] = 16.0*(Dsp*(y[28]-1.0*y[106])+Dsp*(y[106]-2.0*y[21]+y[28]))/Rpp^2
dy[22] = 16.0*(Dsp*(y[29]-1.0*y[107])+Dsp*(y[107]-2.0*y[22]+y[29]))/Rpp^2
dy[23] = 16.0*(Dsp*(y[30]-1.0*y[108])+Dsp*(y[108]-2.0*y[23]+y[30]))/Rpp^2
dy[24] = 16.0*(Dsp*(y[31]-1.0*y[109])+Dsp*(y[109]-2.0*y[24]+y[31]))/Rpp^2
dy[25] = 4.0*(2.0*Dsp*(y[32]-1.0*y[18])+4.0*Dsp*(y[18]-2.0*y[25]+y[32]))/Rpp^2
dy[26] = 4.0*(2.0*Dsp*(y[33]-1.0*y[19])+4.0*Dsp*(y[19]-2.0*y[26]+y[33]))/Rpp^2
dy[27] = 4.0*(2.0*Dsp*(y[34]-1.0*y[20])+4.0*Dsp*(y[20]-2.0*y[27]+y[34]))/Rpp^2
dy[28] = 4.0*(2.0*Dsp*(y[35]-1.0*y[21])+4.0*Dsp*(y[21]-2.0*y[28]+y[35]))/Rpp^2
dy[29] = 4.0*(2.0*Dsp*(y[36]-1.0*y[22])+4.0*Dsp*(y[22]-2.0*y[29]+y[36]))/Rpp^2
dy[30] = 4.0*(2.0*Dsp*(y[37]-1.0*y[23])+4.0*Dsp*(y[23]-2.0*y[30]+y[37]))/Rpp^2
dy[31] = 4.0*(2.0*Dsp*(y[38]-1.0*y[24])+4.0*Dsp*(y[24]-2.0*y[31]+y[38]))/Rpp^2
dy[32] = 1.777777778*(3.0*Dsp*(y[110]-1.0*y[25])+9.0*Dsp*(y[25]-2.0*y[32]+y[110]))/Rpp^2
dy[33] = 1.777777778*(3.0*Dsp*(y[111]-1.0*y[26])+9.0*Dsp*(y[26]-2.0*y[33]+y[111]))/Rpp^2
dy[34] = 1.777777778*(3.0*Dsp*(y[112]-1.0*y[27])+9.0*Dsp*(y[27]-2.0*y[34]+y[112]))/Rpp^2
dy[35] = 1.777777778*(3.0*Dsp*(y[113]-1.0*y[28])+9.0*Dsp*(y[28]-2.0*y[35]+y[113]))/Rpp^2
dy[36] = 1.777777778*(3.0*Dsp*(y[114]-1.0*y[29])+9.0*Dsp*(y[29]-2.0*y[36]+y[114]))/Rpp^2
dy[37] = 1.777777778*(3.0*Dsp*(y[115]-1.0*y[30])+9.0*Dsp*(y[30]-2.0*y[37]+y[115]))/Rpp^2
dy[38] = 1.777777778*(3.0*Dsp*(y[116]-1.0*y[31])+9.0*Dsp*(y[31]-2.0*y[38]+y[116]))/Rpp^2
dy[39] = 16.0*Dsn*(2.0*y[46]-2.0*y[39])/Rpn^2
dy[40] = 16.0*Dsn*(2.0*y[47]-2.0*y[40])/Rpn^2
dy[41] = 16.0*Dsn*(2.0*y[48]-2.0*y[41])/Rpn^2
dy[42] = 16.0*Dsn*(2.0*y[49]-2.0*y[42])/Rpn^2
dy[43] = 16.0*Dsn*(2.0*y[50]-2.0*y[43])/Rpn^2
dy[44] = 16.0*Dsn*(2.0*y[51]-2.0*y[44])/Rpn^2
dy[45] = 16.0*Dsn*(2.0*y[52]-2.0*y[45])/Rpn^2
dy[46] = 4.0*Dsn*(6.0*y[53]+2.0*y[39]-8.0*y[46])/Rpn^2
dy[47] = 4.0*Dsn*(6.0*y[54]+2.0*y[40]-8.0*y[47])/Rpn^2
dy[48] = 4.0*Dsn*(6.0*y[55]+2.0*y[41]-8.0*y[48])/Rpn^2
dy[49] = 4.0*Dsn*(6.0*y[56]+2.0*y[42]-8.0*y[49])/Rpn^2
dy[50] = 4.0*Dsn*(6.0*y[57]+2.0*y[43]-8.0*y[50])/Rpn^2
dy[51] = 4.0*Dsn*(6.0*y[58]+2.0*y[44]-8.0*y[51])/Rpn^2
dy[52] = 4.0*Dsn*(6.0*y[59]+2.0*y[45]-8.0*y[52])/Rpn^2
dy[53] = 1.777777778*Dsn*(12.0*y[124]+6.0*y[46]-18.0*y[53])/Rpn^2
dy[54] = 1.777777778*Dsn*(12.0*y[125]+6.0*y[47]-18.0*y[54])/Rpn^2
dy[55] = 1.777777778*Dsn*(12.0*y[126]+6.0*y[48]-18.0*y[55])/Rpn^2
dy[56] = 1.777777778*Dsn*(12.0*y[127]+6.0*y[49]-18.0*y[56])/Rpn^2
dy[57] = 1.777777778*Dsn*(12.0*y[128]+6.0*y[50]-18.0*y[57])/Rpn^2
dy[58] = 1.777777778*Dsn*(12.0*y[129]+6.0*y[51]-18.0*y[58])/Rpn^2
dy[59] = 1.777777778*Dsn*(12.0*y[130]+6.0*y[52]-18.0*y[59])/Rpn^2
dy[60] = 4.0*D2pos*(-1.0*y[2]-3.0*y[60]+4.0*y[1])/lp
dy[61] = 4.0*D2pos*(y[6]+3.0*y[61]-4.0*y[7])/lp-2.0*D2sep*(-1.0*y[9]-3.0*y[61]+4.0*y[8])/ls
dy[62] = 2.0*D2sep*(y[9]+3.0*y[62]-4.0*y[10])/ls-4.0*D2neg*(-1.0*y[12]-3.0*y[62]+4.0*y[11])/ln1
dy[63] = 4.0*D2neg*(y[16]+3.0*y[63]-4.0*y[17])/ln1
dy[64] = 4.0*(-1.0*y[66]-3.0*y[64]+4.0*y[65])/lp
dy[65] = -4.0*sigmap/lp*(y[87]-1.0*y[85])-4.0*ep^brugp*(.41253e-1+.5007e-3*y[1]-.47212e-6*y[1]^2+.15904e-9*y[1]^3-.16018e-13*y[1]^4)/lp*(y[66]-1.0*y[64])+8.0*ep^brugp*(.41253e-1+.5007e-3*y[1]-.47212e-6*y[1]^2+.15904e-9*y[1]^3-.16018e-13*y[1]^4)*R*Tr/F*(1.0-1.0*t1)/lp*(y[2]-1.0*y[60])/y[1]-iapp
dy[66] = -4.0*sigmap/lp*(y[88]-1.0*y[86])-4.0*ep^brugp*(.41253e-1+.5007e-3*y[2]-.47212e-6*y[2]^2+.15904e-9*y[2]^3-.16018e-13*y[2]^4)/lp*(y[67]-1.0*y[65])+8.0*ep^brugp*(.41253e-1+.5007e-3*y[2]-.47212e-6*y[2]^2+.15904e-9*y[2]^3-.16018e-13*y[2]^4)*R*Tr/F*(1.0-1.0*t1)/lp*(y[3]-1.0*y[1])/y[2]-iapp
dy[67] = -4.0*sigmap/lp*(y[89]-1.0*y[87])-4.0*ep^brugp*(.41253e-1+.5007e-3*y[3]-.47212e-6*y[3]^2+.15904e-9*y[3]^3-.16018e-13*y[3]^4)/lp*(y[68]-1.0*y[66])+8.0*ep^brugp*(.41253e-1+.5007e-3*y[3]-.47212e-6*y[3]^2+.15904e-9*y[3]^3-.16018e-13*y[3]^4)*R*Tr/F*(1.0-1.0*t1)/lp*(y[4]-1.0*y[2])/y[3]-iapp
dy[68] = -4.0*sigmap/lp*(y[90]-1.0*y[88])-4.0*ep^brugp*(.41253e-1+.5007e-3*y[4]-.47212e-6*y[4]^2+.15904e-9*y[4]^3-.16018e-13*y[4]^4)/lp*(y[69]-1.0*y[67])+8.0*ep^brugp*(.41253e-1+.5007e-3*y[4]-.47212e-6*y[4]^2+.15904e-9*y[4]^3-.16018e-13*y[4]^4)*R*Tr/F*(1.0-1.0*t1)/lp*(y[5]-1.0*y[3])/y[4]-iapp
dy[69] = -4.0*sigmap/lp*(y[91]-1.0*y[89])-4.0*ep^brugp*(.41253e-1+.5007e-3*y[5]-.47212e-6*y[5]^2+.15904e-9*y[5]^3-.16018e-13*y[5]^4)/lp*(y[70]-1.0*y[68])+8.0*ep^brugp*(.41253e-1+.5007e-3*y[5]-.47212e-6*y[5]^2+.15904e-9*y[5]^3-.16018e-13*y[5]^4)*R*Tr/F*(1.0-1.0*t1)/lp*(y[6]-1.0*y[4])/y[5]-iapp
dy[70] = -4.0*sigmap/lp*(y[92]-1.0*y[90])-4.0*ep^brugp*(.41253e-1+.5007e-3*y[6]-.47212e-6*y[6]^2+.15904e-9*y[6]^3-.16018e-13*y[6]^4)/lp*(y[71]-1.0*y[69])+8.0*ep^brugp*(.41253e-1+.5007e-3*y[6]-.47212e-6*y[6]^2+.15904e-9*y[6]^3-.16018e-13*y[6]^4)*R*Tr/F*(1.0-1.0*t1)/lp*(y[7]-1.0*y[5])/y[6]-iapp
dy[71] = -4.0*sigmap/lp*(y[93]-1.0*y[91])-4.0*ep^brugp*(.41253e-1+.5007e-3*y[7]-.47212e-6*y[7]^2+.15904e-9*y[7]^3-.16018e-13*y[7]^4)/lp*(y[72]-1.0*y[70])+8.0*ep^brugp*(.41253e-1+.5007e-3*y[7]-.47212e-6*y[7]^2+.15904e-9*y[7]^3-.16018e-13*y[7]^4)*R*Tr/F*(1.0-1.0*t1)/lp*(y[61]-1.0*y[6])/y[7]-iapp
dy[72] = 4.0*ep^brugp*(.41253e-1+.5007e-3*y[61]-.47212e-6*y[61]^2+.15904e-9*y[61]^3-.16018e-13*y[61]^4)*(y[70]+3.0*y[72]-4.0*y[71])/lp-2.0*es^brugs*(.41253e-1+.5007e-3*y[61]-.47212e-6*y[61]^2+.15904e-9*y[61]^3-.16018e-13*y[61]^4)*(-1.0*y[74]-3.0*y[72]+4.0*y[73])/ls
dy[73] = -2.0*es^brugs*(.41253e-1+.5007e-3*y[8]-.47212e-6*y[8]^2+.15904e-9*y[8]^3-.16018e-13*y[8]^4)/ls*(y[74]-1.0*y[72])+4.0*es^brugs*(.41253e-1+.5007e-3*y[8]-.47212e-6*y[8]^2+.15904e-9*y[8]^3-.16018e-13*y[8]^4)*R*Tr/F*(1.0-1.0*t1)/ls*(y[9]-1.0*y[61])/y[8]-iapp
dy[74] = -2.0*es^brugs*(.41253e-1+.5007e-3*y[9]-.47212e-6*y[9]^2+.15904e-9*y[9]^3-.16018e-13*y[9]^4)/ls*(y[75]-1.0*y[73])+4.0*es^brugs*(.41253e-1+.5007e-3*y[9]-.47212e-6*y[9]^2+.15904e-9*y[9]^3-.16018e-13*y[9]^4)*R*Tr/F*(1.0-1.0*t1)/ls*(y[10]-1.0*y[8])/y[9]-iapp
dy[75] = -2.0*es^brugs*(.41253e-1+.5007e-3*y[10]-.47212e-6*y[10]^2+.15904e-9*y[10]^3-.16018e-13*y[10]^4)/ls*(y[76]-1.0*y[74])+4.0*es^brugs*(.41253e-1+.5007e-3*y[10]-.47212e-6*y[10]^2+.15904e-9*y[10]^3-.16018e-13*y[10]^4)*R*Tr/F*(1.0-1.0*t1)/ls*(y[62]-1.0*y[9])/y[10]-iapp
dy[76] = 2.0*es^brugs*(.41253e-1+.5007e-3*y[62]-.47212e-6*y[62]^2+.15904e-9*y[62]^3-.16018e-13*y[62]^4)*(y[74]+3.0*y[76]-4.0*y[75])/ls-4.0*en^brugn*(.41253e-1+.5007e-3*y[62]-.47212e-6*y[62]^2+.15904e-9*y[62]^3-.16018e-13*y[62]^4)*(-1.0*y[78]-3.0*y[76]+4.0*y[77])/ln1
dy[77] = -4.0*sigman/ln1*(y[96]-1.0*y[94])-4.0*en^brugn*(.41253e-1+.5007e-3*y[11]-.47212e-6*y[11]^2+.15904e-9*y[11]^3-.16018e-13*y[11]^4)/ln1*(y[78]-1.0*y[76])+8.0*en^brugn*(.41253e-1+.5007e-3*y[11]-.47212e-6*y[11]^2+.15904e-9*y[11]^3-.16018e-13*y[11]^4)*R*Tr/F*(1.0-1.0*t1)/ln1*(y[12]-1.0*y[62])/y[11]-iapp
dy[78] = -4.0*sigman/ln1*(y[97]-1.0*y[95])-4.0*en^brugn*(.41253e-1+.5007e-3*y[12]-.47212e-6*y[12]^2+.15904e-9*y[12]^3-.16018e-13*y[12]^4)/ln1*(y[79]-1.0*y[77])+8.0*en^brugn*(.41253e-1+.5007e-3*y[12]-.47212e-6*y[12]^2+.15904e-9*y[12]^3-.16018e-13*y[12]^4)*R*Tr/F*(1.0-1.0*t1)/ln1*(y[13]-1.0*y[11])/y[12]-iapp
dy[79] = -4.0*sigman/ln1*(y[98]-1.0*y[96])-4.0*en^brugn*(.41253e-1+.5007e-3*y[13]-.47212e-6*y[13]^2+.15904e-9*y[13]^3-.16018e-13*y[13]^4)/ln1*(y[80]-1.0*y[78])+8.0*en^brugn*(.41253e-1+.5007e-3*y[13]-.47212e-6*y[13]^2+.15904e-9*y[13]^3-.16018e-13*y[13]^4)*R*Tr/F*(1.0-1.0*t1)/ln1*(y[14]-1.0*y[12])/y[13]-iapp
dy[80] = -4.0*sigman/ln1*(y[99]-1.0*y[97])-4.0*en^brugn*(.41253e-1+.5007e-3*y[14]-.47212e-6*y[14]^2+.15904e-9*y[14]^3-.16018e-13*y[14]^4)/ln1*(y[81]-1.0*y[79])+8.0*en^brugn*(.41253e-1+.5007e-3*y[14]-.47212e-6*y[14]^2+.15904e-9*y[14]^3-.16018e-13*y[14]^4)*R*Tr/F*(1.0-1.0*t1)/ln1*(y[15]-1.0*y[13])/y[14]-iapp
dy[81] = -4.0*sigman/ln1*(y[100]-1.0*y[98])-4.0*en^brugn*(.41253e-1+.5007e-3*y[15]-.47212e-6*y[15]^2+.15904e-9*y[15]^3-.16018e-13*y[15]^4)/ln1*(y[82]-1.0*y[80])+8.0*en^brugn*(.41253e-1+.5007e-3*y[15]-.47212e-6*y[15]^2+.15904e-9*y[15]^3-.16018e-13*y[15]^4)*R*Tr/F*(1.0-1.0*t1)/ln1*(y[16]-1.0*y[14])/y[15]-iapp
dy[82] = -4.0*sigman/ln1*(y[101]-1.0*y[99])-4.0*en^brugn*(.41253e-1+.5007e-3*y[16]-.47212e-6*y[16]^2+.15904e-9*y[16]^3-.16018e-13*y[16]^4)/ln1*(y[83]-1.0*y[81])+8.0*en^brugn*(.41253e-1+.5007e-3*y[16]-.47212e-6*y[16]^2+.15904e-9*y[16]^3-.16018e-13*y[16]^4)*R*Tr/F*(1.0-1.0*t1)/ln1*(y[17]-1.0*y[15])/y[16]-iapp
dy[83] = -4.0*sigman/ln1*(y[102]-1.0*y[100])-4.0*en^brugn*(.41253e-1+.5007e-3*y[17]-.47212e-6*y[17]^2+.15904e-9*y[17]^3-.16018e-13*y[17]^4)/ln1*(y[84]-1.0*y[82])+8.0*en^brugn*(.41253e-1+.5007e-3*y[17]-.47212e-6*y[17]^2+.15904e-9*y[17]^3-.16018e-13*y[17]^4)*R*Tr/F*(1.0-1.0*t1)/ln1*(y[63]-1.0*y[16])/y[17]-iapp
dy[84] = y[84]
dy[85] = 4.0*(-1.0*y[87]-3.0*y[85]+4.0*y[86])/lp+1.0*iapp/sigmap
dy[86] = 64.0*sigmap/lp^2*(y[85]-2.0*y[86]+y[87])-2.0*ap*F*kp*(y[1]*c0)^.5*(-1.0*y[110]*ctp+ctp)^.5*(y[110]*ctp)^.5*sinh(.5*F/R/Tr*(y[86]-1.0*y[65]-1.0*ocp_cathode(y[110],up)))
dy[87] = 64.0*sigmap/lp^2*(y[86]-2.0*y[87]+y[88])-2.0*ap*F*kp*(y[2]*c0)^.5*(-1.0*y[111]*ctp+ctp)^.5*(y[111]*ctp)^.5*sinh(.5*F/R/Tr*(y[87]-1.0*y[66]-1.0*ocp_cathode(y[111],up)))
dy[88] = 64.0*sigmap/lp^2*(y[87]-2.0*y[88]+y[89])-2.0*ap*F*kp*(y[3]*c0)^.5*(-1.0*y[112]*ctp+ctp)^.5*(y[112]*ctp)^.5*sinh(.5*F/R/Tr*(y[88]-1.0*y[67]-1.0*ocp_cathode(y[112],up)))
dy[89] = 64.0*sigmap/lp^2*(y[88]-2.0*y[89]+y[90])-2.0*ap*F*kp*(y[4]*c0)^.5*(-1.0*y[113]*ctp+ctp)^.5*(y[113]*ctp)^.5*sinh(.5*F/R/Tr*(y[89]-1.0*y[68]-1.0*ocp_cathode(y[113],up)))
dy[90] = 64.0*sigmap/lp^2*(y[89]-2.0*y[90]+y[91])-2.0*ap*F*kp*(y[5]*c0)^.5*(-1.0*y[114]*ctp+ctp)^.5*(y[114]*ctp)^.5*sinh(.5*F/R/Tr*(y[90]-1.0*y[69]-1.0*ocp_cathode(y[114],up)))
dy[91] = 64.0*sigmap/lp^2*(y[90]-2.0*y[91]+y[92])-2.0*ap*F*kp*(y[6]*c0)^.5*(-1.0*y[115]*ctp+ctp)^.5*(y[115]*ctp)^.5*sinh(.5*F/R/Tr*(y[91]-1.0*y[70]-1.0*ocp_cathode(y[115],up)))
dy[92] = 64.0*sigmap/lp^2*(y[91]-2.0*y[92]+y[93])-2.0*ap*F*kp*(y[7]*c0)^.5*(-1.0*y[116]*ctp+ctp)^.5*(y[116]*ctp)^.5*sinh(.5*F/R/Tr*(y[92]-1.0*y[71]-1.0*ocp_cathode(y[116],up)))
dy[93] = 4.0*(y[91]+3.0*y[93]-4.0*y[92])/lp
dy[94] = 4.0*(-1.0*y[96]-3.0*y[94]+4.0*y[95])/ln1
dy[95] = 64.0*sigman/ln1^2*(y[94]-2.0*y[95]+y[96])-2.0*an*F*kn*(y[11]*c0)^.5*(-1.0*y[124]*ctn+ctn)^.5*(y[124]*ctn)^.5*sinh(.5*F/R/Tr*(y[95]-1.0*y[77]+.5e-1-2.325*exp(-100.0*y[124]^1.15)+.1721*tanh(20.000000*y[124]-21.000000)+.25e-2*tanh(39.347962*y[124]-37.241379)+.34e-1*tanh(89.411765*y[124]-6.0294118)+.2e-2*tanh(7.0422535*y[124]-1.3661972)+.155e-1*tanh(77.519380*y[124]-8.1395349)))
dy[96] = 64.0*sigman/ln1^2*(y[95]-2.0*y[96]+y[97])-2.0*an*F*kn*(y[12]*c0)^.5*(-1.0*y[125]*ctn+ctn)^.5*(y[125]*ctn)^.5*sinh(.5*F/R/Tr*(y[96]-1.0*y[78]+.5e-1-2.325*exp(-100.0*y[125]^1.15)+.1721*tanh(20.000000*y[125]-21.000000)+.25e-2*tanh(39.347962*y[125]-37.241379)+.34e-1*tanh(89.411765*y[125]-6.0294118)+.2e-2*tanh(7.0422535*y[125]-1.3661972)+.155e-1*tanh(77.519380*y[125]-8.1395349)))
dy[97] = 64.0*sigman/ln1^2*(y[96]-2.0*y[97]+y[98])-2.0*an*F*kn*(y[13]*c0)^.5*(-1.0*y[126]*ctn+ctn)^.5*(y[126]*ctn)^.5*sinh(.5*F/R/Tr*(y[97]-1.0*y[79]+.5e-1-2.325*exp(-100.0*y[126]^1.15)+.1721*tanh(20.000000*y[126]-21.000000)+.25e-2*tanh(39.347962*y[126]-37.241379)+.34e-1*tanh(89.411765*y[126]-6.0294118)+.2e-2*tanh(7.0422535*y[126]-1.3661972)+.155e-1*tanh(77.519380*y[126]-8.1395349)))
dy[98] = 64.0*sigman/ln1^2*(y[97]-2.0*y[98]+y[99])-2.0*an*F*kn*(y[14]*c0)^.5*(-1.0*y[127]*ctn+ctn)^.5*(y[127]*ctn)^.5*sinh(.5*F/R/Tr*(y[98]-1.0*y[80]+.5e-1-2.325*exp(-100.0*y[127]^1.15)+.1721*tanh(20.000000*y[127]-21.000000)+.25e-2*tanh(39.347962*y[127]-37.241379)+.34e-1*tanh(89.411765*y[127]-6.0294118)+.2e-2*tanh(7.0422535*y[127]-1.3661972)+.155e-1*tanh(77.519380*y[127]-8.1395349)))
dy[99] = 64.0*sigman/ln1^2*(y[98]-2.0*y[99]+y[100])-2.0*an*F*kn*(y[15]*c0)^.5*(-1.0*y[128]*ctn+ctn)^.5*(y[128]*ctn)^.5*sinh(.5*F/R/Tr*(y[99]-1.0*y[81]+.5e-1-2.325*exp(-100.0*y[128]^1.15)+.1721*tanh(20.000000*y[128]-21.000000)+.25e-2*tanh(39.347962*y[128]-37.241379)+.34e-1*tanh(89.411765*y[128]-6.0294118)+.2e-2*tanh(7.0422535*y[128]-1.3661972)+.155e-1*tanh(77.519380*y[128]-8.1395349)))
dy[100] = 64.0*sigman/ln1^2*(y[99]-2.0*y[100]+y[101])-2.0*an*F*kn*(y[16]*c0)^.5*(-1.0*y[129]*ctn+ctn)^.5*(y[129]*ctn)^.5*sinh(.5*F/R/Tr*(y[100]-1.0*y[82]+.5e-1-2.325*exp(-100.0*y[129]^1.15)+.1721*tanh(20.000000*y[129]-21.000000)+.25e-2*tanh(39.347962*y[129]-37.241379)+.34e-1*tanh(89.411765*y[129]-6.0294118)+.2e-2*tanh(7.0422535*y[129]-1.3661972)+.155e-1*tanh(77.519380*y[129]-8.1395349)))
dy[101] = 64.0*sigman/ln1^2*(y[100]-2.0*y[101]+y[102])-2.0*an*F*kn*(y[17]*c0)^.5*(-1.0*y[130]*ctn+ctn)^.5*(y[130]*ctn)^.5*sinh(.5*F/R/Tr*(y[101]-1.0*y[83]+.5e-1-2.325*exp(-100.0*y[130]^1.15)+.1721*tanh(20.000000*y[130]-21.000000)+.25e-2*tanh(39.347962*y[130]-37.241379)+.34e-1*tanh(89.411765*y[130]-6.0294118)+.2e-2*tanh(7.0422535*y[130]-1.3661972)+.155e-1*tanh(77.519380*y[130]-8.1395349)))
dy[102] = 4.0*(y[100]+3.0*y[102]-4.0*y[101])/ln1+1.0*iapp/sigman
dy[103] = -2.0*(-1.0*y[25]-3.0*y[103]+4.0*y[18])/Rpp
dy[104] = -2.0*(-1.0*y[26]-3.0*y[104]+4.0*y[19])/Rpp
dy[105] = -2.0*(-1.0*y[27]-3.0*y[105]+4.0*y[20])/Rpp
dy[106] = -2.0*(-1.0*y[28]-3.0*y[106]+4.0*y[21])/Rpp
dy[107] = -2.0*(-1.0*y[29]-3.0*y[107]+4.0*y[22])/Rpp
dy[108] = -2.0*(-1.0*y[30]-3.0*y[108]+4.0*y[23])/Rpp
dy[109] = -2.0*(-1.0*y[31]-3.0*y[109]+4.0*y[24])/Rpp
dy[110] = 2.0*Dsp*(y[25]+3.0*y[110]-4.0*y[32])/Rpp+2.0*kp*(y[1]*c0)^.5*(-1.0*y[110]*ctp+ctp)^.5*(y[110]*ctp)^.5*sinh(.5*F/R/Tr*(y[86]-1.0*y[65]-1.0*ocp_cathode(y[110],up)))/ctp
dy[111] = 2.0*Dsp*(y[26]+3.0*y[111]-4.0*y[33])/Rpp+2.0*kp*(y[2]*c0)^.5*(-1.0*y[111]*ctp+ctp)^.5*(y[111]*ctp)^.5*sinh(.5*F/R/Tr*(y[87]-1.0*y[66]-1.0*ocp_cathode(y[111],up)))/ctp
dy[112] = 2.0*Dsp*(y[27]+3.0*y[112]-4.0*y[34])/Rpp+2.0*kp*(y[3]*c0)^.5*(-1.0*y[112]*ctp+ctp)^.5*(y[112]*ctp)^.5*sinh(.5*F/R/Tr*(y[88]-1.0*y[67]-1.0*ocp_cathode(y[112],up)))/ctp
dy[113] = 2.0*Dsp*(y[28]+3.0*y[113]-4.0*y[35])/Rpp+2.0*kp*(y[4]*c0)^.5*(-1.0*y[113]*ctp+ctp)^.5*(y[113]*ctp)^.5*sinh(.5*F/R/Tr*(y[89]-1.0*y[68]-1.0*ocp_cathode(y[113],up)))/ctp
dy[114] = 2.0*Dsp*(y[29]+3.0*y[114]-4.0*y[36])/Rpp+2.0*kp*(y[5]*c0)^.5*(-1.0*y[114]*ctp+ctp)^.5*(y[114]*ctp)^.5*sinh(.5*F/R/Tr*(y[90]-1.0*y[69]-1.0*ocp_cathode(y[114],up)))/ctp
dy[115] = 2.0*Dsp*(y[30]+3.0*y[115]-4.0*y[37])/Rpp+2.0*kp*(y[6]*c0)^.5*(-1.0*y[115]*ctp+ctp)^.5*(y[115]*ctp)^.5*sinh(.5*F/R/Tr*(y[91]-1.0*y[70]-1.0*ocp_cathode(y[115],up)))/ctp
dy[116] = 2.0*Dsp*(y[31]+3.0*y[116]-4.0*y[38])/Rpp+2.0*kp*(y[7]*c0)^.5*(-1.0*y[116]*ctp+ctp)^.5*(y[116]*ctp)^.5*sinh(.5*F/R/Tr*(y[92]-1.0*y[71]-1.0*ocp_cathode(y[116],up)))/ctp
dy[117] = -2.0*(-1.0*y[46]-3.0*y[117]+4.0*y[39])/Rpn
dy[118] = -2.0*(-1.0*y[47]-3.0*y[118]+4.0*y[40])/Rpn
dy[119] = -2.0*(-1.0*y[48]-3.0*y[119]+4.0*y[41])/Rpn
dy[120] = -2.0*(-1.0*y[49]-3.0*y[120]+4.0*y[42])/Rpn
dy[121] = -2.0*(-1.0*y[50]-3.0*y[121]+4.0*y[43])/Rpn
dy[122] = -2.0*(-1.0*y[51]-3.0*y[122]+4.0*y[44])/Rpn
dy[123] = -2.0*(-1.0*y[52]-3.0*y[123]+4.0*y[45])/Rpn
dy[124] = 2.0*Dsn*(y[46]+3.0*y[124]-4.0*y[53])/Rpn+2.0*kn*(y[11]*c0)^.5*(-1.0*y[124]*ctn+ctn)^.5*(y[124]*ctn)^.5*sinh(.5*F/R/Tr*(y[95]-1.0*y[77]+.5e-1-2.325*exp(-100.0*y[124]^1.15)+.1721*tanh(20.000000*y[124]-21.000000)+.25e-2*tanh(39.347962*y[124]-37.241379)+.34e-1*tanh(89.411765*y[124]-6.0294118)+.2e-2*tanh(7.0422535*y[124]-1.3661972)+.155e-1*tanh(77.519380*y[124]-8.1395349)))/ctn
dy[125] = 2.0*Dsn*(y[47]+3.0*y[125]-4.0*y[54])/Rpn+2.0*kn*(y[12]*c0)^.5*(-1.0*y[125]*ctn+ctn)^.5*(y[125]*ctn)^.5*sinh(.5*F/R/Tr*(y[96]-1.0*y[78]+.5e-1-2.325*exp(-100.0*y[125]^1.15)+.1721*tanh(20.000000*y[125]-21.000000)+.25e-2*tanh(39.347962*y[125]-37.241379)+.34e-1*tanh(89.411765*y[125]-6.0294118)+.2e-2*tanh(7.0422535*y[125]-1.3661972)+.155e-1*tanh(77.519380*y[125]-8.1395349)))/ctn
dy[126] = 2.0*Dsn*(y[48]+3.0*y[126]-4.0*y[55])/Rpn+2.0*kn*(y[13]*c0)^.5*(-1.0*y[126]*ctn+ctn)^.5*(y[126]*ctn)^.5*sinh(.5*F/R/Tr*(y[97]-1.0*y[79]+.5e-1-2.325*exp(-100.0*y[126]^1.15)+.1721*tanh(20.000000*y[126]-21.000000)+.25e-2*tanh(39.347962*y[126]-37.241379)+.34e-1*tanh(89.411765*y[126]-6.0294118)+.2e-2*tanh(7.0422535*y[126]-1.3661972)+.155e-1*tanh(77.519380*y[126]-8.1395349)))/ctn
dy[127] = 2.0*Dsn*(y[49]+3.0*y[127]-4.0*y[56])/Rpn+2.0*kn*(y[14]*c0)^.5*(-1.0*y[127]*ctn+ctn)^.5*(y[127]*ctn)^.5*sinh(.5*F/R/Tr*(y[98]-1.0*y[80]+.5e-1-2.325*exp(-100.0*y[127]^1.15)+.1721*tanh(20.000000*y[127]-21.000000)+.25e-2*tanh(39.347962*y[127]-37.241379)+.34e-1*tanh(89.411765*y[127]-6.0294118)+.2e-2*tanh(7.0422535*y[127]-1.3661972)+.155e-1*tanh(77.519380*y[127]-8.1395349)))/ctn
dy[128] = 2.0*Dsn*(y[50]+3.0*y[128]-4.0*y[57])/Rpn+2.0*kn*(y[15]*c0)^.5*(-1.0*y[128]*ctn+ctn)^.5*(y[128]*ctn)^.5*sinh(.5*F/R/Tr*(y[99]-1.0*y[81]+.5e-1-2.325*exp(-100.0*y[128]^1.15)+.1721*tanh(20.000000*y[128]-21.000000)+.25e-2*tanh(39.347962*y[128]-37.241379)+.34e-1*tanh(89.411765*y[128]-6.0294118)+.2e-2*tanh(7.0422535*y[128]-1.3661972)+.155e-1*tanh(77.519380*y[128]-8.1395349)))/ctn
dy[129] = 2.0*Dsn*(y[51]+3.0*y[129]-4.0*y[58])/Rpn+2.0*kn*(y[16]*c0)^.5*(-1.0*y[129]*ctn+ctn)^.5*(y[129]*ctn)^.5*sinh(.5*F/R/Tr*(y[100]-1.0*y[82]+.5e-1-2.325*exp(-100.0*y[129]^1.15)+.1721*tanh(20.000000*y[129]-21.000000)+.25e-2*tanh(39.347962*y[129]-37.241379)+.34e-1*tanh(89.411765*y[129]-6.0294118)+.2e-2*tanh(7.0422535*y[129]-1.3661972)+.155e-1*tanh(77.519380*y[129]-8.1395349)))/ctn
dy[130] = 2.0*Dsn*(y[52]+3.0*y[130]-4.0*y[59])/Rpn+2.0*kn*(y[17]*c0)^.5*(-1.0*y[130]*ctn+ctn)^.5*(y[130]*ctn)^.5*sinh(.5*F/R/Tr*(y[101]-1.0*y[83]+.5e-1-2.325*exp(-100.0*y[130]^1.15)+.1721*tanh(20.000000*y[130]-21.000000)+.25e-2*tanh(39.347962*y[130]-37.241379)+.34e-1*tanh(89.411765*y[130]-6.0294118)+.2e-2*tanh(7.0422535*y[130]-1.3661972)+.155e-1*tanh(77.519380*y[130]-8.1395349)))/ctn
dy[131] = y[131]-y[85]+y[102]+ (current*Kr)
end
function ocp_cathode(theta_p,up)
Up = (-5.057-13.6*theta_p^2+121.6*theta_p^4-185.1*theta_p^6-45.43*theta_p^8+127.7*theta_p^10)/(-1-5.04*theta_p^2+32.3*theta_p^4-40.95*theta_p^6-25.94*theta_p^8+40.63*theta_p^10)
end
yyy = function(pars::Array{Float64})
socn=pars[24]
socp=pars[25]
y0 = [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,socp,socp,socp,socp,socp,socp,socp,socp,socp,socp,socp,socp,socp,socp,socp,socp,socp,socp,socp,socp,socp,socn,socn,socn,socn,socn,socn,socn,socn,socn,socn,socn,socn,socn,socn,socn,socn,socn,socn,socn,socn,socn,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,(-0.5057e1 - 0.136e2 * socp ^ 2 + 0.1216e3 * socp ^ 4 - 0.1851e3 * socp ^ 6 - 0.4543e2 * socp ^ 8 + 0.1277e3 * socp ^ 10) / (-0.1e1 - 0.504e1 * socp ^ 2 + 0.323e2 * socp ^ 4 - 0.4095e2 * socp ^ 6 - 0.2594e2 * socp ^ 8 + 0.4063e2 * socp ^ 10),(-0.5057e1 - 0.136e2 * socp ^ 2 + 0.1216e3 * socp ^ 4 - 0.1851e3 * socp ^ 6 - 0.4543e2 * socp ^ 8 + 0.1277e3 * socp ^ 10) / (-0.1e1 - 0.504e1 * socp ^ 2 + 0.323e2 * socp ^ 4 - 0.4095e2 * socp ^ 6 - 0.2594e2 * socp ^ 8 + 0.4063e2 * socp ^ 10),(-0.5057e1 - 0.136e2 * socp ^ 2 + 0.1216e3 * socp ^ 4 - 0.1851e3 * socp ^ 6 - 0.4543e2 * socp ^ 8 + 0.1277e3 * socp ^ 10) / (-0.1e1 - 0.504e1 * socp ^ 2 + 0.323e2 * socp ^ 4 - 0.4095e2 * socp ^ 6 - 0.2594e2 * socp ^ 8 + 0.4063e2 * socp ^ 10),(-0.5057e1 - 0.136e2 * socp ^ 2 + 0.1216e3 * socp ^ 4 - 0.1851e3 * socp ^ 6 - 0.4543e2 * socp ^ 8 + 0.1277e3 * socp ^ 10) / (-0.1e1 - 0.504e1 * socp ^ 2 + 0.323e2 * socp ^ 4 - 0.4095e2 * socp ^ 6 - 0.2594e2 * socp ^ 8 + 0.4063e2 * socp ^ 10),(-0.5057e1 - 0.136e2 * socp ^ 2 + 0.1216e3 * socp ^ 4 - 0.1851e3 * socp ^ 6 - 0.4543e2 * socp ^ 8 + 0.1277e3 * socp ^ 10) / (-0.1e1 - 0.504e1 * socp ^ 2 + 0.323e2 * socp ^ 4 - 0.4095e2 * socp ^ 6 - 0.2594e2 * socp ^ 8 + 0.4063e2 * socp ^ 10),(-0.5057e1 - 0.136e2 * socp ^ 2 + 0.1216e3 * socp ^ 4 - 0.1851e3 * socp ^ 6 - 0.4543e2 * socp ^ 8 + 0.1277e3 * socp ^ 10) / (-0.1e1 - 0.504e1 * socp ^ 2 + 0.323e2 * socp ^ 4 - 0.4095e2 * socp ^ 6 - 0.2594e2 * socp ^ 8 + 0.4063e2 * socp ^ 10),(-0.5057e1 - 0.136e2 * socp ^ 2 + 0.1216e3 * socp ^ 4 - 0.1851e3 * socp ^ 6 - 0.4543e2 * socp ^ 8 + 0.1277e3 * socp ^ 10) / (-0.1e1 - 0.504e1 * socp ^ 2 + 0.323e2 * socp ^ 4 - 0.4095e2 * socp ^ 6 - 0.2594e2 * socp ^ 8 + 0.4063e2 * socp ^ 10),(-0.5057e1 - 0.136e2 * socp ^ 2 + 0.1216e3 * socp ^ 4 - 0.1851e3 * socp ^ 6 - 0.4543e2 * socp ^ 8 + 0.1277e3 * socp ^ 10) / (-0.1e1 - 0.504e1 * socp ^ 2 + 0.323e2 * socp ^ 4 - 0.4095e2 * socp ^ 6 - 0.2594e2 * socp ^ 8 + 0.4063e2 * socp ^ 10),(-0.5057e1 - 0.136e2 * socp ^ 2 + 0.1216e3 * socp ^ 4 - 0.1851e3 * socp ^ 6 - 0.4543e2 * socp ^ 8 + 0.1277e3 * socp ^ 10) / (-0.1e1 - 0.504e1 * socp ^ 2 + 0.323e2 * socp ^ 4 - 0.4095e2 * socp ^ 6 - 0.2594e2 * socp ^ 8 + 0.4063e2 * socp ^ 10),-0.5e-1 + 0.2325e1 * exp(-0.100e3 * socn ^ 0.115e1) - 0.1721e0 * tanh(0.20000000e2 * socn - 0.21000000e2) - 0.25e-2 * tanh(0.39347962e2 * socn - 0.37241379e2) - 0.34e-1 * tanh(0.89411765e2 * socn - 0.60294118e1) - 0.2e-2 * tanh(0.70422535e1 * socn - 0.13661972e1) - 0.155e-1 * tanh(0.77519380e2 * socn - 0.81395349e1),-0.5e-1 + 0.2325e1 * exp(-0.100e3 * socn ^ 0.115e1) - 0.1721e0 * tanh(0.20000000e2 * socn - 0.21000000e2) - 0.25e-2 * tanh(0.39347962e2 * socn - 0.37241379e2) - 0.34e-1 * tanh(0.89411765e2 * socn - 0.60294118e1) - 0.2e-2 * tanh(0.70422535e1 * socn - 0.13661972e1) - 0.155e-1 * tanh(0.77519380e2 * socn - 0.81395349e1),-0.5e-1 + 0.2325e1 * exp(-0.100e3 * socn ^ 0.115e1) - 0.1721e0 * tanh(0.20000000e2 * socn - 0.21000000e2) - 0.25e-2 * tanh(0.39347962e2 * socn - 0.37241379e2) - 0.34e-1 * tanh(0.89411765e2 * socn - 0.60294118e1) - 0.2e-2 * tanh(0.70422535e1 * socn - 0.13661972e1) - 0.155e-1 * tanh(0.77519380e2 * socn - 0.81395349e1),-0.5e-1 + 0.2325e1 * exp(-0.100e3 * socn ^ 0.115e1) - 0.1721e0 * tanh(0.20000000e2 * socn - 0.21000000e2) - 0.25e-2 * tanh(0.39347962e2 * socn - 0.37241379e2) - 0.34e-1 * tanh(0.89411765e2 * socn - 0.60294118e1) - 0.2e-2 * tanh(0.70422535e1 * socn - 0.13661972e1) - 0.155e-1 * tanh(0.77519380e2 * socn - 0.81395349e1),-0.5e-1 + 0.2325e1 * exp(-0.100e3 * socn ^ 0.115e1) - 0.1721e0 * tanh(0.20000000e2 * socn - 0.21000000e2) - 0.25e-2 * tanh(0.39347962e2 * socn - 0.37241379e2) - 0.34e-1 * tanh(0.89411765e2 * socn - 0.60294118e1) - 0.2e-2 * tanh(0.70422535e1 * socn - 0.13661972e1) - 0.155e-1 * tanh(0.77519380e2 * socn - 0.81395349e1),-0.5e-1 + 0.2325e1 * exp(-0.100e3 * socn ^ 0.115e1) - 0.1721e0 * tanh(0.20000000e2 * socn - 0.21000000e2) - 0.25e-2 * tanh(0.39347962e2 * socn - 0.37241379e2) - 0.34e-1 * tanh(0.89411765e2 * socn - 0.60294118e1) - 0.2e-2 * tanh(0.70422535e1 * socn - 0.13661972e1) - 0.155e-1 * tanh(0.77519380e2 * socn - 0.81395349e1),-0.5e-1 + 0.2325e1 * exp(-0.100e3 * socn ^ 0.115e1) - 0.1721e0 * tanh(0.20000000e2 * socn - 0.21000000e2) - 0.25e-2 * tanh(0.39347962e2 * socn - 0.37241379e2) - 0.34e-1 * tanh(0.89411765e2 * socn - 0.60294118e1) - 0.2e-2 * tanh(0.70422535e1 * socn - 0.13661972e1) - 0.155e-1 * tanh(0.77519380e2 * socn - 0.81395349e1),-0.5e-1 + 0.2325e1 * exp(-0.100e3 * socn ^ 0.115e1) - 0.1721e0 * tanh(0.20000000e2 * socn - 0.21000000e2) - 0.25e-2 * tanh(0.39347962e2 * socn - 0.37241379e2) - 0.34e-1 * tanh(0.89411765e2 * socn - 0.60294118e1) - 0.2e-2 * tanh(0.70422535e1 * socn - 0.13661972e1) - 0.155e-1 * tanh(0.77519380e2 * socn - 0.81395349e1),-0.5e-1 + 0.2325e1 * exp(-0.100e3 * socn ^ 0.115e1) - 0.1721e0 * tanh(0.20000000e2 * socn - 0.21000000e2) - 0.25e-2 * tanh(0.39347962e2 * socn - 0.37241379e2) - 0.34e-1 * tanh(0.89411765e2 * socn - 0.60294118e1) - 0.2e-2 * tanh(0.70422535e1 * socn - 0.13661972e1) - 0.155e-1 * tanh(0.77519380e2 * socn - 0.81395349e1),socp,socp,socp,socp,socp,socp,socp,socp,socp,socp,socp,socp,socp,socp,socn,socn,socn,socn,socn,socn,socn,socn,socn,socn,socn,socn,socn,socn,(-0.5057e1 - 0.136e2 * socp ^ 2 + 0.1216e3 * socp ^ 4 - 0.1851e3 * socp ^ 6 - 0.4543e2 * socp ^ 8 + 0.1277e3 * socp ^ 10) / (-0.1e1 - 0.504e1 * socp ^ 2 + 0.323e2 * socp ^ 4 - 0.4095e2 * socp ^ 6 - 0.2594e2 * socp ^ 8 + 0.4063e2 * socp ^ 10) + 0.5e-1 - 0.2325e1 * exp(-0.100e3 * socn ^ 0.115e1) + 0.1721e0 * tanh(0.20000000e2 * socn - 0.21000000e2) + 0.25e-2 * tanh(0.39347962e2 * socn - 0.37241379e2) + 0.34e-1 * tanh(0.89411765e2 * socn - 0.60294118e1) + 0.2e-2 * tanh(0.70422535e1 * socn - 0.13661972e1) + 0.155e-1 * tanh(0.77519380e2 * socn - 0.81395349e1)]
end
y0=zeros(131)
n=length(y0)
MM=zeros((n,n))
node = 59
for i in 1:node
MM[i,i]=1.
end
pars = [7.33234842e-08,6.40146313e-12,1.56954148e-13,1.65301662e-06,
1.28277079e-05,1.39115369e+00,1.70413971e+00,1.93944836e+00,
1.15245117e+03,3.46144727e+04,4.02951641e+04,3.24118063e-02,
2.70097144e-02,5.20209372e-01,3.74625146e-01,4.16417062e-01,
5.51529622e-09,6.95793778e-09,4.91284482e-05,4.88584483e-05,
1.28054144e-05,1.38146591e+02,9.15019798e+00,9.78830695e-01,
4.58723247e-01,4.32648301e-01,0.02]
y0 = yyy(pars)
dy = zeros(n)
f = (dy,y) -> P2D2c(zero(y),y,pars,0.)
using ForwardDiff
# Really long compile time
ForwardDiff.jacobian(f,dy,y0)
This comes up quite a lot in differential equation work since writing (or generating) complex nonlinear functions is not that uncommon. @MartinOtter mentioned that Modia.jl sometimes hits this issue.
ChrisRackauckas
on 25 Aug 2018
julia> @time f(rand(50))
0.017597 seconds (18.85 k allocations: 824.408 KiB)
0.5138156603313624
julia> @time f(rand(50))
0.000462 seconds (5.20 k allocations: 119.625 KiB)
0.5487568472487016
vtjnash
on 15 Mar 2020
👍1
Should it be added to some kind of compiler benchmarks collection so that it doesn't regress?
StefanKarpinski
on 16 Mar 2020
Was this page helpful?
0 / 5 - 0 ratings
Related issues
helgee
·
3Comments
i-apellaniz
·
3Comments
TotalVerb
·
3Comments
sbromberger
·
3Comments
Keno
·
3Comments
Most helpful comment
Now that the new optimizer is in, I gave this another whirl (script above slightly updated):
julia 0.6.3:
julia 0.7.0-beta.129: still waiting after 20 minutes...
I trimmed the file size down to about 10% of the full size and collected this profile. The biggest offender by far seems to be
construct_domtree!, specifically its usage ofreduce. Hope that's useful.