%{
# Free-surface-gravity: displaying total pressure

In this code, we take the code that display the animation of the free surface waves [free_surface_gravity_particles.m](/sandbox/easystab/free_surface_gravity_particles.m) and we show total pressure instead of only the pressure disturbance (it means that we add hydrostatic pressure to the pressure field) to see what changes.

In a fluid, hydrostatic pressure is determined by the formula below:
$$
p=\rho.g.y
$$
where:

- p is the hydrostatic pressure (Pa)

- $\rho$ is the fluid density ($kg/m^{3}$)

- $g$ is the gravational acceleration ($m/s^{2}$)

- $y$ is the height of the point ($m$)

So we will add the hydrostatic pressure $p$ to our pressure field to obtain the total pressure in our case.

Dependency:

* [chebdif.m](/sandbox/easystab/chebdif.m) for the Chebychev differentiation matrices
%}

clear all; clf;

n=100;      % number of gridpoints
alpha=1;    % wavenumber in x
L=2;        % Fluid height in y
rho=1;      % fluid density
mu=0.0001;    % fuid viscosity
g=1;        % gravity

% differentiation matrices
scale=-2/L;
[y,DM] = chebdif(n,2); 
D=DM(:,:,1)*scale;    
DD=DM(:,:,2)*scale^2;    
y=(y-1)/scale; 
I=eye(n); Z=zeros(n,n);

% renaming the matrices
dy=D; dyy=DD;
dx=i*alpha*I; dxx=-alpha^2*I;
Delta=dxx+dyy;

% System matrices
A=[mu*Delta, Z, -dx, Z(:,1); ...
   Z, mu*Delta, -dy, Z(:,1); ...
   dx, dy, Z, Z(:,1); ...
   Z(1,:),I(n,:),Z(1,:),0];

E=blkdiag(rho*I,rho*I,Z,1);

% boundary conditions
loc=[1,n,n+1,2*n];
C=[I(1,:),Z(1,:),Z(1,:),0; ... 
   Z(1,:),I(1,:),Z(1,:),0; ...
   Z(1,:),Z(1,:),-I(n,:),rho*g; ...
   dy(n,:),dx(n,:),Z(1,:),0]; 

E(loc,:)=0;  
A(loc,:)=C;

% compute eigenmodes
[U,S]=eig(A,E);
s=diag(S);  [t,o]=sort(-real(s)); s=s(o); U=U(:,o);
rem=abs(s)>1000; s(rem)=[]; U(:,rem)=[];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% particle animation of the eigenmodes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% parameters
nper=1;      % number of periods of oscillation
nt=30;       % number of time steps per period
nx=40;       % number of points in x
modesel=1;  % which more do animate


% select the eigenmode
u=1:n; v=u+n; p=v+n; eta=3*n+1;
q=U(:,modesel); 
lambda=s(modesel);

% The time and x extent
tvec=linspace(0,nper*2*pi/abs(lambda),nper*nt); 
dt=tvec(2)-tvec(1);
Lx=2*pi/alpha; 
x=linspace(0,Lx,nx);

% scale mode amplitude 
q=0.05*q/abs(q(eta)); 

% initialize tracer particles
[px,py]=meshgrid(linspace(0,Lx,60),linspace(0,L,30)); 
py=py.*(1+2*real(exp(lambda*tvec(1))*q(eta)*exp(i*alpha*px))/L);
px=px(:);py=py(:);  

figure(1)
% time loop
for ind=1:nper*nt
    
    % expand mode to physical space 
    qq=2*real(exp(lambda*tvec(ind))*q*exp(i*alpha*x));

    % plot pressure
    surf(x,y,qq(p,:)-10+kron(ones(1,nx),rho*g*y),'facealpha',0.3); view(2); shading interp; hold on
    
    % plot free surface
     plot(x,L+qq(eta,:),'k-',x,0*x+L,'k--'); hold on
    
    %%%% plot the particles
    plot(mod(px,Lx),py,'k.');
    
    pu=interp1(y,q(u),py);
    pv=interp1(y,q(v),py);

    % For particles above L, use Taylor expansion for velocity
    pu(py>L)=q(u(n))+q(eta)*dy(n,:)*q(u);
    pv(py>L)=q(v(n))+q(eta)*dy(n,:)*q(v);
    
    % expand to physical space
    puu=2*real(exp(lambda*tvec(ind))*pu.*exp(i*alpha*px));
    pvv=2*real(exp(lambda*tvec(ind))*pv.*exp(i*alpha*px));
 
    % advect particles
    px=px+puu*dt;
    py=py+pvv*dt;
    
    xlabel('x');    
    ylabel('y');
    title('Free surface gravity with total pressure - Mode 1');
    axis equal; axis([0,Lx,0,1.3*L]); 
    grid off
    hold off
    drawnow

end

set(gcf,'paperpositionmode','auto');
print('-dpng','-r100','free_surface_gravity_total_pressure.png');

figure(2)

% initialize tracer particles
[px,py]=meshgrid(linspace(0,Lx,60),linspace(0,L,30)); 
py=py.*(1+2*real(exp(lambda*tvec(1))*q(eta)*exp(i*alpha*px))/L);
px=px(:);py=py(:); 

for ind=1:nper*nt
    
    % expand mode to physical space 
    qq=2*real(exp(lambda*tvec(ind))*q*exp(i*alpha*x));
    
    % plot pressure
    surf(x,y,qq(p,:)-10,'facealpha',0.3); view(2); shading interp; hold on
    
    % plot free surface
     plot(x,L+qq(eta,:),'k-',x,0*x+L,'k--'); hold on
    
    %%%% plot the particles
    plot(mod(px,Lx),py,'k.');
    
    pu=interp1(y,q(u),py);
    pv=interp1(y,q(v),py);

    % For particles above L, use Taylor expansion for velocity
    pu(py>L)=q(u(n))+q(eta)*dy(n,:)*q(u);
    pv(py>L)=q(v(n))+q(eta)*dy(n,:)*q(v);
    
    % expand to physical space
    puu=2*real(exp(lambda*tvec(ind))*pu.*exp(i*alpha*px));
    pvv=2*real(exp(lambda*tvec(ind))*pv.*exp(i*alpha*px));
 
    % advect particles
    px=px+puu*dt;
    py=py+pvv*dt;
    
    xlabel('x');    
    ylabel('y');
    title('Free surface gravity with pressure disturbance - Mode 1');
    axis equal; axis([0,Lx,0,1.3*L]); 
    grid off
    hold off
    drawnow   
    
end

set(gcf,'paperpositionmode','auto');
print('-dpng','-r100','free_surface_gravity_pressure_disturbance.png');

%{
#Figures

![Pressure disturbance](/sandbox/easystab/stab2014/free_surface_gravity_pressure_disturbance.png)

![Total pressure](/sandbox/easystab/stab2014/free_surface_gravity_total_pressure.png)

So we see that because the value of hydrostatic pressure is much higher than the pressure disturbance so the total pressure is almost the same as the hydrostatic pressure, it only depends on the height of the points.
%}