%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This script computes several bi-fluidic Green-Naghdi models  %
% for the propagation of internal waves                        %
% The derivation and presentation of the models is detailed in %
% "A new class of two-layer Green-Naghdi systems with improved %
% frequency dispersion" by V. Duchêne, S. Israwi, R. Taalhouk  %
%                                                              %
% The author of this script is V. Duchêne                      %
% please let me know if you use it, found a bug or improvement %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% This script should run with odetpbar.m and textprogressbar.m
% in the same folder (or delete "'OutputFcn',@odetpbar" occurences)
% odetbar and textprogressbar are available at
% http://www.mathworks.com/matlabcentral/fileexchange/28509
% Copyright (c) 2010, David Romero-Antequera
% Copyright (c) 2010, Paul Proteus 
% All rights reserved.



clear all
close all

% % % PREPARATION % % %

%Parameters of the model
mu=.1; % shallowness parameter (mu=0 for Saint-Venant system)
eps=.5; % nonlinearity parameter
d=.5; % depth ratio, delta=d_1/d_2
g=.95; % density ratio, gamma=rho_1/rho_2 (g=0,d=1 yields one-layer case)
iBo=mu/100; % inverse of Bond number (iBo=0 to neglect surface tension)

%Parameters of the numerical scheme
L=4; % size of the spatial domain
N=2^9; % number of discretization points
T=2; % final time of integration
tol=1e-12; % tolerance in the ode45 method
dx=2*L/N;x=(-L:dx:L-dx)'; % spatial mesh
time=0:T/100:T; % recorded time values

%Initial data
% the first half component is the interface deformation, z=zeta
% the second half component is the flux, w
U0=[-exp(-4*(x+1/3).^2);exp(-4*(x-1/3).^2)];

