% [ShortOutput,LongOutput] = fulldynV(reft)%% PURPOSE:	2D Spectral BIEM for fault dynamics%		Antiplane mode, VELOCITY formulation%		References: Morrisey and Geubelle 1997%			    Lapusta et al. 2000%%		+ second order time evolution scheme (see Lapusta et al. 2000)%		+ second order semi-open integration rule for stress functional%% INPUTS: 	reft	[integer] grid refinement ratio.%			The default (1) uses the grid as defined in the script%% SETTINGS:	Some parameters need to be set inside this script,%		see section "SIMULATION PARAMETERS"%% OUTPUTS: 	see section "STORE RESULTS"%		ShortOutput	at every time step, for selected fault nodes%		LongOutput	at every fault node, for selected times %% USES: 	C3.m, W3.m, friction.m%% Jean-Paul Ampuero	ampuero@erdw.ethz.ch 	% Last modified: August 2005% % TO DO: adaptive timestep%	 optimize the storage scheme of time-truncated fields%        eliminate replication%function [ShortOutput,LongOutput] = fulldynV(reft)%-----------------------------------------------% UNITSMPa = 1e6;year = 356*24*60*60 ;%-----------------------------------------------% SIMULATION PARAMETERS%-- physical properties --L = 20e3 ; 	% fault length (simulated period)CS = 3464.; 	% S wave velocityRHO = 2670.;	% densityMU = RHO*CS^2; 	% shear modulus%-- background stress --TAUINIT1 = 70*MPa ;	% initial shear stress TAUINIT2 = 120 *MPa ;	% normal stress (+ compressive)TAURATE1 = 0; %5e6*MPa/year; % tectonic shear loading rate%-- heterogeneity of initial shear stress --HET.L = 3e3 ;	% length scale [3e3]HET.A = 11.6*MPa;	% amplitude [11.6]%-- numerical settings --if ~exist('reft','var'), reft = 1; end	% refinement ratioN = 128*reft; 	% number of space elements, grid size [128]TMAX = 130*reft; 	% number of time steps [130]%N = 32*reft; TMAX = 32*reft;%-- friction --FRIC.Dc = 0.4;	% friction critical slip DcFRIC.MUs = 0.677;	% static friction coefficientFRIC.MUd = 0.525;	% dynamic friction coefficient%-- output settings --TOUT = 1;	% time stride for large outputs (snapshots)XOUT = (N/2:reft:N);	% nodes for large outputs%XOUT = (N/2:N);	% nodes for large outputsXSIS = [N/2 3*N/8];	% nodes for full time outputsECHO = 2;	% output info%-- expert numerical settings --RTCUT = 2 ;	% kernel_time-cut / fault_travel_time for mode_1, typically = 2QW = 4 ; 	% kernel_time-cut Nyquist_mode / mode_1 /(N/2), typically = 4                % N/2 gives a non-truncated kernelCFL = 0.5 ;	% stability number (Courant), typically = 0.5%-----------------------------------------------% INITIALIZEdx = L/N ;x = (-N/2+1:N/2).' *dx ;dt = CFL*dx/CS ; LONG_OUTPUT= nargout>1;SHORT_OUTPUT= nargout>0;if LONG_OUTPUT  LongNT = floor(TMAX/TOUT);  LongNX = length(XOUT);  LongOutput.Slip = zeros(LongNX,LongNT);  LongOutput.SlipRate = zeros(LongNX,LongNT);  LongOutput.Stress = zeros(LongNX,LongNT);  LongOutput.Strength = zeros(LongNX,LongNT);  LongOutput.X = x(XOUT);  LongOutput.Time = (1:LongNT)*TOUT*dt;endif SHORT_OUTPUT  ShortNX = length(XSIS);  ShortOutput.SlipRate = zeros(TMAX,ShortNX);  ShortOutput.Slip = zeros(TMAX,ShortNX);  ShortOutput.Stress = zeros(TMAX,ShortNX);  ShortOutput.MeanSlip = zeros(TMAX,1);  ShortOutput.MaxSlip = zeros(TMAX,1);  ShortOutput.MeanSlipRate = zeros(TMAX,1);  ShortOutput.MaxSlipRate = zeros(TMAX,1);  ShortOutput.CrackLength = zeros(TMAX,1);  ShortOutput.MeanStress = zeros(TMAX,1);  ShortOutput.Time = (1:TMAX)*dt;  ShortOutput.X = x(XSIS);endimpedance = MU/(2*CS) ;dload1 = TAURATE1*dt ;%load1 = TAUINIT1 + HET.A*exp(-x.^2/HET.L^2);load1 = TAUINIT1 + HET.A*(abs(x)<=HET.L/2);stress2 = repmat(TAUINIT2,N,1);% barrier:NWEAK = N/2 ; % size of the central weak zoneRBAR = 1 ; %10 ; % relative strength of the barrierstress2(1:N/2-NWEAK/2) = RBAR*TAUINIT2 ;stress2(N/2+NWEAK/2+1:end) = RBAR*TAUINIT2 ;% frictionFRIC.DMU = FRIC.MUs-FRIC.MUd ;FRIC.W = FRIC.DMU ./ FRIC.Dc ;W = TAUINIT2*FRIC.W ; 	% slip weakening rateSm = W / impedance ; 	% growth rateLc = 2*0.57888694*MU/W ;	% nucleation sizeV0 = TAUINIT2*FRIC.DMU/impedance ;	% typical earthquake slip ratedisp(sprintf( 'Nucleation length Lc = %0.4g',Lc));disp(sprintf( 'Space resolution Lc/dx = %0.4g',Lc/dx));disp(sprintf( 'Nucleation time scale 1/Sm = %0.4g',1/Sm) );disp(sprintf( 'Time resolution 1/Sm*dt = %0.4g',1/(Sm*dt)) );%--% Initialize kernels%% Convolution integrated with semi-open second order rule% (Numerical Recipes, combined eqs. 4.1.7 and 4.1.11)% Use of a semi-open rule simplifies the second pass of the % time evolution scheme.%% Mode 0 has null contribution, don't need to store it,% only strictly positive wavenumbers.% For negative wavenumbers we will use symetries of real FFT%% Also implemented: mode-dependent KERNEL CUT-OFF (TCUT) as in % Lapusta et al (2000)NK=N/2; % number of strictly positive wavenumbersNKnyq = N/2+1 ; % Nyquist wavenumber position in classical FFT storagek = (1:NK)';TCUT1 = floor(RTCUT*N/CFL);TCUT = ceil(  TCUT1*( 1+(QW-1)*(k-1)/(NK-1) ) ./k ); dk = 2*pi/L ;k = k*dk;pathstr = fileparts(mfilename('fullpath'));kernel_dir = 'kernels';kernel_name = fullfile(pathstr,kernel_dir);if exist(kernel_name,'file') ~= 7, mkdir(pathstr,kernel_dir); endkernel_name = fullfile(kernel_name, sprintf('kernelV_%u',N) );Compute_Kernel = 1;if exist([kernel_name '.mat'],'file')  if ECHO > 1, disp('Reading the kernel'), end  load(kernel_name); % --> kernel, KER_NX, KER_RTCUT, KER_QW, KER_CFL  Compute_Kernel = any([KER_NX,KER_RTCUT,KER_QW,KER_CFL] ~= [N,RTCUT,QW,CFL]);endif Compute_Kernel  if ECHO > 1, disp('Computing the kernel'), end  for ik=1:NK,    arg = k(ik)*CS* (1:TCUT(ik))'*dt; % = k*CS*t    factor = 0.5*MU*k(ik) *dt ;    kernel(ik).f = factor* W3(arg) ;   % factors for semi-open rule (Numerical Recipes eq. 4.1.20)   % it=1, t=0  : ff1 = 0   % it=2, t=dt : ff1 = K_1*V_1   % it=3       : ff1 = 3/2*K_1*V_2 + 1/2*K_2*V_1   % it=4       : ff1 = 3/2*K_1*V_3 + K_2*V_2 + 1/2*K_3*V_1   % it=n       : ff1 = 3/2*K_1*V_n + K_2*V_n-1 + ... + 1/2*K_n*V_1    kernel(ik).f(1) = 3/2*kernel(ik).f(1) ;  end  clear arg  KER_NX = N;  KER_RTCUT = RTCUT;  KER_QW = QW;  KER_CFL = CFL;  save(kernel_name,'kernel','KER_NX','KER_RTCUT','KER_QW','KER_CFL')endstkernel = - 0.5*MU*k ; % Static kernel%--% initialize fields :u = zeros(N,1);v = zeros(N,1);for ik=1:NK,  fv(ik).f = zeros(TCUT(ik),1);end  ff1 = zeros(N,1);ff2 = zeros(N,1);% countersito = 1;techo = ceil(TMAX/10);ite = 1;itcyc = ones(NK,1);%-----------------------------------------------% BEGIN TIME LOOP% ABOUT THE TIME EVOLUTION SCHEME:% Equations are of the form%   impedance*v = -Ffriction(u) + Fstatic(u) + Fdynamic(v) % where, owing to the use of a semi-open integration rule,% Fdynamic(v) depends stricly on previous values of v.% A model equation is:%   v = s*u% Discrete version:%   v(n+1) = s*u(n+1)%   u(n+1) = u(n) + 0.5*dt*[ v(n)+v(n+1) ]% Implemented as a two-pass explicit scheme (SECOND ORDER):%   1. u_1(n+1) = u(n) + dt*v(n)%   2. v_1(n+1) = s*u_1(n+1)%   3. u_2(n+1) = u(n) + 0.5*dt*[ v(n)+v_1(n+1) ]%               = u_1(n+1) +0.5*dt*[ v_1(n+1)-v(n) ]%   4. v_2(n+1) = s*u_2(n+1)if ECHO > 1, disp('Begin time loop'), endfor it = 1:TMAX ,   if it==1  % initial conditions, t=0  % to guarantee second order when initial stress drop is abrupt,   % the initial velocity must be set here, analytically  % impedance*v(0) = max( initial_stress - initial_strength, 0)  stress1 = load1;  strength = stress2.*friction(u,FRIC);  v = max( stress1-strength, 0 )/impedance; else  u_old = u;  v_old = v; %Update external load    load1 = load1 + dload1 ;  for ipass=1:2,   %Update slip    u = u_old + 0.5*dt*(v+v_old);       %Update static stress term    fftu = fft(u,N);     ff2(1) = 0 ;    ff2(2:NK+1) = stkernel.*fftu(2:NK+1) ;    ff2(NKnyq+1:N) = conj( ff2(NK:-1:2) );       %Trial state: assuming no further slip   % = load(n+1) +Fstatic(u) +Fdynamic(n+1)    stress1 = real(ifft(ff1+ff2,N)) ;    stress1 = stress1 + load1 ;     %Update strength Ffriction(u)    frimu = friction(u,FRIC);    strength = stress2.*frimu ; % normal stress >0 is compressive       %Solve friction -> update slip rate    v = max( stress1-strength ,0)/impedance ;    stress1 = min(stress1,strength) ;  end   end%--- Begin STORE RESULTS --- % (x,t) fields, coarsely sampled  if LONG_OUTPUT & it == ito*TOUT    LongOutput.Slip(:,ito) = u(XOUT) ;    LongOutput.SlipRate(:,ito) = v(XOUT) ;    LongOutput.Stress(:,ito) = stress1(XOUT) ;    LongOutput.Strength(:,ito) = strength(XOUT) ;    ito = ito+1;  end   % full time sampled outputs  if SHORT_OUTPUT   % fields(t) at selected locations    ShortOutput.SlipRate(it,:) = v(XSIS)';    ShortOutput.Slip(it,:) = u(XSIS)';    ShortOutput.Stress(it,:) = stress1(XSIS)';   % macroscopic fields    ShortOutput.MeanSlip(it) = mean(u) ;    ShortOutput.MaxSlip(it) = max(u) ;    ShortOutput.MeanSlipRate(it) = mean(v) ;    ShortOutput.MaxSlipRate(it) = max(v) ;    incrack = find(v>0);    ShortOutput.CrackLength(it) = length(incrack)*dx;     ShortOutput.MeanStress(it) = mean(stress1) ;%    ShortOutput.smcrack(it) = sum( stress1(incrack) )*dx;  end%--- End STORE RESULTS ---% Some progress info    if ECHO > 1 & ite == techo    ite = 1;    disp(sprintf('Step %i/%i, Vmax = %0.4g, Umax = %0.4g' ...                ,it,TMAX,max(v),max(u)))  else    ite = ite+1;  end    if it==TMAX, break, end  fftv = fft(v,N);% this loop is the most time-consuming section of the code  for ik=1:NK ,    % Update slip rate spectrum for NEXT time step   % "fv" stores the slip rate spectrum history: fv(MODE).f(TIME)   % The storage is CYCLIC in the time dimension:   % don't need to keep slip history longer than TCUT(mode).   % "itcyc" points, on entry, to the position of the current time step    %         in "fv(mode).f(:)"   = mod(it,TCUT(mode))+1    itc = itcyc(ik);    fv(ik).f(itc) = fftv(ik+1);   % itclast => earliest velocity stored    tcut = TCUT(ik);    itclast = itc+1;    if itclast==tcut, itclast=1; end   % Update stress functional for NEXT time step    if it==1     % fix weights of semi-open quadrature rule      ff1(ik+1) = 2/3*kernel(ik).f(1)* fv(ik).f(1);      fv(ik).f(1) = 0.5 *fv(ik).f(1);    elseif it <= tcut      ff1(ik+1) = sum( kernel(ik).f(1:it) .* fv(ik).f(it:-1:1) );    else     % "ittab" points to the slip history from it+1 back to it-TCUT+2      fv(ik).f(itclast) = 0.5*fv(ik).f(itclast); % fix quadrature weight      ittab = [ (itc:-1:1) (tcut:-1:itc+1) ];      ff1(ik+1) = sum( kernel(ik).f .* fv(ik).f(ittab) );    end   % Update "itcyc": point to next time step    itcyc(ik) = itclast;  end  ff1(1) = 0 ; % mode 0 is null  ff1(NKnyq+1:N) = conj( ff1(NK:-1:2) ); % real signal: complete the spectrumend %----- END TIME LOOP