% MF_ADAPT -- Pipe Cooling using Seawat for Viscosity and heat transport
%  TO 100405

clear variables
close all

%You can swicht between different cases by changing iCase [1,2,3,4]
iCase=5;

basename='Diffusion';

% This exampple shows pure diffusion without even groundwater flow. It
% shows a steady-state solution of stagnant groundwater in a half-infinite
% "rod" with a mass pulse at x=0.
% It also shows as case 2, the same rod but a fixed concentration at x=0
% It then continuous with a heat pulse at x=0 and, finally uses
% a fixed temperature at x=0.
% All thsese cases are easily verified with analytical solutions.
% Diffusion and conduction are computed using seawat (also possible
% directly in MT3DMS) without flow. They could alternatively be computed
% using the flow model so that head stands for concentration.
%
% Cases 2,3,4 are cases with heat flow, they are all situations without
% groundwater flow and are compared with analytical solutions
% case 2: heat flow from a constant temperature at x=0;
% case 3: heat dissipaation from an iniital temperature 40C over the
% width of the aquifer W (only the half rod is modelled. The analytical
% solution is the sum of two erfc funcions.
% case 4: heat dissipation from an initial pulse in the first cell
% such that the total heat in this cell equals that of the previous case
% defined by T0 over the width of the aquifer. Analytical solution is a
% pulse

% there are two stress periods, a very short one 0.1 day and a long one,
% 1000 years (365000 days). The two stress periods are not needed for the
% four cases of this example. They were meant to be used in conjunction
% with mass loading boundary condtion and the analytical solution of a
% pulse.

%% Parameters
W     =25;      % [m] thickness of aquifer
L     =2000;    % [m]  half width of model
peff  = 0.35;   % [ ]  effective porosity
rhos  =2650;    % [kg/m3] density of solids
rhow  =1000;    % [kg/m3] density of water
rhob  =(1-peff)*rhos; % [kg/m3] bulk density of porous medium
k     = 0.1;    % [m/d] hydraulic conductivity
ds    =1.5e-10;  % [m2/s] Diffusion coefficient
Ds    =ds*86400; % [m2/d] Diffusion coefficient
%Kd    =1/1700;   % [m3/kg]distribution coefficient for linear sorption
%R     =1+rhob*Kd/peff;  % [-] retardation
%% Specifics for heat
Lw    =0.6*86400;     % [J/d/m/K] heat conductance of water
Ls    =3.0*86400;     % [J/d/m/K] heat conductance of solides
cw    =4200;          % [J/kg/K]  heat capacity of the water
cs    = 800;          % [J/kg/K]  heat capacity of the solids
Lambda=peff*Lw+(1-peff)*Ls ;            % [J/d/m/K] bulk heat conductance
rhoc  =peff*rhow*cw +(1-peff)*rhob*cs;  % [J/m3/K]  bulk heat capacity
DH    =Lambda/(peff*rhow*cw);     % [m2/d]  diffusivity for heat
Kd=cs/(rhow*cw);        % [m3/kg] distribution coefficient for heat
DTEMP=40;               % initial temperature

%% Different cases with the same model

switch iCase
    case 1, L=2000/100;  % mass transport fixed conc at x=0
            D=Ds;
            C0=DTEMP;
            C0max=C0;
            R=1+rhob*Kd/peff;
            iPer=2;
            iType=-1;
    case 2, L=2000;      % heat transport fixed temp at x=0
            D=DH;
            C0=DTEMP;
            C0max=C0;
            R=1+rhob*Kd/peff;
            iPer=2;    % boundary in stress period 2
            iType=-1;  % constant conc
    case 3, L=2000;      % heat transport, initial temp in aquifer, die out of initial temp
            D=DH;
            R=1+rhob*Kd/peff;
            C0=DTEMP;  % to suddenly inject a mass M use C0=M/dt
            C0max=C0;
    case 4, L=2000;      % heat transport, heat pulse at x=0, die out of initial temp
            D=DH;
            R=1+rhob*Kd/peff;
            C0=DTEMP;
            C0max=C0;
    case 5, L=2000;      % heat transport, heat pulse at x=0, die out after puse (mass loading)
            D=DH;
            R=1+rhob*Kd/peff;
            iType=15;    % mass loading source
            iPer=1;    % boundary in stress period 1
            dt=0.1; % [d] check that this is the first stress period length
            C0=peff*R*DTEMP*(W/2)/dt;
            C0max=DTEMP;
end

%% Mesh
xGr=sinespace(0,L,2000,pi/20,pi/2);
zGr=[-1 1]/2;
yGr=[-1 1]/2;

[xGr,yGr,zGr,xm,ym,zm,Dx,Dy,Dz,NX,NY,NZ]=modelsize3(xGr,yGr,zGr);

Z=zGr(ones(NY,1),ones(1,NX),1:NZ+1);

[XM,YM,ZM]=meshgrid(xm,ym,zm);
[DX,DY,DZ]=meshgrid(Dx,Dy,Dz);

%% Generate all other matrices
IBOUND=ones(NY,NX,NZ);
ICBUND=ones(NY,NX,NZ);
STRTHD=ones(NY,NX,NZ) * 0;
STCONC=ones(NY,NX,NZ)*0;
HK    =ones(NY,NX,NZ)  *k;
VK    =HK;
PEFF  =ones(NY,NX,NZ) *peff;
DMCOEF=ones(size(IBOUND))*D;

%% Adaptation to boundaries

switch iCase
    case {1,2}, I=find(XM<xGr(2));  % fixed conc boundary condition
                
                %These fixed concentration or fixed temperature boundary condition
                %is only used in cases 1 (conc) and 2 (temp)
                %is sets the pointsrc in stress period 2 to the desired
                %concentration or temperature. Stress period 1 is only used
                %ot set an arbitrary initial concentration or temperature.
                %Therefore, this stress period is very short.
                
                % Set point source boundary for stress period 2 (ignore
                % period 1, whic is very short
                LRC =cellIndices(I ,size(IBOUND),'LRC');
                u =ones(size(LRC( :,1)));
                PNTSRC= [iPer*u  LRC  u.*C0   u*iType u.*C0];
                
    case 3,     STCONC(:,XM<W/2,:)=C0;          % no boundary only initial temp
    case 4,     STCONC(:,1,:) =C0*(W/2)/Dx(1);  % no boundary only initial temp
    case 5,     I=find(XM<xGr(2));              % pulse at t=0, by mass loading
                LRC =cellIndices(I ,size(IBOUND),'LRC');
                u =ones(size(LRC( :,1)));
                PNTSRC= [iPer*u  LRC  u.*C0   u*iType u.*C0];

end

%% Save variable that we need in mf_analyze.m
save underneath.mat D C0 R iCase C0max peff W

