%% Clear version of Pressure Diffusion in 2D - Italy 
%% Thanushika Gunatilake & Stephen A. Miller 

clear all
close all

%% Load Data
Italy_AVN_EQData=csvread('Italy_AVN_EQData.csv');
save('Italy_AVN_EQData.mat','Italy_AVN_EQData');
load('Italy_AVN_EQData.mat')
Lat = Italy_AVN_EQData(1:4369,2);
Long= Italy_AVN_EQData(1:4369,3);
M=Italy_AVN_EQData(1:4369,11);
depth=Italy_AVN_EQData(1:4369,4);
yr = datetime(Italy_AVN_EQData(1:4369,5),Italy_AVN_EQData(1:4369,6),Italy_AVN_EQData(1:4369,7),Italy_AVN_EQData(1:4369,8),Italy_AVN_EQData(1:4369,9),Italy_AVN_EQData(1:4369,10));
absoluteDays = days(yr-yr(1));
relativeDays = diff(yr);
relativeDays=[0;days(relativeDays)];
day = relativeDays;
data_days_check_plot=absoluteDays;

%% Domain parameters
L                       = 40*1000;      %   Length of domain  [m]
H                       = 15*1000;      %   Height of domain  [m]
nx                      = 300;          %   Number of gridpoints in x-direction
ny                      = 300;          %   Number of gridpoints in y-direction
dy                      = H/ny;         %   Spacing between gridpoints in y-direction
dx                      = L/nx;         %   Spacing between gridpoints in x-direction
[x2d,y2d]               = meshgrid(0:dx:(nx-1)*dx,(ny-1)*dy:-dy:0); %Computational domain
x2d                     = x2d'./1000;
y2d                     = y2d'./1000;
material                = ones(nx,ny);  %   Model geometry
beta_fluid              = 1e-10;        %   Fluid compressibility       [1/Pa]
beta_phi                = 1e-8;         %   Pore compressibility        [1/Pa]
porosity_rock           = 0.03;         %   Rock porosity
rho_fluid               = 1000;         %   Density of fluid            [kg/m3]
rho_rock                = 2700;         %   Density of rocks            [kg/m3]
g                       = 10;           %   Gravitational acceleration  [m2/s]
Geotherm_grad           = 20;           %   Geothermal gradient         [K/km]
Stress_state            = 2;            %   Stress state: 1-Extension, 2-Compression, 3-Transtension
num_faults              = 6;            %   Specify the number of discrete faults in the model (current maximum 5)
Dir_name                = '2D_output_1';%   Directory where the data should be saved
Day                     = 3600*24;      %   Number of seconds per day
dt                      = Day;          %   Timestep in seconds
tend                    = 351;           %   Overall sumulated time in days (%tend=1000,Save_timestep=2000;)
tot_time                = (tend*Day);   %   Total nunmber of timesteps
Num_timestep            = 0;            %   number of time steps simulated
Save_timestep           = 5;            % 0.01*2000=20Bilder, How many timesteps are there between 2 savings?  0.05*20 = 1Bild pro Tag  
Save_data               = true;         %   Save the data or not?
Save_tiff_pictures      = true;         %   Save pictures as tiff?
error_level             = 1e-2;
error_iterations        = realmax;
n_iter                  = 0;
%% Faults
Faults = repmat(struct('x0', 10000, 'y0', -4000, 'yend', -11000, 'thickness', 150, 'dip', 60), num_faults, 6);
Faults(1).x0            = 28*1000;      %   x-coordinate of faults upper limit
Faults(1).y0            = -2200;        %   y-coordinate of faults upper limit
Faults(1).yend          = -12500;        %   y-coordinate of faults lower limit (x-coordinate is calculated)
Faults(1).thickness     = 0.1*1000;     %   Thickness of the fault [m]
Faults(1).dip           = 135;          %   Dip of the fault. Measured anti-clock wise from the x-axis [deg]

Faults(2).x0            = 14.78*1000;       %   x-coordinate of faults upper limit
Faults(2).y0            = -4569;        %   y-coordinate of faults upper limit
Faults(2).yend          = -7300;        %   y-coordinate of faults lower limit (x-coordinate is calculated)
Faults(2).thickness     = 0.1*1000;     %   Thickness of the fault [m]
Faults(2).dip           = 60;           %   Dip of the fault. Measured anti-clock wise from the x-axis [deg]

Faults(3).x0            = 11.1*1000;       %   x-coordinate of faults upper limit
Faults(3).y0            = -6326;        %   y-coordinate of faults upper limit
Faults(3).yend          = -12600;       %   y-coordinate of faults lower limit (x-coordinate is calculated)
Faults(3).thickness     = 0.1*1000;     %   Thickness of the fault [m]
Faults(3).dip           = 45;           %   Dip of the fault. Measured anti-clock wise from the x-axis [deg]

