

clear all
format long g
tic
%Cavity
%parameters
A=1; %aspect ratio
Ra=10;
RaI=0;
RaE=1;
Q=RaI/RaE;
Ha=0;
g=0;
%discretization and initialization%%%%%%%%%%
%discretization%%%%%%%%%%
[x,y,N]=BLGrid(0.2,15);
for i=1:N-1
    dx(i)=x(i+1)-x(i);
    dy(i)=y(i+1)-y(i);
end
dx(1)
N %%number of nodes in x- direction

%  N=5; %number of nodes in x- and y- direction
%  dx=ones(N-1,1)*1/(N-1);
%  dy=dx;
%  

M=N;%%number of nodes in y- direction

dt=dx(1)/5;

uo=zeros(N,M);un=zeros(N,M);%u_old,u_new
uo(2:N-1,2:M-1)=ones(N-2,M-2);

To=zeros(N,M);Ti=zeros(N,M);Tn=zeros(N,M);%T_old,T_intermediar, T_new
To(1,:)=0.5;To(N,:)=-0.5;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%coefficients in energy discretized equation
%coefficients in momentum and energy discretized equations
for i=2:N-1
    c1(i)=2/dx(i)/dx(i-1);
    c2(i)=2/dy(i)/dy(i-1);
    c3(i)=2/dx(i)/(dx(i)+dx(i-1));
    c4(i)=2/dx(i-1)/(dx(i)+dx(i-1));
    c5(i)=2/dy(i)/(dy(i)+dy(i-1));
    c6(i)=2/dy(i-1)/(dy(i)+dy(i-1));
    c7(i)=Ra/(dx(i)+dx(i-1));
    c8(i)=(dx(i)+dx(i-1));
    c9(i)=(dy(i)+dy(i-1));
end;
a1=c4/2;a3=c1/2;a4=c3/2;a5=A*A*a1;a6=A*A*a3;a7=A*A*a4;

for i=2:N-1
    for j=2:N-1
        Ku1(i,j)=c3(i)*(1+Ha*cos(g)*cos(g))/(c1(i)*(1+Ha*cos(g)*cos(g))+c2(j)*(A*A+Ha*sin(g)*sin(g)));
        Ku2(i,j)=c4(i)*(1+Ha*cos(g)*cos(g))/(c1(i)*(1+Ha*cos(g)*cos(g))+c2(j)*(A*A+Ha*sin(g)*sin(g)));
        Ku3(i,j)=c5(j)*(A*A+Ha*sin(g)*sin(g))/(c1(i)*(1+Ha*cos(g)*cos(g))+c2(j)*(A*A+Ha*sin(g)*sin(g)));
        Ku4(i,j)=c6(j)*(A*A+Ha*sin(g)*sin(g))/(c1(i)*(1+Ha*cos(g)*cos(g))+c2(j)*(A*A+Ha*sin(g)*sin(g)));
        Ku5(i,j)=c7(i)/(c1(i)*(1+Ha*cos(g)*cos(g))+c2(j)*(A*A+Ha*sin(g)*sin(g)));
        Ku6(i,j)=2*A*Ha*sin(g)*cos(g)/c8(i)/c9(j)/(c1(i)*(1+Ha*cos(g)*cos(g))+c2(j)*(A*A+Ha*sin(g)*sin(g)));
        a2(i,j)=0.5*A/(dx(i)+dx(i-1))/(dy(j)+dy(j-1));
    end
end 


kk=0;
t=0;
tf=0.5;
titer=1;
steady=1;
errT=1;
ur=1.8;
i_Nu=0;
while (errT>1e-6)

t=t+dt;
clear Ax bx
%solve temperature
%implicit in x direction
    for j=2:M-1
    %%%%%%%solve system for each fixed j %%%%%%%%%%%%%%%%
        
        Ax=zeros(N,N);bx=zeros(N,1);
        %boundary conditions for x=0
        Ax(1,1)=1;
        %bx(1)=0;
        bx(1)=0.5;
        %%%%%%%%%system%%%%%%%%%%%%%                
        for i=2:N-1
         Ax(i,i-1)=-a1(i)-a2(i,j)*(uo(i,j+1)-uo(i,j-1)); Ax(i,i)=1/dt+a3(i); Ax(i,i+1)=-a4(i)+a2(i,j)*(uo(i,j+1)-uo(i,j-1));
         bx(i)=To(i,j-1)*(a5(j)-a2(i,j)*(uo(i+1,j)-uo(i-1,j)))+To(i,j)*(1/dt-a6(j))+To(i,j+1)*(a7(j)+a2(i,j)*(uo(i+1,j)-uo(i-1,j)))+0.5*Q;
        end
        %boundary conditions for x=1
        Ax(N,N)=1;
        %bx(N)=0;
        bx(N)=-0.5;

        Ti(:,j)=Ax\bx;
    end;
    
    clear Ax bx
