[cig-commits] r19352 - seismo/3D/FAULT_SOURCE/branches/new_fault_db/src
surendra at geodynamics.org
surendra at geodynamics.org
Thu Jan 12 16:20:31 PST 2012
Author: surendra
Date: 2012-01-12 16:20:31 -0800 (Thu, 12 Jan 2012)
New Revision: 19352
Modified:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
Log:
Final update of solver for TPV16/17 : properly benchmarked
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-01-12 23:31:23 UTC (rev 19351)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90 2012-01-13 00:20:31 UTC (rev 19352)
@@ -59,7 +59,7 @@
private
integer :: kind
logical :: healing = .false.
- real(kind=CUSTOM_REAL), dimension(:), pointer :: Dc=>null(), mus=>null(), mud=>null(), theta=>null(), T=>null()
+ real(kind=CUSTOM_REAL), dimension(:), pointer :: Dc=>null(), mus=>null(), mud=>null(), theta=>null(), T=>null(), C=>null()
end type swf_type
type bc_dynflt_type
@@ -228,8 +228,10 @@
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, dimension(:), allocatable :: iLoc
+ integer :: iLoc
+
NAMELIST / INIT_STRESS / S1,S2,S3,n1,n2,n3
NAMELIST / SWF / mus,mud,dc,nmus,nmud,ndc
@@ -323,6 +325,22 @@
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) )
@@ -330,22 +348,45 @@
allocate( bc%swf%Dc(bc%nglob) )
allocate( bc%swf%theta(bc%nglob) )
allocate( bc%swf%T(bc%nglob) )
+ allocate( bc%swf%C(bc%nglob) )
- allocate(iLoc(bc%nglob))
- do ipar=1,bc%nglob
- iLoc = minloc( (loc_str(ipar)-bc%coord(1,:))**2 + (-loc_dip(ipar)-bc%coord(3,:))**2 ) !iloc_dip is negative of Z-coord
-print*,minval(iLoc)
- bc%T0(1,iLoc(1)) = sigma0(ipar)
- bc%T0(2,iLoc(1)) = tau0_str(ipar)
- bc%T0(3,iLoc(1)) = tau0_dip(ipar)
+! 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(1)) = static_fc(ipar)
- bc%swf%mud(iLoc(1)) = dyn_fc(ipar)
- bc%swf%Dc(iLoc(1)) = swcd(ipar)
- bc%swf%T(iLoc(1)) = tim_forcedRup(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
+ 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
+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)
+
! WARNING: if V_HEALING is negative we turn off healing
bc%swf%healing = (V_HEALING > 0e0_CUSTOM_REAL)
@@ -356,7 +397,15 @@
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
!---------------------------------------------------------------------
@@ -695,7 +744,7 @@
! 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,:))
@@ -960,7 +1009,6 @@
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
@@ -1112,7 +1160,7 @@
character :: NTchar*5
- integer*4 today(3), now(3)
+ integer :: today(3), now(3)
call idate(today) ! today(1)=day, (2)=month, (3)=year
call itime(now) ! now(1)=hour, (2)=minute, (3)=second
@@ -1129,10 +1177,10 @@
open(IOUT,file='OUTPUT_FILES/'//trim(dataT%name(i))//'.dat',status='replace')
write(IOUT,*) "# problem=TPV16"
write(IOUT,*) "# author=Surendra Nadh Somala"
- write(IOUT,1000) "# date=",today(2), today(1), today(3), now
+ write(IOUT,1000) today(2), today(1), today(3), now
write(IOUT,*) "# code=SPECFEM3D_SESAME (split nodes)"
write(IOUT,*) "# code_version=1.1"
- write(IOUT,*) "# element_size=75 m (*3 GLL nodes)"
+ write(IOUT,*) "# element_size=75 m (*5 GLL nodes)"
write(IOUT,*) "# time_step=",DT
write(IOUT,*) "# location=",trim(dataT%name(i))
write(IOUT,*) "# Column #1 = Time (s)"
@@ -1145,18 +1193,17 @@
write(IOUT,*) "# Column #8 = normal stress (MPa)"
write(IOUT,*) "#"
write(IOUT,*) "# The line below lists the names of the data fields:"
- write(IOUT,*) "#t h-slip h-slip-rate h-shear-stress v-slip v-slip-rate v-shear-stress n-stress"
+ write(IOUT,*) "t h-slip h-slip-rate h-shear-stress v-slip v-slip-rate v-shear-stress n-stress"
write(IOUT,*) "#"
do k=1,NT
write(IOUT,'(8(E15.7))') k*DT, dataT%d1(k,i), dataT%v1(k,i), dataT%t1(k,i)/1.0e6_CUSTOM_REAL, &
- dataT%d2(k,i), dataT%v2(k,i), dataT%t2(k,i)/1.0e6_CUSTOM_REAL, &
+ -dataT%d2(k,i), -dataT%v2(k,i), -dataT%t2(k,i)/1.0e6_CUSTOM_REAL, &
dataT%t3(k,i)/1.0e6_CUSTOM_REAL
enddo
close(IOUT)
- 1000 format ( 'Date ', i2.2, '/', i2.2, '/', i4.4, '; time ',
- & i2.2, ':', i2.2, ':', i2.2 )
+ 1000 format ( ' # Date = ', i2.2, '/', i2.2, '/', i4.4, '; time = ',i2.2, ':', i2.2, ':', i2.2 )
enddo
@@ -1187,24 +1234,25 @@
open(IOUT,file=trim(filename),status='replace')
write(IOUT,*) "# problem=TPV16"
write(IOUT,*) "# author=Surendra Nadh Somala"
- write(IOUT,1000) "# date=",today(2), today(1), today(3), now
+ write(IOUT,1000) today(2), today(1), today(3), now
write(IOUT,*) "# code=SPECFEM3D_SESAME (split nodes)"
write(IOUT,*) "# code_version=1.1"
- write(IOUT,*) "# element_size=75 m (*3 GLL nodes)"
+ write(IOUT,*) "# element_size=75 m (*5 GLL nodes)"
write(IOUT,*) "# Column #1 = horizontal coordinate, distance along strike (m)"
write(IOUT,*) "# Column #2 = vertical coordinate, distance down-dip (m)"
write(IOUT,*) "# Column #3 = rupture time (s)"
- write(IOUT,*) "# x y z time"
+ write(IOUT,*) "# "
+ write(IOUT,*) "j k t"
do i = 1,size(dataXZ%tRUP)
- write(IOUT,'(4(E15.7))') dataXZ%xcoord(i), dataXZ%ycoord(i), dataXZ%zcoord(i), dataXZ%tRUP(i)
+ write(IOUT,'(3(E15.7))') dataXZ%xcoord(i), -dataXZ%zcoord(i), dataXZ%tRUP(i)
end do
close(IOUT)
- 1000 format ( 'Date ', i2.2, '/', i2.2, '/', i4.4, '; time ',
- & i2.2, ':', i2.2, ':', i2.2 )
+ 1000 format ( ' # Date = ', i2.2, '/', i2.2, '/', i4.4, '; time = ',i2.2, ':', i2.2, ':', i2.2 )
+
end subroutine SCEC_Write_RuptureTime
!-------------------------------------------------------------------------------------------------
More information about the CIG-COMMITS
mailing list