[cig-commits] r19876 - seismo/3D/FAULT_SOURCE/branches/new_fault_db/src
surendra at geodynamics.org
surendra at geodynamics.org
Mon Mar 26 13:02:19 PDT 2012
Author: surendra
Date: 2012-03-26 13:02:18 -0700 (Mon, 26 Mar 2012)
New Revision: 19876
Modified:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
Log:
Move loop over fault nodes outside of NRsearch. Corrected typos in forced Time rupture
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 18:26:58 UTC (rev 19875)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90 2012-03-26 20:02:18 UTC (rev 19876)
@@ -642,7 +642,7 @@
real(kind=CUSTOM_REAL) :: half_dt
real(kind=CUSTOM_REAL), dimension(bc%nglob) :: ta, Vf,Vf1,tau1,Vf2,tau2,Vfavg
- integer, dimension(bc%nglob) :: ierr
+ integer :: ierr,iNode
real(kind=CUSTOM_REAL), dimension(bc%nglob) :: TxExt
real(kind=CUSTOM_REAL) :: TLoad,DTau0,GLoad
real(kind=CUSTOM_REAL) :: time
@@ -691,7 +691,7 @@
! 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
+ where (bc%swf%T <= it*bc%dt) bc%MU = bc%swf%mud
endif
! Update strength
@@ -728,40 +728,44 @@
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)
+ do iNode=1,bc%nglob
+ call NRsearch(Vf1(iNode),tau1(iNode),bc%rsf%f0(iNode),bc%rsf%V0(iNode),bc%rsf%a(iNode),bc%rsf%b(iNode),-T(3,iNode),log(bc%rsf%theta(iNode)),bc%Z(iNode),ta(iNode),tStick(iNode),ierr)
- ! 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)
- tau1 = tStick
- Vf1 = 0._CUSTOM_REAL
- endwhere
+ ! If Vf1 or tau1 is negative, set Vf1 = 0 and tau1 = TauStick
+ if (Vf1(iNode) < 0._CUSTOM_REAL .or. tau1(iNode) < 0._CUSTOM_REAL .or. ierr == 1) then
+ tau1(iNode) = tStick(iNode)
+ Vf1(iNode) = 0._CUSTOM_REAL
+ endif
+
+ if (ierr == 1) then
+ print*, 'NR search 1st loop failed!'
+ print*, 'iteration = ', it
+ !call exit_MPI(myrank,'NR search 1st loop failed!')
+ stop
+ endif
+ enddo
- 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
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
- call NRsearch(Vf2,tau2,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)
+ do iNode=1,bc%nglob
+ call NRsearch(Vf2(iNode),tau2(iNode),bc%rsf%f0(iNode),bc%rsf%V0(iNode),bc%rsf%a(iNode),bc%rsf%b(iNode),-T(3,iNode),log(bc%rsf%theta(iNode)),bc%Z(iNode),ta(iNode),tStick(iNode),ierr)
- where (Vf2 < 0._CUSTOM_REAL .or. tau2 < 0._CUSTOM_REAL .or. ierr == 1)
- tau2 = tStick
- Vf2 = 0._CUSTOM_REAL
- endwhere
+ if (Vf2(iNode) < 0._CUSTOM_REAL .or. tau2(iNode) < 0._CUSTOM_REAL .or. ierr == 1) then
+ tau2(iNode) = tStick(iNode)
+ Vf2(iNode) = 0._CUSTOM_REAL
+ endif
- if (any(ierr == 1) ) then
- print*, 'NR search 2nd loop failed!'
- print*, 'iteration = ', it
- !call exit_MPI(myrank,'NR search 2nd loop failed!')
- stop
- endif
+ if (ierr == 1) then
+ print*, 'NR search 2nd loop failed!'
+ print*, 'iteration = ', it
+ !call exit_MPI(myrank,'NR search 2nd loop failed!')
+ stop
+ endif
+ enddo
tStick = max(tStick,1e0_CUSTOM_REAL)
T(1,:) = tau2 * T(1,:)/tStick
@@ -957,78 +961,70 @@
! OUTPUT VARIABLES Vfout,tauout
!
! Yoshihiro Kaneko: ykaneko at gps.caltech.edu
-! Surendra : Modified to loop over fault nodes
-subroutine NRsearch(Vfout,tauout,fo,Vo,cca,ccb,Seff,psi,FaultZ,tau,TauStick,ierror,nglob)
+subroutine NRsearch(Vfout,tauout,fo,Vo,cca,ccb,Seff,psi,FaultZ,tau,TauStick,ierror)
- !implicit none
- !include "constants.h"
-
- integer :: iglob,nglob
-
- 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
+ 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
- integer, dimension(nglob) :: ierror
+ integer :: ierror
- do iglob=1,nglob
- cA(iglob) = cca(iglob)*Seff(iglob)
- eps(iglob) = 0.001*cA(iglob)
- k=0
- delta(iglob) = HUGEVAL
+ cA = cca*Seff
+ eps = 0.001*cA
+ k=0
+ delta = HUGEVAL
- 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
+ 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(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
+ 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
- END DO
+ 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
+ !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
end subroutine NRsearch
More information about the CIG-COMMITS
mailing list