% Author: Layla S. Mayboudi
% 2021-05-27
% This code creates a transient thermal PDE model and plots the geometry 
% and heat flux flow diagrams.
% Selected nodes and their temperatures are also displayed.

% Note:
% A transient model uses ThemalConductivity, MassDensity and SpecificHeat. 


clear all; close all; clf; clc;

set(groot,'DefaultAxesColorOrder',[0 0 0],...
    'DefaultSurfaceEdgeColor',[0 0 0],...
    'DefaultLineLineWidth',1,...
    'DefaultAxesLineStyleOrder','-|--|:|-.|+|o|*|.|x|s|d|^|v|>|<|p|h');      %may be added to set the line style and color.
ax.FontSize = 20; ax.FontName = 'Cambria'; ax.FontSmoothing = 'on'; ax.FontWeight = 'normal'; %'bold'; 

  %'DefaultAxesLineStyleOrder','-|--|:|-.|+|o|*|.|x|s|d|^|v|>|<|p|h'      %may be added to set the line style and color.
  %set(gcf,'defaultSurfaceEdgeColor','b')

thermalModelT = createpde('thermal','transient'); % create a model object: thermal
geom1 = [3;4;0;60;60;0;0;0;26;26]; % create geometry: rectangle solid
geom2 = [3;4;0;60;60;0;26;26;30;30]; % create geometry: rectangle solid

gd = [geom1 geom2];
sf = 'geom1+geom2';
ns = char('geom1','geom2');
ns = ns';
dl = decsg(gd,sf,ns);
geometryFromEdges(thermalModelT,dl); % associate geometry and thermal model object
figure(1);
pdegplot(thermalModelT,'EdgeLabels','on','FaceLabels','on');
xlim([0 60]); ylim([0 30]); axis equal;
grid on; title 'Geometry with Edge and Face Labels Displayed';

%Water
k1 = 0.6305; % thermal conductivity (W/mK)
rho1 = 988; %density (kg/m3)
cp1 = 4181; % heat capacity (J/kgK)
q1 = 0; %heat generation (W/m^3)
kFunc1 = @(region,state) k1*region.y;
cFunc1 = @(region,state) cp1*region.y;
qFunc1 = @(region,state) q1*region.y;
rhoFunc1 = @(region,state) rho1*region.y;
thermalProperties(thermalModelT,'ThermalConductivity',k1,...
                               'MassDensity',rho1,...
                               'SpecificHeat',cp1,'Face',1);
internalHeatSource(thermalModelT,q1,'Face',1);

%Aluminum
k2 = 210; % thermal conductivity (W/mK)
rho2 = 2698.9; %density (kg/m3)
cp2 = 900; % heat capacity (J/kgK)
q2 = 1000; %heat generation (W/m^3)
kFunc2 = @(region,state) k2*region.y;
cFunc2 = @(region,state) cp2*region.y;
qFunc2 = @(region,state) q2*region.y;
rhoFunc2 = @(region,state) rho2*region.y;
thermalProperties(thermalModelT,'ThermalConductivity',k2,...
                               'MassDensity',rho2,...
                               'SpecificHeat',cp2,'Face',2);
internalHeatSource(thermalModelT,q2,'Face',2);

% Boundary Conditions
T_ambient = 25; %degC
T_water = 50; %degC
outerCC_V = @(region,~) 10*region.y; %horizontal heat convection coefficient (W/m2K)
top_BC_HF = @(region,~) 100*region.y; %heat source distribution (W/m2)
thermalBC(thermalModelT,'Edge',2,'HeatFlux',top_BC_HF); %boundary condition: heat source (W/m2)
thermalBC(thermalModelT,'Edge',1,'Temperature',T_water); %boundary condition: temperature
thermalBC(thermalModelT,'Edge',[4 5 6 7],'ConvectionCoefficient',outerCC_V,'AmbientTemperature',T_ambient); %boundary condition: convection heat transfer

