clear all
close all
clc

% Numero di nodi
Nx = 80;
Ny = 80;

% Parametri della griglia
Lx = 1;
Ly = 1;

Dx = Lx/(Nx-1);
Dy = Ly/(Ny-1);

% Parametri simulazione
Pr = [1     1       1       1       1       0.1     0.5     1       2       5   ];
Ra = [10^3  10^4    10^5    10^6    5*10^2  10^4    10^4    10^4    10^4    10^4];
Dt = 0.001;
TIMEMAX = 100000;

% Temperature delle pareti verticali
TC = -0.5;
TH = 0.5;

% Griglia computazionale
X = zeros(1,Nx);
for i=1:(Nx-1)
    X(i+1)=X(i)+Dx;
end
Y = zeros(1,Ny);
for j=1:(Ny-1)
    Y(j+1)=Y(j)+Dy;
end

% simulazione
for caso = 1:10
       
% Inizializzazione vettori
U = zeros(Nx,Ny);
V = zeros(Nx,Ny);
T = zeros(Nx,Ny);
T2 = zeros(Nx,Ny);
PSI = zeros(Nx,Ny);
VOR = zeros(Nx,Ny);
VOR2 = zeros(Nx,Ny);

% Condizioni al contorno temperatura
for j=1:Ny
    T2(Nx,j) = TH;
    T2(1,j) = TC;
end

% Contatori
EPS = 1.e-6;
TIME = 0;
NODI = Nx*Ny;
NUMERO = 0;     % Nodi convergenti
ITERTOT = 0;

% Solutore

while NUMERO<NODI & TIME<TIMEMAX

    TIME = TIME+Dt
    
    % Risoluzione dell'equazione di Poisson
    W = 1.9;
    EPSILP = 1.e-6;
    ITERMAX = 100000; %1000
    ITER = 0;
    NUM = 0;
    NUMEQ = (Nx-2)*(Ny-2);
    
    while NUM<NUMEQ & ITER<ITERMAX
    
        NUM = 0;

        for i=2:(Nx-1)
            for j=2:(Ny-1)

                PSI2 = ( 1/( 2*(Dx^2+Dy^2) ) ) * ( ( PSI(i+1,j)+PSI(i-1,j) )*Dy^2 + ( PSI(i,j+1)+PSI(i,j-1) )*Dx^2 - VOR2(i,j)*Dx^2*Dy^2 );
                PSI2 = (1-W)*PSI(i,j)+W*PSI2;
                ERROR1 = abs( PSI2-PSI(i,j) );
                PSI(i,j) = PSI2;

                if ERROR1<EPSILP
                    NUM = NUM+1;
                end

            end
        end

        ITER = ITER+1;

    end
    
    % Calcolo vorticità al bordo espandendo in serie di Taylor la funzione di corrente;
    
    % Bordi orizzontali
    for i=1:Nx
        VOR2(i,1) = (0.5/Dy^2)*(8*PSI(i,2)-PSI(i,3));          % bordo basso
        VOR2(i,Ny) = (0.5/Dy^2)*(8*PSI(i,Ny-1)-PSI(i,Ny-2));   % bordo alto
    end
    
    % Bordi verticali
    for j=1:Ny
        VOR2(1,j) = (0.5/Dx^2)*(8*PSI(2,j)-PSI(3,j));         % bordo sinistro
        VOR2(Nx,j) = (0.5/Dx^2)*(8*PSI(Nx-1,j)-PSI(Nx-2,j));  % bordo destro
    end

    % Calcolo campo di moto
    for i=2:(Nx-1)
        for j=2:(Ny-1)
            U(i,j) = (PSI(i,j+1)-PSI(i,j-1))/(2*Dy);
            V(i,j) = -(PSI(i+1,j)-PSI(i-1,j))/(2*Dx);
        end
    end
    
    % Aggiornamento campi per avanzamento temporale
    VOR = VOR2;
    T = T2;
    
    % Calcolo della vorticità nel nucleo del dominio all'istante successivo attraverso l'equazione di trasporto della vorticità
    
    for i=2:(Nx-1)
        for j=2:(Ny-1)
            
            CONV(i,j) = 0.5/Dx * ( U(i+1,j)*VOR(i+1,j)-U(i-1,j)*VOR(i-1,j) ) + 0.5/Dy * ( V(i,j+1)*VOR(i,j+1)-V(i,j-1)*VOR(i,j-1) );
            DIFF(i,j) = ( VOR(i+1,j)-2*VOR(i,j)+VOR(i-1,j) )/Dx^2 + ( VOR(i,j+1)-2*VOR(i,j)+VOR(i,j-1) )/Dy^2;
            BUOY(i,j) = ( T(i+1,j)-T(i-1,j) )/Dx;
            
            RDIFF = Pr(caso)/(Ra(caso)^0.5);
            RBUOY = Pr(caso);
            
            VOR2(i,j) = VOR(i,j) + Dt*( -CONV(i,j) + RDIFF*DIFF(i,j) + RBUOY*BUOY(i,j) );

        end
    end 
    
    % Calcolo del campo termico
    
    % Bordi orizzontali - dt/dy=0
    for i=1:Nx
        T2(i,1) = T(i,2);       % bordo basso
        T2(i,Ny) = T(i,Ny-1);   % bordo alto
    end
    
    % Bordi verticali - T=Tfix
    for j=1:Ny
        T2(1,j) = TC;       % bordo sinistro
        T2(Nx,j) = TH;      % bordo destro
    end
    
    % Calcolo della temperatura nel nucleo del dominio all'istante sucessivo attraverso l'equazione dell'energia

    for i=2:(Nx-1)
        for j=2:(Ny-1)
            
            CONVT(i,j) = 0.5/Dx * ( U(i+1,j)*T(i+1,j)-U(i-1,j)*T(i-1,j) ) + 0.5/Dy * ( V(i,j+1)*T(i,j+1)-V(i,j-1)*T(i,j-1) );
            DIFFT(i,j) = ( T(i+1,j)-2*T(i,j)+T(i-1,j) )/Dx^2 + ( T(i,j+1)-2*T(i,j)+T(i,j-1) )/Dy^2;
            
            T2(i,j) = T(i,j) + Dt*( -CONVT(i,j) + 1/(Ra(caso)^0.5)*DIFFT(i,j) );
 
        end
    end 
    
    % Calcolo dell'errore sulla vorticità
    
    NUMERO = 0;
    ERRTOT = 0;
    
    for i=1:Nx
        for j=1:Ny
    
            ERROR = abs( VOR2(i,j)-VOR(i,j) );
            ERRTOT = ERRTOT+ERROR;
            
            if ERROR<EPS
                NUMERO = NUMERO+1;
            end  
            
        end
    end
    
    MEANERR = ERRTOT/NODI;
    
