[cig-commits] r19863 - seismo/3D/FAULT_SOURCE/branches/new_fault_db/src
surendra at geodynamics.org
surendra at geodynamics.org
Sat Mar 24 10:46:48 PDT 2012
Author: surendra
Date: 2012-03-24 10:46:48 -0700 (Sat, 24 Mar 2012)
New Revision: 19863
Modified:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
Log:
Added slip-law and loop over fault nodes for NRsearch
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-24 12:04:53 UTC (rev 19862)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90 2012-03-24 17:46:48 UTC (rev 19863)
@@ -103,6 +103,8 @@
logical, save :: Rate_AND_State = .false.
+ logical, save :: SlipLaw = .false.
+
real(kind=CUSTOM_REAL), allocatable, save :: Kelvin_Voigt_eta(:)
public :: BC_DYNFLT_init, BC_DYNFLT_set3d_all, Kelvin_Voigt_eta, &
@@ -885,13 +887,10 @@
real(kind=CUSTOM_REAL), dimension(bc%nglob) :: theta_old, Vnorm, Vnorm_old
real(kind=CUSTOM_REAL) :: half_dt
- integer :: ierr
real(kind=CUSTOM_REAL), dimension(bc%nglob) :: tStick,ta,tau_ratio
- real(kind=CUSTOM_REAL), dimension(bc%nglob) :: Vf,Vf1,tau1
- real(kind=CUSTOM_REAL) :: ONE
- logical, dimension(bc%nglob) :: ierr
+ real(kind=CUSTOM_REAL), dimension(bc%nglob) :: Vf,Vf1,tau1,Vf2,tau2,Vfavg
+ integer, dimension(bc%nglob) :: ierr
- ONE = 1.e0_CUSTOM_REAL
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,:))
@@ -959,7 +958,7 @@
tau_ratio = T(2,:)/T(1,:)
tStick = sqrt(T(1,:)**2 + T(2,:)**2)
ta = sqrt(bc%T(1,:)**2 + bc%T(2,:)**2)
- call NRsearch(Vf1,tau1,bc%rs%fo,bc%rs%Vo,bc%rs%a,bc%rs%b,-T(3,:),log(bc%rs%theta),bc%Z,ta,tStick,ierr)
+ call NRsearch(Vf1,tau1,bc%rs%f0,bc%rs%V0,bc%rs%a,bc%rs%b,-T(3,:),log(bc%rs%theta),bc%Z,ta,tStick,ierr,bc%nglob)
! If Vf1 or tau1 is negative, set Vf1 = 0 and tau1 = TauStick
where (Vf1 < 0._CUSTOM_REAL .or. tau1 < 0._CUSTOM_REAL .or. ierr == 1)
@@ -967,9 +966,11 @@
Vf1 = 0._CUSTOM_REAL
endwhere
- if (any(ierr) == 1) then
+ if (any(ierr == 1) ) then
print*, 'NR search 1st loop failed!'
print*, 'iteration = ', it
+ !call exit_MPI(myrank,'NR search 2nd loop failed!')
+ stop
endif
! Updating state variable: 2nd loop
@@ -977,7 +978,7 @@
call rs_update_state(Vfavg,bc%dt,bc%rs)
! NR search 2nd loop
- call NRsearch(Vf2,tau2,bc%rs%fo,bc%rs%Vo,bc%rs%a,bc%rs%b,-T(3,:),log(bc%rs%theta),bc%Z,ta,tStick,ierr)
+ call NRsearch(Vf2,tau2,bc%rs%f0,bc%rs%V0,bc%rs%a,bc%rs%b,-T(3,:),log(bc%rs%theta),bc%Z,ta,tStick,ierr,bc%nglob)
where (Vf2 < 0._CUSTOM_REAL .or. tau2 < 0._CUSTOM_REAL .or. ierr == 1)
tau2 = tStick
@@ -997,10 +998,11 @@
T(2,:) = -tau2*abs(tau_ratio)/sqrt(ONE+tau_ratio**2) - bc%T0(2,:)
endwhere
- if (any(ierr) == 1 ) then
+ if (any(ierr == 1) ) then
print*, 'NR search 2nd loop failed!'
print*, 'iteration = ', it
- call exit_MPI(myrank,'NR search 2nd loop failed!')
+ !call exit_MPI(myrank,'NR search 2nd loop failed!')
+ stop
endif
endif
@@ -1127,16 +1129,22 @@
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
+ type(rs_type), intent(inout) :: f
+ real(kind=CUSTOM_REAL), intent(in) :: dt
+ real(kind=CUSTOM_REAL) :: vDtL(size(V))
+
vDtL = V * dt * f%L
- where(vDtL > 1.e-20_CUSTOM_REAL)
- f%theta = f%theta * exp(vDtL) + f%L/V * (1 - exp(vDtL))
+ if(.not. SlipLaw) 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 - 0.5*V/f%L*dt**2
+ f%theta = f%theta * exp(-vDtL) + dt - HALF*V/f%L*dt**2
endwhere
+ else
+ f%theta = f%L/V * (f%theta*V/f%L)**(exp(-vDtL))
+ endif
end subroutine rs_update_state
@@ -1178,11 +1186,11 @@
function rs_mu(f,V) result(mu)
type(rs_type), intent(in) :: f
- real(kind=CUSTOM_REAL), dimension(:), intent(in) :: mu
real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: V
+ real(kind=CUSTOM_REAL) :: mu(size(V,2))
!-- 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))
+ 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
@@ -1513,73 +1521,78 @@
! OUTPUT VARIABLES Vfout,tauout
!
! Yoshihiro Kaneko: ykaneko at gps.caltech.edu
- !
- subroutine NRsearch(Vfout,tauout,fo,Vo,cca,ccb,Seff,psi,FaultZ,tau,TauStick,ierror)
+ ! Surendra : Modified to loop over fault nodes
+ subroutine NRsearch(Vfout,tauout,fo,Vo,cca,ccb,Seff,psi,FaultZ,tau,TauStick,ierror,nglob)
- implicit none
- include "constants.h"
+ !implicit none
+ !include "constants.h"
- real(kind=CUSTOM_REAL) :: cA,eps,delta,fa,help,help1,help2,Vfprime
- real(kind=CUSTOM_REAL) :: Vfout,tauout,fo,Vo,cca,ccb,Seff,tau,psi,FaultZ,TauStick
- real(kind=CUSTOM_REAL) :: f_bisec,f_bisec_Vflow,Vfbisec,Vfhigh,Vflow
- integer :: j,k,ierror
+ integer :: iglob,nglob
- cA = cca*Seff
- eps = 0.001*cA
- k=0
- delta = HUGEVAL
+ real(kind=CUSTOM_REAL), dimension(nglob) :: cA,eps,delta,fa,help,help1,help2,Vfprime
+ real(kind=CUSTOM_REAL), dimension(nglob) :: Vfout,tauout,fo,Vo,cca,ccb,Seff,tau,psi,FaultZ,TauStick
+ real(kind=CUSTOM_REAL), dimension(nglob) :: f_bisec,f_bisec_Vflow,Vfbisec,Vfhigh,Vflow
+ integer :: j,k
+ integer, dimension(nglob) :: ierror
- DO WHILE ( abs(delta) > eps )
- fa = tau/(Seff*cca)
- help = -(fo+ccb*psi)/cca
- help1 = exp(help+fa)
- help2 = exp(help-fa)
- Vfout = Vo*(help1 - help2)
- Vfprime = (Vo/cA)*(help1 + help2)
- delta = (TauStick - HALF*FaultZ*Vfout - tau)/(ONE + HALF*FaultZ*Vfprime)
- tau = tau + delta
- k = k + 1
- if ( abs(delta) > 1.e10 .OR. k > 100 ) then
- ! try bisection if NR search fails
- Vflow = 0._CUSTOM_REAL ! lower bound of Vf for bisection
- Vfhigh = 5._CUSTOM_REAL ! upper bound of Vf for bisection
- do j = 1,200
- Vfbisec = (Vflow + Vfhigh)/TWO
- f_bisec = TauStick - HALF*FaultZ*Vfbisec - cA*asinh(Vfbisec/(TWO*Vo)*exp(-help))
- if ( abs(f_bisec) < eps ) then
- Vfout = Vfbisec
- tauout = TauStick - HALF*FaultZ*Vfout
- return
- endif
- f_bisec_Vflow = TauStick - HALF*FaultZ*Vflow - cA*asinh(Vflow/(TWO*Vo)*exp(-help))
- if ( f_bisec*f_bisec_Vflow < 0._CUSTOM_REAL ) then
- Vfhigh = Vfbisec
- else
- Vflow = Vfbisec
- endif
- enddo
- ierror = 1
- !compute updated Vf before exiting NR search
- fa = tau/(Seff*cca)
- help = -(fo+ccb*psi)/cca
- help1 = exp(help+fa)
- help2 = exp(help-fa)
- Vfout = Vo*(help1 - help2)
- tauout = tau
- return
- endif
+ do iglob=1,nglob
+ cA(iglob) = cca(iglob)*Seff(iglob)
+ eps(iglob) = 0.001*cA(iglob)
+ k=0
+ delta(iglob) = HUGEVAL
- END DO
+ DO WHILE ( abs(delta(iglob)) > eps(iglob) )
+ fa(iglob) = tau(iglob)/(Seff(iglob)*cca(iglob))
+ help(iglob) = -(fo(iglob)+ccb(iglob)*psi(iglob))/cca(iglob)
+ help1(iglob) = exp(help(iglob)+fa(iglob))
+ help2(iglob) = exp(help(iglob)-fa(iglob))
+ Vfout(iglob) = Vo(iglob)*(help1(iglob) - help2(iglob))
+ Vfprime(iglob) = (Vo(iglob)/cA(iglob))*(help1(iglob) + help2(iglob))
+ delta(iglob) = (TauStick(iglob) - HALF*FaultZ(iglob)*Vfout(iglob) - tau(iglob))/(ONE + HALF*FaultZ(iglob)*Vfprime(iglob))
+ tau(iglob) = tau(iglob) + delta(iglob)
+ k = k + 1
+ if ( abs(delta(iglob)) > 1.e10 .OR. k > 100 ) then
+ ! try bisection if NR search fails
+ Vflow(iglob) = 0._CUSTOM_REAL ! lower bound of Vf for bisection
+ Vfhigh(iglob) = 5._CUSTOM_REAL ! upper bound of Vf for bisection
+ do j = 1,200
+ Vfbisec(iglob) = (Vflow(iglob) + Vfhigh(iglob))/TWO
+ f_bisec(iglob) = TauStick(iglob) - HALF*FaultZ(iglob)*Vfbisec(iglob) - cA(iglob)*asinh(Vfbisec(iglob)/(TWO*Vo(iglob))*exp(-help(iglob)))
+ if ( abs(f_bisec(iglob)) < eps(iglob) ) then
+ Vfout(iglob) = Vfbisec(iglob)
+ tauout(iglob) = TauStick(iglob) - HALF*FaultZ(iglob)*Vfout(iglob)
+ return
+ endif
+ f_bisec_Vflow(iglob) = TauStick(iglob) - HALF*FaultZ(iglob)*Vflow(iglob) - cA(iglob)*asinh(Vflow(iglob)/(TWO*Vo(iglob))*exp(-help(iglob)))
+ if ( f_bisec(iglob)*f_bisec_Vflow(iglob) < 0._CUSTOM_REAL ) then
+ Vfhigh(iglob) = Vfbisec(iglob)
+ else
+ Vflow(iglob) = Vfbisec(iglob)
+ endif
+ enddo
- !compute the updated Vf before exiting NR search
- fa = tau/(Seff*cca)
- help = -(fo+ccb*psi)/cca
- help1 = exp(help+fa)
- help2 = exp(help-fa)
- Vfout = Vo*(help1 - help2)
- tauout = tau
+ ierror(iglob) = 1
+ !compute updated Vf before exiting NR search
+ fa(iglob) = tau(iglob)/(Seff(iglob)*cca(iglob))
+ help(iglob) = -(fo(iglob)+ccb(iglob)*psi(iglob))/cca(iglob)
+ help1(iglob) = exp(help(iglob)+fa(iglob))
+ help2(iglob) = exp(help(iglob)-fa(iglob))
+ Vfout(iglob) = Vo(iglob)*(help1(iglob) - help2(iglob))
+ tauout(iglob) = tau(iglob)
+ return
+ endif
+ END DO
+
+ !compute the updated Vf before exiting NR search
+ fa(iglob) = tau(iglob)/(Seff(iglob)*cca(iglob))
+ help(iglob) = -(fo(iglob)+ccb(iglob)*psi(iglob))/cca(iglob)
+ help1(iglob) = exp(help(iglob)+fa(iglob))
+ help2(iglob) = exp(help(iglob)-fa(iglob))
+ Vfout(iglob) = Vo(iglob)*(help1(iglob) - help2(iglob))
+ tauout(iglob) = tau(iglob)
+ enddo
end subroutine NRsearch
!---------------------------------------------------------------
More information about the CIG-COMMITS
mailing list