[cig-commits] r20829 - in seismo/3D/FAULT_SOURCE/branches/new_fault_db: decompose_mesh_SCOTCH src
ampuero at geodynamics.org
ampuero at geodynamics.org
Fri Oct 12 09:16:34 PDT 2012
Author: ampuero
Date: 2012-10-12 09:16:33 -0700 (Fri, 12 Oct 2012)
New Revision: 20829
Modified:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/fault_scotch.f90
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_object.f90
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
Log:
fixed rate-and-state initialization of theta (made it more general); cleaned up
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/fault_scotch.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/fault_scotch.f90 2012-10-12 15:41:07 UTC (rev 20828)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/decompose_mesh_SCOTCH/fault_scotch.f90 2012-10-12 16:16:33 UTC (rev 20829)
@@ -1,8 +1,9 @@
module fault_scotch
+
implicit none
include "../constants.h"
-
private
+
type fault_type
private
integer :: nspec
@@ -266,7 +267,7 @@
integer, dimension(nspec) :: ind,ninseg,iwork
logical :: ifseg(nspec)
integer :: ispec,i,j
- integer :: ieoff,ilocnum,nseg,ioff,iseg
+ integer :: nseg,ioff,iseg
double precision :: SMALLVALTOL
xp=xyz_c(1,:)
@@ -278,10 +279,7 @@
! establish initial pointers
do ispec=1,nspec
- ieoff=1*(ispec-1)
- do ilocnum=1,1 !NGLLSQUARE
- loc(ilocnum+ieoff)=ilocnum+ieoff
- enddo
+ loc(ispec)=ispec
enddo
ifseg(:)=.false.
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_object.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_object.f90 2012-10-12 15:41:07 UTC (rev 20828)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_object.f90 2012-10-12 16:16:33 UTC (rev 20829)
@@ -308,7 +308,7 @@
double precision :: xp(npointot),yp(npointot),zp(npointot),xmin,xmax
integer :: loc(npointot)
logical :: ifseg(npointot)
- integer :: ispec,i,j,k,igll,ie,je,ke,e
+ integer :: ispec,k,igll,ie,je,ke,e
xmin = minval(nodes_coords_ext_mesh(1,:))
xmax = maxval(nodes_coords_ext_mesh(1,:))
@@ -316,18 +316,14 @@
k = 0
do e = 1,fdb%nspec
ispec = fdb%ispec1(e)
- igll = 0
- do i=1,NGLLX
- do j=1,NGLLX
- igll = igll + 1
- ie=fdb%ijk1(1,igll,e)
- je=fdb%ijk1(2,igll,e)
- ke=fdb%ijk1(3,igll,e)
- k = k+1
- xp(k) = xstore(ie,je,ke,ispec)
- yp(k) = ystore(ie,je,ke,ispec)
- zp(k) = zstore(ie,je,ke,ispec)
- enddo
+ do igll=1,NGLLSQUARE
+ ie=fdb%ijk1(1,igll,e)
+ je=fdb%ijk1(2,igll,e)
+ ke=fdb%ijk1(3,igll,e)
+ k = k+1
+ xp(k) = xstore(ie,je,ke,ispec)
+ yp(k) = ystore(ie,je,ke,ispec)
+ zp(k) = zstore(ie,je,ke,ispec)
enddo
enddo
allocate( fdb%ibool1(NGLLSQUARE,fdb%nspec) )
@@ -340,18 +336,14 @@
k = 0
do e = 1,fdb%nspec
ispec = fdb%ispec2(e)
- igll = 0
- do i=1,NGLLX
- do j=1,NGLLX
- igll = igll + 1
- ie=fdb%ijk2(1,igll,e)
- je=fdb%ijk2(2,igll,e)
- ke=fdb%ijk2(3,igll,e)
- k = k+1
- xp(k) = xstore(ie,je,ke,ispec)
- yp(k) = ystore(ie,je,ke,ispec)
- zp(k) = zstore(ie,je,ke,ispec)
- enddo
+ do igll=1,NGLLSQUARE
+ ie=fdb%ijk2(1,igll,e)
+ je=fdb%ijk2(2,igll,e)
+ ke=fdb%ijk2(3,igll,e)
+ k = k+1
+ xp(k) = xstore(ie,je,ke,ispec)
+ yp(k) = ystore(ie,je,ke,ispec)
+ zp(k) = zstore(ie,je,ke,ispec)
enddo
enddo
allocate( fdb%ibool2(NGLLSQUARE,fdb%nspec) )
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-10-12 15:41:07 UTC (rev 20828)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90 2012-10-12 16:16:33 UTC (rev 20829)
@@ -252,7 +252,6 @@
real(kind=CUSTOM_REAL) :: C,T
integer :: nC,nForcedRup
-
integer, parameter :: IIN_NUC =270
integer :: ipar
integer :: ier
@@ -267,16 +266,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 :: iLoc
real(kind=CUSTOM_REAL) :: minX
-
-
real(kind=CUSTOM_REAL) :: W1,W2,w,hypo_z
real(kind=CUSTOM_REAL) :: x,z
logical :: c1,c2,c3,c4
real(kind=CUSTOM_REAL) :: b11,b12,b21,b22,B1,B2
- integer iglob,nglob_bulk
+ integer :: i,nglob_bulk
double precision, dimension(:), allocatable :: mu1,mu2,mu3
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: init_vel
@@ -326,19 +322,19 @@
if (PARALLEL_FAULT) then
accel=0._CUSTOM_REAL
- if (bc%nspec>0) accel(1,bc%ibulk1) = bc%B(:)
+ if (bc%nspec>0) accel(1,bc%ibulk1) = bc%B
! assembles with other MPI processes
call assemble_MPI_vector_ext_mesh(NPROC,NGLOB_AB,accel, &
num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, &
my_neighbours_ext_mesh)
- if (bc%nspec>0) bc%B(:) = accel(1,bc%ibulk1)
+ if (bc%nspec>0) bc%B = accel(1,bc%ibulk1)
accel=0._CUSTOM_REAL
if (bc%nspec>0) then
- accel(1,bc%ibulk1) = nx(:)
- accel(2,bc%ibulk1) = ny(:)
- accel(3,bc%ibulk1) = nz(:)
+ accel(1,bc%ibulk1) = nx
+ accel(2,bc%ibulk1) = ny
+ accel(3,bc%ibulk1) = nz
endif
! assembles with other MPI processes
call assemble_MPI_vector_ext_mesh(NPROC,NGLOB_AB,accel, &
@@ -346,9 +342,9 @@
nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, &
my_neighbours_ext_mesh)
if (bc%nspec>0) then
- nx(:) = accel(1,bc%ibulk1)
- ny(:) = accel(2,bc%ibulk1)
- nz(:) = accel(3,bc%ibulk1)
+ nx = accel(1,bc%ibulk1)
+ ny = accel(2,bc%ibulk1)
+ nz = accel(3,bc%ibulk1)
endif
endif
@@ -446,18 +442,19 @@
minX = minval(bc%coord(1,:))
- do iLoc=1,bc%nglob
+ do i=1,bc%nglob
- ipar = minloc( (minX+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)
+ !loc_dip is negative of Z-coord
+ ipar = minloc( (minX+loc_str(:)-bc%coord(1,i))**2 + (-loc_dip(:)-bc%coord(3,i))**2 , 1)
+ bc%T0(3,i) = -sigma0(ipar)
+ bc%T0(1,i) = tau0_str(ipar)
+ bc%T0(2,i) = 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)
+ bc%swf%mus(i) = static_fc(ipar)
+ bc%swf%mud(i) = dyn_fc(ipar)
+ bc%swf%Dc(i) = swcd(ipar)
+ bc%swf%C(i) = cohes(ipar)
+ bc%swf%T(i) = tim_forcedRup(ipar)
enddo
endif
@@ -569,13 +566,15 @@
!!$ !init_vel = rotate(bc,init_vel,-1) ! directly assigned in global coordinates here as it is a simplified case
!!$ vel = vel + init_vel
+ !WARNING: below is an ad-hoc setting of a(x,z) for some SCEC benchmark
+ ! This should be instead an option in init_2d_distribution
W1=15000._CUSTOM_REAL
W2=7500._CUSTOM_REAL
w=3000._CUSTOM_REAL
hypo_z = -7500._CUSTOM_REAL
- do iglob=1,bc%nglob
- x=bc%coord(1,iglob)
- z=bc%coord(3,iglob)
+ do i=1,bc%nglob
+ x=bc%coord(1,i)
+ z=bc%coord(3,i)
c1=abs(x)<W1+w
c2=abs(x)>W1
c3=abs(z-hypo_z)<W2+w
@@ -602,30 +601,32 @@
B2 = 0._CUSTOM_REAL
endif
+ bc%rsf%a(i) = 0.008 + 0.008 * (ONE - B1*B2)
- bc%rsf%a(iglob) = 0.008 + 0.008 * (ONE - B1*B2)
elseif( abs(x)<=W1 .and. abs(z-hypo_z)<=W2 ) then
- bc%rsf%a(iglob) = 0.008
+ bc%rsf%a(i) = 0.008
else
- bc%rsf%a(iglob) = 0.016
+ bc%rsf%a(i) = 0.016
endif
- bc%rsf%theta(iglob) = L/V0 * exp( (bc%rsf%a(iglob) * log(2.0*sinh(-S1/S3/bc%rsf%a(iglob))) - f0 - bc%rsf%a(iglob)*log(V_init/V0) )/b )
-
-
enddo
-
-
+ ! WARNING: The line below scratches an earlier initialization of theta through theta_init
+ ! We should implement it as an option for the user
+ bc%rsf%theta = bc%rsf%L/bc%rsf%V0 &
+ * exp( ( bc%rsf%a * log(2.0*sinh(-sqrt(bc%T0(1,:)**2+bc%T0(2,:)**2)/bc%T0(3,:)/bc%rsf%a)) &
+ - bc%rsf%f0 - bc%rsf%a*log(bc%rsf%V_init/bc%rsf%V0) ) &
+ / bc%rsf%b )
+
allocate(bc%MU(bc%nglob))
+ ! WARNING: the line below is only valid for pure strike-slip faulting
bc%V(1,:) = bc%rsf%V_init
+
allocate( bc%asp )
allocate( bc%asp%Fload(bc%nglob) )
-
Fload =0.e0_CUSTOM_REAL
nFload =0
-
read(IIN_PAR, nml=ASP)
bc%asp%Fload =Fload
call init_2d_distribution(bc%asp%Fload,bc%coord,IIN_PAR,nFload)
@@ -776,16 +777,16 @@
! NOTE: On non-split nodes at fault edges, dD=dV=dA=0
! and the net contribution of B*T is =0
!
-subroutine bc_dynflt_set3d_all(F,Vel,Dis)
+subroutine bc_dynflt_set3d_all(F,V,D)
- real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: Vel,Dis
+ real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: V,D
real(kind=CUSTOM_REAL), dimension(:,:), intent(inout) :: F
- integer :: iflt
+ integer :: i
if (.not. allocated(faults)) return
- do iflt=1,size(faults)
- call BC_DYNFLT_set3d(faults(iflt),F,Vel,Dis,iflt)
+ do i=1,size(faults)
+ call BC_DYNFLT_set3d(faults(i),F,V,D,i)
enddo
end subroutine bc_dynflt_set3d_all
@@ -798,26 +799,19 @@
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
+ 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) :: tStick,tnew
- real(kind=CUSTOM_REAL), dimension(3,bc%nglob) :: dD,dV,dA
- real(kind=CUSTOM_REAL), dimension(bc%nglob) :: theta_old, theta_new, Vnorm, Vnorm_old, dc
- real(kind=CUSTOM_REAL) :: half_dt
+ real(kind=CUSTOM_REAL), dimension(3,bc%nglob) :: T,dD,dV,dA
+ real(kind=CUSTOM_REAL), dimension(bc%nglob) :: strength,tStick,tnew, &
+ theta_old, theta_new, dc, &
+ ta,Vf_old,Vf_new,tau1,TxExt
+ real(kind=CUSTOM_REAL) :: half_dt,TLoad,DTau0,GLoad,time
+ integer :: i
- real(kind=CUSTOM_REAL), dimension(bc%nglob) :: ta, Vf,Vf1,tau1,Vf2,tau2,Vfavg
- integer :: ierr,iNode
- real(kind=CUSTOM_REAL), dimension(bc%nglob) :: TxExt
- real(kind=CUSTOM_REAL) :: TLoad,DTau0,GLoad
- real(kind=CUSTOM_REAL) :: time
- real(kind=CUSTOM_REAL) :: psi
-
if (bc%nspec > 0) then !Surendra : for parallel faults
half_dt = 0.5e0_CUSTOM_REAL*bc%dt
- Vnorm_old = sqrt(bc%V(1,:)*bc%V(1,:)+bc%V(2,:)*bc%V(2,:))
+ Vf_old = sqrt(bc%V(1,:)*bc%V(1,:)+bc%V(2,:)*bc%V(2,:))
! get predicted values
dD = get_jump(bc,D) ! dD_predictor
@@ -848,7 +842,25 @@
! Opening implies free stress
if (bc%allow_opening) T(3,:) = min(T(3,:),0.e0_CUSTOM_REAL)
- if(.not. RATE_AND_STATE) then ! Update slip weakening friction:
+ ! smooth loading within nucleation patch
+ !WARNING : ad hoc for SCEC benchmark TPV10x
+ if (RATE_AND_STATE) then
+ TxExt = 0._CUSTOM_REAL
+ TLoad = 1.0_CUSTOM_REAL
+ DTau0 = 25e6_CUSTOM_REAL
+ time = it*bc%dt !time will never be zero. it starts from 1
+ if (time <= TLoad) then
+ GLoad = exp( (time-TLoad)*(time-Tload) / (time*(time-2.0_CUSTOM_REAL*TLoad)) )
+ else
+ GLoad = 1.0_CUSTOM_REAL
+ endif
+ TxExt = DTau0 * bc%asp%Fload * GLoad
+ T(1,:) = T(1,:) + TxExt
+ endif
+
+ tStick = sqrt( T(1,:)*T(1,:) + T(2,:)*T(2,:))
+
+ 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
@@ -867,63 +879,44 @@
strength = -bc%MU * min(T(3,:),0.e0_CUSTOM_REAL) + bc%swf%C
! Solve for shear stress
- 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
else ! Update rate and state friction:
-
- ! smooth loading within nucleation patch
- !WARNING : ad hoc for SCEC benchmark TPV10x
- TxExt = 0._CUSTOM_REAL
- TLoad = 1.0_CUSTOM_REAL
- DTau0 = 25e6_CUSTOM_REAL
- time = it*bc%dt !time will never be zero. it starts from 1
- if (time <= TLoad) then
- GLoad = exp( (time-TLoad)*(time-Tload) / (time*(time-2.0_CUSTOM_REAL*TLoad)) )
- else
- GLoad = 1.0_CUSTOM_REAL
- endif
- TxExt = DTau0 * bc%asp%Fload * GLoad
- T(1,:) = T(1,:) + TxExt
-
!JPA the solver below can be refactored into a loop with two passes
- Vf = Vnorm_old
+
+ ! first pass
theta_old = bc%rsf%theta
- call rsf_update_state(Vf,bc%dt,bc%rsf)
-
- tStick = sqrt(T(1,:)**2 + T(2,:)**2)
- do iNode=1,bc%nglob
- Vf1(iNode)=rtsafe(funcd,0.0,Vf(iNode)+5.0,1e-5,tStick(iNode),-T(3,iNode),bc%Z(iNode),bc%rsf%f0(iNode),bc%rsf%V0(iNode),bc%rsf%a(iNode),bc%rsf%b(iNode),bc%rsf%L(iNode),bc%rsf%theta(iNode))
- tau1(iNode) = tStick(iNode) - bc%Z(iNode)*Vf1(iNode)
+ call rsf_update_state(Vf_old,bc%dt,bc%rsf)
+ do i=1,bc%nglob
+ Vf_new(i)=rtsafe(funcd,0.0,Vf_old(i)+5.0,1e-5,tStick(i),-T(3,i),bc%Z(i),bc%rsf%f0(i), &
+ bc%rsf%V0(i),bc%rsf%a(i),bc%rsf%b(i),bc%rsf%L(i),bc%rsf%theta(i))
enddo
- ! Updating state variable: 2nd loop
- Vfavg = 0.5e0_CUSTOM_REAL*(Vf + Vf1)
+ ! second pass
bc%rsf%theta = theta_old
- call rsf_update_state(Vfavg,bc%dt,bc%rsf)
-
- ! NR search 2nd loop
- do iNode=1,bc%nglob
- Vf2(iNode)=rtsafe(funcd,0.0,Vf(iNode)+5.0,1e-5,tStick(iNode),-T(3,iNode),bc%Z(iNode),bc%rsf%f0(iNode),bc%rsf%V0(iNode),bc%rsf%a(iNode),bc%rsf%b(iNode),bc%rsf%L(iNode),bc%rsf%theta(iNode))
- tau2(iNode) = tStick(iNode) - bc%Z(iNode)*Vf2(iNode)
+ call rsf_update_state(0.5e0_CUSTOM_REAL*(Vf_old + Vf_new),bc%dt,bc%rsf)
+ do i=1,bc%nglob
+ Vf_new(i)=rtsafe(funcd,0.0,Vf_old(i)+5.0,1e-5,tStick(i),-T(3,i),bc%Z(i),bc%rsf%f0(i), &
+ bc%rsf%V0(i),bc%rsf%a(i),bc%rsf%b(i),bc%rsf%L(i),bc%rsf%theta(i))
enddo
- tStick = max(tStick,1e0_CUSTOM_REAL)
- T(1,:) = tau2 * T(1,:)/tStick
- T(2,:) = tau2 * T(2,:)/tStick
-
+ tnew = tStick - bc%Z*Vf_new
+
endif
+ tStick = max(tStick,1e0_CUSTOM_REAL) ! to avoid division by zero
+ T(1,:) = tnew * T(1,:)/tStick
+ T(2,:) = tnew * T(2,:)/tStick
+
! Save total tractions
bc%T = T
! Subtract initial stress
T = T - bc%T0
- if(RATE_AND_STATE) T(1,:) = T(1,:) - TxExt
+ if (RATE_AND_STATE) T(1,:) = T(1,:) - TxExt
+ !JPA: this eliminates the effect of TxExt on the equations of motion. Why is it needed?
+
! 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)
@@ -946,7 +939,7 @@
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,:))
+ Vf_new = sqrt(bc%V(1,:)*bc%V(1,:)+bc%V(2,:)*bc%V(2,:))
if(.not. RATE_AND_STATE) then
theta_new = bc%swf%theta
dc = bc%swf%dc
@@ -956,7 +949,7 @@
endif
call store_dataXZ(bc%dataXZ, strength, theta_old, theta_new, dc, &
- Vnorm_old, Vnorm, it*bc%dt,bc%dt)
+ Vf_old, Vf_new, it*bc%dt,bc%dt)
! call store_dataT(bc%dataT,bc%D,bc%V,bc%T,theta_new,it)
!-- outputs --
More information about the CIG-COMMITS
mailing list