%HEAT1D: Esempio di utilizzo di matrici sparse.
% ---------------------------------------------------------------  
%  In questo script MATLAB, si risolve il transitorio termico
%  monodimensionale di una lastra piana, per mezzo del metodo 
%  dei Volumi Finiti, utilizzando (per semplicità) un algoritmo
%  implicito di Eulero del I ordine.
% ---------------------------------------------------------------
   clear; clc; echo off; help heat1d;
   disp('premere un tasto per continuare...'), disp(' '), pause;
   
disp('Desideriamo studiare come varia la temperatura') 
disp('nel tempo per una parete (lastra) piana che,') 
disp('inizialmente alla temperatura tin = 200 C, viene messa')
disp('in contatto, su ambedue i lati, con aria alla')
disp('temperatura tair = 20 C.')
disp(' ')
%
   tin  = 200;
   tair = 20;
%
   disp('premere un tasto per continuare...'), disp(' '), pause;
%
clc;
disp('Lo spessore della parete è pari a L = 0.25 m, e si')
disp('ipotizza che il coefficiente convettivo alla interfaccia')
disp('aria-parate sia pari a h = 10 [W/m^2 K].')
disp(' ')
%
   L = 0.25;
   h = 10;
%
   disp('premere un tasto per continuare...'), disp(' '), pause;
clc;   
disp('Le proprietà termofisiche della parete, supposta di')
disp('mattoni comuni, sono le seguenti:')
disp(' ') 
disp('rho (massa volumica):     1920    [Kg/m^3]')    
disp('c (calore specifico):      835    [J/Kg K]')
disp('k (conducibilità termica):   0.72 [W/m K')
disp(' ')
%
   rho = 1920;
   c   = 835;
   k   = 0.72;
disp(' ')
disp(' Il valore della diffusività termica risulta quindi:')
disp(' ')
   alpha  = k/(rho*c);          
disp([' alpha = ',num2str(alpha),' [m^2/s]'])
%
   disp('premere un tasto per continuare...'), disp(' '), pause;
%
clc;
disp(' Suddividiamo la parete in M = 20 piastre elementari (nodi),')
disp(' immaginando di associare, a ciascuna parete elementare,')
disp(' un unico valore di temperatura: temp(i), i=1:M')
disp(' ')
%
   M = 20;
%
disp(' Ciascuna delle pareti elementari avrà quindi')
disp(' uno spessore pari a dL = L/(M-1):')
disp(' ')
   dL = L/(M-1);
disp([' dL = ',num2str(dL),' [m]'])
%
   disp('premere un tasto per continuare...'), disp(' '), pause;
%
clc;
disp(' Scegliamo poi il periodo di tempo di osservazione:')
disp(' tau = 3600 secondi (1 ora)')
disp(' ')
%
   tau = 3600;
%
disp(' Suddividiamo quindi il tempo di osservazione in')
disp(' N = 30 intervalli di tempo')
disp(' ')
%
   N = 30;
%
disp(' Ciascuno degli intervalli di tempo avrà quindi')
disp(' una durata pari a dtau = tau/N:')
disp(' ')
   dtau = tau/N;
disp([' dtau = ',num2str(dtau),' [s]'])
%
   disp('premere un tasto per continuare...'), disp(' '), pause;
%
clc;
disp(' Ora calcoliamo alcuni gruppi adimensionali, per')
disp(' operare con la equazione della conduzione')
disp(' in forma adimensionale:')
disp(' ')
disp(' - biot   = (h*L)/(2*k)       [numero di Biot]')
disp(' - bi_loc = h*dL/k            [numero di Biot di griglia]')
disp(' - fo     = alpha*dtau/(dL^2) [numero di Fourier]')
disp(' ')
biot   = (h*L)/(2*k);
bi_loc = h*dL/k;
fo     = alpha*dtau/(dL^2);
disp(' Si ottiene:')
%
disp([' biot   = ',num2str(biot)])
disp([' bi_loc = ',num2str(bi_loc)])
disp([' fo     = ',num2str(fo)])
disp(' ')
%
   disp('premere un tasto per continuare...'), disp(' '), pause;
%
clc;
disp(' Ora pre-allochiamo ed inizializziamo')
disp(' alcuni vettori che servono per il calcolo:')
echo on;
%Termine noto:
   rhs([1:M],1)          = zeros(M,1);
% Valore all'istante precedente:
   temp_old([1:M],1)     = tin*ones(M,1);
% Valore all'istante attuale:  
   temp([1:M],1)         = tin*ones(M,1);
% Storia della temperatura all'asse: 
   central_temp([1:N],1) = zeros(N,1);
% Storia della temperature sulla interfaccia parete-fluido:
lateral_temp([1:N],1)    = zeros(N,1);                        
echo off;
disp(' ')
%
   disp('premere un tasto per continuare...'), disp(' '), pause;
%
clc;
disp(' Ora assegnamo i valori ai vettori contenenti')
disp(' i termini non nulli della matrice')
echo on;

% Diagonale principale:
   m_diag([1:M],1) = 1.+(2.*fo.*ones(M,1));
   m_diag(1,1)     = 1.+2*fo + 2*bi_loc*fo;
   m_diag(M,1)     = 1.+2*fo + 2*bi_loc*fo;
% Sotto-diagonale:
   l_diag([1:M],1) = -fo*ones(M,1);
   l_diag(M-1,1)   = -2.*fo;
   l_diag(M,1)     = 0.;
% Sopra-diagonale:
   u_diag([1:M],1) = -fo*ones(M,1);
   u_diag(2,1)     = -2.*fo;
   u_diag(1,1)     = 0.;
echo off;
%
   disp('premere un tasto per continuare...'), disp(' '), pause;
%
clc;
echo on;
%
% Costruiamo ora la matrice sparsa K:
%
   K = spdiags([l_diag([1:M],1) m_diag([1:M],1) u_diag([1:M],1)], ...
       [-1 0 1],M,M)
%
% Proviamo a vederne la struttura:
%
   clf
   spy(K);
%
echo off;
   disp('premere un tasto per continuare...'), disp(' '), pause;
%
clc;
disp(' ')
disp(' Si tratta, infatti, di una matrice tridiagonale !')
disp(' ')
%
   disp('premere un tasto per continuare...'), disp(' '), pause;
%
clc;
disp(' ')
disp(' Iniziamo il calcolo - loop principale')
disp(' memorizzando le temperature')
disp(' in una serie di frames, in modo da poter')
disp(' poi visualizzare la dinamica del fenomeno')
echo off ;
   clf;
   elapsed = 0.;
   x = [0:dL:L]';
   MOV = moviein(N);
   grid
   xlabel('node');
   ylabel('Temperature');
   axis([0. L tair tin]);
   h_axis=gca;
   set(h_axis,'Xtick',[0 L/4 L/2 3*L/4 L],'Ytick',[tair tin/4 tin/2 tin*3/4 tin],...
       'Xtickmode','manual','Ytickmode','manual');
   set(h_axis,'XLimMode','manual','YLimMode','manual')
   hold on;
   for l = 1:N;       % Per tutti gli istanti
        elapsed   = elapsed + dtau;
        rhs([2:M-1],1)  = temp_old([2:M-1],1);
        rhs(1,1)        = temp_old(1,1) + tair*(2.*fo*bi_loc);
        rhs(M,1)        = temp_old(M,1) + tair*(2.*fo*bi_loc);
        temp = K\rhs; 
        plot(x,temp);
        MOV(:,l) = getframe;
        temp_old([1:M],1) = temp([1:M],1);
   end
   movie(MOV,5)
echo off;