% Initial Condition
T_init = 50; % set initial temperature (degC)
thermalIC(thermalModelT,T_init); % associate initial condition with thermal model object

figure(2);
grid_size=0.5; %identify grid size
mesh1 = generateMesh(thermalModelT,'Hmax',grid_size); %generate mesh *****
pdeplot(mesh1); %plot PDE thermal model object with node ID's displayed.
%axis equal; xlim([0 15]); ylim([0 7.5]); % set axis conditions
grid on; title 'Triangular Elements';
% or pdeplot(mesh.Nodes,mesh.Elements); %plot PDE thermal model object with node ID's displayed.

figure(3);
pdeplot(mesh1,'NodeLabels','on','ElementLabels','on'); %plot PDE thermal model object with node ID's displayed.
axis equal; xlim([0 0.5]); ylim([0 0.5]); % set axis conditions
grid on; title 'Node IDs on';

%tlist = logspace(-2,-1,10);
%thermalresults = solve(thermalModelT,tlist)
%T = thermalresults.Temperature;
tfinal = 340*60; % set final time (s)
tlist = 0.1:100:tfinal; % set time array
% thermalModelT.SolverOptions.PDESolverOptions Properties (or the following)
% AbsoluteTolerance, RelativeTolerance, ResidualTolerance, MaxIterations, MinStep, ResidualNorm, ReportStatistics
thermalModelT.SolverOptions.ReportStatistics = 'on';
results = solve(thermalModelT,tlist); % identify PDESolverOptions Properties
%pde.TransientThermalResults.Temperature (or the following)
% Temperature, XGradients, YGradients, and ZGradients 
T = results.Temperature; % identify data to process and plot *****
Tx = results.XGradients; % identify data gradient along the x-coordinate *****
Ty = results.YGradients; % identify data gradient along the y-coordinate *****
T_st = results.SolutionTimes; % identify data gradient along the y-coordinate *****
 
figure(4);
pdeplot(thermalModelT,'XYData',T(:,end),'ColorMap','hot','Contour','on','Levels',10); %,'ColorMap','hot' ,'Contour','on'
axis equal; xlim([0 30]); ylim([0 30]);
grid on; title(sprintf('Temperature Contours at the End of Heating (%g min)',tfinal/60));
% pdeplot(thermalModelT,'XYdata',T(:,end),'Contour','on','LineWidth',2);

figure(5);
pdeplot(thermalModelT,'FlowData',[Tx,Ty]);
axis equal; xlim([0 1]); ylim([0 30]);
grid on; title(sprintf('Temperature Contours at the End of Heating (%g min)',tfinal/60));

h = pdeplot(thermalModelT,'XYData',T,'ZData',T, ...
                   'FaceAlpha',0.5, ...
                   'FlowData',[Tx,Ty], ...
                   'ColorMap','jet', ...
                   'Mesh','on')
 
pp = thermalModelT.Mesh.Nodes; % recall node locations from thermal model object *****
nodesLeftEnd = find(pp(1,:) < eps); % identify far left boundary nodes
nodeCenter = nodesLeftEnd(pp(2,nodesLeftEnd) < eps+5); % identify central nodes
nodeOuter = nodesLeftEnd(pp(2,nodesLeftEnd) < 4-eps); % identify upper boundary nodes
 
figure(6);
plot(tlist/60,T(nodeCenter,:),'MarkerSize',2); %plot central nodes
hold all % multiple diagrams are plotted in one figure
grid on; title 'Transient Temperature at the Center'; legend ('Center Axis'); 
xlabel 'time (min)'; ylabel 'Temperature (degC)';
 
figure(7);
plot(tlist/60,Tx(nodeCenter,:),'MarkerSize',2); %plot central nodes
hold all % multiple diagrams are plotted in one figure
grid on; title 'Transient Temperature Gradient along the x-Coordinate at the Center'; legend ('Center Axis'); 
xlabel 'time (min)'; ylabel 'Temperature (degC)';
 