%Discrete differentiation operators and Fourier multipliers
% We use spectral methods to construct spatial Fourier multipliers
% (including the discrete differentiation matrix, D1)
% Essentially, this amounts to (fast) Fourier transform, multiply in 
% Fourier variables, then inverse Fourier transform
% F,Fm1,Fm2,Fi1,Fi2 are the functions defining the Fourier multipliers
% D1,DFm1,DFm2,DFi1,DFi2 are the discrete operators (N by N matrices)
F=@(D) 1i*D;
mF=F((pi/L)*[(0:N/2-1) 0 -N/2+1:-1]');
col=real(ifft(mF.*fft([1;zeros(N-1,1)])));
D1=toeplitz(col,-col); 

Fm1=@(D) 1i*D./sqrt(1+mu/15*D.*D);
mF=Fm1((pi/L)*[(0:N/2-1) 0 -N/2+1:-1]');
col=real(ifft(mF.*fft([1;zeros(N-1,1)])));
DFm1=toeplitz(col,-col);
    
Fm2=@(D) 1i*D./sqrt(1+mu/15/d/d*D.*D);
mF=Fm2((pi/L)*[(0:N/2-1) 0 -N/2+1:-1]');
col=real(ifft(mF.*fft([1;zeros(N-1,1)])));
DFm2=toeplitz(col,-col);
    
Fi1=@(D) 1i*D.*sqrt(3./(sqrt(mu)*D.*tanh(sqrt(mu)*D))-3./(mu*D.*D));
mF=Fi1((pi/L)*[(0:N/2-1) 0 -N/2+1:-1]');
mF(end/2+1)=0;mF(1)=0;
col=real(ifft(mF.*fft([1;zeros(N-1,1)])));
DFi1=toeplitz(col,-col);
    
Fi2=@(D) 1i*D.*sqrt(3./(sqrt(mu)*D.*tanh(sqrt(mu)*D/d)/d)-3./(mu*D.*D/d/d));
mF=Fi2((pi/L)*[(0:N/2-1) 0 -N/2+1:-1]');
mF(end/2+1)=0;mF(1)=0;
col=real(ifft(mF.*fft([1;zeros(N-1,1)])));
DFi2=toeplitz(col,-col);


% Construction of the functions used for the time-evolution scheme (ODEs)
D=@(u) spdiags(u,0,length(u),length(u)); % produces a diagonal matrix so 
% that D(u)*v=u.*v for any vectors u,v

h1=@(z) 1-eps*z;h2=@(z) 1/d+eps*z; % depth of the upper and lower layers
G=@(z,w) eps/2*(h1(z).^2-g*h2(z).^2)./((h1(z).*h2(z)).^2).*(w.^2);
H=@(z) (h1(z)+g*h2(z))./(h1(z).*h2(z));
DH=@(z) D(H(z));

% operator Q[epsilon zeta] in the original model (DQ)
DQ1=@(z) D(1./h1(z))*D1*D(h1(z).^3)*D1*D(1./h1(z));
DQ2=@(z) D(1./h2(z))*D1*D(h2(z).^3)*D1*D(1./h2(z));
DQ=@(z) -1/3*DQ2(z)-g/3*DQ1(z);

% operator Q[epsilon zeta] in the regularized model (DQm)
DQm1=@(z) D(1./h1(z))*DFm1*D(h1(z).^3)*DFm1*D(1./h1(z));
DQm2=@(z) D(1./h2(z))*DFm2*D(h2(z).^3)*DFm2*D(1./h2(z));
DQm=@(z) -1/3*DQm2(z)-g/3*DQm1(z);

% operator Q[epsilon zeta] in the improved model (DQi)
DQi1=@(z) D(1./h1(z))*DFi1*D(h1(z).^3)*DFi1*D(1./h1(z));
DQi2=@(z) D(1./h2(z))*DFi2*D(h2(z).^3)*DFi2*D(1./h2(z));
DQi=@(z) -1/3*DQi2(z)-g/3*DQi1(z);

% operator R[epsilon zeta,epsilon w] in the original model (R)
R1=@(z,w) h1(z).*(D1*(w./h1(z)));
R2=@(z,w) h2(z).*(D1*(w./h2(z)));
Q1=@(z,w) 1./h1(z).*(D1*((h1(z).^3).*(D1*(w./h1(z)))));
Q2=@(z,w) 1./h2(z).*(D1*((h2(z).^3).*(D1*(w./h2(z)))));
R=@(z,w) 1/3*w./h2(z).*Q2(z,w)-g/3*w./h1(z).*Q1(z,w)+1/2*R2(z,w).^2-g/2*R1(z,w).^2;

% operator R[epsilon zeta,epsilon w] in the regularized model (Rm)
Rm1=@(z,w) h1(z).*(DFm1*(w./h1(z)));
Rm2=@(z,w) h2(z).*(DFm2*(w./h2(z)));
Qm1=@(z,w) 1./h1(z).*(DFm1*((h1(z).^3).*(DFm1*(w./h1(z)))));
Qm2=@(z,w) 1./h2(z).*(DFm2*((h2(z).^3).*(DFm2*(w./h2(z)))));
Rm=@(z,w) 1/3*w./h2(z).*Qm2(z,w)-g/3*w./h1(z).*Qm1(z,w)+1/2*Rm2(z,w).^2-g/2*Rm1(z,w).^2;

% operator R[epsilon zeta,epsilon w] in the improved model (Ri)
Ri1=@(z,w) h1(z).*(DFi1*(w./h1(z)));
Ri2=@(z,w) h2(z).*(DFi2*(w./h2(z)));
Qi1=@(z,w) 1./h1(z).*(DFi1*((h1(z).^3).*(DFi1*(w./h1(z)))));
Qi2=@(z,w) 1./h2(z).*(DFi2*((h2(z).^3).*(DFi2*(w./h2(z)))));
Ri=@(z,w) 1/3*w./h2(z).*Qi2(z,w)-g/3*w./h1(z).*Qi1(z,w)+1/2*Ri2(z,w).^2-g/2*Ri1(z,w).^2;

% surface tension operator
S=@(z) D1*(1./(sqrt(1+mu*eps^2*(D1*z).^2)).*(D1*z));


% The following functions solve (H(epsilon zeta)+mu*Q[epsilon zeta]) w = v
% and conversely
w2vGN=@(z,w) (DH(z)+mu*DQ(z))*w; %(z,w)->v (original Green-Naghdi model)
v2wGN=@(z,v) (DH(z)+mu*DQ(z))\v; %(z,v)->w (original Green-Naghdi model)

w2vGNm=@(z,w) (DH(z)+mu*DQm(z))*w; %(z,w)->v (regularized Green-Naghdi model)
v2wGNm=@(z,v) (DH(z)+mu*DQm(z))\v; %(z,v)->w (regularized Green-Naghdi model)

w2vGNi=@(z,w) (DH(z)+mu*DQi(z))*w; %(z,w)->v (improved Green-Naghdi model)
v2wGNi=@(z,v) (DH(z)+mu*DQi(z))\v; %(z,v)->w (improved Green-Naghdi model)


% At last, the functions used for the time-evolution scheme
% The models are solved using unknowns z=zeta and v, so that the system is
%    dt (z,v) + g(z,w) = 0
% thus w needs to be reconstructed at each time step
gGN=@(z,w) [D1*w;D1*((g+d)*z+G(z,w)-mu*eps*R(z,w)-iBo*(g+d)*S(z))];
gGNm=@(z,w) [D1*w;D1*((g+d)*z+G(z,w)-mu*eps*Rm(z,w)-iBo*(g+d)*S(z))];
gGNi=@(z,w) [D1*w;D1*((g+d)*z+G(z,w)-mu*eps*Ri(z,w)-iBo*(g+d)*S(z))];
fGN=@(t,V) -gGN(V(1:end/2),v2wGN(V(1:end/2),V(end/2+1:end)));
fGNm=@(t,V) -gGNm(V(1:end/2),v2wGNm(V(1:end/2),V(end/2+1:end)));
fGNi=@(t,V) -gGNi(V(1:end/2),v2wGNi(V(1:end/2),V(end/2+1:end)));




% % % COMPUTATION % % %

% you should have odetpbar.m and textprogressbar.m in your folder
% or delete "'OutputFcn',@odetpbar" occurences
if exist('odetpbar','file')==0
    message=sprintf(...
    ['This script should run with odetpbar.m and textprogressbar.m\n'...
    'available at http://www.mathworks.com/matlabcentral/fileexchange/28509\n'...
    'Copyright (c) 2010, David Romero-Antequera\n'...
    'Copyright (c) 2010, Paul Proteus ']);
    error(message);
end

disp('*** Original Green-Naghdi model ***')
V0=[U0(1:end/2);w2vGN(U0(1:end/2),U0(end/2+1:end))]; %initial data (z,v)
options = odeset('RelTol',tol*100,'AbsTol',tol,'OutputFcn',@odetpbar);
[time,UGN]=ode45(fGN,time,V0,options); %Runge-Kutta scheme

disp('*** Regularized Green-Naghdi model ***')
V0=[U0(1:end/2);w2vGNm(U0(1:end/2),U0(end/2+1:end))]; %initial data (z,v)
options = odeset('RelTol',tol*100,'AbsTol',tol,'OutputFcn',@odetpbar);
[time,UGNm]=ode45(fGNm,time,V0,options); %Runge-Kutta scheme

disp('*** Improved Green-Naghdi model ***')
V0=[U0(1:end/2);w2vGNi(U0(1:end/2),U0(end/2+1:end))]; %initial data (z,v)
options = odeset('RelTol',tol*100,'AbsTol',tol,'OutputFcn',@odetpbar);
[time,UGNi]=ode45(fGNi,time,V0,options); %Runge-Kutta scheme



% % % VISUALISATION % % %

% Energy (Hamiltonian) of each system
T1=@(z,w) h1(z).*(h1(z).*(D1*(w./h1(z)))).^2;
T2=@(z,w) h2(z).*(h2(z).*(D1*(w./h2(z)))).^2;
T1i=@(z,w) h1(z).*(h1(z).*(DFi1*(w./h1(z)))).^2;
T2i=@(z,w) h2(z).*(h2(z).*(DFi2*(w./h2(z)))).^2;
T1m=@(z,w) h1(z).*(h1(z).*(DFm1*(w./h1(z)))).^2;
T2m=@(z,w) h2(z).*(h2(z).*(DFm2*(w./h2(z)))).^2;
iS=@(z) sqrt(1+mu*eps^2*((D1*z).^2))-1;
eGN=@(z,w) (g+d)*z.^2+2*(g+d)/mu/eps/eps*iBo*iS(z)+H(z).*w.^2+mu*g/3*T1(z,w)+mu/3*T2(z,w);
eGNm=@(z,w) (g+d)*z.^2+2*(g+d)/mu/eps/eps*iBo*iS(z)+H(z).*w.^2+mu*g/3*T1m(z,w)+mu/3*T2m(z,w);
eGNi=@(z,w) (g+d)*z.^2+2*(g+d)/mu/eps/eps*iBo*iS(z)+H(z).*w.^2+mu*g/3*T1i(z,w)+mu/3*T2i(z,w);


figure('units','normalized','outerposition',[0 0 1 1])
Col=[0 0 1;0 .5 0;1 0 0];
set(gcf,'DefaultAxesColorOrder',Col);

I=zeros(length(time),3);V=I;E=I;Z=I;
for i=1:length(time)
% Construction of the data

% deformation of the interface, zeta
zGN=UGN(i,1:end/2)';
zGNm=UGNm(i,1:end/2)';
zGNi=UGNi(i,1:end/2)';
% horizontal velocity, v
vGN=UGN(i,1+end/2:end)';
vGNm=UGNm(i,1+end/2:end)';
vGNi=UGNi(i,1+end/2:end)';
% flux, w
wGN=v2wGN(zGN,vGN);
wGNi=v2wGNi(zGNi,vGNi);
wGNm=v2wGNm(zGNm,vGNm);

% amplitude of the Fourier transform
k=(0:(length(x)-1)/2)/2/L*2*pi;
hatzGN=(abs(fft(zGN))*dx/L)';
hatzGNi=(abs(fft(zGNi))*dx/L)';
hatzGNm=(abs(fft(zGNm))*dx/L)';
hatwGN=(abs(fft(wGN))*dx/L)';
hatwGNi=(abs(fft(wGNi))*dx/L)';
hatwGNm=(abs(fft(wGNm))*dx/L)';
hatvGN=(abs(fft(vGN))*dx/L)';
hatvGNi=(abs(fft(vGNi))*dx/L)';
hatvGNm=(abs(fft(vGNm))*dx/L)';

% conserved quantities
Z(i,:)=[sum(zGN);sum(zGNm);sum(zGNi)]*dx;
V(i,:)=[sum(vGN);sum(vGNm);sum(vGNi)]*dx;
I(i,:)=[sum(zGN.*vGN);sum(zGNm.*vGNm);sum(zGNi.*vGNi)]*dx;
E(i,:)=[sum(eGN(zGN,wGN));sum(eGNm(zGNm,wGNm));sum(eGNi(zGNi,wGNi))]*dx;



subplot(2,2,1)  % upper-left panel : surface deformation
plot(x,zGN,x,zGNi,x,zGNm)
axis([-L L -2 2])
xlabel('$x$','interpreter','latex')
ylabel('$\zeta(x)$','interpreter','latex')
legend('original GN','improved GN','regularized GN','Location','SouthEast')
title(['surface deformation, zeta, at t= ' num2str(time(i))])

subplot(2,2,2) % upper-right panel : flux, w
plot(x,wGN,x,wGNi,x,wGNm)
axis([-L L -2 2])
xlabel('$x$','interpreter','latex')
ylabel('$w(x)$','interpreter','latex')
legend('original GN','improved GN','regularized GN','Location','SouthEast')
title(['flux, w, at time t= ' num2str(time(i))])

subplot(2,2,3) % lower-left panel : Fourier transform of surface deformation (log10 scale)
plot(k,log10(hatzGN(1:end/2)),k,log10(hatzGNi(1:end/2)),k,log10(hatzGNm(1:end/2)))
axis([0 k(end) -18 0])
xlabel('$k$','interpreter','latex')
ylabel('$\log(|\widehat{\zeta}(k)|)$','interpreter','latex')
title('amplitude of the Fourier transform of zeta (in log10 scale)')
set(gca,'YGrid','on');

subplot(2,2,4) % lower-right panel : Fourier transform of flux (log10 scale)
plot(k,log10(hatwGN(1:end/2)),k,log10(hatwGNi(1:end/2)),k,log10(hatwGNm(1:end/2)))
axis([0 k(end) -18 0])
xlabel('$k$','interpreter','latex')
ylabel('$\log(|\widehat{w}(k)|)$','interpreter','latex')
title('amplitude of the Fourier transform of w (in log10 scale)')
set(gca,'YGrid','on');

drawnow
end

figure('units','normalized','outerposition',[0 0 1 1])
set(gcf,'DefaultAxesColorOrder',Col, 'defaultAxesLineStyleOrder','-|--|:|-.');

% Evolution in time of the conserved quantities
plot(time,[Z V I E])
legend(...
    'mass (original GN)','mass (improved GN)','mass (regularized GN)',...
    'velocity (original GN)','velocity (improved GN)','velocity (regularized GN)',...
    'impulse (original GN)','impulse (improved GN)','impulse (regularized GN)',...
    'energy (original GN)','energy (improved GN)','energy (regularized GN)',...
    'Location','BestOutside')
xlabel('t')
title('Conserved quantities versus time')