Starting from MATLAB code, how do I use ode45 to solve Navier-Stokes in Julia?

New to Julia
, ,
1

I’m a new Julia user but matlab expert(I am tired of paying them).

I want to call back certain values on my ode during solving. I have watched the tutorial from Chris https://www.youtube.com/watch?v=KPEqYtEd-zY but not clear.

In matlab, I would do it as follows:

function [eqn,Qgap,Force,Fr]=navier(t,y)

%parameters
rho=810;
beta=5e8;
r1=0.01;
r2=0.004;
Ap=pi*(0.01^2-0.004^2);
L=0.06;
a=0.01;
f=10;
l=0.01;
Aeff=pi*2*r1*l;
delta=0.007;
N=20;
dx=delta/N;
rad=linspace(r1,r1+delta,N);

vel=y(1:N);
p1=y(N+1);
p2=y(N+2);

x=a*sin(2*pi*f*t);
u=a*2*pi*f*cos(2*pi*f*t);

eqn=zeros(N+2,1);
vel(1)=u;
vel(end)=0;
dudx=[(vel(2)-vel(1))/dx;(vel(3:end) -vel(1:end-2))/(2*dx);(vel(end)-vel(end-1))/dx];
dudx2=(vel(3:end)-2*vel(2:end-1)+vel(1:end-2))/(dx*dx);
dudx2 = [(dudx(1)-dudx(2))/dx;dudx2;(dudx(end-1)-dudx(end))/dx];
mu=[630 * (1+(0.05*(dudx)).^2).^(-0.23)];
pf=12*a*0.707*u*cos(2*pi*f*t)/(delta^2)*630*(1+(0.05*(0.707*u/delta)^2))^(-0.23);

Qgap =simpsons(2*pi*transpose(rad).*vel,r1,r1+delta);
% Qgap=trapz(rad,2*pi*transpose(rad).*vel);

eqn(1)=(p1-p2+pf)/l + ((mu(1))/rho)*dudx2(1);
eqn(2:N-1)=(p1-p2+pf)/l + ((mu(2:N-1))./rho).*dudx2(2:N-1);
eqn(N)=(p1-p2+pf)/l + ((mu(N))/rho)*dudx2(N);
eqn(N+1)=(beta/(Ap*(L-l-x)))*(-Qgap+Ap*u);
eqn(N+2)=(beta/(Ap*(L-l+x)))*(Qgap-Ap*u);

Force=(p1-p2+pf)*Ap;
Fr=mu(1)*Aeff*dudx(1);

    
end

In matlab, we just have to put the variables as outputs after the eqn, then write a callback function (a bit complicated version of my own)

