[cig-commits] r19861 - seismo/3D/FAULT_SOURCE/branches/new_fault_db/src
surendra at geodynamics.org
surendra at geodynamics.org
Fri Mar 23 19:37:29 PDT 2012
Author: surendra
Date: 2012-03-23 19:37:29 -0700 (Fri, 23 Mar 2012)
New Revision: 19861
Modified:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
Log:
Started working on rate & state friction implementation
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-03-23 22:52:25 UTC (rev 19860)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90 2012-03-24 02:37:29 UTC (rev 19861)
@@ -25,6 +25,8 @@
! This module was written by:
! Percy Galvez, Jean-Paul Ampuero and Tarje Nissen-Meyer
+! Surendra Nadh Somala : Added Heterogenous initial stress capabilities (based on TPV16)
+! Surendra Nadh Somala : Added Rate and State Friction
module fault_solver
@@ -62,6 +64,12 @@
real(kind=CUSTOM_REAL), dimension(:), pointer :: Dc=>null(), mus=>null(), mud=>null(), theta=>null(), T=>null(), C=>null()
end type swf_type
+ type rs_type
+ private
+ logical :: healing = .false.
+ real(kind=CUSTOM_REAL), dimension(:), pointer :: V0=>null(), f0=>null(), L=>null(), theta_init=>null(), V_init=>null(), a=>null(), b=>null(), theta=>null(), T=>null(), C=>null()
+ end type rs_type
+
type bc_dynflt_type
private
integer :: nspec,nglob
@@ -72,6 +80,7 @@
real(kind=CUSTOM_REAL) :: dt
integer, dimension(:), pointer :: ibulk1, ibulk2
type(swf_type), pointer :: swf => null()
+ type(rs_type), pointer :: rs => null()
logical :: allow_opening = .false. ! default : do not allow opening
type(dataT_type) :: dataT
type(dataXZ_type) :: dataXZ
@@ -93,6 +102,8 @@
logical, save :: TPV16 = .false.
+ logical, save :: Rate_AND_State = .false.
+
real(kind=CUSTOM_REAL), allocatable, save :: Kelvin_Voigt_eta(:)
public :: BC_DYNFLT_init, BC_DYNFLT_set3d_all, Kelvin_Voigt_eta, &
@@ -432,10 +443,12 @@
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
-
NAMELIST / INIT_STRESS / S1,S2,S3,n1,n2,n3
NAMELIST / SWF / mus,mud,dc,nmus,nmud,ndc
+ NAMELIST / RS / V0,f0,a,b,L,V_init,theta_init,nV0,nf0,na,nb,nL,nV_init,ntheta_init
read(IIN_BIN) bc%nspec,bc%nglob
if (bc%nspec==0) return
@@ -523,6 +536,7 @@
! if (abs(bc%zcoord(k)-0.e0_CUSTOM_REAL) <= SMALLVAL) bc%T0(2,k) = 0
! end do
+if(.not. Rate_AND_State) then
! Set friction parameters and initialize friction variables
allocate( bc%swf )
allocate( bc%swf%mus(bc%nglob) )
@@ -551,6 +565,56 @@
allocate(bc%MU(bc%nglob))
bc%MU = swf_mu(bc%swf)
+else !Rate AND State (Ageing Law)
+! Set friction parameters and initialize friction variables
+ allocate( bc%rs )
+ allocate( bc%rs%V0(bc%nglob) )
+ allocate( bc%rs%f0(bc%nglob) )
+ allocate( bc%rs%a(bc%nglob) )
+ allocate( bc%rs%b(bc%nglob) )
+ allocate( bc%rs%L(bc%nglob) )
+ allocate( bc%rs%V_init(bc%nglob) )
+ allocate( bc%rs%theta_init(bc%nglob) )
+ ! WARNING: if V_HEALING is negative we turn off healing
+ bc%rs%healing = (V_HEALING > 0e0_CUSTOM_REAL)
+
+V0 =1.e-6_CUSTOM_REAL
+f0 =0.6_CUSTOM_REAL
+a =0.0080_CUSTOM_REAL
+b =0.0120_CUSTOM_REAL
+L =0.0135_CUSTOM_REAL
+V_init =1.e-12_CUSTOM_REAL
+theta_init =1.084207680000000e+09_CUSTOM_REAL
+
+nV0 =0
+nf0 =0
+na =0
+nb =0
+nL =0
+nV_init =0
+ntheta_init =0
+
+ read(IIN_PAR, nml=RS)
+ bc%rs%V0 = V0
+ bc%rs%f0 = f0
+ bc%rs%a = a
+ bc%rs%b = b
+ bc%rs%L = L
+ bc%rs%V_init = V_init
+ bc%rs%theta_init = theta_init
+ call init_2d_distribution(bc%rs%V0,bc%coord,IIN_PAR,nV0)
+ call init_2d_distribution(bc%rs%f0,bc%coord,IIN_PAR,nf0)
+ call init_2d_distribution(bc%rs%a,bc%coord,IIN_PAR,na)
+ call init_2d_distribution(bc%rs%b,bc%coord,IIN_PAR,nb)
+ call init_2d_distribution(bc%rs%L,bc%coord,IIN_PAR,nL)
+ call init_2d_distribution(bc%rs%V_init,bc%coord,IIN_PAR,nV_init)
+ call init_2d_distribution(bc%rs%theta_init,bc%coord,IIN_PAR,ntheta_init)
+
+ bc%rs%theta = 0e0_CUSTOM_REAL
+ allocate(bc%MU(bc%nglob))
+ bc%MU = rs_mu(bc%rs)
+endif
+
call init_dataT(bc%dataT,bc%coord,bc%nglob,NT,iflt)
call init_dataXZ(bc%dataXZ,bc,bc%nglob)
@@ -856,14 +920,22 @@
! Opening implies free stress
if (bc%allow_opening) T(3,:) = min(T(3,:),0.e0_CUSTOM_REAL)
- ! Update slip weakening friction:
+if(.not. Rate_AND_State) then ! Update slip weakening friction:
! Update slip state variable
! WARNING: during opening the friction state variable should not evolve
call swf_update_state(bc%D,dD,bc%V,bc%swf)
! Update friction coeficient
bc%MU = swf_mu(bc%swf)
+else ! Update rate and state friction:
+ ! Update slip state variable
+ ! WARNING: during opening the friction state variable should not evolve
+ call rs_update_state(bc%V,bc%dt,bc%rs)
+ ! Update friction coeficient
+ bc%MU = rs_mu(bc%rs,bc%V)
+endif
+
! combined with time-weakening for nucleation
! if (associated(bc%twf)) bc%MU = min( bc%MU, twf_mu(bc%twf,bc%coord,time) )
@@ -999,6 +1071,28 @@
end subroutine swf_update_state
!=====================================================================
+subroutine rs_update_state(V,dt,f)
+
+ real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: V
+ type(swf_type), intent(inout) :: f
+ double precision, intent(in) :: dt
+
+vDtL = sqrt(V(1,:)**2 + V(2,:)**2) * dt * L
+
+where(vDtL > 1.e-20_CUSTOM_REAL)
+ f%theta = f%theta * exp(vDtL) + L/V * (1 - exp(vDtL))
+elsewhere
+ f%theta = f%theta * exp(vDtL) + dt - 0.5*V/L*dt**2
+endwhere
+
+ if (f%healing) then
+ where(sqrt(V(1,:)**2 + V(2,:)**2) < V_HEALING)
+ f%theta = 0e0_CUSTOM_REAL
+ endwhere
+ endif
+end subroutine rs_update_state
+
+!=====================================================================
! Friction coefficient
function swf_mu_TPV16(f,time,nglob) result(mu)
@@ -1031,6 +1125,20 @@
end function swf_mu
+!=====================================================================
+! Friction coefficient
+function rs_mu(f,V) result(mu)
+
+ type(rs_type), intent(in) :: f
+ real(kind=CUSTOM_REAL) :: mu(size(f%theta))
+ real(kind=CUSTOM_REAL) :: V(size(f%theta))
+ real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: V
+
+ !-- rate and state friction (ageing law):
+ mu = f%f0 + (f%a * log(sqrt(V(1,:)**2 + V(2,:)**2))/f%V0)) + (f%b * log(f%V0*f%theta/f%L))
+
+end function rs_mu
+
!===============================================================
! OUTPUTS
subroutine init_dataT(DataT,coord,nglob,NT,iflt)
More information about the CIG-COMMITS
mailing list