Example input file

The following example is for the simulation of recurring earthquakes along a single one-dimensional fault embedded in a two-dimensional full elastic space within the in-plane strain condition. The Dieterich-Ruina-Rice number Ru = 13.13 and the Rb number is 0.28.

#!/bin/bash -e

# single asperity on a single fault, Ru=13.13, Rb=0.28

selfdir=$(dirname $0)

WDIR=$selfdir/output1

if [ ! -e $WDIR ]; then
	echo adding directory $WDIR
	mkdir $WDIR
fi

# number of fault patches
N2=1024
# patch size
DX=20

# use --verbose=2 to output all parameters
# use --verbose=1 to output only parameters different from previous patch
OMP_NUM_THREADS=2 motorcycle-ps-ratestate-serial \
	--verbose 1 \
	--epsilon 1e-6 \
	--export-state \
	--export-netcdf \
	--export-netcdf-rate 20 \
	--export-netcdf-step 1 \
	--maximum-step 3.15e7 \
	--maximum-iterations 1000000 \
	--friction-law 1 <<EOF
# output directory
$WDIR
# Lame parameters (lambda, mu)
30e3 30e3
# time interval
1e10
# number of faults
1
# grid dimension (N2)
$N2
# sampling (dx2)
$DX
#   n  tau0   mu0   sig   a   b   L   Vo   G/(2Vs)   Vl Dirichlet
$(echo "" | awk -v n2=$N2 -v dx=$DX '
	function abs(x){return (x>0)?x:-x};
	function max(x,y){return (x>y)?x:y};
	function boxcar(x){return (x>=-0.5 && x<=0.5)?1:0};
	function heavi(x){return (x>0)?1:0};
	function ramp(x){return x*boxcar(x-0.5)+heavi(x-1)};
	function asinh(x){return log(x+sqrt(1+x^2))};
	function sinh(x){return (exp(x)-exp(-x))/2};
	BEGIN{
	c=1;
	tau0_p=-1;mu0_p =-1;sig_p=-1;
	a_p=-1;b_p=-1;L_p=-1;
	Vo_p=-1;damping_p=-1;
	Vl_p=-1;dirichlet_p="T";
	}{
	for (i2=0;i2<n2;i2++){
		x2=(i2-n2/2)*dx; 
		# normal of initial shear traction (MPa); use -1 to use steady-state value
		tau0=-1;
		# characteristic weakening distance (m)
		L=1e-2;
		# rate-dependent parameter (unitless)
		a=1e-2; 
		# state-dependent parameter (unitless)
		b=(abs(x2)<2.5e3)?a+4.0e-3:a-4e-3;
		# reference coefficient of friction (unitless)
		mu0=0.6; 
		# effective normal stress (MPa)
		sig=1e2; 
		# reference slip-rate (m/s)
		Vo=1e-6; 
		# loading rate (m/s)
		Vl=1e-9;
		# radiation damping coefficient ( G/(2 Vs) in units of MPa/m*s )
		damping=5;
		# whether to apply Dirichlet boundary condition
		dirichlet=((i2<10) || (i2>=(n2-10)))?"T":"F";
		if (1==boxcar((x2-0.5e3)/1e3)){
			# initial stress triggers nucleation (MPa)
			tau0=mu0*sig*(1.1e-9/Vo)^(a/mu0)*(L/Vl)^(b/mu0);
		}
		# check if all parameters are identical to that of the previous patch
		if ((tau0_p == tau0) && (mu0_p==mu0) && (sig_p==sig) && (a_p==a) && (b_p==b) && (L_p==L) && 
		    (Vo_p==Vo) && (damping_p==damping) && (Vl_p==Vl) && (dirichlet_p==dirichlet)){
			# print minus line number to save space (value of previous patch is used in code)
			printf "%5d\n",-c;
		} else {
			# print new set of parameters
			printf "%5d %10.2e %10.2e %10.2e %10.2e %10.2e %10.2e %10.2e %d %10.2e %s\n", 
	        			c,  tau0,   mu0,   sig,     a,     b,     L,    Vo, damping, Vl, dirichlet;
		}
		c++;
		# save current parameters for reuse
		tau0_p=tau0;mu0_p=mu0;sig_p=sig;
		a_p=a;b_p=b;L_p=L;
		Vo_p=Vo;damping_p=damping;
		Vl_p=Vl;dirichlet_p=dirichlet;
	}
}')
# number of observation patches
1
# n fault i2 rate
  1     1 $(echo "" | awk -v n2=$N2 '{print n2/2+1,1}')
# number of events
0
EOF