end

Nconv = NUMERO/(40*40)*100

% Nusselt TH
for j = 1:Ny
    NH(j) = ( T(Nx,j)-T(Nx-1,j) )/( X(Nx)-X(Nx-1) );
end

% Nusselt TC
for j = 1:Ny
    NC(j) = ( T(2,j)-T(1,j) )/( X(2)-X(1) );
end

% Salvataggio risultati
save( strcat(num2str(caso),'_PSI.txt') , 'PSI' , '-ascii' );
save( strcat(num2str(caso),'_T.txt') , 'T' , '-ascii' );
save( strcat(num2str(caso),'_U.txt') , 'U' , '-ascii' );
save( strcat(num2str(caso),'_V.txt') , 'V' , '-ascii' );
save( strcat(num2str(caso),'_VOR.txt') , 'VOR' , '-ascii' );
save( strcat(num2str(caso),'_NH.txt') , 'NH' , '-ascii' );
save( strcat(num2str(caso),'_NC.txt') , 'NC' , '-ascii' );
nomefile = fopen(strcat(num2str(caso),'_REPORT.txt'),'wt');
    fprintf(nomefile,'Ra = %d \n',Ra(caso));
    fprintf(nomefile,'Pr = %d \n',Pr(caso));
    fprintf(nomefile,'Griglia = %d \n',Nx);
    fprintf(nomefile,'Sovrarilassamento = %d \n',W);
    fprintf(nomefile,'Dt = %d \n',Dt);
    fprintf(nomefile,'tempo tot = %d \n',TIME);
    fprintf(nomefile,'percentuale di nodi convergenti = %d \n',Nconv);
fclose(nomefile);

end