%% DYKE COOLING WITH LATENT HEAT
% Calculates the cooling time for a dyke of uniform initial temperature.
% Uses MATLAB's ODE solver to solve the heat equation.
% Accounts for latent heat and finds position of solidification front.

%% DEFINE MODEL PARAMETERS:

Hd=0.5;     % Dyke half-width (m)
H=Hd*10;    % Location of boundary condition, total width of model (m)

TI=1200+273.15;     % Initial temperature (K)
T0=0+273.15;        % Boundary temperature (K)
Tc=1000+273.15;     % Solidification temperature (K)

cpS=1200;   % Specific heat capacity for solid basalt (J kg-1 K-1)
cpL=1200;   % Specific heat capacity for liquid basalt (J kg-1 K-1)

kS=1.3;     % Thermal conductivity for solid basalt (W m-1 K-1)
kL=1.3;     % Thermal conductivity for liquid basalt (W m-1 K-1)

rho=2700;   % Density of basalt (kg m-3) (assume constant)
L=400000;   % Latent heat of fusion (J kg-1)

tMAX=3*Hd^2 / (10^-6); % Simulation duration (s)

%% SET UP MODEL

% Define spatial grid:
N=2000;                         % Number of cells
n=500;                          % Number of time intervals evaluated
xB=linspace(0,H,N+1)';          % Cell boundaries
x=(xB(1:N,1)+xB(2:N+1,1))/2;    % Mid-points of cells

% Define initial temperature vector:
TIvec=zeros(size(x))+T0;        % Set all cells to T0.
TIvec(x<Hd)=TI;                 % Set all dyke cells to TI.

% Determine initial enthalpy vector:
hIvec=cpL*TIvec+(cpS-cpL)*Tc+L;
hIvec(TIvec<Tc)=cpS*TIvec(TIvec<Tc);

% Plot initial temperature profile:
figure
plot(x,TIvec-273.15);
title("Initial temperature profile");
xline(Hd);
xlabel('Distance (m)');
ylabel('Temperature (deg C)');

%% SOLVE PARTIAL DIFFERENTIAL EQUATION

% Define and set the Jacobian pattern:
JPat=spdiags(ones(N,3),[-1 0 1],N,N);
options=odeset('JPattern',JPat);

% Apply ode15s to obtain the finite difference solution:
sol=ode15s(@myODE,[0 tMAX],hIvec,options,x,xB,cpS,cpL,L,Tc,kS,kL,rho,T0,H);

%% EVALUATE SOLUTION AT TIME INTERVALS n

% Time intervals from 0 to max (s):
t = linspace(0,tMAX,n)';
% Evaluate solution at these intervals:
h = deval(sol,t);
for j=1:numel(t)
    [~,T1(:,j),xc(j)]=myODE(t(j),h(:,j),x,xB,cpS,cpL,L,Tc,kS,kL,rho,T0,H);
end

%% CENTRE TEMPERATURE:
% Find the time taken for the dyke centre to solidify.
% Track the temperature of the dyke centre.

% Isolate the thermal evolution of the dyke centre:
Tcntr = T1(1,:);

% Find the time when centre temp is first < Tc.
t_end = t(find(Tcntr<Tc,1,'first'));

disp("Solidification time (s)");
disp(t_end);

figure
plot(t./60^2,Tcntr-273.15);
xlabel('Time (hr)');
ylabel('Centre temp (deg C)');
xline(t_end);
xlim([0 max(t./60^2)]);
title("Centre temp. over time");
% t is in seconds: when plotting, can convert to minutes, hours, days, etc.

%% SOLIDIFICATION FRONT
% The boundary between T=Tc and T<Tc is the solidification front.

% Find location of solidification front at each timestep:
solPos = zeros(length(t),1);
for i=1:length(solPos)
    % Get the thermal profile for each timestep:
    Ti = T1(:,i);
    % Find first position where T < Tc:
    solPos(i) = x(find(Ti<Tc,1,'first'));
end

% Plot solidification front over time:
figure
plot(t/60^2,solPos);
xlabel('Time (hr)');
ylabel('Distance from dyke centre (m)');
xlim([0 t_end/60^2])
ylim([0 Hd]);
title("Position of solidification front");

%% FUNCTIONS
function [dhdt,T,xc] = myODE(~,h,x,xB,cpS,cpL,L,Tc,kS,kL,rho,T0,H)
% Determine temperature:
T=(h-(cpS-cpL)*Tc-L)/cpL;
T(h<(cpS*Tc+L))=Tc;
T(h<cpS*Tc)=h(h<cpS*Tc)/cpS;

% Interpolate temperature to block boundaries:
TB=interp1(x,T,xB);
% Ignore first node for zero-flux boundary:
TB(1)=[];
% Set last node to boundary temperature:
TB(end)=T0;

% Determine thermal conductivity:
k=zeros(size(TB))+kL;
k(TB<Tc)=kS;

% Calculate derivatives and include T=T0 at x=H boundary condition:
dTdx=diff([T;T0],1,1)./diff([x;H],1,1);
% Calculate heat fluxes and include J=0 at x=0 boundary condition:
J=[0;-k.*dTdx];
% Calculate the flux divergence:
dJdx=diff(J,1,1)./diff(xB,1,1);
% Solve the energy conservation statement:
dhdt=-dJdx/rho;
% Determine location of phase change:
xc=max(x(T>Tc));

if isempty(xc)
    % If no phase change is found:
    xc=NaN;
end

end