# DifferentialEquations -- what goes wrong?

I have solved a PDE successfully in OpenModelica, but I have problems with solving it in Julia. The Julia code seems to work for up to 8 slices (control volumes), but mysteriously collapses for 9 and more slices…

I’m considering a model of the temperature in a heated water tank, possibly distributed in the z direction. Water may flow through the tank from the bottom (advection), there is convective heat input from the ambient, heater at a specific height (“mathematical point”), and come quasi buoyancy effect.

I’ve already solved the model in OpenModelica, and the results make sense – at least up to 120 slices. Now, I’m trying to solve the model in Julia/DifferentialEquations. In the code (below), I can vary the number of discretization slices from 1 to whatever. The results from OpenModelica and Julia seem to be identical for N=5 slices. For OpenModelica:

and for Julia (sorry about different colors in lines, as well as different number of curves):

However, with N=10 slices, the Julia model collapses. First, OpenModelica results:

followed by the collapsed result from Julia:

OpenModelica gives sensible results up to at least N=120 (i.e., that is how far I have tested it). Julia appears to give ok results up to N=8, but collapses at N=9.

I have probably done a mistake in my code…??? If anyone sees the problem, I appreciate suggestions. Also, my code is probably not super efficient wrt. memory allocation and speed. The relevant code is included below.

