%
% 2D SGSN testing code
% Assumes fault is symmetry plane
%
% Slip Weakening Version
%
%  S.M. Day, 02 Jan 2009
%
%    To Get Rough Run time estimate for 
%      ? MacBook Pro 4 GB, 
%      ? 2.53 GHz Intel Core 2 Duo
%      ? OS 10.5.6
%      ? Matlab 7.7.0.471 (R2008b)
%
%     time in sec = ~Nx*Nz*Nt*Nplot/700,000
%.........................................................................
%%Problem Inputs..........................................................
%
%
%    computational inputs
Nx=400;Nz=200;dx=100;        % x runs from 0 to (Nx-1)*dx; likewise z
CFL=0.5;                     % Courant number for max P speed
Nt=100;Nplot=10;             % Nt*Nplot   = total number of time steps
Nsmooth=50;                  % Nsmooth*dt = nucleation smoothing time
gamma=0.1;                   % numerical damping parameter
%    volume inputs
Vp=6000;Vs=3464;rho=2700;    % P speed, S speed, density
%    fault plane inputs
T0=0.5e8;                    % initial shear stress
SigmaN=1e8;                  % initial normal stress
coefs=0.7;                   % static friction coefficient
coefd=0.4;                   % dynamic friction coefficient
Dc=0.5;                      % slip weakening displacement
Xhypo=dx*Nx/2;               % coordinate of nucleation point
Vrup=1000;                   % nucleation-zone minimum rupture velocity
Rcrit=2000;                  % nucleation zone maximum half-length
%
%.........................................................................
%%Set constants...........................................................
%
%
%     General constants
dt=CFL*dx/Vp;
Tsmooth=Nsmooth*dt;
lam=rho*(Vp^2-2*Vs^2);
lam2mu=rho*Vp^2;
mu=rho*Vs^2;
mod=4*mu*(lam+mu)/lam2mu;
%     Fault plane constants
M=0.5*rho*dx^3;
%     Grid coordinates
xgu=(0:Nx-1)*dx;
zgu=(0:Nz-1)*dx;
%
%%Initialize Main Arrays..................................................
%
%
%     volume arrays
u=zeros(Nx,Nz); w=zeros(Nx,Nz);
sxx=zeros(Nx,Nz); szz=zeros(Nx,Nz); 
sxz=zeros(Nx,Nz);
rxx=zeros(Nx,Nz); rzz=zeros(Nx,Nz); 
rxz=zeros(Nx,Nz);
dudx=zeros(Nx-1,Nz); dudz=zeros(Nx,Nz-1);
dwdx=zeros(Nx-1,Nz); dwdz=zeros(Nx,Nz-1);
%     fault-plane arrays
slip=zeros(Nx,1);
V=zeros(Nx,1);
tauc=zeros(Nx,1);
Ttrial=zeros(Nx,1);
T=zeros(Nx,1);
R=zeros(Nx,1);
%
%
%%Time Step (explicit central differences)................................
%
%
for k=1:Nplot
    tic;
for j=1:Nt
    time=((k-1)*Nt+(j-1))*dt;
%
%    calculate derivatives at interior points
%
dudx(1:Nx-1,1:Nz)=(1/dx)*(u(2:Nx,:)-u(1:Nx-1,:));
dwdx(1:Nx-1,1:Nz)=(1/dx)*(w(2:Nx,:)-w(1:Nx-1,:));
dudz(1:Nx,1:Nz-1)=(1/dx)*(u(:,2:Nz)-u(:,1:Nz-1));
dwdz(1:Nx,1:Nz-1)=(1/dx)*(w(:,2:Nz)-w(:,1:Nz-1));
%
%    calculate stresses at interior points
%
rxx(1:Nx-1,2:Nz)=lam2mu*dudx(1:Nx-1,2:Nz)+lam*dwdz(1:Nx-1,1:Nz-1);
rzz(1:Nx-1,2:Nz)=lam*dudx(1:Nx-1,2:Nz)+lam2mu*dwdz(1:Nx-1,1:Nz-1);
rxz(2:Nx,1:Nz-1)=mu*dudz(2:Nx,1:Nz-1)+mu*dwdx(1:Nx-1,1:Nz-1);
sxx(1:Nx-1,2:Nz)=sxx(1:Nx-1,2:Nz)+dt*rxx(1:Nx-1,2:Nz);
szz(1:Nx-1,2:Nz)=szz(1:Nx-1,2:Nz)+dt*rzz(1:Nx-1,2:Nz);
sxz(2:Nx,1:Nz-1)=sxz(2:Nx,1:Nz-1)+dt*rxz(2:Nx,1:Nz-1);
%
%    calculate fault-plane sxx (leave szz zero, by symmetry)
%
PulsePeriod=20*dt;
szz(1:Nx-1,1)=zeros(Nx-1,1);
rzz(1:Nx-1,1)=zeros(Nx-1,1);
rxx(1:Nx-1,1)=mod*dudx(1:Nx-1,1)+szz(1:Nx-1,1)*lam/lam2mu;
sxx(1:Nx-1,1)=sxx(1:Nx-1,1)+dt*rxx(1:Nx-1,1);

