%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright Ville Vuorinen (2019)
% No guarantee of code functionality
% FILE: heat2dfins.m 
% PURPOSE: a code that aims to solve heat transport
% in 2d with a provided (fixed) velocity distribution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
set(0,'DefaultAxesFontName','Times New Roman')
set(0,'DefaultAxesFontSize',15)


clear *
close all

% two data sets available: single fin (VeldataSingle.mat) 
% and multiple (Veldata.mat) fin

Single = 2; % 1=single fin; 2 = many fins

if(Single == 1)
    load VeldataSingle.mat
    MyFile = 'FinsSingle.mat'; % give name of the file to save data in
elseif(Single == 2) 
    load Veldata.mat
    MyFile = 'FinsMany.mat'; % give name of the file to save data in
end





% The chosen values for T are here made to match with classroom exp
Tleft = 273 + 20; Thot = 273+5;

% Fin width Lx = 5cm and channel gap/height Ly = 3mm
Nx = 1024; Ny = 128; Lx = S.Lx; Ly = S.Ly;

% reference tables for partial derivative calculations
inx = 2:(Nx-1); iny = 1:(Ny); 
north = iny+1; north(Ny)=1;
east = inx+1;  
south = iny-1; south(1)=Ny;
west = inx-1;

% grid spacing dx and dy - divide by Nx-2 etc due to 2 ghost cell layers
dx = Lx/(Nx-2); dy = Ly/(Ny-2); 

% pre-step to determine the physical location of cell centroids
x=linspace(0,Lx,Nx-2); 
y=linspace(0,Ly,Ny-2); 

% create the "grid" of (x_i, y_j) pairs
% not needed for solution but may be useful in visualization/analysis
[X,Y] = meshgrid(x,y); 

% initialize temperature
T = Tleft*ones(Ny,Nx);
CFLmax = 0.1;
% The values below are for air at 1atm 300K taken from 
% Incropera Appendix
Pr = 0.707;        % Prandtl number of air
nu = (1.589e-5);   % nu = kinematic viscosity of air [nu] = m^2/s
alpha = nu/Pr;     % alpha = thermal conductivity of air [alpha] = [nu]
rho = 1.1614; cp = 1.007*1000; 
k   = alpha*rho*cp;    % thermal conductivity of air [k] = W/(mK)

% Set U and V constant based on a pre-cursor CFD simulation
U = S.U(:,inx); V=S.V(:,inx);

% timestep, make sure Co and CFL are small enough (please see slides)
dt = 0.5e-5;

dxmin = min(dx,dy);

giveInfo; 
% into giveInfo.m file you can add any info to be printed on the screen

setTBCs;
for(t=1:8000)
setTBCs;
solveTemperature; 
visualizeResults;
end

save(MyFile, 'T', 'U', 'X', 'Y', 'Lx', 'Ly', 'alpha', 'k', 'rho', 'cp', 'iny','inx','dx','dy');