Faults(4).x0            = 22.4*1000;   %   x-coordinate of faults upper limit
Faults(4).y0            = -6150;        %   y-coordinate of faults upper limit
Faults(4).yend          = -9850;        %   y-coordinate of faults lower limit (x-coordinate is calculated)
Faults(4).thickness     = 0.05*1000;     %   Thickness of the fault [m]
Faults(4).dip           = 109;          %   Dip of the fault. Measured anti-clock wise from the x-axis [deg]

Faults(5).x0            = 19*1000;   %   x-coordinate of faults upper limit
Faults(5).y0            = -6250;        %   y-coordinate of faults upper limit
Faults(5).yend          = -7250;        %   y-coordinate of faults lower limit (x-coordinate is calculated)
Faults(5).thickness     = 0.08*1000;     %   Thickness of the fault [m]
Faults(5).dip           = 80;           %   Dip of the fault. Measured anti-clock wise from the x-axis [deg]

Faults(6).x0            = 18.3*1000;     %   x-coordinate of faults upper limit
Faults(6).y0            = -8552;        %   y-coordinate of faults upper limit
Faults(6).yend          = -10350;        %   y-coordinate of faults lower limit (x-coordinate is calculated)
Faults(6).thickness     = 0.05*1000;     %   Thickness of the fault [m]
Faults(6).dip           = 75;           %   Dip of the fault. Measured anti-clock wise from the x-axis [deg]

%% Boundary Conditions
Bound_type_hor          =   {'Neumann'  ,'Neumann'  };
Bound_val_hor           =   [0   , 0];
Bound_type_ver          =   {'Neumann'  ,'Dirichlet'    };
Bound_val_ver           =   [0              0];

%% Material
Material_mat(1,:)       =  [   60        0      0.01    1e-14        35      ];     
Material_mat(2,:)       =  [   60        0      0.05    1e-12        100     ];     
Material_mat(3,:)       =  [   60        0      0.05    1e-12        100     ];     
Material_mat(4,:)       =  [   60        0      0.05    1e-12        100     ];     
Material_mat(5,:)       =  [   60        1e-6   0.05    1e-12        100     ];     
Material_mat(6,:)       =  [   60        1e-6   0.05    1e-12        100     ];     
Material_mat(7,:)       =  [   60        1e-6   0.05    1e-12        100     ];     

mat_matrix=1;
mat_fault2=2;
mat_fault3=3;
mat_fault4=4;
mat_fault5=5;
mat_fault6=6;
mat_fault7=7;

%% Create Domain
for i=1:nx
    for j=1:ny
        material(i,j)=mat_matrix;
    end
end

%% Set fault indentifier
for k=1:num_faults
  if Faults(k).y0 <= Faults(k).yend
   error(['The end of fault %i should be deeper than the start!',k])
  end
  x0 = Faults(k).x0;
  y0 = Faults(k).y0;
  xend = x0-abs(Faults(k).yend+Faults(k).y0) / tan(-Faults(k).dip*pi/180.0);
  yend = Faults(k).yend;
  for i=1:nx
   for j=1:ny  
       xpos = i*dx;
       ypos = -(ny-j)*dy;
       dist_to_fault = abs((yend-y0)*xpos-(xend-x0)*ypos+xend*y0-yend*x0)/sqrt((yend-y0)^2+(xend-x0)^2);
       if (dist_to_fault <= Faults(k).thickness && ypos >= Faults(k).yend && ypos <= Faults(k).y0)
      material(i,j)=1+k;
       end
   end
  end
end

P_fluid_hydro   =   rho_fluid*g*(y2d+.1)*1000;  % Hydrostatic fluid pressure
Sigma_1         =   rho_rock*g*(y2d+.1)*1000;   % Lithostatic rock pressure


if Stress_state==1                              
    Sigma_3 = Sigma_1;
    Sigma_1 = 2.*Sigma_3;                      
elseif Stress_state==2                         
    Sigma_3 = 0.62*Sigma_1;                    
    Sigma_1 = Sigma_1;  
elseif Stress_state==3                        
    Sigma_3 = Sigma_1;
    Sigma_1 = 1.5*Sigma_3;                                        
else
    error('Stress state unknown')
end