# -- Importing packages
using Plots; pyplot()
using LaTeXStrings
using DifferentialEquations
# -- Defining constants
LW1 = 2.5
LW2 = 1.5
LW3 = 1.0
LS1 = :solid
LS2 = :dot
LS3 = :dashdot;
LA1 = 1.0;
LA2 = 0.6;
LA3 = 0.2;
#
# -- model function
function watertank(dT,T,p,t,u)
# Parameters
N = length(T)
# -- mass flow through tank
md = p.rho*u.VdL(t)*u.uv(t)
#
if N == 1
As = 2*p.A + p.per*p.h
V = p.A*p.h
m = p.rho*V
Cp = m*p.chp
Qd_el = p.Po*u.uP(t)
Qd_a = p.Ux*As*(u.Ta(t)-T)
dT = (md*p.chp*(u.Ti(t)-T) + Qd_el + Qd_a)/Cp
else
# N > 1
dz = p.h/N
V = p.A*dz
m = p.rho*V
Cp = m*p.chp
#
# -- heater
Qd_el = zeros(N)
#
jP = Int(ceil(p.zP/dz))
Qd_el[jP] = p.Po*u.uP(t)
#
# -- heat from ambient
Qd_a = zeros(N)
#
Qd_a = p.Ux*(p.per*dz+p.A)*(u.Ta(t)-T)
Qd_a[2:N-1] = p.Ux*p.per*dz*(u.Ta(t).-T[2:N-1])
Qd_a[N] = p.Ux*(p.per*dz+p.A)*(u.Ta(t)-T[N])
#
# -- heat flow by diffusion
kb = zeros(N-1)
for j in 1:(N-1)
kb[j] = p.rho*p.chp*p.kappa^2*p.d^2*sqrt(p.g*p.alpha_p*max(-(T[j+1]-T[j]),0)/dz)
end
#
Qd_d = zeros(N)
#
Qd_d = p.A*(p.kt+p.cb*kb)*(T-T)/dz
Qd_d[2:N-1] = p.A*(p.kt.+p.cb*kb[2:N-1]).*(T[3:N]-T[2:N-1])/dz
-p.A*(p.kt.+p.cb*kb[1:N-2]).*(T[2:N-1]-T[1:N-2])/dz
Qd_d[N] = -p.A*(p.kt+p.cb*kb[N-1])*(T[N]-T[N-1])/dz
#
# -- heat total
Qd = Qd_el + Qd_a + Qd_d
#
# -- temperature equations
dT = (md*p.chp*(u.Ti(t)-T) + Qd)/Cp
dT[2:N] = (md*p.chp*(T[1:N-1]-T[2:N]) + Qd[2:N])/Cp
end
#
end
#
# -- function for computing output
function Ti_loop(Te,Ti,uv)
return uv*Te + (1-uv)*Ti
end
;
# Model parameters
# -- constants
g = 9.81y
kappa = 0.41
# -- geometry
h = 1.5
d = 0.5
per = pi*d
A = pi*d^2/4
# -- equipment location
zP = 1.15
z1s = 1.3
z2s = 0.23
# -- thermodynamics
rho = 1.0e3
chp = 4.19e3y
alpha_p = 303e-6
# -- heat transport
kt = 0.6
ha = 3
hw = 50
ki_di = 0.5
Ux = 1/(1/ha + 1/hw + 1/ki_di)
# -- heat source
Po = 15e3
#
p = (g=g,kappa=kappa,h=h,d=d,per=per,A=A,zP=zP,rho=rho,chp=chp,alpha_p=alpha_p,kt=kt,cb=cb,Ux=Ux,Po=Po);
#
# Input functions
#
t_stop = 3600*4.;
#
uP_i = t->  t < t_stop/5 ? 1.0 : 0
uv_i = t->  begin
if t < t_stop/6
return 0.75
elseif t > t_stop/2
return 0.25
else
return 0
end
end
VdL_i = t->  t < t_stop/8 ? 8e-3/60 : 4e-3/60
#
Ta_i = t->  25.
Ti_i = t->  t < t_stop/3 ? 28. : 20.
#
u = (uP=uP_i,uv=uv_i,Ti=Ti_i,Ta=Ta_i,VdL=VdL_i);
# -- creating function into form suitable for ODEProblem
wt_i = (dx,x,p,t) -> watertank(dx,x,p,t,u);
#
# Setting up Simulations
#
# -- time span
tspan = (0.0,t_stop)
#
# --number of slices
nx = 8
#
# -- initial value of states
T0 = fill(25.,nx)
#
# -- defining and solving problem
prob = ODEProblem(wt_i,T0,tspan,p)
sol = solve(prob,Tsit5(),saveat=60)
#
# -- post processing
Sol = reduce(hcat,sol.u)
Tm_ii_L = u.Ti.(sol.t).*(1.0.-u.uv.(sol.t)) .+ Sol[nx,:].*u.uv.(sol.t)
#
dz = h/nx # slice thickness
j1s = Int(ceil(z1s/dz))
j2s = Int(ceil(z2s/dz))
jP = Int(ceil(p.zP/dz))
#
# PLOTTING RESULTS
#
col = range(RGB(0.6,0.6,1),stop=RGB(1,0.7,0.7),length=nx)
#
plot()
for j in 1:nx
plot!(sol.t/3600,Sol[j,:],lw=LW2,lc=col[j],la=0.6,label="")
end
# Model temperature range
plot!(sol.t[1:5:end]/3600,Sol[nx,1:5:end],st=:scatter,m=:x,mc=col[nx],msc=col[nx],ms=4,label=L"\hat{T}_\mathrm{e}")
plot!(sol.t[1:5:end]/3600,Sol[1,1:5:end],st=:scatter,m=:o,mc=col,msc=:blue,ms=4,label=L"\hat{T}_1")
plot!(sol.t[1:5:end]/3600,Sol[jP,1:5:end],st=:scatter,m=:v,mc=col[jP],msc=col[jP],ms=4,label=L"\hat{T}(z_\mathrm{P})")
# Sensor temperatures
plot!(sol.t/3600,Sol[j1s,:],lw=LW0,ls=LS2,lc=:maroon,label=L"\hat{T}_1^\mathrm{s}")
plot!(sol.t/3600,Sol[j2s,:],lw=LW0,ls=LS2,lc=:orange,label=L"\hat{T}_2^\mathrm{s}")
plot!(sol.t/3600,u.Ti.(sol.t),lc=:blue,ls=LS1,lw=LW1,la=0.4,label=L"T_\mathrm{i}")
plot!(sol.t/3600,Tm_ii_L,lw=LW0,lc=:coral,ls=LS2,label=L"\hat{T}_\mathrm{i}^\ell")
#
plot!(ylabel=L"temperature $T$ [$^\circ$C]")
plot!(xlim=tspan./3600,legend=:topright,box=true)

