% *************************************************
%      resolution of the cavity flow problem      *
%           schema explicite centre d'ordre 2     *
%           pour la vorticite                     *
%           resolution directe du laplacien       *
%           pour la fonction courant.             *
% *************************************************
% * M. Garbey 00                                  *
% *************************************************
N=Ncell+1;
% space grid:
if cas_test==1,
   Lx=2.*pi; Ly=2.*pi;
else
   Lx=1.; Ly=1.;
end

hx=Lx/(N-1);hy=Ly/(N-1);

% time step:
deltaT_Diffusion=0.9*1/4*min(hx,hy)^2*Re;
deltaT=deltaT_Diffusion;
for i=1:N,
   x(i)=(i-1)*hx; y(i)=(i-1)*hy;
end
% initialisation:
newomega=omega;

% boundary condition on the sliding wall:
if cas_test==1, % exact analytical solution
   for j=1:N,
      g(j)=sin(y(j));
   end
end
if cas_test==2, % bell function for sliding wall velocity
	for j=1:N,
 		g(j)=-16*x(j)^2*(1.-x(j))^2;
   end
end
if cas_test==3, % uniform sliding wall velocity
	for j=1:N,
 		g(j)=-1.;
   end
end

% getting the right hand side corresponding to the BC:
if cas_test==1, % exact analytical solution.
   for i=1:N,
      for j=1:N,
         rhsomega(i,j)=4./Re*sin(x(j))*sin(y(i));
      end
   end
end
if cas_test>1, % numerical solution for sliding wall.
   rhsomega=zeros(N,N);
end
if cas_test==1,
   for i=1:N,
      for j=1:N,
         omega_exact(i,j)=2.*sin(x(j))*sin(y(i));
      end
   end
   psi_exact=0.5.*omega_exact;
end

temps=0.;
% operators for the Poisson problem:
M=N-2;
M2=(N-2)*(N-2);
DC=spalloc(M2,M2,5*M2);

for i=1:M2,
   DC(i,i)=(-2./hx/hx-2./hy/hy);
   if i>1, DC(i,i-1)=1./hy/hy; end;
   if i<M2, DC(i,i+1)=1./hy/hy; end;
   if (i+M) < M2+1, DC(i,i+M)=1./hx/hx; end;
   if (i-M) > 0, DC(i,i-M)=1./hx/hx; end;
end;
for i=1:M-1,
   DC(i*M,i*M+1)=0;
   DC(i*M+1,i*M)=0;
end;
if solveur==1,
   [LP,UP,PP]=lu(DC);
end
if solveur==2,
   [LP,UP]=luinc(DC,1e-6);
