[cig-commits] r19875 - seismo/3D/FAULT_SOURCE/branches/new_fault_db/src
ampuero at geodynamics.org
ampuero at geodynamics.org
Mon Mar 26 11:26:58 PDT 2012
Author: ampuero
Date: 2012-03-26 11:26:58 -0700 (Mon, 26 Mar 2012)
New Revision: 19875
Modified:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
Log:
fixed bugs in vDtL, Vfavg, tnorm
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-26 17:24:01 UTC (rev 19874)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90 2012-03-26 18:26:58 UTC (rev 19875)
@@ -343,61 +343,62 @@
! end do
- if(TPV16) then
+ if (TPV16) then
- allocate(inp_nx(bc%nglob))
- allocate(inp_nz(bc%nglob))
- allocate(loc_str(bc%nglob))
- allocate(loc_dip(bc%nglob))
- allocate(sigma0(bc%nglob))
- allocate(tau0_str(bc%nglob))
- allocate(tau0_dip(bc%nglob))
- allocate(Rstress_str(bc%nglob))
- allocate(Rstress_dip(bc%nglob))
- allocate(static_fc(bc%nglob))
- allocate(dyn_fc(bc%nglob))
- allocate(swcd(bc%nglob))
- allocate(cohes(bc%nglob))
- allocate(tim_forcedRup(bc%nglob))
+ allocate(inp_nx(bc%nglob))
+ allocate(inp_nz(bc%nglob))
+ allocate(loc_str(bc%nglob))
+ allocate(loc_dip(bc%nglob))
+ allocate(sigma0(bc%nglob))
+ allocate(tau0_str(bc%nglob))
+ allocate(tau0_dip(bc%nglob))
+ allocate(Rstress_str(bc%nglob))
+ allocate(Rstress_dip(bc%nglob))
+ allocate(static_fc(bc%nglob))
+ allocate(dyn_fc(bc%nglob))
+ allocate(swcd(bc%nglob))
+ allocate(cohes(bc%nglob))
+ allocate(tim_forcedRup(bc%nglob))
+
+ open(unit=IIN_NUC,file='DATA/FAULT/input_file.txt',status='old',iostat=ier)
+ read(IIN_NUC,*) relz_num,sub_relz_num
+ read(IIN_NUC,*) num_cell_str,num_cell_dip,siz_str,siz_dip
+ read(IIN_NUC,*) hypo_cell_str,hypo_cell_dip,hypo_loc_str,hypo_loc_dip,rad_T_str,rad_T_dip
+ do ipar=1,bc%nglob
+ read(IIN_NUC,*) inp_nx(ipar),inp_nz(ipar),loc_str(ipar),loc_dip(ipar),sigma0(ipar),tau0_str(ipar),tau0_dip(ipar), &
+ Rstress_str(ipar),Rstress_dip(ipar),static_fc(ipar),dyn_fc(ipar),swcd(ipar),cohes(ipar),tim_forcedRup(ipar)
+ enddo
+ close(IIN_NUC)
+
+ allocate(bc%T0(3,bc%nglob))
+ allocate( bc%swf )
+ allocate( bc%swf%mus(bc%nglob) )
+ allocate( bc%swf%mud(bc%nglob) )
+ allocate( bc%swf%Dc(bc%nglob) )
+ allocate( bc%swf%theta(bc%nglob) )
+ allocate( bc%swf%T(bc%nglob) )
+ allocate( bc%swf%C(bc%nglob) )
+
+ do iLoc=1,bc%nglob
+
+ ipar = minloc( (-24000.0+loc_str(:)-bc%coord(1,iLoc))**2 + (-loc_dip(:)-bc%coord(3,iLoc))**2 , 1) !iloc_dip is negative of Z-coord
+ bc%T0(3,iLoc) = -sigma0(ipar)
+ bc%T0(1,iLoc) = tau0_str(ipar)
+ bc%T0(2,iLoc) = tau0_dip(ipar)
+
+ bc%swf%mus(iLoc) = static_fc(ipar)
+ bc%swf%mud(iLoc) = dyn_fc(ipar)
+ bc%swf%Dc(iLoc) = swcd(ipar)
+ bc%swf%C(iLoc) = cohes(ipar)
+ bc%swf%T(iLoc) = tim_forcedRup(ipar)
+ enddo
- open(unit=IIN_NUC,file='DATA/FAULT/input_file.txt',status='old',iostat=ier)
- read(IIN_NUC,*) relz_num,sub_relz_num
- read(IIN_NUC,*) num_cell_str,num_cell_dip,siz_str,siz_dip
- read(IIN_NUC,*) hypo_cell_str,hypo_cell_dip,hypo_loc_str,hypo_loc_dip,rad_T_str,rad_T_dip
- do ipar=1,bc%nglob
- read(IIN_NUC,*) inp_nx(ipar),inp_nz(ipar),loc_str(ipar),loc_dip(ipar),sigma0(ipar),tau0_str(ipar),tau0_dip(ipar), &
- Rstress_str(ipar),Rstress_dip(ipar),static_fc(ipar),dyn_fc(ipar),swcd(ipar),cohes(ipar),tim_forcedRup(ipar)
- enddo
- close(IIN_NUC)
-
- allocate(bc%T0(3,bc%nglob))
- allocate( bc%swf )
- allocate( bc%swf%mus(bc%nglob) )
- allocate( bc%swf%mud(bc%nglob) )
- allocate( bc%swf%Dc(bc%nglob) )
- allocate( bc%swf%theta(bc%nglob) )
- allocate( bc%swf%T(bc%nglob) )
- allocate( bc%swf%C(bc%nglob) )
-
- do iLoc=1,bc%nglob
-
- ipar = minloc( (-24000.0+loc_str(:)-bc%coord(1,iLoc))**2 + (-loc_dip(:)-bc%coord(3,iLoc))**2 , 1) !iloc_dip is negative of Z-coord
- bc%T0(3,iLoc) = -sigma0(ipar)
- bc%T0(1,iLoc) = tau0_str(ipar)
- bc%T0(2,iLoc) = tau0_dip(ipar)
-
- bc%swf%mus(iLoc) = static_fc(ipar)
- bc%swf%mud(iLoc) = dyn_fc(ipar)
- bc%swf%Dc(iLoc) = swcd(ipar)
- bc%swf%C(iLoc) = cohes(ipar)
- bc%swf%T(iLoc) = tim_forcedRup(ipar)
- enddo
-
endif
+ ! Set friction parameters and initialize friction variables
+ ! Slip weakening friction
if(.not. Rate_AND_State) then
- ! Set friction parameters and initialize friction variables
allocate( bc%swf )
allocate( bc%swf%mus(bc%nglob) )
allocate( bc%swf%mud(bc%nglob) )
@@ -425,8 +426,8 @@
allocate(bc%MU(bc%nglob))
bc%MU = swf_mu(bc%swf)
- else !Rate AND State
- ! Set friction parameters and initialize friction variables
+ ! Rate and state friction
+ else
allocate( bc%rsf )
allocate( bc%rsf%V0(bc%nglob) )
allocate( bc%rsf%f0(bc%nglob) )
@@ -435,7 +436,6 @@
allocate( bc%rsf%L(bc%nglob) )
allocate( bc%rsf%V_init(bc%nglob) )
allocate( bc%rsf%theta(bc%nglob) )
- ! WARNING: if V_HEALING is negative we turn off healing
V0 =1.e-6_CUSTOM_REAL
f0 =0.6_CUSTOM_REAL
@@ -636,21 +636,18 @@
real(kind=CUSTOM_REAL), dimension(bc%nglob) :: strength
real(kind=CUSTOM_REAL), dimension(3,bc%nglob) :: T
- real(kind=CUSTOM_REAL), dimension(bc%nglob) :: t1,t2,tnorm,tnew
+ real(kind=CUSTOM_REAL), dimension(bc%nglob) :: tStick,tnew
real(kind=CUSTOM_REAL), dimension(3,bc%nglob) :: dD,dV,dA
real(kind=CUSTOM_REAL), dimension(bc%nglob) :: theta_old, Vnorm, Vnorm_old
real(kind=CUSTOM_REAL) :: half_dt
- real(kind=CUSTOM_REAL), dimension(bc%nglob) :: tStick,ta
- real(kind=CUSTOM_REAL), dimension(bc%nglob) :: Vf,Vf1,tau1,Vf2,tau2,Vfavg
+ real(kind=CUSTOM_REAL), dimension(bc%nglob) :: ta, Vf,Vf1,tau1,Vf2,tau2,Vfavg
integer, dimension(bc%nglob) :: ierr
real(kind=CUSTOM_REAL), dimension(bc%nglob) :: TxExt
real(kind=CUSTOM_REAL) :: TLoad,DTau0,GLoad
real(kind=CUSTOM_REAL) :: time
- real(kind=CUSTOM_REAL), dimension(bc%nglob) :: thetan
half_dt = 0.5e0_CUSTOM_REAL*bc%dt
- theta_old = bc%swf%theta
Vnorm_old = sqrt(bc%V(1,:)*bc%V(1,:)+bc%V(2,:)*bc%V(2,:))
! get predicted values
@@ -685,35 +682,32 @@
if(.not. Rate_AND_State) then ! Update slip weakening friction:
! Update slip state variable
! WARNING: during opening the friction state variable should not evolve
+ theta_old = bc%swf%theta
call swf_update_state(bc%D,dD,bc%V,bc%swf)
- if(.not. TPV16) then
! Update friction coeficient
- bc%MU = swf_mu(bc%swf)
- else
- bc%MU = swf_mu_TPV16(bc%swf,it*bc%dt,bc%nglob)
- endif
+ bc%MU = swf_mu(bc%swf)
! combined with time-weakening for nucleation
! if (associated(bc%twf)) bc%MU = min( bc%MU, twf_mu(bc%twf,bc%coord,time) )
+ if (TPV16) then
+ where (f%T <= it*bc%dt) bc%MU = f%mud
+ endif
! Update strength
strength = -bc%MU * min(T(3,:),0.e0_CUSTOM_REAL) + bc%swf%C
! Solve for shear stress
- tnorm = sqrt( T(1,:)*T(1,:) + T(2,:)*T(2,:))
- tnorm = max(tnorm,1e0_CUSTOM_REAL)
- t1 = T(1,:)/tnorm
- t2 = T(2,:)/tnorm
+ tStick = sqrt( T(1,:)*T(1,:) + T(2,:)*T(2,:))
+ tnew = min(tStick,strength)
+ tStick = max(tStick,1e0_CUSTOM_REAL)
+ T(1,:) = tnew * T(1,:)/tStick
+ T(2,:) = tnew * T(2,:)/tStick
- tnew = min(tnorm,strength)
- T(1,:) = tnew * t1
- T(2,:) = tnew * t2
-
else ! Update rate and state friction:
! smooth loading within nucleation patch
- !WARNING : ad hoc for SCEC becnhmark TPV10x
+ !WARNING : ad hoc for SCEC benchmark TPV10x
TxExt = 0._CUSTOM_REAL
TLoad = 1.0_CUSTOM_REAL
DTau0 = 25e6_CUSTOM_REAL
@@ -728,13 +722,12 @@
TxExt = DTau0 * bc%asp%Fload * GLoad
T(1,:) = T(1,:) + TxExt
- tStick = sqrt(T(1,:)**2 + T(2,:)**2)
-
- Vf = sqrt(bc%V(1,:)**2 + bc%V(2,:)**2)
- thetan = bc%rsf%theta
+ Vf = Vnorm_old
+ theta_old = bc%rsf%theta
call rsf_update_state(Vf,bc%dt,bc%rsf)
ta = sqrt(bc%T(1,:)**2 + bc%T(2,:)**2)
+ tStick = sqrt(T(1,:)**2 + T(2,:)**2)
call NRsearch(Vf1,tau1,bc%rsf%f0,bc%rsf%V0,bc%rsf%a,bc%rsf%b,-T(3,:),log(bc%rsf%theta),bc%Z,ta,tStick,ierr,bc%nglob)
! If Vf1 or tau1 is negative, set Vf1 = 0 and tau1 = TauStick
@@ -749,9 +742,10 @@
!call exit_MPI(myrank,'NR search 2nd loop failed!')
stop
endif
+
! Updating state variable: 2nd loop
- Vfavg = 0.5e0_CUSTOM_REAL**abs(Vf + Vf1)
- bc%rsf%theta = thetan
+ Vfavg = 0.5e0_CUSTOM_REAL*(Vf + Vf1)
+ bc%rsf%theta = theta_old
call rsf_update_state(Vfavg,bc%dt,bc%rsf)
! NR search 2nd loop
@@ -770,10 +764,8 @@
endif
tStick = max(tStick,1e0_CUSTOM_REAL)
- t1 = T(1,:)/tnorm
- t2 = T(2,:)/tnorm
- T(1,:) = tau2 * t1
- T(2,:) = tau2 * t2
+ T(1,:) = tau2 * T(1,:)/tStick
+ T(2,:) = tau2 * T(2,:)/tStick
endif
@@ -895,60 +887,12 @@
endif
end subroutine swf_update_state
-!=====================================================================
-subroutine rsf_update_state(V,dt,f)
-
- real(kind=CUSTOM_REAL), dimension(:), intent(in) :: V
- type(rsf_type), intent(inout) :: f
- real(kind=CUSTOM_REAL), intent(in) :: dt
-
- real(kind=CUSTOM_REAL) :: vDtL(size(V))
-
- vDtL = V * dt * f%L
-
- ! ageing law
- if (f%StateLaw == 1) then
- where(vDtL > 1.e-5_CUSTOM_REAL)
- f%theta = f%theta * exp(-vDtL) + f%L/V * (ONE - exp(-vDtL))
- elsewhere
- f%theta = f%theta * exp(-vDtL) + dt - HALF*V/f%L*dt**2
- endwhere
-
- ! slip law
- else
- f%theta = f%L/V * (f%theta*V/f%L)**(exp(-vDtL))
- endif
-
-end subroutine rsf_update_state
-
-!=====================================================================
-! Friction coefficient
-function swf_mu_TPV16(f,time,nglob) result(mu)
-
- type(swf_type), intent(in) :: f
- real(kind=CUSTOM_REAL) :: mu(size(f%theta))
- real(kind=CUSTOM_REAL) :: time
-
- integer :: i,nglob
-
- !-- linear slip weakening:
- mu = f%mus -(f%mus-f%mud)/f%dc *f%theta
- do i=1,nglob
- if ( (f%theta(i) >= f%Dc(i)) .or. (time >= f%T(i)) ) then
- mu(i) = f%mud(i)
- endif
- enddo
-
-end function swf_mu_TPV16
-
-!=====================================================================
-! Friction coefficient
+!---------------------------------------------------------------------
function swf_mu(f) result(mu)
type(swf_type), intent(in) :: f
real(kind=CUSTOM_REAL) :: mu(size(f%theta))
- !-- linear slip weakening:
mu = f%mus -(f%mus-f%mud)/f%dc *f%theta
mu = max( mu, f%mud)
@@ -969,7 +913,33 @@
end function rsf_mu
-!===============================================================
+!---------------------------------------------------------------------
+subroutine rsf_update_state(V,dt,f)
+
+ real(kind=CUSTOM_REAL), dimension(:), intent(in) :: V
+ type(rsf_type), intent(inout) :: f
+ real(kind=CUSTOM_REAL), intent(in) :: dt
+
+ real(kind=CUSTOM_REAL) :: vDtL(size(V))
+
+ vDtL = V * dt / f%L
+
+ ! ageing law
+ if (f%StateLaw == 1) then
+ where(vDtL > 1.e-5_CUSTOM_REAL)
+ f%theta = f%theta * exp(-vDtL) + f%L/V * (ONE - exp(-vDtL))
+ elsewhere
+ f%theta = f%theta * exp(-vDtL) + dt*( ONE - HALF*vDtL )
+ endwhere
+
+ ! slip law
+ else
+ f%theta = f%L/V * (f%theta*V/f%L)**(exp(-vDtL))
+ endif
+
+end subroutine rsf_update_state
+
+!---------------------------------------------------------------
! Newton-Raphson search for rate-and-state solver
!
! NRsearch(Vfout,tauout,fo,Vo,cca,ccb,Seff,psi,FaultZ,tau,TauStick,ierror)
@@ -1059,6 +1029,7 @@
Vfout(iglob) = Vo(iglob)*(help1(iglob) - help2(iglob))
tauout(iglob) = tau(iglob)
enddo
+
end subroutine NRsearch
!===============================================================
More information about the CIG-COMMITS
mailing list