%implicit in y direction
    for i=2:N-1
    %%%%%%%solve system for each fixed i %%%%%%%%%%%%%%%%
        
        Ax=zeros(M,M);bx=zeros(M,1);
        %boundary conditions for y=0
        Ax(1,1)=1;Ax(1,2)=-1;bx(1)=0;
        %%%%%%%%%system%%%%%%%%%%%%%                
        for j=2:M-1
         Ax(j,j-1)=-a5(j)+a2(i,j)*(uo(i+1,j)-uo(i-1,j));Ax(j,j)=1/dt+a6(j); Ax(j,j+1)=-a7(j)-a2(i,j)*(uo(i+1,j)-uo(i-1,j));
         bx(j)=Ti(i-1,j)*(a1(i)+a2(i,j)*(uo(i,j+1)-uo(i,j-1)))+Ti(i,j)*(1/dt-a3(i))+Ti(i+1,j)*(a4(i)-a2(i,j)*(uo(i,j+1)-uo(i,j-1)))+0.5*Q;
        end
        %boundary conditions for x=0
        Ax(M,M-1)=-1;Ax(M,M)=1;bx(M)=0;
        Tn(i,:)=(Ax\bx)';
    end;
   Tn(1,:)=0.5;Tn(N,:)=-0.5;
  errT=norm(To-Tn);
%   steady=abs(max(max(To-Tn))/max(max(Tn)));
   To=Tn;

     %solve stream
     nr_it=0;err_u=1;
   while (err_u>1e-6)
   nr_it=nr_it+1;   
   err_u=0;  
    for i=2:N-1
        for j=2:M-1
           un(i,j)=Ku1(i,j)*uo(i+1,j)+Ku2(i,j)*uo(i-1,j)+Ku3(i,j)*uo(i,j+1)+Ku4(i,j)*uo(i,j-1)+Ku5(i,j)*(To(i+1,j)-To(i-1,j))...
                 +Ku6(i,j)*(uo(i+1,j+1)-uo(i+1,j-1)-uo(i-1,j+1)+uo(i-1,j-1));
            un(i,j)=uo(i,j)+ur*(un(i,j)-uo(i,j));
         
          if abs(uo(i,j)-un(i,j))>err_u 
                err_u=abs(uo(i,j)-un(i,j));
            end%if
         
            uo(i,j)=un(i,j);
        end%for j
    end;%for i
   
   uo=un;
    
%    if mod(nr_it,10)==0
%     fprintf('nr_it=%g err_u=%g err_T=%g\n', nr_it, err_u, errT);
%    end
   
  end %while stream
   
   titer=titer+1;
   
    if mod(titer,10)==0
        kk=kk+1;
       i_Nu=i_Nu+1;
       timp(i_Nu)=t;
       Nuloc=(Tn(2,:)-Tn(1,:))/dx(1);
       Nu(i_Nu)=trapz(y,Nuloc);
        fprintf('\nt=%g  errT=%g \n', t, errT);
%        figure(1)
%        contourf(x,y,To',20)
%        axis equal
%        m1(kk)=getframe;
%        figure(2)
%  contourf(x,y,uo',20)
%  axis equal
%  m2(kk)=getframe;
        
    end

end%while main    

%x=0:dx(1):(N-1)*dx(1);y=x;
% plot(x, Tn(:,floor(M/2)+1),'k ')    

Nuloc=(Tn(2,:)-Tn(1,:))/dx(1);
i_Nu=i_Nu+1;
timp(i_Nu)=t;
Nuloc=(Tn(2,:)-Tn(1,:))/dx(1);
Nu(i_Nu)=trapz(y,Nuloc); %using trapezoidal rule
% figure(3)
% hold on;
% plot(timp(2:i_Nu),-Nu(2:i_Nu))
% save Ra1000_87x87

 
 
 figure(1)
 contour(x,y,To',20)
 axis equal
 axis([0,1,0,1])
 
 figure(2)
 contour(x,y,uo',20)
 axis equal
 axis([0,1,0,1])
toc