end
cpt=0; % ATTENTION
test=1e-9;
flops(0);
test_iteration=tol+1; % initialisation of test.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
kt=0;
while test_iteration>tol, % test on steady solution.
   temps=temps+deltaT; % update time
   kt=kt+1;
   % Boundary conditions:
   if cas_test==1,
      if CL==2,
         omega(1,2:N-1)=1./(2.*hy^2)*(psi(3,2:N-1)-8.*psi(2,2:N-1))+3./hy*g(2:N-1);
 	      omega(N,2:N-1)=1./(2.*hy^2)*(psi(N-2,2:N-1)-8.*psi(N-1,2:N-1))-3./hy*g(2:N-1);
 	      omega(2:N-1,1)=1./(2.*hx^2)*(psi(2:N-1,3)-8.*psi(2:N-1,2))+3./hx*g(2:N-1)';
         omega(2:N-1,N)=1./(2.*hx^2)*(psi(2:N-1,N-2)-8.*psi(2:N-1,N-1))-3./hx*g(2:N-1)';
      end 
      if CL==1,
         omega(1,2:N-1)=-2./hy^2*psi(2,2:N-1)+2./hy*g(2:N-1);
         omega(N,2:N-1)=-2./hy^2*psi(N-1,2:N-1)-2./hy*g(2:N-1);
         omega(2:N-1,1)=-2./hx^2*psi(2:N-1,2)+2./hy*g(2:N-1)';
         omega(2:N-1,N)=-2./hx^2*psi(2:N-1,N-1)-2./hy*g(2:N-1)';
      end
      if CL==3,
         omega(1,2:N-1)=-3./hy^2*psi(2,2:N-1)-0.5*omega(2,2:N-1)+3./hy*g(2:N-1)+hy/2.*g(2:N-1);
         omega(N,2:N-1)=-3./hy^2*psi(N-1,2:N-1)-0.5*omega(N-1,2:N-1)-3./hy*g(2:N-1)-hy/2.*g(2:N-1);
         omega(2:N-1,1)=-3./hx^2*psi(2:N-1,2)-0.5*omega(2:N-1,2)+3./hy*g(2:N-1)'+hy/2.*g(2:N-1)';
         omega(2:N-1,N)=-3./hx^2*psi(2:N-1,N-1)-0.5*omega(2:N-1,N-1)-3./hy*g(2:N-1)'-hy/2.*g(2:N-1)';
      end
   end
   if cas_test>1,
      if CL==1,
         omega(1,2:N-1)=-2./hy^2*psi(2,2:N-1)+2./hy*g(2:N-1);
         omega(N,2:N-1)=-2./hy^2*psi(N-1,2:N-1);
         omega(2:N-1,1)=-2./hx^2*psi(2:N-1,2);
         omega(2:N-1,N)=-2./hx^2*psi(2:N-1,N-1);
      end
      if CL==2,
 	      omega(1,2:N-1)=1./(2.*hy^2)*(psi(3,2:N-1)-8.*psi(2,2:N-1))+3./hy*g(2:N-1);
 	      omega(N,2:N-1)=1./(2.*hy^2)*(psi(N-2,2:N-1)-8.*psi(N-1,2:N-1));
 	      omega(2:N-1,1)=1./(2.*hy^2)*(psi(2:N-1,3)-8.*psi(2:N-1,2));
         omega(2:N-1,N)=1./(2.*hy^2)*(psi(2:N-1,N-2)-8.*psi(2:N-1,N-1));
      end 
   end

   % finite difference computation of derivatives:
   for i=2:N-1,
      for j=2:N-1,
         psiy(i,j)=(psi(i+1,j)-psi(i-1,j))/(2*hy);
         psix(i,j)=(psi(i,j+1)-psi(i,j-1))/(2*hx);
         omegay(i,j)=(omega(i+1,j)-omega(i-1,j))/(2*hy);
         omegax(i,j)=(omega(i,j+1)-omega(i,j-1))/(2*hx);
         
         % check time step for stability condition:
         CFL=max(max(max(abs(psix))),max(max(abs(psiy))));
         if CFL<=0.
            deltaT=0.01*min(hx,hy);
         else
            deltaT=0.25*min(hx,hy)/CFL;
         end
         if deltaT>1/4*min(hx,hy)^2*Re,
            deltaT=deltaT_Diffusion;
         end
         convect(i,j)=psiy(i,j)*omegax(i,j)-psix(i,j)*omegay(i,j);
         omegayy(i,j)=(omega(i+1,j)+omega(i-1,j)-2*omega(i,j))/hy^2;
         omegaxx(i,j)=(omega(i,j+1)+omega(i,j-1)-2*omega(i,j))/hx^2;
         diffu(i,j)=1./Re*(omegaxx(i,j)+omegayy(i,j));
      end
   end
   % Time stepping for vorticity:
   if verif ==1,
      newomega(2:N-1,2:N-1)=omega(2:N-1,2:N-1)-deltaT*convect(2:N-1,2:N-1)+deltaT*(diffu(2:N-1,2:N-1)+rhsomega(2:N-1,2:N-1));
   else
      omega(2:N-1,2:N-1)=omega(2:N-1,2:N-1)-deltaT*convect(2:N-1,2:N-1)+deltaT*(diffu(2:N-1,2:N-1)+rhsomega(2:N-1,2:N-1));
   end
   % verification of the vorticity solution:
   if verif ==1,
      for i=2:N-1,
         for j=2:N-1,
            omegayy(i,j)=(omega(i+1,j)+omega(i-1,j)-2*omega(i,j))/hy^2;
            omegaxx(i,j)=(omega(i,j+1)+omega(i,j-1)-2*omega(i,j))/hx^2;
            argu=(newomega(i,j)-omega(i,j))/deltaT;
            argu=argu-(omegaxx(i,j)+omegayy(i,j))/Re-rhsomega(i,j);
            argu=argu+psiy(i,j)*omegax(i,j)-psix(i,j)*omegay(i,j);
            residus1(i,j)=argu;
         end
      end
      omega(2:N-1,2:N-1)=newomega(2:N-1,2:N-1);
   end
   histoire(kt,1)=max(max(abs(convect(2:N-1,2:N-1)-diffu(2:N-1,2:N-1)-rhsomega(2:N-1,2:N-1))));
   test_iteration=histoire(kt,1);
   
   histoire(kt,2)=deltaT;
   histoire(kt,1)
   if cas_test==1,
      histoire(kt,3)=max(max(max(abs(omega-omega_exact))),max(max(abs(psi-psi_exact))));
   end

   % resolution of the Poisson problem:
   RHSC=reshape(omega(2:N-1,2:N-1),M2,1);
   % choice is for direct solver or iterative solver:
   if solveur==1, % LU decomposition solver.
      % newpsi=DC\RHSC;
      newpsi=UP\(LP\(PP*RHSC));
   end
   if solveur==2, % preconjugate gradient solver.
      if choix==0,
         test=1e-9;
      else
         test=min(histoire(kt,1)*min(hx,hy)^(1.5),1e-6);
      end
      newpsi=pcg(DC,RHSC,test,10,LP,UP);
   end
   psi=zeros(N);
   psi(2:N-1,2:N-1)=-reshape(newpsi,M,M);
   % verification of the Poisson solution:
   if verif ==1,
      for i=2:N-1,
         for j=2:N-1,
            psiyy(i,j)=(psi(i+1,j)+psi(i-1,j)-2*psi(i,j))/hy^2;
            psixx(i,j)=(psi(i,j+1)+psi(i,j-1)-2*psi(i,j))/hx^2;
            residus(i,j)=psixx(i,j)+psiyy(i,j)+omega(i,j);
         end
      end
   end
end % end time stepping.
cpt=cpt+flops;
cpt=round(cpt/10^6)
figure(1)
subplot(2,2,1)
contour(omega,10);
title('vorticity');
subplot(2,2,2)
contour(psi,10);
title('courant');
subplot(2,2,3);
plot(log10(histoire));
title('steady solution?');
if verif==1,
   subplot(2,2,4);
   mesh(residus);
   title('verification?');
   pause
   subplot(2,2,4);
   mesh(residus1);
   title('verification?');
else
   subplot(2,2,4);
   mesh(convect-diffu);
   title('stationary solution?');
end

% computation of the numerical error when the analytical solution is available: 
if cas_test==1,
   AA=max(max(abs(omega-omega_exact)))
   BB=max(max(abs(psi-psi_exact)))
end