%{
# Free-surface: displaying the eigenmodes using particles

In this code, we take the code that computes the eigenmodes for free surface waves [free_surface_gravity.m]() and we visualize the flow motion of the first eigenmode by expanding the solution to physical space and advecting tracer particles, like we did in a simple case in [particles.m]().
%}

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)=[];
%{
# The animation

We show the position of the free surface stored in *q(eta)*, and the particles are initialy set as a mesh, and stretched in *y* to fit with the initial position of the free surface for the top prticles.
%}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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=2;  % 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(:);  

% time loop
for ind=1:nper*nt
    
    
%{
# Expand the solution to physical space

We have assumed a wavelike behaviour of the flow solution in the $x$ direction, thus we have writen for all the flow variables, for instance for the $u$ velocity
$$
u(x,y,t)=\hat{u(y)}\exp(\alpha x+s t)+CC
$$
where $CC$ is the complex conjugate of the other term. This way we have transformed a variable $u$ depending on three coordinates, to a variable $\hat{u}$ depending on one single coordinate, thus transforming everything into a 1D problem. What we do now is to transform back the 1D problem into a full original problem by using the exponential. 

A complex quantity plus its complex conjugate gives twice the real part of it, this is what we code here.

In the code, we do the expansion to hysical space for all the variables $(u,v,p,\eta)$ at once in a single command.

%}
    % 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.');

%{
# Interpolation of the particle positions

The particles are not on a cartesian grid, they can be wherever the field has brought them, the only thing that matter here for the interpolation is where they are on the $y$ grid, and then we can do the expansion in $x$ as a second step, this is less computations.
%}
    % compute velocity at position of particles
    pu=interp1(y,q(u),py);
    pv=interp1(y,q(v),py);

%{
# The Taylor expansion

When we wrote the boundary conditionat the free surface for setting up the linear model, we have used Taylor expansions to express the velocity and pressure at the position of the moving interface as a function of the values at the fixed position $L$. here we do again the same thing for the particles that are above $L$.

%}

    % 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'); 
    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_particles.png');

%{
![The figure at final iteration of the time loop, shwing the free surface at the top, the particles in black, and the disturbance pressure field in color ](/sandbox/easystab/free_surface_gravity_particles.png)

# Exercices/Contributions

* Please make a movie instead of just an image to display on the page --> [free_surface_gravity_particles_gif.m]()
* Please play with the parameters of the particles, and the amplitude of the eigenmode
* Please check that the first and second eigenmodes correspond to waves propagating in opposite directions -->[free_surface_gravity_first_modes.m](/sandbox/easystab/stab2014/free_surface_gravity_first_modes)
* Please use this animation to check the theory for the wave speed of gravity waves
* Please plot also the vector field using quiver in a way that confirms the motion of the particles that we observe -->[free_surface_gravity_vector_field.m](/sandbox/easystab/stab2014/free_surface_gravity_vector_field)
* Please investigate the motion of the other eigenmodes (which correspond to damped stationnary waves, for now, there is a bug that makes things strange for the non propagative modes...)
* Please consider the unstable case with *g* negative and find the nice scaling of the mode amplitude and the time of the animation so that we nicely see the flow in this instability
* Please investigate how the Taylor expansino is not accurate when the surface deformation is large (you can see this with the particles that are close to the interface)
* Please do the animation showing the total pressure instead of only the pressure disturbance (that is, add the hydrostatic pressure field as well to see what this changes)[free_surface_gravity_total_pressure](/sandbox/easystab/stab2014/free_surface_gravity_total_pressure)
* Please do the same visualisation with particle for the other flow instabilities
%}