%% Variables
P_source = zeros(nx,ny);         % Fluid source pressure field (impermeable source is handled in initial conditions)
P_source2 = zeros(nx,ny); 
P_source3 = zeros(nx,ny); 
P_source4 = zeros(nx,ny); 
parameter
P_old = zeros(nx,ny);            % Solution field
kc0 = ones(nx,ny);               % Permeability matrix (m^2)
sig_star=zeros(nx,ny);           % Stress norm factor matrix
Temp = -y2d./1e3*Geotherm_grad + 0; % Temperature, assuming a constant gradient and 0 degrees on top
Fi = porosity_rock.*ones(nx,ny);  %Porosity matrix
Beta = Fi.*(beta_phi + beta_fluid);% beta coefficient used for the diffusion equation.
Visc = ones(nx,ny)*10^(-3)-0.9*(y2d./max(max(y2d))).*ones(nx,ny)*10^(-3); %visc in Pa-s (20 dec C=1e-3 Pa-s)
theal = zeros(nx,ny);             % Healing time. Basically stores the number of timesteps since healing became active.
trig_tog = zeros(nx,ny);          % Stores if the point has already been triggered.
plot_trig = zeros(nx,ny);         %Stores the points that were triggered in the current timestep
parameter;

%% Mechanic (Plot initial failure condition)
Sigma_normal = (Sigma_1 + Sigma_3 - 2.*(P_old + P_fluid_hydro))./2+((Sigma_1 - Sigma_3)./2).*(cos(2.*Dip_angle.*pi/180));
Tau = ((Sigma_1-Sigma_3)./2).*sin(2.*Dip_angle.*pi./180);
fric_fail=0.48; 

%% Simulation Begin
time_count=1;           
time = 0;               

k_mean = NaN(floor(tot_time/dt)+1,1);
p_mean = NaN(floor(tot_time/dt)+1,1);
Q = NaN(floor(tot_time/dt)+1,1);
total_trig = zeros(floor(tot_time/dt)+1,1);
time_trig = zeros(floor(tot_time/dt)+1,1);
P_source0 = P_source;


while time<tot_time
  count_trig=1;
  time =time+dt;
  n_iter= 0;

  if time_count==int16(66/dt*Day) || time_count==int16(136/dt*Day) || time_count==int16(198/dt*Day)
    theal = 0.*theal;
    theal2 = 0.*theal2;
  end
  

  sigma_n=(Sigma_1+Sigma_3-2.*(P_old + P_fluid_hydro))./2+((Sigma_1-Sigma_3)./2).*(cos(2.*Dip_angle.*pi/180));
  tau_a=((Sigma_1-Sigma_3)./2).*sin(2.*Dip_angle.*pi./180);

  for i=1:nx
    for j=1:ny
        if (material(i,j)==1 || material(i,j)==2 ||material(i,j)==3 ||material(i,j)==4 ...
                || material(i,j)==5 || material(i,j)==6 || material(i,j)==7)
            if ((tau_a(i,j)./sigma_n(i,j)>fric_fail||tau_a(i,j)./sigma_n(i,j)<-fric_fail)&&trig_tog(i,j)==0)
                trig_tog(i,j)=1;
                plot_trig(i,j)= plot_trig(i,j)+1;
                count_trig=count_trig+1;
                theal(i,j)=0;
            end
        
                  theal(i,j)=theal(i,j)+dt;
                  theal2(i,j)=theal(i,j)/24/3600;

        if (trig_tog(i,j)==1)

        if time_count<=int16(65/dt*Day)
            alpha = 0.03;
            P_source(i,j)=P_source0(i,j).*exp(-alpha.*theal2(i,j));
            kc0(i,j)=1e-14+kc0(i,j).*exp(-alpha.*theal2(i,j));
            P_source2(i,j) = P_source(i,j);
        elseif time_count>int16(65/dt*Day) && time_count<=int16(135/dt*Day)
            alpha2 = 0.0567;
             P_source(i,j)=P_source2(i,j).*exp(-alpha2.*theal2(i,j));
            kc0(i,j)=5e-14+kc0(i,j).*exp(-alpha2.*theal2(i,j));
            P_source3(i,j) = P_source(i,j);
        else
            alpha4 = 0.01;
            P_source(i,j)=P_source3(i,j).*exp(-alpha4.*theal2(i,j));
            kc0(i,j)=1e-14+kc0(i,j).*exp(-alpha4*theal2(i,j));
        end
        end
end
        
    end
  end
    K                   =   kc0.*exp(-Sigma_normal./1e6./sig_star);         
    k_mu                =   (1./Visc).*K;
    [P_new]             =   p2d(P_old, Beta, k_mu, P_source, nx, ny, dx, dy, dt, Bound_type_hor, Bound_val_hor, Bound_type_ver, Bound_val_ver);
    error_iterations    =   max(max(abs(P_old-P_new)))./max(max(abs(P_new)));
    dP_dt               =  (P_new-P_old)/(dt/Day);
    P_old               =   P_new;
    n_iter              =   n_iter + 1;
    Q                   =   (P_source./Beta);
    disp(['The error of iterations # ', num2str(n_iter),' = ',num2str(error_iterations)]);



  time_count=time_count+1;

end
            

