[cig-commits] r19869 - seismo/3D/FAULT_SOURCE/branches/new_fault_db/src
surendra at geodynamics.org
surendra at geodynamics.org
Mon Mar 26 06:45:16 PDT 2012
Author: surendra
Date: 2012-03-26 06:45:16 -0700 (Mon, 26 Mar 2012)
New Revision: 19869
Modified:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
Log:
Merged TPV16 implementation into the flow of code
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 13:13:17 UTC (rev 19868)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90 2012-03-26 13:45:16 UTC (rev 19869)
@@ -187,11 +187,7 @@
allocate( faults(nbfaults) )
do iflt=1,nbfaults
read(IIN_PAR,nml=BEGIN_FAULT,end=100)
- if(.not.TPV16) then
- call init_one_fault(faults(iflt),IIN_BIN,IIN_PAR,Minv,dt,nt,iflt)
- else
- call init_one_fault_TPV16(faults(iflt),IIN_BIN,IIN_PAR,Minv,dt,nt,iflt)
- endif
+ call init_one_fault(faults(iflt),IIN_BIN,IIN_PAR,Minv,dt,nt,iflt)
enddo
close(IIN_BIN)
close(IIN_PAR)
@@ -214,10 +210,9 @@
end subroutine BC_DYNFLT_init
-
!---------------------------------------------------------------------
-subroutine init_one_fault_TPV16(bc,IIN_BIN,IIN_PAR,Minv,dt,NT,iflt)
+subroutine init_one_fault(bc,IIN_BIN,IIN_PAR,Minv,dt,NT,iflt)
type(bc_dynflt_type), intent(inout) :: bc
real(kind=CUSTOM_REAL), intent(in) :: Minv(:)
@@ -233,7 +228,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
+ real(kind=CUSTOM_REAL) :: AspTx,Fload
+ integer :: nAsp,nFload
+
integer, parameter :: IIN_NUC =270
integer :: ipar
integer :: ier
@@ -248,12 +248,13 @@
real(kind=CUSTOM_REAL), dimension(:), allocatable :: loc_str,loc_dip,sigma0,tau0_str,tau0_dip,Rstress_str,Rstress_dip,static_fc, &
dyn_fc,swcd,cohes,tim_forcedRup
- !integer, dimension(:), allocatable :: iLoc
-
integer :: iLoc
+
NAMELIST / INIT_STRESS / S1,S2,S3,n1,n2,n3
NAMELIST / SWF / mus,mud,dc,nmus,nmud,ndc
+ NAMELIST / RSF / V0,f0,a,b,L,V_init,theta_init,nV0,nf0,na,nb,nL,nV_init,ntheta_init
+ NAMELIST / ASP / AspTx,Fload,nAsp,nFload
read(IIN_BIN) bc%nspec,bc%nglob
if (bc%nspec==0) return
@@ -274,7 +275,7 @@
read(IIN_BIN) bc%coord(3,:)
bc%dt = dt
- allocate( bc%B(bc%nglob) )
+ allocate( bc%B(bc%nglob) )
bc%B = 0e0_CUSTOM_REAL
allocate( nx(bc%nglob),ny(bc%nglob),nz(bc%nglob) )
nx = 0e0_CUSTOM_REAL
@@ -289,17 +290,11 @@
bc%B(k) = bc%B(k) + jacobian2Dw(ij,e)
enddo
enddo
-
- !JPA assemble bc%B across processors
- !JPA 1. put bc%B in the vector 'accel', which is defined in every GLL node including in the bulk (contains zeros off the fault)
- !JPA 2. call assemble_MPI_vector_* (or maybe assemble_MPI_scalar_*)
- !JPA do the same for bc%n
-
do k=1,bc%nglob
norm = sqrt( nx(k)*nx(k) + ny(k)*ny(k) + nz(k)*nz(k) )
nx(k) = nx(k) / norm
- ny(k) = ny(k) / norm
- nz(k) = nz(k) / norm
+ ny(k) = ny(k) / norm
+ nz(k) = nz(k) / norm
enddo
allocate( bc%R(3,3,bc%nglob) )
@@ -325,7 +320,31 @@
bc%D = 0e0_CUSTOM_REAL
bc%V = 0e0_CUSTOM_REAL
+ ! Set initial fault stresses
+ allocate(bc%T0(3,bc%nglob))
+ S1 = 0e0_CUSTOM_REAL
+ S2 = 0e0_CUSTOM_REAL
+ S3 = 0e0_CUSTOM_REAL
+ n1=0
+ n2=0
+ n3=0
+ read(IIN_PAR, nml=INIT_STRESS)
+ bc%T0(1,:) = S1
+ bc%T0(2,:) = S2
+ bc%T0(3,:) = S3
+ call init_2d_distribution(bc%T0(1,:),bc%coord,IIN_PAR,n1)
+ call init_2d_distribution(bc%T0(2,:),bc%coord,IIN_PAR,n2)
+ call init_2d_distribution(bc%T0(3,:),bc%coord,IIN_PAR,n3)
+
+ !WARNING : Quick and dirty free surface condition at z=0
+ ! do k=1,bc%nglob
+ ! if (abs(bc%zcoord(k)-0.e0_CUSTOM_REAL) <= SMALLVAL) bc%T0(2,k) = 0
+ ! end do
+
+
+ if(TPV16) then
+
allocate(inp_nx(bc%nglob))
allocate(inp_nz(bc%nglob))
allocate(loc_str(bc%nglob))
@@ -351,22 +370,6 @@
enddo
close(IIN_NUC)
-
- print*,minval(inp_nx),maxval(inp_nx)
- print*,minval(inp_nz),maxval(inp_nz)
- print*,minval(loc_str),maxval(loc_str)
- print*,minval(loc_dip),maxval(loc_dip)
- print*,minval(sigma0),maxval(sigma0)
- print*,minval(tau0_str),maxval(tau0_str)
- print*,minval(tau0_dip),maxval(tau0_dip)
- print*,minval(Rstress_str),maxval(Rstress_str)
- print*,minval(Rstress_dip),maxval(Rstress_dip)
- print*,minval(static_fc),maxval(static_fc)
- print*,minval(dyn_fc),maxval(dyn_fc)
- print*,minval(swcd),maxval(swcd)
- print*,minval(cohes),maxval(cohes)
- print*,minval(tim_forcedRup),maxval(tim_forcedRup)
-
allocate(bc%T0(3,bc%nglob))
allocate( bc%swf )
allocate( bc%swf%mus(bc%nglob) )
@@ -376,19 +379,6 @@
allocate( bc%swf%T(bc%nglob) )
allocate( bc%swf%C(bc%nglob) )
- ! do ipar=1,bc%nglob
- ! iLoc = minloc( (loc_str(ipar)-bc%coord(1,:))**2 + (-loc_dip(ipar)-bc%coord(3,:))**2 , 1) !iloc_dip is negative of Z-coord
- ! bc%T0(1,iLoc) = sigma0(ipar)
- ! bc%T0(2,iLoc) = tau0_str(ipar)
- ! bc%T0(3,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
-
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
@@ -403,151 +393,9 @@
bc%swf%T(iLoc) = tim_forcedRup(ipar)
enddo
- print*,' '
- print*,minval(bc%T0(1,:)),maxval(bc%T0(1,:))
- print*,minval(bc%T0(2,:)),maxval(bc%T0(2,:))
- print*,minval(bc%T0(3,:)),maxval(bc%T0(3,:))
- print*,minval(bc%swf%mus),maxval(bc%swf%mus)
- print*,minval(bc%swf%mud),maxval(bc%swf%mud)
- print*,minval(bc%swf%Dc),maxval(bc%swf%Dc)
- print*,minval(bc%swf%C),maxval(bc%swf%C)
- print*,minval(bc%swf%T),maxval(bc%swf%T)
+ endif
- ! WARNING: if V_HEALING is negative we turn off healing
- bc%swf%healing = (V_HEALING > 0e0_CUSTOM_REAL)
- bc%swf%theta = 0e0_CUSTOM_REAL
- allocate(bc%MU(bc%nglob))
- bc%MU = swf_mu_TPV16(bc%swf,0.0,bc%nglob)
-
- call init_dataT(bc%dataT,bc%coord,bc%nglob,NT,iflt)
- call init_dataXZ(bc%dataXZ,bc,bc%nglob)
-
- bc%DataXZ%t1 => bc%T0(1,:)
- bc%DataXZ%t2 => bc%T0(2,:)
- bc%DataXZ%t3 => bc%T0(3,:)
- call write_dataXZ(bc%dataXZ,0,iflt)
- bc%DataXZ%t1 => bc%t(1,:)
- bc%DataXZ%t2 => bc%t(2,:)
- bc%DataXZ%t3 => bc%t(3,:)
-
-
-end subroutine init_one_fault_TPV16
-
-!---------------------------------------------------------------------
-
-subroutine init_one_fault(bc,IIN_BIN,IIN_PAR,Minv,dt,NT,iflt)
-
- type(bc_dynflt_type), intent(inout) :: bc
- real(kind=CUSTOM_REAL), intent(in) :: Minv(:)
- integer, intent(in) :: IIN_BIN,IIN_PAR,NT,iflt
- real(kind=CUSTOM_REAL), intent(in) :: dt
-
- real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: jacobian2Dw
- real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: normal
- integer, dimension(:,:), allocatable :: ibool1
- real(kind=CUSTOM_REAL) :: norm
- real(kind=CUSTOM_REAL) :: S1,S2,S3
- integer :: n1,n2,n3
- 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) :: AspTx,Fload
- integer :: nAsp,nFload
-
- NAMELIST / INIT_STRESS / S1,S2,S3,n1,n2,n3
- NAMELIST / SWF / mus,mud,dc,nmus,nmud,ndc
- NAMELIST / RSF / V0,f0,a,b,L,V_init,theta_init,nV0,nf0,na,nb,nL,nV_init,ntheta_init
- NAMELIST / ASP / AspTx,Fload,nAsp,nFload
-
- read(IIN_BIN) bc%nspec,bc%nglob
- if (bc%nspec==0) return
- allocate( bc%ibulk1(bc%nglob) )
- allocate( bc%ibulk2(bc%nglob) )
- allocate( ibool1(NGLLSQUARE,bc%nspec) )
- allocate(normal(NDIM,NGLLSQUARE,bc%nspec))
- allocate(jacobian2Dw(NGLLSQUARE,bc%nspec))
-
- allocate(bc%coord(3,(bc%nglob)))
- read(IIN_BIN) ibool1
- read(IIN_BIN) jacobian2Dw
- read(IIN_BIN) normal
- read(IIN_BIN) bc%ibulk1
- read(IIN_BIN) bc%ibulk2
- read(IIN_BIN) bc%coord(1,:)
- read(IIN_BIN) bc%coord(2,:)
- read(IIN_BIN) bc%coord(3,:)
- bc%dt = dt
-
- allocate( bc%B(bc%nglob) )
- bc%B = 0e0_CUSTOM_REAL
- allocate( nx(bc%nglob),ny(bc%nglob),nz(bc%nglob) )
- nx = 0e0_CUSTOM_REAL
- ny = 0e0_CUSTOM_REAL
- nz = 0e0_CUSTOM_REAL
- do e=1,bc%nspec
- do ij = 1,NGLLSQUARE
- k = ibool1(ij,e)
- nx(k) = nx(k) + normal(1,ij,e)
- ny(k) = ny(k) + normal(2,ij,e)
- nz(k) = nz(k) + normal(3,ij,e)
- bc%B(k) = bc%B(k) + jacobian2Dw(ij,e)
- enddo
- enddo
- do k=1,bc%nglob
- norm = sqrt( nx(k)*nx(k) + ny(k)*ny(k) + nz(k)*nz(k) )
- nx(k) = nx(k) / norm
- ny(k) = ny(k) / norm
- nz(k) = nz(k) / norm
- enddo
-
- allocate( bc%R(3,3,bc%nglob) )
- call compute_R(bc%R,bc%nglob,nx,ny,nz)
-
- ! Needed in dA_Free = -K2*d2/M2 + K1*d1/M1
- allocate(bc%invM1(bc%nglob))
- allocate(bc%invM2(bc%nglob))
- bc%invM1 = Minv(bc%ibulk1)
- bc%invM2 = Minv(bc%ibulk2)
-
- ! Fault impedance, Z in : Trac=T_Stick-Z*dV
- ! Z = 1/( B1/M1 + B2/M2 ) / (0.5*dt)
- ! T_stick = Z*Vfree traction as if the fault was stuck (no displ discontinuity)
- ! NOTE: same Bi on both sides, see note above
- allocate(bc%Z(bc%nglob))
- bc%Z = 1.e0_CUSTOM_REAL/(0.5e0_CUSTOM_REAL*dt * bc%B *( bc%invM1 + bc%invM2 ))
-
- allocate(bc%T(3,bc%nglob))
- allocate(bc%D(3,bc%nglob))
- allocate(bc%V(3,bc%nglob))
- bc%T = 0e0_CUSTOM_REAL
- bc%D = 0e0_CUSTOM_REAL
- bc%V = 0e0_CUSTOM_REAL
-
- ! Set initial fault stresses
- allocate(bc%T0(3,bc%nglob))
- S1 = 0e0_CUSTOM_REAL
- S2 = 0e0_CUSTOM_REAL
- S3 = 0e0_CUSTOM_REAL
- n1=0
- n2=0
- n3=0
- read(IIN_PAR, nml=INIT_STRESS)
- bc%T0(1,:) = S1
- bc%T0(2,:) = S2
- bc%T0(3,:) = S3
-
- call init_2d_distribution(bc%T0(1,:),bc%coord,IIN_PAR,n1)
- call init_2d_distribution(bc%T0(2,:),bc%coord,IIN_PAR,n2)
- call init_2d_distribution(bc%T0(3,:),bc%coord,IIN_PAR,n3)
-
- !WARNING : Quick and dirty free surface condition at z=0
- ! do k=1,bc%nglob
- ! 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 )
@@ -771,136 +619,12 @@
if (.not. allocated(faults)) return
do iflt=1,size(faults)
- if(.not.TPV16) then
- if (faults(iflt)%nspec>0) call BC_DYNFLT_set3d(faults(iflt),F,Vel,Dis,iflt)
- else
- if (faults(iflt)%nspec>0) call BC_DYNFLT_set3d_TPV16(faults(iflt),F,Vel,Dis,iflt)
- endif
+ if (faults(iflt)%nspec>0) call BC_DYNFLT_set3d(faults(iflt),F,Vel,Dis,iflt)
enddo
end subroutine bc_dynflt_set3d_all
!---------------------------------------------------------------------
-subroutine BC_DYNFLT_set3d_TPV16(bc,MxA,V,D,iflt)
-
- use specfem_par, only:it,NSTEP
-
- real(kind=CUSTOM_REAL), intent(inout) :: MxA(:,:)
- type(bc_dynflt_type), intent(inout) :: bc
- real(kind=CUSTOM_REAL), intent(in) :: V(:,:),D(:,:)
- integer,intent(in) :: iflt
-
- 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(3,bc%nglob) :: dD,dV,dA
- real(kind=CUSTOM_REAL), dimension(bc%nglob) :: theta_old, Vnorm, Vnorm_old
- real(kind=CUSTOM_REAL) :: half_dt
- ! integer :: k
-
- 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
- dD = get_jump(bc,D) ! dD_predictor
- dV = get_jump(bc,V) ! dV_predictor
- dA = get_weighted_jump(bc,MxA) ! dA_free
-
- ! rotate to fault frame (tangent,normal)
- ! component 3 is normal to the fault
- dD = rotate(bc,dD,1)
- dV = rotate(bc,dV,1)
- dA = rotate(bc,dA,1)
-
-
- ! T_stick
- T(1,:) = bc%Z * ( dV(1,:) + half_dt*dA(1,:) )
- T(2,:) = bc%Z * ( dV(2,:) + half_dt*dA(2,:) )
- T(3,:) = bc%Z * ( dV(3,:) + half_dt*dA(3,:) )
-
- !Warning : dirty particular free surface condition z = 0.
- ! where (bc%zcoord(:) > - SMALLVAL) T(2,:) = 0
- ! do k=1,bc%nglob
- ! if (abs(bc%zcoord(k)-0.e0_CUSTOM_REAL) < SMALLVAL) T(2,k) = 0.e0_CUSTOM_REAL
- ! end do
-
- ! add initial stress
- T = T + bc%T0
-
- ! Solve for normal stress (negative is compressive)
- ! Opening implies free stress
- if (bc%allow_opening) T(3,:) = min(T(3,:),0.e0_CUSTOM_REAL)
-
- ! 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_TPV16(bc%swf,it*bc%dt,bc%nglob)
-
- ! combined with time-weakening for nucleation
- ! if (associated(bc%twf)) bc%MU = min( bc%MU, twf_mu(bc%twf,bc%coord,time) )
-
- ! 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
-
- tnew = min(tnorm,strength)
- T(1,:) = tnew * t1
- T(2,:) = tnew * t2
-
- ! Save total tractions
- bc%T = T
-
- ! Subtract initial stress
- T = T - bc%T0
-
- ! Update slip acceleration da=da_free-T/(0.5*dt*Z)
- dA(1,:) = dA(1,:) - T(1,:)/(bc%Z*half_dt)
- dA(2,:) = dA(2,:) - T(2,:)/(bc%Z*half_dt)
- dA(3,:) = dA(3,:) - T(3,:)/(bc%Z*half_dt)
-
- ! Update slip and slip rate, in fault frame
- bc%D = dD
- bc%V = dV + half_dt*dA
-
- ! Rotate tractions back to (x,y,z) frame
- T = rotate(bc,T,-1)
-
- ! Add boundary term B*T to M*a
- MxA(1,bc%ibulk1) = MxA(1,bc%ibulk1) + bc%B*T(1,:)
- MxA(2,bc%ibulk1) = MxA(2,bc%ibulk1) + bc%B*T(2,:)
- MxA(3,bc%ibulk1) = MxA(3,bc%ibulk1) + bc%B*T(3,:)
-
- MxA(1,bc%ibulk2) = MxA(1,bc%ibulk2) - bc%B*T(1,:)
- MxA(2,bc%ibulk2) = MxA(2,bc%ibulk2) - bc%B*T(2,:)
- MxA(3,bc%ibulk2) = MxA(3,bc%ibulk2) - bc%B*T(3,:)
-
-
- !-- intermediate storage of outputs --
- Vnorm = sqrt(bc%V(1,:)*bc%V(1,:)+bc%V(2,:)*bc%V(2,:))
- call store_dataXZ(bc%dataXZ, strength, theta_old, bc%swf%theta, bc%swf%dc, &
- Vnorm_old, Vnorm, it*bc%dt,bc%dt)
- call store_dataT(bc%dataT,bc%D,bc%V,bc%T,it)
-
-
- !-- 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)
- ! write dataXZ every NSNAP time step
- if ( mod(it,NSNAP) == 0) call write_dataXZ(bc%dataXZ,it,iflt)
- if ( it == NSTEP) call SCEC_Write_RuptureTime(bc%dataXZ,bc%dt,NSTEP,iflt)
-
-end subroutine BC_DYNFLT_set3d_TPV16
-
-!---------------------------------------------------------------------
subroutine BC_DYNFLT_set3d(bc,MxA,V,D,iflt)
use specfem_par, only:it,NSTEP
@@ -963,14 +687,18 @@
! WARNING: during opening the friction state variable should not evolve
call swf_update_state(bc%D,dD,bc%V,bc%swf)
+ if(TPV16) then
! Update friction coeficient
- bc%MU = swf_mu(bc%swf)
+ bc%MU = swf_mu(bc%swf)
+ else
+ bc%MU = swf_mu_TPV16(bc%swf,it*bc%dt,bc%nglob)
+ endif
! combined with time-weakening for nucleation
! if (associated(bc%twf)) bc%MU = min( bc%MU, twf_mu(bc%twf,bc%coord,time) )
! Update strength
- strength = -bc%MU * min(T(3,:),0.e0_CUSTOM_REAL)
+ 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,:))
More information about the CIG-COMMITS
mailing list