[cig-commits] r20831 - seismo/3D/FAULT_SOURCE/branches/new_fault_db/src
surendra at geodynamics.org
surendra at geodynamics.org
Fri Oct 12 13:56:58 PDT 2012
Author: surendra
Date: 2012-10-12 13:56:58 -0700 (Fri, 12 Oct 2012)
New Revision: 20831
Modified:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
Log:
Modified slip law of rate-and-state to use strong rate-weakening by default. Benchmarked with TPV103 & TPV104
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90 2012-10-12 16:23:53 UTC (rev 20830)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90 2012-10-12 20:56:58 UTC (rev 20831)
@@ -78,7 +78,8 @@
real(kind=CUSTOM_REAL), dimension(:), pointer :: V0=>null(), f0=>null(), L=>null(), &
V_init=>null(), &
a=>null(), b=>null(), theta=>null(), &
- T=>null(), C=>null()
+ T=>null(), C=>null(), &
+ fw=>null(), Vw=>null()
end type rsf_type
type asperity_type
@@ -245,8 +246,8 @@
real(kind=CUSTOM_REAL) :: mus,mud,dc
integer :: nmus,nmud,ndc,ij,k,e
real(kind=CUSTOM_REAL), dimension(:), allocatable :: nx,ny,nz
- real(kind=CUSTOM_REAL) :: V0,f0,a,b,L,theta,theta_init,V_init
- integer :: nV0,nf0,na,nb,nL,nV_init,ntheta_init
+ real(kind=CUSTOM_REAL) :: V0,f0,a,b,L,theta,theta_init,V_init,fw,Vw
+ integer :: nV0,nf0,na,nb,nL,nV_init,ntheta_init,nfw,nVw
real(kind=CUSTOM_REAL) :: AspTx,Fload
integer :: nAsp,nFload
real(kind=CUSTOM_REAL) :: C,T
@@ -280,7 +281,7 @@
NAMELIST / INIT_STRESS / S1,S2,S3,n1,n2,n3
NAMELIST / SWF / mus,mud,dc,nmus,nmud,ndc,C,T,nC,nForcedRup
- NAMELIST / RSF / V0,f0,a,b,L,V_init,theta_init,nV0,nf0,na,nb,nL,nV_init,ntheta_init,C,T,nC,nForcedRup
+ NAMELIST / RSF / V0,f0,a,b,L,V_init,theta_init,nV0,nf0,na,nb,nL,nV_init,ntheta_init,C,T,nC,nForcedRup,Vw,fw,nVw,nfw
NAMELIST / ASP / Fload,nFload
read(IIN_BIN) bc%nspec,bc%nglob
@@ -513,6 +514,8 @@
allocate( bc%rsf%theta(bc%nglob) )
allocate( bc%rsf%C(bc%nglob) )
allocate( bc%rsf%T(bc%nglob) )
+ allocate( bc%rsf%fw(bc%nglob) )
+ allocate( bc%rsf%Vw(bc%nglob) )
V0 =1.e-6_CUSTOM_REAL
f0 =0.6_CUSTOM_REAL
@@ -535,6 +538,11 @@
nC = 0
nForcedRup = 0
+ fw = 0.2_CUSTOM_REAL
+ Vw = 0.1_CUSTOM_REAL
+ nfw =0
+ nVw =0
+
read(IIN_PAR, nml=RSF)
bc%rsf%V0 = V0
bc%rsf%f0 = f0
@@ -545,6 +553,8 @@
bc%rsf%theta = theta_init
bc%rsf%C = C
bc%rsf%T = T
+ bc%rsf%fw = fw
+ bc%rsf%Vw = Vw
call init_2d_distribution(bc%rsf%V0,bc%coord,IIN_PAR,nV0)
call init_2d_distribution(bc%rsf%f0,bc%coord,IIN_PAR,nf0)
call init_2d_distribution(bc%rsf%a,bc%coord,IIN_PAR,na)
@@ -554,6 +564,8 @@
call init_2d_distribution(bc%rsf%theta,bc%coord,IIN_PAR,ntheta_init)
call init_2d_distribution(bc%rsf%C ,bc%coord,IIN_PAR,nC)
call init_2d_distribution(bc%rsf%T ,bc%coord,IIN_PAR,nForcedRup)
+ call init_2d_distribution(bc%rsf%fw ,bc%coord,IIN_PAR,nfw)
+ call init_2d_distribution(bc%rsf%Vw ,bc%coord,IIN_PAR,nVw)
!!$ ! WARNING : Not general enough
!!$ nglob_bulk = size(vel,2)
@@ -584,7 +596,7 @@
if (c1 .and. c2) then
b11 = w/(abs(x)-W1-w)
b12 = w/(abs(x)-W1)
- B1 = 0.5 * (ONE + tanh(b11 + b12))
+ B1 = HALF * (ONE + tanh(b11 + b12))
elseif(abs(x)<=W1) then
B1 = 1._CUSTOM_REAL
else
@@ -594,7 +606,7 @@
if (c3 .and. c4) then
b21 = w/(abs(z-hypo_z)-W2-w)
b22 = w/(abs(z-hypo_z)-W2)
- B2 = 0.5 * (ONE + tanh(b21 + b22))
+ B2 = HALF * (ONE + tanh(b21 + b22))
elseif(abs(z-hypo_z)<=W2) then
B2 = 1._CUSTOM_REAL
else
@@ -602,22 +614,29 @@
endif
bc%rsf%a(i) = 0.008 + 0.008 * (ONE - B1*B2)
+ bc%rsf%Vw(iglob) = 0.1 + 0.9 * (ONE - B1*B2)
elseif( abs(x)<=W1 .and. abs(z-hypo_z)<=W2 ) then
bc%rsf%a(i) = 0.008
+ bc%rsf%Vw(iglob) = 0.1_CUSTOM_REAL
else
bc%rsf%a(i) = 0.016
+ bc%rsf%Vw(iglob) = 1.0_CUSTOM_REAL
endif
enddo
! WARNING: The line below scratches an earlier initialization of theta through theta_init
! We should implement it as an option for the user
- bc%rsf%theta = bc%rsf%L/bc%rsf%V0 &
- * exp( ( bc%rsf%a * log(2.0*sinh(-sqrt(bc%T0(1,:)**2+bc%T0(2,:)**2)/bc%T0(3,:)/bc%rsf%a)) &
+ if(bc%rsf%stateLaw == 1) then
+ bc%rsf%theta = bc%rsf%L/bc%rsf%V0 &
+ * exp( ( bc%rsf%a * log(TWO*sinh(-sqrt(bc%T0(1,:)**2+bc%T0(2,:)**2)/bc%T0(3,:)/bc%rsf%a)) &
- bc%rsf%f0 - bc%rsf%a*log(bc%rsf%V_init/bc%rsf%V0) ) &
/ bc%rsf%b )
-
+ else
+ bc%rsf%theta = bc%rsf%a * log(TWO*V0/V_init * sinh(-sqrt(bc%T0(1,:)**2+bc%T0(2,:)**2)/bc%T0(3,:)/bc%rsf%a)
+ endif
+
allocate(bc%MU(bc%nglob))
! WARNING: the line below is only valid for pure strike-slip faulting
@@ -954,7 +973,7 @@
!-- outputs --
! write dataT every NTOUT time step or at the end of simulation
-! if ( mod(it,NTOUT) == 0 .or. it==NSTEP) call SCEC_write_dataT(bc%dataT,bc%dt,it)
+! if ( mod(it,NTOUT) == 0 .or. it==NSTEP) call SCEC_write_dataT(bc%dataT,bc%dt,it,bc%rsf%StateLaw)
if ( it == NSTEP) call SCEC_Write_RuptureTime(bc%dataXZ,bc%dt,NSTEP,iflt)
endif
@@ -1073,18 +1092,22 @@
end function rsf_mu
!---------------------------------------------------------------------
-subroutine funcd(x,fn,df,tStick,Seff,Z,f0,V0,a,b,L,theta)
+subroutine funcd(x,fn,df,tStick,Seff,Z,f0,V0,a,b,L,theta,statelaw)
real(kind=CUSTOM_REAL) :: tStick,Seff,Z,f0,V0,a,b,L,theta
double precision :: arg,fn,df,x
-
- arg = exp((f0+dble(b)*log(V0*theta/L))/a)/2._CUSTOM_REAL/V0
+
+ if(statelaw == 1) then
+ arg = exp((f0+dble(b)*log(V0*theta/L))/a)/TWO/V0
+ else
+ arg = exp(theta/a)/TWO/V0
+ endif
fn = tStick - Z*x - a*Seff*asinh(x*arg)
df = -Z - a*Seff/sqrt(ONE + (x*arg)**2)*arg
end subroutine funcd
!---------------------------------------------------------------------
-function rtsafe(funcd,x1,x2,xacc,tStick,Seff,Z,f0,V0,a,b,L,theta)
+function rtsafe(funcd,x1,x2,xacc,tStick,Seff,Z,f0,V0,a,b,L,theta,statelaw)
integer, parameter :: MAXIT=200
real(kind=CUSTOM_REAL) :: x1,x2,xacc
@@ -1093,6 +1116,7 @@
!real(kind=CUSTOM_REAL) :: df,dx,dxold,f,fh,fl,temp,xh,xl
double precision :: df,dx,dxold,f,fh,fl,temp,xh,xl,rtsafe
real(kind=CUSTOM_REAL) :: tStick,Seff,Z,f0,V0,a,b,L,theta
+ integer :: statelaw
call funcd(dble(x1),fl,df,tStick,Seff,Z,f0,V0,a,b,L,theta)
call funcd(dble(x2),fh,df,tStick,Seff,Z,f0,V0,a,b,L,theta)
@@ -1149,6 +1173,7 @@
real(kind=CUSTOM_REAL), intent(in) :: dt
real(kind=CUSTOM_REAL) :: vDtL(size(V))
+ real(kind=CUSTOM_REAL) :: f_ss(size(V)),xi_ss(size(V)),fLV(size(V))
vDtL = V * dt / f%L
@@ -1160,9 +1185,17 @@
f%theta = f%theta * exp(-vDtL) + dt*( ONE - HALF*vDtL )
endwhere
- ! slip law
+ ! slip law : by default use strong rate-weakening
else
- f%theta = f%L/V * (f%theta*V/f%L)**(exp(-vDtL))
+! f%theta = f%L/V * (f%theta*V/f%L)**(exp(-vDtL))
+ where(V /= 0._CUSTOM_REAL)
+ fLV = f%f0 - (f%b - f%a)*log(V/f%V0)
+ f_ss = f%fw + (fLV - f%fw)/(ONE + (V/f%Vw)**8)**0.125
+ xi_ss = f%a * log( TWO*f%V0/V * sinh(f_ss/f%a) )
+ f%theta = xi_ss + (f%theta - xi_ss) * exp(-vDtL)
+ elsewhere
+ f%theta = f%theta
+ endwhere
endif
end subroutine rsf_update_state
@@ -1314,26 +1347,28 @@
end subroutine store_dataT
!-----------------------------------------------------------------
-subroutine write_dataT_all(nt)
+subroutine write_dataT_all(nt,statelaw)
integer, intent(in) :: nt
+ integer, intent(in) :: statelaw
integer :: i
if (.not.allocated(faults)) return
do i = 1,size(faults)
- call SCEC_write_dataT(faults(i)%dataT,faults(i)%dt,nt)
+ call SCEC_write_dataT(faults(i)%dataT,faults(i)%dt,nt,statelaw)
enddo
end subroutine write_dataT_all
!------------------------------------------------------------------------
-subroutine SCEC_write_dataT(dataT,DT,NT)
+subroutine SCEC_write_dataT(dataT,DT,NT,statelaw)
type(dataT_type), intent(in) :: dataT
real(kind=CUSTOM_REAL), intent(in) :: DT
integer, intent(in) :: NT
-
+ integer, intent(in) :: statelaw
+
integer :: i,k,IOUT,ipoin
character :: NTchar*5
integer :: today(3), now(3)
@@ -1355,7 +1390,7 @@
endif
open(IOUT,file='OUTPUT_FILES/'//trim(dataT%name(i))//'.dat',status='replace')
- write(IOUT,*) "# problem=TPV102"
+ write(IOUT,*) "# problem=TPV104"
write(IOUT,*) "# author=Surendra Nadh Somala"
write(IOUT,1000) today(2), today(1), today(3), now
write(IOUT,*) "# code=SPECFEM3D_SESAME (split nodes)"
@@ -1383,6 +1418,7 @@
dataT%t3(k,i)/1.0e6_CUSTOM_REAL
enddo
else
+ if(statelaw == 1) then
write(IOUT,*) "t h-slip h-slip-rate h-shear-stress v-slip v-slip-rate v-shear-stress n-stress log-theta"
write(IOUT,*) "#"
do k=1,NT
@@ -1390,6 +1426,15 @@
-dataT%d2(k,i), -dataT%v2(k,i), -dataT%t2(k,i)/1.0e6_CUSTOM_REAL, &
dataT%t3(k,i)/1.0e6_CUSTOM_REAL, log10(dataT%theta(k,i))
enddo
+ else
+ write(IOUT,*) "t h-slip h-slip-rate h-shear-stress v-slip v-slip-rate v-shear-stress n-stress psi"
+ write(IOUT,*) "#"
+ do k=1,NT
+ write(IOUT,'(9(E15.7))') k*DT, dataT%d1(k,i), dataT%v1(k,i), dataT%t1(k,i)/1.0e6_CUSTOM_REAL, &
+ -dataT%d2(k,i), -dataT%v2(k,i), -dataT%t2(k,i)/1.0e6_CUSTOM_REAL, &
+ dataT%t3(k,i)/1.0e6_CUSTOM_REAL, dataT%theta(k,i)
+ enddo
+ endif
endif
close(IOUT)
@@ -1421,7 +1466,7 @@
IOUT = 121 !WARNING: not very robust. Could instead look for an available ID
open(IOUT,file=trim(filename),status='replace')
- write(IOUT,*) "# problem=TPV102"
+ write(IOUT,*) "# problem=TPV104"
write(IOUT,*) "# author=Surendra Nadh Somala"
write(IOUT,1000) today(2), today(1), today(3), now
write(IOUT,*) "# code=SPECFEM3D_SESAME (split nodes)"
More information about the CIG-COMMITS
mailing list