1 Like

To isolate that it is your implementation, make sure to try multiple solver libraries at lower tolerances. If everything gives the same answer, you can be sure it’s the implementation of the model. Exactly why is something that would be hard to quickly parse from the code you have there.

1 Like

Move the - to the previous line. It won’t properly parse as one expression like this. See the 7 Julia Gotchas as to why.

1 Like

Good idea. I had already isolated the problem to be related to the diffusive heat flow by in turn setting terms to zero and comparing with OpenModelica.

Argh… yes, of course. Thanks!!

What bugs me is that it worked for N<=8… or maybe the diffusion terms were small for N <= 8…

Now, I need to try to make the code more efficient. With N=120 (or nx = 120…), it took quite some while to solve the problem. I’m now trying to solve it for N=300.

Possible reasons why it becomes slower:

• massive memory usage due to inefficient code
• the code gets stiff when N increases. From the method of separation of variables, I seem to recall that the smallest eigenvalue is roughly proportional to N^2…

Thus, at some stage, I probably should switch to a stiff solver…
(Yup. With N=300, I get

 Warning: Interrupted. Larger maxiters is needed.


Yes, there’s a lot of unnecessary allocations in there. You can just broadcast those 2:N lines.

The other thing you might want to do is set tstops to handle those discontinuities.

1 Like

Should I also…

• transfer a temporary array as a parameter, and use that as a working array to avoid allocating new memory inside the function?
• if so, can I use “view” to give that working array a name that is easy to work with?

tstops" and discontinuities… Are you referring to the step changes in the inputs? Modelica has very good event detection, and figures out when there are step changes to avoid the discretization step stepping over the changes --> inaccuracies.

With N=200, a solution is found. But some of the solutions start to get jittery:

So I should try to use some stiff solver, perhaps.

Modelica probably isn’t using event detection there. It’s probably looking at the way you defined your problem and giving the solver a tstop. That’s an advantage of a DSL, but if you’re calling the solver directly you’ll need to do what Modelica is doing in the backend.

Yes

1 Like

I tried several of the stiff solvers… all of them crashed:

MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardD...


It looks like these methods use AD/ForwardDiff. If so, that may be a problem, since my model contains square root functions, the max function, etc. Is there a method that doesn’t use AD?
[I even tried the Kvaerno methods… I suspect the developer of these methods were a class mate of mine :-). ]

On tstops… I found a workbook of yours at http://tutorials.juliadiffeq.org/html/exercises/02-workshop_solutions.html, but I don’t quite understand it. Here is what I tried:

tstops = [t_stop/8,t_stop/6,t_stop/5,t_stop/3,t_stop/2]
condition(u,t,integrator) = t in tstops
affect!(integrator) = (integrator.u += 100)
cb = DiscreteCallback(condition,affect!)
#
sol = @btime solve(prob,Tsit5(),callback=cb, tstops=tstops,saveat=60,reltol=1e-6);


This doesn’t work properly, obviously because I don’t quite understand it.

I have assumed that vector tstops contains the time points where there are jumps in my inputs.

What I don’t understand is the affect! function. Does this function inject jumps of 100 in state no. 1? That is not what I want. I just want the jumps to come from the input functions…

Just add the tstops. Don’t add events with it…

KenCarp4(autodiff=false) etc., though you can also just fix the caches.

1 Like

Thanks. I have now (i) made the code more efficient by reducing memory allocation to ca. 10% of the original use, (ii) set up tstops [I found the documentation], and (iii) setting autodiff=false for a couple of stiff solvers.

For N=20 slices, the KenCarp4 algorithm was slightly faster than Tsit5 and used some 10% of the memory allocation. But with N=50, the memory allocation for the two was more or less the same. For N=200, Tsit5 finds the solution in some 30 s. But KenCarp4` takes “forever”. This was on my 16GB RAM laptop.

OK – I’ll test it tomorrow with my 32GB RAM desktop, and test more stiff solvers.

Your problem is probably sparse? Check the docs on how to handle that. A big dense Jacobian is slow to factorize

1 Like