figure(8);
plot(tlist/60,Ty(nodeCenter,:),'MarkerSize',2); %plot central nodes
hold all % multiple diagrams are plotted in one figure
grid on; title 'Transient Temperature Gradient along the y-Coordinate at the Center'; legend ('Center Axis'); 
xlabel 'time (min)'; ylabel 'Temperature (degC)';
 
%tlist2 = 0:1000:tfinal; % perform a new solution for shorter time.
%resultT2 = solve(thermalModelT,tlist2)

figure(9)
v = linspace(0,30,15);
[X,Y] = meshgrid(v);
[qx,qy] = evaluateHeatFlux(results,X,Y,1:length(tlist));
tlist(end)
quiver(X(:),Y(:),qx(:,end),qy(:,end),'MarkerSize',2);
axis equal; xlim([0 30]); ylim([25 50]);
grid on; title 'Heat Flux at the End of Heating Time'); 
xlabel 'X'; ylabel 'Heat flux (W/m^2)';

figure(10)
Qh = evaluateHeatRate(results,'Edge',1);
plot(tlist,Qh,'MarkerSize',2);
xlabel('Time')
ylabel('Heat Flow Rate across edge 1')
ax.FontSize = 20; ax.FontName = 'Cambria'; ax.FontWeight = 'bold'; ax.FontSmoothing = 'on'
grid on; title 'Heat Flow at the End of Heating Time'); 
xlabel 'time (min)'; ylabel 'Heat flow rate (W/m^2)';

figure(11)
v = linspace(0,30,15);
[X1,Y1] = meshgrid(v);
[Tx,Ty] = evaluateTemperatureGradient(results,X1,Y1,1:length(tlist));
tlist(end);
axis equal; xlim([0 30]); ylim([25 50]);
grid on; title 'Heat Flux for the 20th Solution Step'); 
xlabel 'time (min)'; ylabel 'Heat flux (W/m^2)';
plot(tlist,Qh,'MarkerSize',2);
%quiver(X,Y,gradTx,gradTy)
grid on; title 'Transient Heat Rate'; 

X = 0:0.1:1;
Y = ones(size(X));
[gradTx,gradTy] = evaluateTemperatureGradient(results,X,Y,1:length(tlist));
[qx,qy] = evaluateTemperatureGradient(results,X,Y,1:length(tlist));

%set(groot,'defaultAxesColorOrder',[1 0 0;1 0 1;0 1 0;0 0 1],...
      %'defaultAxesLineStyleOrder','+-|o-|v-|s-|p-|h-|d-|^-') %this line
      %may be added to set the line style and color.
%%Define a line at convection boundary and compute temperature gradients across that line.
%%Plot the interpolated gradient component gradTx along the x axis for the following values from the time interval tlist.  
X = 0:1:60;
Y = ones(size(X));
[gradTx,gradTy] = evaluateTemperatureGradient(results,X,Y,1:length(tlist));

figure(12)
legend; xlabel('x'); ylabel('Temperature gradient (degC/m)');
t = [51:500:201];
%t = [1:100:204];
for i = t
legend; xlabel('x'); ylabel('Temperature gradient (degC/m)');
    p(i) = plot(X,gradTx(:,i),'DisplayName', strcat('t=', num2str(tlist(i)))); %,'MarkerSize',2);
  hold on
end
legend; xlabel('x'); ylabel('Temperature gradient (degC/m)');
grid on; title 'Spatial Temperature Gradient'; 

figure(13)
v = linspace(0,30,15);
[X1,Y1] = meshgrid(v);
[Tx,Ty] = evaluateTemperatureGradient(results,X1,Y1,1:length(tlist));
size(Tx); size(Ty);
X=[100:1:303]/2000;
Y=[0:1:224]/10000;
grid on; title 'Temperature Gradient for the Selected Grid'); 
xlabel 'x (m)'; ylabel 'Temperature gradient (degC/m)';
quiver(X,Y,Tx,Ty,'MarkerSize',2); 
grid on; title 'Temperature Gradient Flow'; 