There are 3 very simple improvements here. The first is @views function navier!(dy,y,p,t) which will make array slices lazy rather than copies. The second is Qgap = 2*pi*rad*simps(vel,r₁,r₁+δ) which will save a bunch of multiplies. The third is using a better solver. as you increase N, this problem will become increasingly stiff, so you should probably be using a solver like FBDF.
1 Like
Thanks for the tips.
Q_{gap}=\int2\pi r*u(r)dr mathematically incorrect with what you’ve suggested.
oh right. you can’t move the r out of the integral.
although the code you wrote I think is equally wrong.
I have tried to use CVODE_BDF from Sundials. Still loading and time consuming.
using Sundials
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,CVODE_BDF(linear_solver=:GMRES),abstol=1e-6,reltol=1e-6,saveat=0.3:0.001:0.4);
Got this scary memory allocation report 
17.961 s (32365094 allocations: 13.89 GiB)
And I get this as solution report
retcode: MaxIters
Interpolation: 1st order linear
t: Float64
u: Vector{Float64}
yeah. your code allocates way more than it should. if you rewrote this to loop over dy, it would be way faster.
What exactly do you mean by that?
Thank you Chris. I have read the article multiple times. That is the least I could get out of it.
Please help me if possible.
I wished it wasn’t physically feasible, but I did manage to solve it with matlab and got a valid solution.
Are the MATLAB and Julia codes exactly the same, i.e. if you put in a random u,p,t do you get the exact same derivative?
Here is the matlab code. I’m quite sure have done every conversion procedure correctly.
f=10;
a=0.01;
time=linspace(0,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','LineWidth',3);
plot(tsol,ysol(:,N+2),'b-','LineWidth',3);
legend('ComCha','RebCha')
yyaxis right
ylabel('dist[m]')
set(gcf,'Position', [2 2 1000 1000])
plot(tsol,xglob,'HandleVisibility','off');
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);
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-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
What is your evidence? On a randomly generated u,p,t, do you have a script that compares the output of the derivative calculation and shows that they are exactly the same? That’s the standard starting point to show such a claim.
1 Like
Yes.
I have tried two set of parameters. y=y0 and t=1e-7
I have got the same dudx and dudx2 for both. (first values of 8975.979010079373 and 1.28…e8 respectively).
y=ysol(9) from matlab at t=0.0162
(0.32668, 0.33158, 0.33626, 0.34074, 0.34501, 0.34909, 0.35298, 0.35667, 0.36018, 0.36351, 0.36666, 0.36963, 0.37243, 0.37507, 0.37754, 0.37985, 0.382, 0.38401, 0.38586, 0.38757, 0.38913, 0.39056, 0.39185, 0.39301, 0.39404, 0.39492, 0.39572, 0.39633, 0.3969, 0.39723, 0.39758, 0.39763, 0.39777, 0.39757, 0.39745, 0.39704, 0.39664, 0.39601, 0.39532, 0.39446, 0.39351, 0.3924, 0.39118, 0.38982, 0.38832, 0.38667, 0.38489, 0.38296, 0.38087, 0.37863, 0.37624, 0.37368, 0.37095, 0.36806, 0.365, 0.36175, 0.35833, 0.35472, 0.35093, 0.34694, 0.34275, 0.33837, 0.33378, 0.32899, 0.32398, 0.31876, 0.31332, 0.30765, 0.30176, 0.29564, 0.28928, 0.28268, 0.27585, 0.26876, 0.26143, 0.25384, 0.24599, 0.23788, 0.22951, 0.22087, 0.21195, 0.20276, 0.19329, 0.18353, 0.17349, 0.16315, 0.15252, 0.14159, 0.13035, 0.11881, 0.10696, 0.09479, 0.082305, 0.069498, 0.056366, 0.042906, 0.029112, 0.014982, 0.00051189, -0.014302, 768650.0, 843650.0)
I get this error in julia for dudx;
MethodError: no method matching *(::Int64, ::NTuple{98, Float64})
Closest candidates are:
*(::Any, ::Any, ::Any, ::Any…) at C:\Users\Simon Luzara\AppData\Local\Programs\Julia-1.7.3\share\julia\base\operators.jl:655
*(::T, ::T) where T<:Union{Int128, Int16, Int32, Int64, Int8, UInt128, UInt16, UInt32, UInt64, UInt8} at C:\Users\Simon Luzara\AppData\Local\Programs\Julia-1.7.3\share\julia\base\int.jl:88
*(::Union{Real, Complex}, ::Union{LinearAlgebra.Adjoint{var"#s861", var"#s8611"}, LinearAlgebra.Transpose{var"#s861", var"#s8611"}} where {var"#s861"<:Union{Real, Complex}, var"#s8611"<:(AbstractVector)}, ::AbstractMatrix{<:Union{Real, Complex}}, ::AbstractMatrix{<:Union{Real, Complex}}) at C:\Users\Simon Luzara\AppData\Local\Programs\Julia-1.7.3\share\julia\stdlib\v1.7\LinearAlgebra\src\matmul.jl:1172
…Stacktrace:
[1] top-level scope
@ In[29]:1
[2] eval
@ .\boot.jl:373 [inlined]
[3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base .\loading.jl:1196
Here is a quick snip of the solution from matlab. vel(t) along the gap.
A tuple of 98 things will kill compile and run time performance. I would avoid that.
Can you share the test script? It should be a simple copy-paste that just runs the two (MATLAB.jl helps here).
Here is the julia code I created.
using DifferentialEquations
using ModelingToolkit
using PyPlot
using BenchmarkTools
using NumericalIntegration
N=100
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,LinRange(0.01,0.017,N))
@views 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₁+δ)
Qgap = integrate(rad, 2*pi*rad.*vel)
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)
nothing
end
using Sundials
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);
soln=solve(ode_prob,CVODE_BDF(linear_solver=:GMRES),abstol=1e-3,reltol=1e-5,saveat=0.3:0.005:0.4)
Can I share the test script? What I’ve done to test the value is run the matlab function above itself by saying y=ysol(9) (the array I posted with t=0.0162).
Keep things simple. I don’t have MATLAB installed, so could you run:
using NumericalIntegration
N=100
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,LinRange(0.01,0.017,N))
y=rand(N+2)
t=5rand()
@views 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₁+δ)
Qgap = integrate(rad, 2*pi*rad.*vel)
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)
nothing
end
dy = similar(y0)
navier!(dy,y,p,t)
dy
using MATLAB
mat"""
f=10;
a=0.01;
time=linspace(0,2/f,100);
%
xglob=a*sin(2*pi*f*time);
uglob=a*2*pi*f*cos(2*pi*f*time);
N=100;
y = $y
$mdy = navier($t,y.')
"""
@show mdy - dy
and confirm that shows all zeros? I think I got something of the call to navier incorrect: it’s not enclosing all of the parameters you have here.
2 Likes
Nothing good is coming out of it. I have created the function code on the same directory as julia.
I have runned the code and am getting this error.
Error using vertcat
Dimensions of arrays being concatenated are not consistent.Error in navier (line 28)
dudx=[(u-vel(1))/dx;(vel(3:end) -vel(1:end-2))/(2*dx);(vel(end)-0)/dx]';MATLAB.MEngineError(“failed to get variable matlab_jl_4 from MATLAB session”)
Stacktrace:
[1] get_mvariable(session::MSession, name::Symbol)
@ MATLAB C:\Users\Simon Luzara.julia\packages\MATLAB\SVjnA\src\engine.jl:164
[2] get_mvariable
@ C:\Users\Simon Luzara.julia\packages\MATLAB\SVjnA\src\engine.jl:168 [inlined]
[3] get_variable(name::Symbol)
@ MATLAB C:\Users\Simon Luzara.julia\packages\MATLAB\SVjnA\src\engine.jl:170
[4] top-level scope
@ C:\Users\Simon Luzara.julia\packages\MATLAB\SVjnA\src\matstr.jl:162
[5] eval
@ .\boot.jl:373 [inlined]
[6] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
@ Base .\loading.jl:1196
I think I would stick to matlab for unsteady flow, and use Julia for steady flow analysis with no du/dt. I came up with an analytical formulae.
Big thanks to you Chris and everyone. Now I think am quite ready to be a julia user 
!