%
%    update interior velocities
%

u(2:Nx,2:Nz)=u(2:Nx,2:Nz)+...
    (1/rho)*dt*(1/dx)*(sxx(2:Nx,2:Nz)-sxx(1:Nx-1,2:Nz)+...
    sxz(2:Nx,2:Nz)-sxz(2:Nx,1:Nz-1)+...
    gamma*dt*(rxx(2:Nx,2:Nz)-rxx(1:Nx-1,2:Nz)+...
    rxz(2:Nx,2:Nz)-rxz(2:Nx,1:Nz-1)));
w(1:Nx-1,1:Nz-1)=w(1:Nx-1,1:Nz-1)+...
    (1/rho)*dt*(1/dx)*(sxz(2:Nx,1:Nz-1)-sxz(1:Nx-1,1:Nz-1)+...
    szz(1:Nx-1,2:Nz)-szz(1:Nx-1,1:Nz-1)+...
    gamma*dt*(rxz(2:Nx,1:Nz-1)-rxz(1:Nx-1,1:Nz-1)+...
    rzz(1:Nx-1,2:Nz)-rzz(1:Nx-1,1:Nz-1)));
%
%    calculate trial shear traction
%
R(1,1)=0; 
R(2:Nx,1)=dx^2*( 0.5*(  sxx(2:Nx,1)-sxx(1:Nx-1,1)  )+sxz(2:Nx,1)  );
R(2:Nx,1)=dx^2*(   0.5*(  sxx(2:Nx,1)-sxx(1:Nx-1,1)  )+sxz(2:Nx,1)+...
    dt*gamma*(  0.5*(  rxx(2:Nx,1)-rxx(1:Nx-1,1)  )+rxz(2:Nx,1)));
Ttrial(1:Nx,1)=(1/dx^2)*((1/dt)*M*u(1:Nx,1)+R(1:Nx,1))+T0;
%
%    calculate fault strength
%
tauc=(SigmaN*(coefs-(coefs-coefd)*min([abs(slip),Dc*ones(size(slip))]')./Dc))';
%
%    artificial nucleation
%
indxNuc=find(abs(xgu-Xhypo)<=min([Rcrit,time*Vrup])');
tauc(indxNuc)=min([abs(Ttrial(indxNuc))-(abs(Ttrial(indxNuc))-SigmaN*coefd)...
    .*min([  (abs(time-abs(xgu(indxNuc)-Xhypo)./Vrup)./Tsmooth)' ,...
    ones(size(Ttrial(indxNuc))) ]')',tauc(indxNuc)]')';
%
%    set fault shear traction T
%
T=Ttrial.*(min([tauc./abs(Ttrial),ones(size(tauc))]'))';
%
%    update fault-plane slip velocity
% 
u(1:Nx,1)=u(1:Nx,1)+(dt/M)*(R-dx^2*(T-T0));
slip=slip+2*dt*u(1:Nx,1);
V=abs(u(1:Nx,1));
end
duration=toc;
PlotNumber=k;
[num2str(PlotNumber),  ' ..... ', num2str(Nt), ...
    ' steps in ', num2str(duration), ' sec']
%
%    make plots
%
figure(k)
subplot(4,1,1)
imagesc(xgu,-zgu,u');
axis xy
colorbar
xlabel('Fault-parallel Distance (m)')
ylabel ('Fault-normal Distance (m)')
title ('Fault-parallel velocity (m/s)')
subplot(4,1,2)
imagesc(xgu,-(zgu+dx/2),w');
axis xy
colorbar
xlabel('Fault-parallel Distance (m)')
ylabel ('Fault-normal Distance (m)')
title ('Fault-normal velocity (m/s)')
subplot(4,1,3)
plot(xgu,2*u(:,1))
xlabel ('Fault-parallel Distance (m)')
ylabel ('Slip rate (m/s)')
subplot(4,1,4)
plot(xgu,T,'r',xgu,tauc,'b')
legend ( 'stress' , 'strength')
xlabel ('Fault-parallel Distance (m)')
ylabel ('Shear Stress (Pa)')
end