%% COOLING MULTIPLE PULSES WITH LATENT HEAT
% Same as original Cooling_latentHeat, but repeated for multiple pulses.
% Array is shunted aside to allow injection of new magma.
% Uses MATLAB's ODE solver to solve the heat equation.
% Accounts for latent heat and finds position of solidification front.

%% DEFINE MODEL PARAMETERS:

% Marginal band widths/pulse half-widths (m):
Hd=[0.01051,0.00442,0.00754,0.01170,0.01028,0.03431];
% Location of boundary condition, total width of model (m):
H=sum(Hd)*10;

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)

Tlim = 700+273.15;  % Temperature reached before next injection.
tMod = 100;         % Controls duration of model.

%% SET UP MODEL
% These values will remain constant on each iteration.

% Output: [band no., time to solidify, time to reach Tlim]
coolTimes = zeros(length(Hd),3);

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

% Determine initial tMAX:
tMAX = tMod*Hd(1)^2/10^-6;
dt = tMAX/n;

% Define initial temperature vector:
TIvec=zeros(size(x))+T0;        % Set all cells to T0.
TIvec(x<Hd(1))=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);

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

%% RUN MODEL FOR EACH MAGMA PULSE
% Determined by length of array of band widths, Hd.

solPos = zeros(1,2);            % Solidus position and time.
coolt = zeros(length(Hd),1);    % Cooling times for each pulse.
t_tot = 0;

for i=1:length(Hd)
    % SOLVE PARTIAL DIFFERENTIAL EQUATION
    % Using MATLAB's PDE solver.

    % 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);

    % Time intervals:
    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.

    % Isolate the thermal evolution of the dyke centre (K):
    Tcntr = T1(1,:);
    % Get the time when it solidified:
    t_s = find(Tcntr<Tc,1,'first');
    % Time when centre temp < Tlim:
    til = find(Tcntr<Tlim,1,'first');
    
    if isempty(til)
        disp("Centre temp. did not reach Tlim within tMAX.");
        disp("Try using a higher value of tMAX.");
    end
    
    coolt(i) = t_tot + t(til);
    % Display cooling times:
    disp([i, coolt(i)]);
    % Record data:
    coolTimes(i,:) = [i,t_tot+t(t_s),t_tot+t(til)];

    % POSITION OF SOLIDIFICATION FRONT:
    % Should move from x=Hd to x=0 in time t_s.
    % Array T1 contains temperature profiles through time.
    
    % For each timestep, find location of solidification front:
    hPos = zeros(til,1);
    for j=1:til
        % Get the thermal profile for this time interval:
        Tj = T1(:,j);
        % Find solidification front:
        hPos(j) = x(find(Tj<Tc,1,'first'));
    end

    % Give solidus position relative to outermost margin.
    m = sum(Hd(1:i));   % Margin position in m.
    hPosm = m-hPos;     % Sol. pos. relative to margin.
    % Only record positions up to t=t_s:
    solPos=[solPos;([hPosm,t_tot+t(1:til)])];
    
    % Update t_tot:
    t_tot = coolt(i);

    % SET PARAMETERS FOR NEXT ITERATION:
    if (i < length(Hd))
        % The temp. profile needed for the next iteration:
        Tnxt = T1(:,til);
        
        % The model width is set so that H=sum(Hd)*10.
        % Even if all pulses were emplaced as one, would not heat boundary.
        % Number of cells is constant between iterations.
        % Unheated wall cells are shunted off beyond the boundary.

        % Calculate the number of cells in the next pulse:
        dNxt = round(Hd(i+1)/dx);
        % Create new temperature vector using Tnxt:
        TIvec=zeros(size(x))+TI;        % Set all cells to TI.
        for j=dNxt+1:length(TIvec)
            TIvec(j)=Tnxt(j-dNxt);
        end
        % Ensure boundary remains at T0:
        TIvec(end)=T0;

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

        % Determine new tMAX:
        tMAX = tMod*Hd(i+1)^2/10^-6;
        dt = tMAX/n;

    end
end

%% FIGURES

% Solidification front over time:
figure
plot(solPos(:,2)./60^2,solPos(:,1));
xlabel('Time (hr)')
ylabel('Solidus position (m)');
ylim([0 sum(Hd)]);
hold on
for i=1:length(coolt)
    xline(coolt(i)./60^2)
end
hold off

%% FUNCTIONS

function [dhdt,T,xc] = myODE(t,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