[~,Qgap,Force,Fr] = cellfun(@(tsol,ysol)  navier(tsol,ysol.'), num2cell(tsol), num2cell(ysol,2),'uni',0);
Qgap=cell2mat(Qgap);
Force=cell2mat(Force);
Fr=cell2mat(Fr);

where as it gives me back a cell of the values with the time points, and then convert to a matlab simple array.

Can you help me achieve this in Julia?

2

You’re just trying to save a few extra variables of interest, right? You should be able to handle that with a SavingCallback. Do you have a working Julia version of the code without callbacks?

3

No I don’t have it. I was looking for every option of conversion from Matlab to Julia so I can do it smoothly.
If you have got free time and alot of coding interest, please fell free.

4

It would be easier to help you if you could attempt some of the conversion yourself. For example, I think you could have ported navier for us and maybe asked how to do the numerical quadrature.

Nonetheless, I’m feeling generous. I’m attempted to port navier. I did not think too much about it, so there is likely an error:

using NumericalIntegration

function navier(t,y)

#parameters
rho=810
beta=5e8
r1=0.01
r2=0.004
Ap=pi*(0.01^2-0.004^2)
L=0.06
a=0.01
f=10
l=0.01
Aeff=pi*2*r1*l
delta=0.007
N=20
dx=delta/N
rad=range(r1,r1+delta,length=N)

vel=y[1:N]
p1=y[N+1]
p2=y[N+2]

x=a*sin(2*pi*f*t)
u=a*2*pi*f*cos(2*pi*f*t)

eqn=zeros(N+2,1)
vel[1]=u
vel[end]=0
dudx=[(vel[2]-vel[1])/dx;(vel[3:end] -vel[1:end-2])/(2*dx);(vel[end]-vel[end-1])/dx]
dudx2=(vel[3:end]-2*vel[2:end-1].+vel[1:end-2])/(dx*dx)
dudx2 = [(dudx[1]-dudx[2])/dx;dudx2;(dudx[end-1]-dudx[end])/dx]
mu=630 * (1 .+ (0.05*(dudx)).^2).^(-0.23)
pf=12*a*0.707*u*cos(2*pi*f*t)/(delta^2)*630*(1+(0.05*(0.707*u/delta)^2))^(-0.23)

#Qgap =simpsons(2*pi*transpose(rad).*vel,r1,r1+delta)
# Qgap=trapz(rad,2*pi*transpose(rad).*vel)
Qgap = integrate(rad, 2*pi*rad.*vel)

eqn[1]=(p1-p2+pf)/l .+ ((mu[1])/rho)*dudx2[1]
eqn[2:N-1]=(p1-p2+pf)/l .+ ((mu[2:N-1])./rho).*dudx2[2:N-1]
eqn[N]=(p1-p2+pf)/l .+ ((mu[N])/rho)*dudx2[N]
eqn[N+1]=(beta/(Ap*(L-l-x)))*(-Qgap+Ap*u)
eqn[N+2]=(beta/(Ap*(L-l+x)))*(Qgap-Ap*u)

Force=(p1-p2+pf)*Ap
Fr=mu[1]*Aeff*dudx[1]

return (
    eqn=eqn,
    Qgap=Qgap,
    Force=Force,
    Fr=Fr
)

end

It would probably be best to use a for loop to aggregate the results, but I’ll show you how to do with map.

julia> results = map((t,y)->navier(t,y)[(:Qgap, :Force, :Fr)], 1:10, eachcol(rand(22,10)))
10-element Vector{NamedTuple{(:Qgap, :Force, :Fr), Tuple{Float64, Float64, Float64}}}:
 (Qgap = 0.00030102609426749455, Force = 53.32666898640483, Fr = 24.619494874431233)
 (Qgap = 0.00029967156523574053, Force = 53.32663257434115, Fr = -49.498430150318434)
 (Qgap = 0.0002801507212840063, Force = 53.32680302845422, Fr = 61.985357870666995)
 (Qgap = 0.00033119472944960455, Force = 53.327010098387575, Fr = 37.68745577297903)
 (Qgap = 0.00029727868343977574, Force = 53.326764122592515, Fr = -50.54111347534392)
 (Qgap = 0.0003496617278965436, Force = 53.3268011948002, Fr = 8.655337852075665)
 (Qgap = 0.0003150886996319404, Force = 53.32660908124239, Fr = -81.82836381327539)
 (Qgap = 0.00028876439037434895, Force = 53.32686035546042, Fr = 61.154294238471735)
 (Qgap = 0.0002664868528073511, Force = 53.326850589478646, Fr = -68.03440757714097)
 (Qgap = 0.00024699380480928317, Force = 53.326924487620666, Fr = 5.137304480117839)

julia> Qgap = getproperty.(results, :Qgap)
10-element Vector{Float64}:
 0.00030102609426749455
 0.00029967156523574053
 0.0002801507212840063
 0.00033119472944960455
 0.00029727868343977574
 0.0003496617278965436
 0.0003150886996319404
 0.00028876439037434895
 0.0002664868528073511
 0.00024699380480928317

julia> Force = getproperty.(results, :Force)
10-element Vector{Float64}:
 53.32666898640483
 53.32663257434115
 53.32680302845422
 53.327010098387575
 53.326764122592515
 53.3268011948002
 53.32660908124239
 53.32686035546042
 53.326850589478646
 53.326924487620666

julia> Fr = getproperty.(results, :Fr)
10-element Vector{Float64}:
  24.619494874431233
 -49.498430150318434
  61.985357870666995
  37.68745577297903
 -50.54111347534392
   8.655337852075665
 -81.82836381327539
  61.154294238471735
 -68.03440757714097
   5.137304480117839
3 Likes
5

Thank you very much for your generosity. I was looking forward to write it tomorrow morning, fresh start with a cup of coffee (Numerical Modelling is Numerically tiresome).

That’s all. No ODEProblem(navier,u0,tspan…), no choosing methods to solve?

6

I only parrotted back the code you gave us. If you want more, we’ll need more executable MATLAB code.

1 Like
7

Here is the full code. I have altered quite few parameteres (prototype of the lab).

f=10;
a=0.01;
time=linspace(0.25/f,2/f,100);

% 
xglob=a*sin(2*pi*f*time);
uglob=a*2*pi*f*cos(2*pi*f*time);

N=100;
y0(1:N,1)=0;
y0(N+1)=8e5;
y0(N+2)=8e5;
tic
[tsol,ysol]=ode45(@navier,time,y0);
toc

[~,Qgap,Force,Fr] = cellfun(@(tsol,ysol)  navier(tsol,ysol.'), num2cell(tsol), num2cell(ysol,2),'uni',0);
Qgap=cell2mat(Qgap);

figure(1)
hold on
yyaxis left
title('Qgap and Velocity')
xlabel('time')
ylabel('V(m/s)')
plot(tsol,ysol(:,1),'b--',tsol,ysol(:,N/2),'r--',tsol,ysol(:,N),'m--',tsol,uglob);            
yyaxis right
ylabel('m^3/s')
plot(tsol,Qgap,'ko-');

hold off



Force=cell2mat(Force);
Fr=cell2mat(Fr);

figure(2)
hold on
yyaxis left
title('Chamber Pressures')
xlabel('time[s]')
ylabel('Pressure[Pa]')
plot(tsol,ysol(:,N+1),'k',tsol,ysol(:,N+2));            
yyaxis right
ylabel('dist[m]')
legend('ComCha','RebCha')
set(gcf,'Position', [2 2 1000 1000])
plot(tsol,xglob);

hold off



function [eqn,Qgap,Force,Fr]=navier(t,y)

%parameters
rho=810;
beta=5e8;
r1=0.01;
r2=0.004;
Ap=pi*(0.01^2-0.004^2);
L=0.06;
a=0.01;
f=10;
l=0.01;
Aeff=pi*2*r1*l;
delta=0.007;
N=100;
dx=delta/N;
rad=linspace(r1,r1+delta,N);

vel=y(1:N);
p1=y(N+1);
p2=y(N+2);

x=a*sin(2*pi*f*t);
u=a*2*pi*f*cos(2*pi*f*t);

eqn=zeros(N+2,1);

dudx=[(u-vel(1))/dx;(vel(3:end) -vel(1:end-2))/(2*dx);(vel(end)-0)/dx];
dudx2=(vel(3:end)-2*vel(2:end-1)+vel(1:end-2))/(dx*dx);
dudx2 = [(dudx(1)-dudx(2))/dx;dudx2;(dudx(end-1)-dudx(end))/dx];
mu=[630 * (1+((0.05*(dudx)).^2)).^(-0.23)];
pf=12*a*0.707*u*cos(2*pi*f*t)/(delta^2)*630*(1+(0.05*(0.707*u/delta)^2))^(-0.23);


eqn(1)=(p1-p2+pf)/l + ((mu(1))/rho)*dudx2(1);
eqn(2:N-1)=(p1-p2+pf)/l + ((mu(2:N-1))./rho).*dudx2(2:N-1);
eqn(N)=(p1-p2+pf)/l + ((mu(N))/rho)*dudx2(N);

Qgap =simpsons(2*pi*transpose(rad).*vel,r1,r1+delta);
% Qgap=trapz(rad,2*pi*transpose(rad).*vel);

eqn(N+1)=(beta/(Ap*(L-x)))*(-Qgap+Ap*u);
eqn(N+2)=(beta/(Ap*(L+x)))*(Qgap-Ap*u);

Force=(p1-p2+pf)*Ap;
Fr=mu(1)*Aeff*dudx(1);

%eqn = transpose([temp(1,:) temp(2,:) temp(3,:)]);
    




end
8

I suggest posting the example as a new post asking specifically about ODE solvers. Perhaps you could even combine the video from Chris with the Julia function I ported for you to make an attempt and then ask a more specific question.

1 Like
9

That makes it more complicated don’t you think. I wish we worked on this, that it makes it easier for the other users new and experts out there.

10

If you want to learn how to solve ODEs with callbacks in DifferentialEquations.jl I suggest carefully going through the docs at Event Handling and Callback Functions · SciML and at Callback Library · SciML. These have examples showing how to solve differential equations using various kinds of callbacks.

Once you’ve read them, if you want to put together a simple illustrative Julia example and post it here as a new issue I’m sure you’ll be able to get help with any problems that you may encounter solving ODEs in Julia. Getting a simple example running with the features you need (i.e. callbacks and such) would allow you to understand the DifferentialEquations.jl workflow, and presumably make it easier for you to port your more complicated codes over.

2 Likes
11

I changed the title. Maybe @ChrisRackauckas or someone else is feeling generous enough to port the code further for you.

12

Per https://diffeq.sciml.ai/stable/solvers/ode_solve/

For users familiar with MATLAB/Python/R, good translations of the standard library methods are as follows:

  • ode45/dopri5 –> DP5(), though in most cases Tsit5() is more efficient
3 Likes
13

Just make it navier(dy,y,p,t) from @mkitti 's code and you’re done. Did you try that?

2 Likes
15

Hi Chris. Been tough with company meetings and holidays. Came back fresh and attacked the problem. But still I have got a bug.

I have followed your tutorial in setting up constant parameters and calculated ones inside the function from https://tutorials.sciml.ai/html/introduction/03-optimizing_diffeq_code.html

p = (1.0,1.0,1.0,10.0,0.001,100.0,Ayu,uAx,Du,Ayv,vAx,Dv) # a,α,ubar,β,D1,D2
function gm4!(dr,r,p,t)
  a,α,ubar,β,D1,D2,Ayu,uAx,Du,Ayv,vAx,Dv = p
  u = @view r[:,:,1]
  v = @view r[:,:,2]
  du = @view dr[:,:,1]
  dv = @view dr[:,:,2]
  mul!(Ayu,Ay,u)
  mul!(uAx,u,Ax)
  mul!(Ayv,Ay,v)
  mul!(vAx,v,Ax)
  @. Du = D1*(Ayu + uAx)
  @. Dv = D2*(Ayv + vAx)
  @. du = Du + a*u*u./v + ubar - α*u
  @. dv = Dv + a*u*u - β*v
end
prob = ODEProblem(gm4!,r0,(0.0,0.1),p)
@benchmark solve(prob,Tsit5())

So I have done the same with my code, but I get an error;

#parameters
p=(810,5e8,0.01,0.004,pi*(0.01^2-0.004^2),0.06,0.01,10,0.01,pi*2*0.01*0.01,0.007,100,0.0007/100,
    range(0.01,0.017,length=100),x,u,μ,Pf,Qg)

function navier!(dy,y,p,t)
    #parameters
    ρ,β,r₁,r₂,Aₚ,L₀,amp,f,Lₚ,Aₑ,δ,N,Δx,rad,x,u,μ,Pf,Qg=p
    vel=y[1:N]
    p1=y[N+1]
    p2=y[N+2]

    x=amp*sin(2*pi*f*t)
    u=amp*2*pi*f*cos(2*pi*f*t)

    dudx=[(u-vel(1))/Δx;(vel(3:end) -vel(1:end-2))/(2*Δx);(vel(end)-0)/Δx]
    dudx2=(vel[3:end]-2*vel[2:end-1].+vel[1:end-2])/(Δx*Δx)
    dudx2 = [(dudx[1]-dudx[2])/Δx;dudx2;(dudx[end-1]-dudx[end])/Δx]
    μ=630 * (1 .+ (0.05*(dudx)).^2).^(-0.23)
    Pf=12*amp*0.707*u*cos(2*pi*f*t)/(δ^2)*630*(1+(0.05*(0.707*u/δ)^2))^(-0.23)

    #Qgap =simpsons(2*pi*transpose(rad).*vel,r1,r1+delta)
    # Qgap=trapz(rad,2*pi*transpose(rad).*vel)
    Qgap = simps(2*pi*rad.*vel,r₁,r₁+δ)

    dy[1:end-2]=(p1-p2+pf)/Lₚ. + ((μ[1:N])./rho).*dudx2[1:N]
    dy[end-1]=(β/(Aₚ*(L₀-Lₚ/2-x)))*(-Qgap+Aₚ*u)
    dy[end]=(β/(Aₚ*(L₀-Lₚ/2+x)))*(Qgap-Aₚ*u)

    Force=(p1-p2+Pf)*Aₚ
    Fr=μ[1]*Aₑ*dudx[1]

    return (Qgap=Qgap,Force=Force,Fr=Fr)

end

UndefVarError: x not defined

Stacktrace:
[1] top-level scope
@ In[9]:2
[2] eval
@ .\boot.jl:373 [inlined]
[3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base .\loading.jl:1196

Can you or anyone help me please?

16

I haven’t read this thread and am not sure whether the code snippet you posted is supposed to be a self contained MWE, but if it is the error is obvious:

#parameters
p=(810,5e8,0.01,0.004,pi*(0.01^2-0.004^2),0.06,0.01,10,0.01,pi*2*0.01*0.01,0.007,100,0.0007/100,
    range(0.01,0.017,length=100),x,u,μ,Pf,Qg)

here you are referencing a variable x which is apparently undefined. In the notebook you linked, Ayu, uAx etc are all defined in previous cells, so you would similarly have to define x, u, μ, Pf, and Qg before using them.

17

Yes I thought of that too. Since x is dependent in time, it will be an iterative parameter. I have removed it from the list of p with u,Pf,Qg.

Now I get another error which I think comes from me not knowing to create non-allocation arrays inside the function. Here is the modified code.

#parameters
p=(810,5e8,0.01,0.004,pi*(0.01^2-0.004^2),0.06,0.01,10,0.01,pi*2*0.01*0.01,0.007,100,0.007/100,0.01:0.007/100:0.017)

function navier!(dy,y,p,t)
    #parameters
    ρ,β,r₁,r₂,Aₚ,L₀,amp,f,Lₚ,Aₑ,δ,N,Δx,rad=p
    vel=y[1:N]
    p1=y[N+1]
    p2=y[N+2]

    x=amp*sin(2*pi*f*t)
    u=amp*2*pi*f*cos(2*pi*f*t)

    dudx=[(u-vel(1))/Δx;(vel(3:end)-vel(1:end-2))/(2*Δx);(vel(end)-0)/Δx]
    dudx2=(vel[3:end]-2*vel[2:end-1].+vel[1:end-2])/(Δx*Δx)
    dudx2 = [(dudx[1]-dudx[2])/Δx;dudx2;(dudx[end-1]-dudx[end])/Δx]
    μ=630 * (1 .+ (0.05*(dudx)).^2).^(-0.23)
    Pf=12*amp*0.707*u*cos(2*pi*f*t)/(δ^2)*630*(1+(0.05*(0.707*u/δ)^2))^(-0.23)

    #Qgap =simpsons(2*pi*transpose(rad).*vel,r1,r1+delta)
    # Qgap=trapz(rad,2*pi*transpose(rad).*vel)
    Qgap = simps(2*pi*rad.*vel,r₁,r₁+δ)

    dy[1:end-2]=(p1-p2+pf)/Lₚ. + ((μ[1:N])./rho).*dudx2[1:N]
    dy[end-1]=(β/(Aₚ*(L₀-Lₚ/2-x)))*(-Qgap+Aₚ*u)
    dy[end]=(β/(Aₚ*(L₀-Lₚ/2+x)))*(Qgap-Aₚ*u)

    Force=(p1-p2+Pf)*Aₚ
    Fr=μ[1]*Aₑ*dudx[1]

    return (Qgap=Qgap,Force=Force,Fr=Fr)

end

I am getting this error;

syntax: missing last argument in “3:” range expression

Stacktrace:
[1] top-level scope
@ In[3]:14
[2] eval
@ .\boot.jl:373 [inlined]
[3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base .\loading.jl:1196

18

Julia uses square brackets for array indexes [ and ] go through and change all instances of using ( and )

19

I feel dumb as a primary school kid not seeing this :grinning:.

I think now the generic function has been created with this code.
Two things;

  1. I want later on to get to the origin of this post (calling the variables Qgap,Force and Fr after the solution).
  2. Any suggestion to optimise the code, as I want to increase the FDM spacing to N=1000 for stability reasons and am afraid it will take ages.
20

the main thing is to stop rolling your own integration. Solving Integro Differential Equations · NeuralPDE.jl describes how to deal with integrodifferential equations more efficiently.

3 Likes
21

Please be generous enough to help me through it. Here is the finest version.

#parameters
N=101
p=(810,5e8,0.01,0.004,pi*(0.01^2-0.004^2),0.06,0.01,10,0.01,pi*2*0.01*0.01,0.007,N,0.007/N,range(0.01,0.017,length=N))

function navier!(dy,y,p,t)
    #parameters
    ρ,β,r₁,r₂,Aₚ,L₀,amp,f,Lₚ,Aₑ,δ,N,Δx,rad=p
    vel=y[1:N]
    p1=y[N+1]
    p2=y[N+2]

    x=amp*sin(2*pi*f*t)
    u=amp*2*pi*f*cos(2*pi*f*t)

    dudx=[(u-vel[1])/Δx;(vel[3:end]-vel[1:end-2])/(2*Δx);(vel[end]-0)/Δx]
    dudx2=(vel[3:end]-2*vel[2:end-1].+vel[1:end-2])/(Δx*Δx)
    dudx2 = [(dudx[1]-dudx[2])/Δx;dudx2;(dudx[end-1]-dudx[end])/Δx]
    μ=630 * (1 .+ (0.05*(dudx)).^2).^(-0.23)
    Pf=12*amp*0.707*u*cos(2*pi*f*t)/(δ^2)*630*(1+(0.05*(0.707*u/δ)^2))^(-0.23)

    #Qgap =simpsons(2*pi*transpose(rad).*vel,r1,r1+delta)
    # Qgap=trapz(rad,2*pi*transpose(rad).*vel)
    Qgap = simps(2*pi*rad.*vel,r₁,r₁+δ)

    dy[1:N]=(p1-p2+Pf)/Lₚ.+((μ[1:N])./ρ).*dudx2[1:N]
    dy[N+1]=(β/(Aₚ*(L₀-Lₚ/2-x)))*(-Qgap+Aₚ*u)
    dy[N+2]=(β/(Aₚ*(L₀-Lₚ/2+x)))*(Qgap-Aₚ*u)

    Force=(p1-p2+Pf)*Aₚ
    Fr=μ[1]*Aₑ*dudx[1]

    return (Qgap=Qgap,Force=Force,Fr=Fr)

end

with the solving method chosen as follows;

y0=zeros(Float64,N+2)
y0[1:N].=0
y0[N+1]=8e5
y0[N+2]=8e5
time=(0.0,0.5)
ode_prob=ODEProblem(navier!,y0,time,p);
@btime soln=solve(ode_prob,Rodas5(),abstol=1e-6,reltol=1e-6,saveat=0.3:0.001:0.4);