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

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?

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?

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.

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

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?

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

1 Like

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

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

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.

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.

1 Like

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

Per ODE Solvers · DifferentialEquations.jl

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

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

2 Likes

I am on it for the weekend…