function [w_tilde, w_tilde_evo, t, t_evo] = RK4_FGM2D(f,x,y,kkx,kx,kky,ky,kk_2,k_2,t,ts,te,dt,Re,N,w0_tilde,w0,autosaveInterval,monitorOrNot)
w_tilde_evo=[];
t_evo=[];
if monitorOrNot
[~,H]=contourf(x(1:4:end,1:4:end),y(1:4:end,1:4:end),w0(1:4:end,1:4:end),80,'edgecolor','none');
axis equal
colorbar;
title("Time= "+ts+"s")
xlabel("Dimensionless $$x$$","Interpreter","latex")
ylabel("Dimensionless $$y$$","Interpreter","latex")
drawnow
end
range=N/4+1:N/4+N;
M=N*3/2;
displaySampling=1:6:M;
if isempty(dt)
t(1)=ts;
w_tilde(:,:)=w0_tilde;
w_tilde(:,1)=0;
w_tilde(1,:)=0;
counter_frame=0;
counter_autosave=0;
counter_evo=1;
frameInterval=floor(N/32);
if autosaveInterval
% Variable dt with autosave
n=1;
while t(n)<te
[w_tilde , ww, h]=RK4_sub1(f,w_tilde,kkx,kx,kky,ky,kk_2,k_2,Re,N,M,range,t(n),te);
t(n+1)=t(n)+h;
if counter_autosave==autosaveInterval
[w_tilde_evo(:,:,counter_evo), t_evo(counter_evo)]=saveEvo(w_tilde,t(n+1)); %#ok<AGROW>
counter_autosave=0;
counter_evo=counter_evo+1;
end
if counter_frame==frameInterval
if monitorOrNot
monitor(H,real(ww(displaySampling,displaySampling)),t(n+1))
counter_frame=0;
end
end
counter_frame=counter_frame+1;
counter_autosave=counter_autosave+1;
n=n+1;
end
if t_evo(counter_evo-1)~=te
[w_tilde_evo(:,:,counter_evo), t_evo(counter_evo)]=saveEvo(w_tilde,te);
end
else
% Variable dt without autosave
n=1;
while t(n)<te
[w_tilde , ww, h]=RK4_sub1(f,w_tilde,kkx,kx,kky,ky,kk_2,k_2,Re,N,M,range,t(n),te);
t(n+1)=t(n)+h;
if counter_frame==frameInterval
if monitorOrNot
monitor(H,real(ww(displaySampling,displaySampling)),t(n+1))
counter_frame=0;
end
end
counter_frame=counter_frame+1;
n=n+1;
end
end
else
if autosaveInterval
w_tilde_evo=zeros(N,N,ceil(length(t)/autosaveInterval)-1);
t_evo=zeros(1,ceil(length(t)/autosaveInterval)-1);
end
w_tilde(:,:,1)=w0_tilde;
w_tilde(:,1,1)=0;
w_tilde(1,:,1)=0;
counter_frame=0;
counter_autosave=0;
counter_evo=1;
frameInterval=floor(length(t)/50);
if autosaveInterval
% Fixed dt with autosave
for n=1:length(t)-1
[w_tilde , ww]=RK4_sub2(f,w_tilde,dt,kkx,kx,kky,ky,kk_2,k_2,Re,M,range);
if counter_autosave==autosaveInterval
[w_tilde_evo(:,:,counter_evo), t_evo(counter_evo)]=saveEvo(w_tilde,t(n+1)); %#ok<AGROW>
counter_autosave=0;
counter_evo=counter_evo+1;
end
if counter_frame==frameInterval
if monitorOrNot
monitor(H,real(ww(displaySampling,displaySampling)),t(n+1))
counter_frame=0;
end
end
counter_frame=counter_frame+1;
counter_autosave=counter_autosave+1;
end
if t_evo(counter_evo-1)~=te
[w_tilde_evo(:,:,counter_evo), t_evo(counter_evo)]=saveEvo(w_tilde,te);
end
else
% Fixed dt without autsave
for n=1:length(t)-1
[w_tilde , ww]=RK4_sub2(f,w_tilde,dt,kkx,kx,kky,ky,kk_2,k_2,Re,M,range);
if counter_frame==frameInterval
if monitorOrNot
monitor(H,real(ww(displaySampling,displaySampling)),t(n+1))
counter_frame=0;
end
end
counter_frame=counter_frame+1;
end
end
end
return
function [w_tilde , ww, h] = RK4_sub1(f,w_tilde,kkx,kx,kky,ky,kk_2,k_2,Re,N,M,range,t,te)
[f_tpm, ww, uu, vv]=f(w_tilde,kkx,kx,kky,ky,kk_2,k_2,Re,M,range);
% Determine dt;
a=(N-2)/2*(max(abs(real(uu(:))))+max(abs(real(vv(:)))));
b=N^2/2/Re;
h=min([2.82/a, 2.8/b]);
s=-a*1i*h-b*h;
ss=a*1i*h-b*h;
while abs(abs(1+s+s^2/2+s^3/6+s^4/24)-1)>10^-3
F=(abs(1+s+s^2/2+s^3/6+s^4/24)-1);
dF=((s^4/(6*h)+s^3/(2*h)+s^2/h+s/h)*(ss^4/24+ss^3/6+ss^2/2+ss+1)+(ss^4/(6*h)+ss^3/(2*h)+ss^2/h+ss/h)*(s^4/24+s^3/6+s^2/2+s+1))/(2*((s^4/24+s^3/6+s^2/2+s+1)*(ss^4/24+ss^3/6+ss^2/2+ss+1))^(1/2));
h=h-F/dF;
s=-a*1i*h-b*h;
ss=a*1i*h-b*h;
end
if te-t<h
h=te-t;
end
F1=h.*f_tpm;
F2=h.*f(w_tilde+F1/2,kkx,kx,kky,ky,kk_2,k_2,Re,M,range);
F3=h.*f(w_tilde+F2/2,kkx,kx,kky,ky,kk_2,k_2,Re,M,range);
F4=h.*f(w_tilde+F3,kkx,kx,kky,ky,kk_2,k_2,Re,M,range);
w_tilde=w_tilde+(F1+2*F2+2*F3+F4)./6;
w_tilde(:,1)=0;
w_tilde(1,:)=0;
return
function [w_tilde , ww] = RK4_sub2(f,w_tilde,h,kkx,kx,kky,ky,kk_2,k_2,Re,M,range)
[f_tpm, ww, ~, ~]=f(w_tilde,kkx,kx,kky,ky,kk_2,k_2,Re,M,range);
F1=h.*f_tpm;
F2=h.*f(w_tilde+F1/2,kkx,kx,kky,ky,kk_2,k_2,Re,M,range);
F3=h.*f(w_tilde+F2/2,kkx,kx,kky,ky,kk_2,k_2,Re,M,range);
F4=h.*f(w_tilde+F3,kkx,kx,kky,ky,kk_2,k_2,Re,M,range);
w_tilde=w_tilde+(F1+2*F2+2*F3+F4)./6;
w_tilde(:,1)=0;
w_tilde(1,:)=0;
return
function monitor(H,ww,t)
H.ZData=ww;
H.LevelList=linspace(min(ww(:)),max(ww(:)),80);
title("Time= "+t+"s")
pause(0)
return
function [w_tilde, t] = saveEvo(w_tilde,t)