[cig-commits] r20907 - seismo/3D/FAULT_SOURCE/branches/new_fault_db/src
surendra at geodynamics.org
surendra at geodynamics.org
Wed Oct 24 11:35:57 PDT 2012
Author: surendra
Date: 2012-10-24 11:35:56 -0700 (Wed, 24 Oct 2012)
New Revision: 20907
Modified:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_common.f90
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_kinematic.f90
Log:
More work on re-factoring fault routines
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-24 17:14:27 UTC (rev 20906)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90 2012-10-24 18:35:56 UTC (rev 20907)
@@ -37,7 +37,33 @@
private
+ ! outputs on selected fault nodes at every time step:
+ ! slip, slip velocity, fault stresses
+ type dataT_type
+ integer :: npoin=0,npoin_local=0
+ integer, dimension(:), pointer :: iglob=>null() ! on-fault global index of output nodes
+ integer, dimension(:,:), pointer :: iglob_all=>null()
+ integer, dimension(:), pointer :: islice=>null(),glob_indx=>null()
+ real(kind=CUSTOM_REAL), dimension(:), pointer :: dist=>null()
+ real(kind=CUSTOM_REAL), dimension(:,:), pointer :: dist_all=>null()
+ real(kind=CUSTOM_REAL), dimension(:,:), pointer :: d1=>null(),v1=>null(),t1=>null(), &
+ d2=>null(),v2=>null(),t2=>null(), &
+ t3=>null(),theta=>null()
+ character(len=70), dimension(:), pointer :: name=>null()
+ end type dataT_type
+
+ ! outputs(dyn) /inputs (kind) at selected times for all fault nodes:
+ ! strength, state, slip, slip velocity, fault stresses, rupture time, process zone time
+ ! rupture time = first time when slip velocity = threshold V_RUPT (defined below)
+ ! process zone time = first time when slip = Dc
+ type dataXZ_type
+ real(kind=CUSTOM_REAL), dimension(:), pointer :: stg=>null(), sta=>null(), d1=>null(), d2=>null(), v1=>null(), v2=>null(), &
+ t1=>null(), t2=>null(), t3=>null(), tRUP=>null(), tPZ=>null()
+ real(kind=CUSTOM_REAL), dimension(:), pointer :: xcoord=>null(), ycoord=>null(), zcoord=>null()
+ integer :: npoin=0
+ end type dataXZ_type
+
type swf_type
private
integer :: kind
@@ -63,6 +89,10 @@
type, EXTENDS (fault_type) :: bc_dynflt_type
private
+ real(kind=CUSTOM_REAL), dimension(:,:), pointer :: T0=>null()
+ real(kind=CUSTOM_REAL), dimension(:), pointer :: MU=>null()
+ type(dataT_type) :: dataT
+ type(dataXZ_type) :: dataXZ,dataXZ_all
type(swf_type), pointer :: swf => null()
type(rsf_type), pointer :: rsf => null()
type(asperity_type), pointer :: asp => null()
@@ -78,6 +108,9 @@
!slip velocity threshold for definition of rupture front
real(kind=CUSTOM_REAL), save :: V_RUPT
+ !Number of time steps defined by the user : NTOUT
+ integer, save :: NTOUT,NSNAP
+
logical, save :: SIMULATION_TYPE_DYN = .false.
logical, save :: TPV16 = .false.
@@ -244,99 +277,10 @@
NAMELIST / RSF / V0,f0,a,b,L,V_init,theta_init,nV0,nf0,na,nb,nL,nV_init,ntheta_init,C,T,nC,nForcedRup,Vw,fw,nVw,nfw
NAMELIST / ASP / Fload,nFload
- read(IIN_BIN) bc%nspec,bc%nglob
- if (.NOT.PARALLEL_FAULT .and. bc%nspec==0) return
- if (bc%nspec>0) then
- 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_tmp
+ call initialize_simulation(bc,IIN_BIN,Minv,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
- endif
-
- if (PARALLEL_FAULT) then
-
- accel=0._CUSTOM_REAL
- 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)
-
- accel=0._CUSTOM_REAL
- if (bc%nspec>0) then
- 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, &
- 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) then
- nx = accel(1,bc%ibulk1)
- ny = accel(2,bc%ibulk1)
- nz = accel(3,bc%ibulk1)
- endif
-
- endif
-
if (bc%nspec>0) then
- 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 ))
- ! WARNING: In non-split nodes at fault edges M is assembled across the fault.
- ! hence invM1+invM2=2/(M1+M2) instead of 1/M1+1/M2
- ! In a symmetric mesh (M1=M2) Z will be twice its intended value
-
allocate(bc%T(3,bc%nglob))
allocate(bc%D(3,bc%nglob))
allocate(bc%V(3,bc%nglob))
@@ -619,6 +563,92 @@
end subroutine init_one_fault
+!---------------------------------------------------------------------
+! adds a value to a fault parameter inside an area with prescribed shape
+subroutine init_2d_distribution(a,coord,iin,n)
+!JPA refactor: background value shuld be an argument
+
+ real(kind=CUSTOM_REAL), intent(inout) :: a(:)
+ real(kind=CUSTOM_REAL), intent(in) :: coord(:,:)
+ integer, intent(in) :: iin,n
+
+ real(kind=CUSTOM_REAL) :: b(size(a))
+ character(len=20) :: shape
+ real(kind=CUSTOM_REAL) :: val,valh, xc, yc, zc, r, l, lx,ly,lz
+ real(kind=CUSTOM_REAL) :: r1(size(a))
+ integer :: i
+
+ NAMELIST / DIST2D / shape, val,valh, xc, yc, zc, r, l, lx,ly,lz
+
+ if (n==0) return
+
+ do i=1,n
+ shape = ''
+ xc = 0e0_CUSTOM_REAL
+ yc = 0e0_CUSTOM_REAL
+ zc = 0e0_CUSTOM_REAL
+ r = 0e0_CUSTOM_REAL
+ l = 0e0_CUSTOM_REAL
+ lx = 0e0_CUSTOM_REAL
+ ly = 0e0_CUSTOM_REAL
+ lz = 0e0_CUSTOM_REAL
+ valh = 0e0_CUSTOM_REAL
+
+ read(iin,DIST2D)
+ select case(shape)
+ case ('circle')
+ b = heaviside( r - sqrt((coord(1,:)-xc)**2 + (coord(2,:)-yc)**2 + (coord(3,:)-zc)**2 ) )
+ case ('circle-exp')
+ r1 = sqrt((coord(1,:)-xc)**2 + (coord(2,:)-yc)**2 + (coord(3,:)-zc)**2 )
+ where(r1<r)
+ b =exp(r1**2/(r1**2 - r**2) )
+ elsewhere
+ b =0._CUSTOM_REAL
+ endwhere
+ case ('ellipse')
+ b = heaviside( 1e0_CUSTOM_REAL - sqrt( (coord(1,:)-xc)**2/lx**2 + (coord(2,:)-yc)**2/ly**2 + (coord(3,:)-zc)**2/lz**2 ) )
+ case ('square')
+ b = heaviside((l/2._CUSTOM_REAL)-abs(coord(1,:)-xc)+SMALLVAL) * &
+ heaviside((l/2._CUSTOM_REAL)-abs(coord(2,:)-yc)+SMALLVAL) * &
+ heaviside((l/2._CUSTOM_REAL)-abs(coord(3,:)-zc)+SMALLVAL) * &
+ val
+ case ('cylinder')
+ b = heaviside(r - sqrt((coord(1,:)-xc)**2 + (coord(2,:)-yc)**2)) * &
+ heaviside((lz/2._CUSTOM_REAL)-abs(coord(3,:)-zc)+SMALLVAL) * &
+ val
+ case ('rectangle')
+ b = heaviside((lx/2._CUSTOM_REAL)-abs(coord(1,:)-xc)+SMALLVAL) * &
+ heaviside((ly/2._CUSTOM_REAL)-abs(coord(2,:)-yc)+SMALLVAL) * &
+ heaviside((lz/2._CUSTOM_REAL)-abs(coord(3,:)-zc)+SMALLVAL) * &
+ val
+ case ('rectangle-taper')
+ b = heaviside((lx/2._CUSTOM_REAL)-abs(coord(1,:)-xc)+SMALLVAL) * &
+ heaviside((ly/2._CUSTOM_REAL)-abs(coord(2,:)-yc)+SMALLVAL) * &
+ heaviside((lz/2._CUSTOM_REAL)-abs(coord(3,:)-zc)+SMALLVAL) * &
+ (val + ( coord(3,:) - zc + lz/2._CUSTOM_REAL ) * ((valh-val)/lz))
+ case default
+ stop 'bc_dynflt_3d::init_2d_distribution:: unknown shape'
+ end select
+
+ where (b /= 0e0_CUSTOM_REAL) a = b
+ enddo
+
+end subroutine init_2d_distribution
+
+!---------------------------------------------------------------------
+elemental function heaviside(x)
+
+ real(kind=CUSTOM_REAL), intent(in) :: x
+ real(kind=CUSTOM_REAL) :: heaviside
+
+ if (x>=0e0_CUSTOM_REAL) then
+ heaviside = 1e0_CUSTOM_REAL
+ else
+ heaviside = 0e0_CUSTOM_REAL
+ endif
+
+end function heaviside
+
!=====================================================================
! adds boundary term Bt into Force array for each fault.
!
@@ -1011,6 +1041,255 @@
end subroutine rsf_update_state
+
+!===============================================================
+! OUTPUTS
+subroutine init_dataT(DataT,coord,nglob,NT,iflt)
+
+ use specfem_par, only : NPROC,myrank
+ ! NT = total number of time steps
+
+ integer, intent(in) :: nglob,NT,iflt
+ real(kind=CUSTOM_REAL), intent(in) :: coord(3,nglob)
+ type (dataT_type), intent(out) :: DataT
+
+ real(kind=CUSTOM_REAL) :: xtarget,ytarget,ztarget,dist,distkeep
+ integer :: i, iglob , IIN, ier, jflt, np, k
+ character(len=70) :: tmpname
+
+ integer :: ipoin, ipoin_local, iproc
+
+ ! 1. read fault output coordinates from user file,
+ ! 2. define iglob: the fault global index of the node nearest to user
+ ! requested coordinate
+
+ IIN = 251
+ open(IIN,file='DATA/FAULT/FAULT_STATIONS',status='old',action='read',iostat=ier)
+ read(IIN,*) np
+ DataT%npoin =0
+ do i=1,np
+ read(IIN,*) xtarget,ytarget,ztarget,tmpname,jflt
+ if (jflt==iflt) DataT%npoin = DataT%npoin +1
+ enddo
+ close(IIN)
+
+ if (DataT%npoin == 0) return
+
+ allocate(DataT%iglob(DataT%npoin))
+ allocate(DataT%name(DataT%npoin))
+ allocate(DataT%dist(DataT%npoin)) !Surendra : for parallel fault
+
+ open(IIN,file='DATA/FAULT/FAULT_STATIONS',status='old',action='read',iostat=ier)
+ if( ier /= 0 ) stop 'error opening FAULT_STATIONS file'
+ read(IIN,*) np
+ k = 0
+ do i=1,np
+ read(IIN,*) xtarget,ytarget,ztarget,tmpname,jflt
+ if (jflt/=iflt) cycle
+ k = k+1
+ DataT%name(k) = tmpname
+ !search nearest node
+ distkeep = huge(distkeep)
+
+ do iglob=1,nglob
+ dist = sqrt((coord(1,iglob)-xtarget)**2 &
+ + (coord(2,iglob)-ytarget)**2 &
+ + (coord(3,iglob)-ztarget)**2)
+ if (dist < distkeep) then
+ distkeep = dist
+ DataT%iglob(k) = iglob
+ endif
+ enddo
+ DataT%dist(k) = distkeep !Surendra : for parallel fault
+ enddo
+
+ !Surendra : for parallel fault
+ if (PARALLEL_FAULT) then
+ allocate(DataT%islice(DataT%npoin))
+ allocate(DataT%iglob_all(DataT%npoin,0:NPROC-1))
+ allocate(DataT%dist_all(DataT%npoin,0:NPROC-1))
+ call gather_all_i(DataT%iglob,DataT%npoin,DataT%iglob_all,DataT%npoin,NPROC)
+ call gather_all_cr(DataT%dist,DataT%npoin,DataT%dist_all,DataT%npoin,NPROC)
+ if(myrank==0) then
+ do ipoin = 1,DataT%npoin
+ iproc = minloc(DataT%dist_all(ipoin,:), 1) - 1
+ DataT%islice(ipoin) = iproc
+ DataT%iglob(ipoin) = DataT%iglob_all(ipoin,iproc)
+ enddo
+ endif
+
+ call bcast_all_i(DataT%islice,DataT%npoin)
+ call bcast_all_i(DataT%iglob,DataT%npoin)
+
+ DataT%npoin_local = 0
+ do ipoin = 1,DataT%npoin
+ if(myrank == DataT%islice(ipoin)) DataT%npoin_local = DataT%npoin_local + 1
+ enddo
+ allocate(DataT%glob_indx(DataT%npoin_local))
+ do ipoin = 1,DataT%npoin
+ if(myrank == DataT%islice(ipoin)) then
+ ipoin_local = ipoin_local + 1
+ DataT%glob_indx(ipoin_local) = ipoin
+ endif
+ enddo
+ else
+ DataT%npoin_local = DataT%npoin
+ endif !Parallel_fault
+
+
+ ! 3. allocate arrays and set to zero
+ allocate(DataT%d1(NT,DataT%npoin))
+ allocate(DataT%v1(NT,DataT%npoin))
+ allocate(DataT%t1(NT,DataT%npoin))
+ allocate(DataT%d2(NT,DataT%npoin))
+ allocate(DataT%v2(NT,DataT%npoin))
+ allocate(DataT%t2(NT,DataT%npoin))
+ allocate(DataT%t3(NT,DataT%npoin))
+ allocate(DataT%theta(NT,DataT%npoin))
+ DataT%d1 = 0e0_CUSTOM_REAL
+ DataT%v1 = 0e0_CUSTOM_REAL
+ DataT%t1 = 0e0_CUSTOM_REAL
+ DataT%d2 = 0e0_CUSTOM_REAL
+ DataT%v2 = 0e0_CUSTOM_REAL
+ DataT%t2 = 0e0_CUSTOM_REAL
+ DataT%t3 = 0e0_CUSTOM_REAL
+ DataT%theta = 0e0_CUSTOM_REAL
+
+ close(IIN)
+
+end subroutine init_dataT
+
+!---------------------------------------------------------------
+subroutine store_dataT(dataT,d,v,t,theta,itime)
+
+ type(dataT_type), intent(inout) :: dataT
+ real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: d,v,t
+ real(kind=CUSTOM_REAL), dimension(:), intent(in) :: theta
+ integer, intent(in) :: itime
+
+ integer :: i,k,ipoin
+
+ do ipoin=1,dataT%npoin_local
+ if (PARALLEL_FAULT) then
+ i = DataT%glob_indx(ipoin)
+ else
+ i = ipoin
+ endif
+ k = dataT%iglob(i)
+ dataT%d1(itime,i) = d(1,k)
+ dataT%d2(itime,i) = d(2,k)
+ dataT%v1(itime,i) = v(1,k)
+ dataT%v2(itime,i) = v(2,k)
+ dataT%t1(itime,i) = t(1,k)
+ dataT%t2(itime,i) = t(2,k)
+ dataT%t3(itime,i) = t(3,k)
+ dataT%theta(itime,i) = theta(k)
+ enddo
+
+end subroutine store_dataT
+
+!-----------------------------------------------------------------
+subroutine write_dataT_all(nt,statelaw)
+
+ integer, intent(in) :: nt
+ integer, intent(in) :: statelaw
+
+ integer :: i
+
+ if (.not.allocated(faults)) return
+ do i = 1,size(faults)
+ call SCEC_write_dataT(faults(i)%dataT,faults(i)%dt,nt,statelaw)
+ enddo
+
+end subroutine write_dataT_all
+
+!------------------------------------------------------------------------
+subroutine SCEC_write_dataT(dataT,DT,NT,statelaw)
+
+ type(dataT_type), intent(in) :: dataT
+ real(kind=CUSTOM_REAL), intent(in) :: DT
+ integer, intent(in) :: NT
+ integer, intent(in) :: statelaw
+
+ integer :: i,k,IOUT,ipoin
+ character :: NTchar*5
+ 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
+
+ IOUT = 121 !WARNING: not very robust. Could instead look for an available ID
+
+ write(NTchar,1) NT
+ NTchar = adjustl(NTchar)
+
+1 format(I5)
+ do ipoin=1,dataT%npoin_local
+ if (PARALLEL_FAULT) then
+ i = DataT%glob_indx(ipoin)
+ else
+ i = ipoin
+ endif
+
+ open(IOUT,file='OUTPUT_FILES/'//trim(dataT%name(i))//'.dat',status='replace')
+ write(IOUT,*) "# problem=TPV104"
+ write(IOUT,*) "# author=Surendra Nadh Somala"
+ 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=100 m (*5 GLL nodes)"
+ write(IOUT,*) "# time_step=",DT
+ write(IOUT,*) "# location=",trim(dataT%name(i))
+ write(IOUT,*) "# Column #1 = Time (s)"
+ write(IOUT,*) "# Column #2 = horizontal right-lateral slip (m)"
+ write(IOUT,*) "# Column #3 = horizontal right-lateral slip rate (m/s)"
+ write(IOUT,*) "# Column #4 = horizontal right-lateral shear stress (MPa)"
+ write(IOUT,*) "# Column #5 = vertical up-dip slip (m)"
+ write(IOUT,*) "# Column #6 = vertical up-dip slip rate (m/s)"
+ write(IOUT,*) "# Column #7 = vertical up-dip shear stress (MPa)"
+ write(IOUT,*) "# Column #8 = normal stress (MPa)"
+ if (RATE_AND_STATE) write(IOUT,*) "# Column #9 = log10 of state variable (log-seconds)"
+ write(IOUT,*) "#"
+ write(IOUT,*) "# The line below lists the names of the data fields:"
+ if (.not. RATE_AND_STATE) then
+ 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%t3(k,i)/1.0e6_CUSTOM_REAL
+ enddo
+ else
+ if(statelaw == 1) then
+ write(IOUT,*) "t h-slip h-slip-rate h-shear-stress v-slip v-slip-rate v-shear-stress n-stress log-theta"
+ write(IOUT,*) "#"
+ do k=1,NT
+ write(IOUT,'(9(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%t3(k,i)/1.0e6_CUSTOM_REAL, log10(dataT%theta(k,i))
+ enddo
+ else
+ write(IOUT,*) "t h-slip h-slip-rate h-shear-stress v-slip v-slip-rate v-shear-stress n-stress psi"
+ write(IOUT,*) "#"
+ do k=1,NT
+ write(IOUT,'(9(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%t3(k,i)/1.0e6_CUSTOM_REAL, dataT%theta(k,i)
+ enddo
+ endif
+ endif
+ close(IOUT)
+
+
+1000 format ( ' # Date = ', i2.2, '/', i2.2, '/', i4.4, '; time = ',i2.2, ':', i2.2, ':', i2.2 )
+
+
+ enddo
+
+end subroutine SCEC_write_dataT
+
+!-------------------------------------------------------------------------------------------------
+
subroutine SCEC_Write_RuptureTime(dataXZ,DT,NT,iflt)
type(dataXZ_type), intent(in) :: dataXZ
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_common.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_common.f90 2012-10-24 17:14:27 UTC (rev 20906)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_common.f90 2012-10-24 18:35:56 UTC (rev 20907)
@@ -7,53 +7,20 @@
implicit none
private
- ! outputs on selected fault nodes at every time step:
- ! slip, slip velocity, fault stresses
- type dataT_type
- integer :: npoin=0,npoin_local=0
- integer, dimension(:), pointer :: iglob=>null() ! on-fault global index of output nodes
- integer, dimension(:,:), pointer :: iglob_all=>null()
- integer, dimension(:), pointer :: islice=>null(),glob_indx=>null()
- real(kind=CUSTOM_REAL), dimension(:), pointer :: dist=>null()
- real(kind=CUSTOM_REAL), dimension(:,:), pointer :: dist_all=>null()
- real(kind=CUSTOM_REAL), dimension(:,:), pointer :: d1=>null(),v1=>null(),t1=>null(), &
- d2=>null(),v2=>null(),t2=>null(), &
- t3=>null(),theta=>null()
- character(len=70), dimension(:), pointer :: name=>null()
- end type dataT_type
-
-
- ! outputs(dyn) /inputs (kind) at selected times for all fault nodes:
- ! strength, state, slip, slip velocity, fault stresses, rupture time, process zone time
- ! rupture time = first time when slip velocity = threshold V_RUPT (defined below)
- ! process zone time = first time when slip = Dc
- type dataXZ_type
- real(kind=CUSTOM_REAL), dimension(:), pointer :: stg=>null(), sta=>null(), d1=>null(), d2=>null(), v1=>null(), v2=>null(), &
- t1=>null(), t2=>null(), t3=>null(), tRUP=>null(), tPZ=>null()
- real(kind=CUSTOM_REAL), dimension(:), pointer :: xcoord=>null(), ycoord=>null(), zcoord=>null()
- integer :: npoin=0
- end type dataXZ_type
-
-
type fault_type
private
integer :: nspec=0, nglob=0
- real(kind=CUSTOM_REAL), dimension(:,:), pointer :: T0=>null(),T=>null(),V=>null(),D=>null()
+ real(kind=CUSTOM_REAL), dimension(:,:), pointer :: T=>null(),V=>null(),D=>null()
real(kind=CUSTOM_REAL), dimension(:,:), pointer :: coord=>null()
real(kind=CUSTOM_REAL), dimension(:,:,:), pointer :: R=>null()
- real(kind=CUSTOM_REAL), dimension(:), pointer :: MU=>null(),B=>null(),invM1=>null(),invM2=>null(),Z=>null()
+ real(kind=CUSTOM_REAL), dimension(:), pointer :: B=>null(),invM1=>null(),invM2=>null(),Z=>null()
real(kind=CUSTOM_REAL) :: dt
integer, dimension(:), pointer :: ibulk1=>null(), ibulk2=>null()
- type(dataT_type) :: dataT
- type(dataXZ_type) :: dataXZ,dataXZ_all
integer, dimension(:), pointer :: poin_offset=>null(),npoin_perproc=>null()
end type fault_type
-!Number of time steps defined by the user : NTOUT
- integer, save :: NTOUT,NSNAP
-
logical, save :: PARALLEL_FAULT = .false.
public :: initialize_fault,PARALLEL_FAULT
@@ -66,7 +33,7 @@
use specfem_par
- type(bc_dynflt_type), intent(inout) :: bc
+ type(fault_type), intent(inout) :: bc
real(kind=CUSTOM_REAL), intent(in) :: Minv(:)
integer, intent(in) :: IIN_BIN
real(kind=CUSTOM_REAL), intent(in) :: dt
@@ -175,6 +142,8 @@
! hence invM1+invM2=2/(M1+M2) instead of 1/M1+1/M2
! In a symmetric mesh (M1=M2) Z will be twice its intended value
+ endif
+
end subroutine initialize
!---------------------------------------------------------------------
@@ -216,7 +185,7 @@
!===============================================================
function get_jump (bc,v) result(dv)
- type(bc_dynflt_type), intent(in) :: bc
+ type(fault_type), intent(in) :: bc
real(kind=CUSTOM_REAL), intent(in) :: v(:,:)
real(kind=CUSTOM_REAL) :: dv(3,bc%nglob)
@@ -230,7 +199,7 @@
!---------------------------------------------------------------------
function get_weighted_jump (bc,f) result(da)
- type(bc_dynflt_type), intent(in) :: bc
+ type(fault_type), intent(in) :: bc
real(kind=CUSTOM_REAL), intent(in) :: f(:,:)
real(kind=CUSTOM_REAL) :: da(3,bc%nglob)
@@ -248,7 +217,7 @@
!----------------------------------------------------------------------
function rotate(bc,v,fb) result(vr)
- type(bc_dynflt_type), intent(in) :: bc
+ type(fault_type), intent(in) :: bc
real(kind=CUSTOM_REAL), intent(in) :: v(3,bc%nglob)
integer, intent(in) :: fb
real(kind=CUSTOM_REAL) :: vr(3,bc%nglob)
@@ -271,340 +240,7 @@
end function rotate
-!---------------------------------------------------------------------
-! adds a value to a fault parameter inside an area with prescribed shape
-subroutine init_2d_distribution(a,coord,iin,n)
-!JPA refactor: background value shuld be an argument
+!----------------------------------------------------------------------
- real(kind=CUSTOM_REAL), intent(inout) :: a(:)
- real(kind=CUSTOM_REAL), intent(in) :: coord(:,:)
- integer, intent(in) :: iin,n
- real(kind=CUSTOM_REAL) :: b(size(a))
- character(len=20) :: shape
- real(kind=CUSTOM_REAL) :: val,valh, xc, yc, zc, r, l, lx,ly,lz
- real(kind=CUSTOM_REAL) :: r1(size(a))
- integer :: i
-
- NAMELIST / DIST2D / shape, val,valh, xc, yc, zc, r, l, lx,ly,lz
-
- if (n==0) return
-
- do i=1,n
- shape = ''
- xc = 0e0_CUSTOM_REAL
- yc = 0e0_CUSTOM_REAL
- zc = 0e0_CUSTOM_REAL
- r = 0e0_CUSTOM_REAL
- l = 0e0_CUSTOM_REAL
- lx = 0e0_CUSTOM_REAL
- ly = 0e0_CUSTOM_REAL
- lz = 0e0_CUSTOM_REAL
- valh = 0e0_CUSTOM_REAL
-
- read(iin,DIST2D)
- select case(shape)
- case ('circle')
- b = heaviside( r - sqrt((coord(1,:)-xc)**2 + (coord(2,:)-yc)**2 + (coord(3,:)-zc)**2 ) )
- case ('circle-exp')
- r1 = sqrt((coord(1,:)-xc)**2 + (coord(2,:)-yc)**2 + (coord(3,:)-zc)**2 )
- where(r1<r)
- b =exp(r1**2/(r1**2 - r**2) )
- elsewhere
- b =0._CUSTOM_REAL
- endwhere
- case ('ellipse')
- b = heaviside( 1e0_CUSTOM_REAL - sqrt( (coord(1,:)-xc)**2/lx**2 + (coord(2,:)-yc)**2/ly**2 + (coord(3,:)-zc)**2/lz**2 ) )
- case ('square')
- b = heaviside((l/2._CUSTOM_REAL)-abs(coord(1,:)-xc)+SMALLVAL) * &
- heaviside((l/2._CUSTOM_REAL)-abs(coord(2,:)-yc)+SMALLVAL) * &
- heaviside((l/2._CUSTOM_REAL)-abs(coord(3,:)-zc)+SMALLVAL) * &
- val
- case ('cylinder')
- b = heaviside(r - sqrt((coord(1,:)-xc)**2 + (coord(2,:)-yc)**2)) * &
- heaviside((lz/2._CUSTOM_REAL)-abs(coord(3,:)-zc)+SMALLVAL) * &
- val
- case ('rectangle')
- b = heaviside((lx/2._CUSTOM_REAL)-abs(coord(1,:)-xc)+SMALLVAL) * &
- heaviside((ly/2._CUSTOM_REAL)-abs(coord(2,:)-yc)+SMALLVAL) * &
- heaviside((lz/2._CUSTOM_REAL)-abs(coord(3,:)-zc)+SMALLVAL) * &
- val
- case ('rectangle-taper')
- b = heaviside((lx/2._CUSTOM_REAL)-abs(coord(1,:)-xc)+SMALLVAL) * &
- heaviside((ly/2._CUSTOM_REAL)-abs(coord(2,:)-yc)+SMALLVAL) * &
- heaviside((lz/2._CUSTOM_REAL)-abs(coord(3,:)-zc)+SMALLVAL) * &
- (val + ( coord(3,:) - zc + lz/2._CUSTOM_REAL ) * ((valh-val)/lz))
- case default
- stop 'bc_dynflt_3d::init_2d_distribution:: unknown shape'
- end select
-
- where (b /= 0e0_CUSTOM_REAL) a = b
- enddo
-
-end subroutine init_2d_distribution
-
-!---------------------------------------------------------------------
-elemental function heaviside(x)
-
- real(kind=CUSTOM_REAL), intent(in) :: x
- real(kind=CUSTOM_REAL) :: heaviside
-
- if (x>=0e0_CUSTOM_REAL) then
- heaviside = 1e0_CUSTOM_REAL
- else
- heaviside = 0e0_CUSTOM_REAL
- endif
-
-end function heaviside
-
-
-
-!===============================================================
-! OUTPUTS
-subroutine init_dataT(DataT,coord,nglob,NT,iflt)
-
- use specfem_par, only : NPROC,myrank
- ! NT = total number of time steps
-
- integer, intent(in) :: nglob,NT,iflt
- real(kind=CUSTOM_REAL), intent(in) :: coord(3,nglob)
- type (dataT_type), intent(out) :: DataT
-
- real(kind=CUSTOM_REAL) :: xtarget,ytarget,ztarget,dist,distkeep
- integer :: i, iglob , IIN, ier, jflt, np, k
- character(len=70) :: tmpname
-
- integer :: ipoin, ipoin_local, iproc
-
- ! 1. read fault output coordinates from user file,
- ! 2. define iglob: the fault global index of the node nearest to user
- ! requested coordinate
-
- IIN = 251
- open(IIN,file='DATA/FAULT/FAULT_STATIONS',status='old',action='read',iostat=ier)
- read(IIN,*) np
- DataT%npoin =0
- do i=1,np
- read(IIN,*) xtarget,ytarget,ztarget,tmpname,jflt
- if (jflt==iflt) DataT%npoin = DataT%npoin +1
- enddo
- close(IIN)
-
- if (DataT%npoin == 0) return
-
- allocate(DataT%iglob(DataT%npoin))
- allocate(DataT%name(DataT%npoin))
- allocate(DataT%dist(DataT%npoin)) !Surendra : for parallel fault
-
- open(IIN,file='DATA/FAULT/FAULT_STATIONS',status='old',action='read',iostat=ier)
- if( ier /= 0 ) stop 'error opening FAULT_STATIONS file'
- read(IIN,*) np
- k = 0
- do i=1,np
- read(IIN,*) xtarget,ytarget,ztarget,tmpname,jflt
- if (jflt/=iflt) cycle
- k = k+1
- DataT%name(k) = tmpname
- !search nearest node
- distkeep = huge(distkeep)
-
- do iglob=1,nglob
- dist = sqrt((coord(1,iglob)-xtarget)**2 &
- + (coord(2,iglob)-ytarget)**2 &
- + (coord(3,iglob)-ztarget)**2)
- if (dist < distkeep) then
- distkeep = dist
- DataT%iglob(k) = iglob
- endif
- enddo
- DataT%dist(k) = distkeep !Surendra : for parallel fault
- enddo
-
- !Surendra : for parallel fault
- if (PARALLEL_FAULT) then
- allocate(DataT%islice(DataT%npoin))
- allocate(DataT%iglob_all(DataT%npoin,0:NPROC-1))
- allocate(DataT%dist_all(DataT%npoin,0:NPROC-1))
- call gather_all_i(DataT%iglob,DataT%npoin,DataT%iglob_all,DataT%npoin,NPROC)
- call gather_all_cr(DataT%dist,DataT%npoin,DataT%dist_all,DataT%npoin,NPROC)
- if(myrank==0) then
- do ipoin = 1,DataT%npoin
- iproc = minloc(DataT%dist_all(ipoin,:), 1) - 1
- DataT%islice(ipoin) = iproc
- DataT%iglob(ipoin) = DataT%iglob_all(ipoin,iproc)
- enddo
- endif
-
- call bcast_all_i(DataT%islice,DataT%npoin)
- call bcast_all_i(DataT%iglob,DataT%npoin)
-
- DataT%npoin_local = 0
- do ipoin = 1,DataT%npoin
- if(myrank == DataT%islice(ipoin)) DataT%npoin_local = DataT%npoin_local + 1
- enddo
- allocate(DataT%glob_indx(DataT%npoin_local))
- do ipoin = 1,DataT%npoin
- if(myrank == DataT%islice(ipoin)) then
- ipoin_local = ipoin_local + 1
- DataT%glob_indx(ipoin_local) = ipoin
- endif
- enddo
- else
- DataT%npoin_local = DataT%npoin
- endif !Parallel_fault
-
-
- ! 3. allocate arrays and set to zero
- allocate(DataT%d1(NT,DataT%npoin))
- allocate(DataT%v1(NT,DataT%npoin))
- allocate(DataT%t1(NT,DataT%npoin))
- allocate(DataT%d2(NT,DataT%npoin))
- allocate(DataT%v2(NT,DataT%npoin))
- allocate(DataT%t2(NT,DataT%npoin))
- allocate(DataT%t3(NT,DataT%npoin))
- allocate(DataT%theta(NT,DataT%npoin))
- DataT%d1 = 0e0_CUSTOM_REAL
- DataT%v1 = 0e0_CUSTOM_REAL
- DataT%t1 = 0e0_CUSTOM_REAL
- DataT%d2 = 0e0_CUSTOM_REAL
- DataT%v2 = 0e0_CUSTOM_REAL
- DataT%t2 = 0e0_CUSTOM_REAL
- DataT%t3 = 0e0_CUSTOM_REAL
- DataT%theta = 0e0_CUSTOM_REAL
-
- close(IIN)
-
-end subroutine init_dataT
-
-!---------------------------------------------------------------
-subroutine store_dataT(dataT,d,v,t,theta,itime)
-
- type(dataT_type), intent(inout) :: dataT
- real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: d,v,t
- real(kind=CUSTOM_REAL), dimension(:), intent(in) :: theta
- integer, intent(in) :: itime
-
- integer :: i,k,ipoin
-
- do ipoin=1,dataT%npoin_local
- if (PARALLEL_FAULT) then
- i = DataT%glob_indx(ipoin)
- else
- i = ipoin
- endif
- k = dataT%iglob(i)
- dataT%d1(itime,i) = d(1,k)
- dataT%d2(itime,i) = d(2,k)
- dataT%v1(itime,i) = v(1,k)
- dataT%v2(itime,i) = v(2,k)
- dataT%t1(itime,i) = t(1,k)
- dataT%t2(itime,i) = t(2,k)
- dataT%t3(itime,i) = t(3,k)
- dataT%theta(itime,i) = theta(k)
- enddo
-
-end subroutine store_dataT
-
-!-----------------------------------------------------------------
-subroutine write_dataT_all(nt,statelaw)
-
- integer, intent(in) :: nt
- integer, intent(in) :: statelaw
-
- integer :: i
-
- if (.not.allocated(faults)) return
- do i = 1,size(faults)
- call SCEC_write_dataT(faults(i)%dataT,faults(i)%dt,nt,statelaw)
- enddo
-
-end subroutine write_dataT_all
-
-!------------------------------------------------------------------------
-subroutine SCEC_write_dataT(dataT,DT,NT,statelaw)
-
- type(dataT_type), intent(in) :: dataT
- real(kind=CUSTOM_REAL), intent(in) :: DT
- integer, intent(in) :: NT
- integer, intent(in) :: statelaw
-
- integer :: i,k,IOUT,ipoin
- character :: NTchar*5
- 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
-
- IOUT = 121 !WARNING: not very robust. Could instead look for an available ID
-
- write(NTchar,1) NT
- NTchar = adjustl(NTchar)
-
-1 format(I5)
- do ipoin=1,dataT%npoin_local
- if (PARALLEL_FAULT) then
- i = DataT%glob_indx(ipoin)
- else
- i = ipoin
- endif
-
- open(IOUT,file='OUTPUT_FILES/'//trim(dataT%name(i))//'.dat',status='replace')
- write(IOUT,*) "# problem=TPV104"
- write(IOUT,*) "# author=Surendra Nadh Somala"
- 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=100 m (*5 GLL nodes)"
- write(IOUT,*) "# time_step=",DT
- write(IOUT,*) "# location=",trim(dataT%name(i))
- write(IOUT,*) "# Column #1 = Time (s)"
- write(IOUT,*) "# Column #2 = horizontal right-lateral slip (m)"
- write(IOUT,*) "# Column #3 = horizontal right-lateral slip rate (m/s)"
- write(IOUT,*) "# Column #4 = horizontal right-lateral shear stress (MPa)"
- write(IOUT,*) "# Column #5 = vertical up-dip slip (m)"
- write(IOUT,*) "# Column #6 = vertical up-dip slip rate (m/s)"
- write(IOUT,*) "# Column #7 = vertical up-dip shear stress (MPa)"
- write(IOUT,*) "# Column #8 = normal stress (MPa)"
- if (RATE_AND_STATE) write(IOUT,*) "# Column #9 = log10 of state variable (log-seconds)"
- write(IOUT,*) "#"
- write(IOUT,*) "# The line below lists the names of the data fields:"
- if (.not. RATE_AND_STATE) then
- 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%t3(k,i)/1.0e6_CUSTOM_REAL
- enddo
- else
- if(statelaw == 1) then
- write(IOUT,*) "t h-slip h-slip-rate h-shear-stress v-slip v-slip-rate v-shear-stress n-stress log-theta"
- write(IOUT,*) "#"
- do k=1,NT
- write(IOUT,'(9(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%t3(k,i)/1.0e6_CUSTOM_REAL, log10(dataT%theta(k,i))
- enddo
- else
- write(IOUT,*) "t h-slip h-slip-rate h-shear-stress v-slip v-slip-rate v-shear-stress n-stress psi"
- write(IOUT,*) "#"
- do k=1,NT
- write(IOUT,'(9(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%t3(k,i)/1.0e6_CUSTOM_REAL, dataT%theta(k,i)
- enddo
- endif
- endif
- close(IOUT)
-
-
-1000 format ( ' # Date = ', i2.2, '/', i2.2, '/', i4.4, '; time = ',i2.2, ':', i2.2, ':', i2.2 )
-
-
- enddo
-
-end subroutine SCEC_write_dataT
-
-!-------------------------------------------------------------------------------------------------
-
end module fault_solver_common
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_kinematic.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_kinematic.f90 2012-10-24 17:14:27 UTC (rev 20906)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_kinematic.f90 2012-10-24 18:35:56 UTC (rev 20907)
@@ -35,8 +35,38 @@
private
+ ! outputs on selected fault nodes at every time step:
+ ! slip, slip velocity, fault stresses
+ type dataT_type
+ integer :: npoin=0,npoin_local=0
+ integer, dimension(:), pointer :: iglob=>null() ! on-fault global index of output nodes
+ integer, dimension(:,:), pointer :: iglob_all=>null()
+ integer, dimension(:), pointer :: islice=>null(),glob_indx=>null()
+ real(kind=CUSTOM_REAL), dimension(:), pointer :: dist=>null()
+ real(kind=CUSTOM_REAL), dimension(:,:), pointer :: dist_all=>null()
+ real(kind=CUSTOM_REAL), dimension(:,:), pointer :: d1=>null(),v1=>null(),t1=>null(), &
+ d2=>null(),v2=>null(),t2=>null(), &
+ t3=>null()
+ character(len=70), dimension(:), pointer :: name=>null()
+ end type dataT_type
+
+
+ ! outputs(dyn) /inputs (kind) at selected times for all fault nodes:
+ ! strength, state, slip, slip velocity, fault stresses, rupture time, process zone time
+ ! rupture time = first time when slip velocity = threshold V_RUPT (defined below)
+ ! process zone time = first time when slip = Dc
+ type dataXZ_type
+ real(kind=CUSTOM_REAL), dimension(:), pointer :: d1=>null(), d2=>null(), v1=>null(), v2=>null(), &
+ t1=>null(), t2=>null(), t3=>null()
+ real(kind=CUSTOM_REAL), dimension(:), pointer :: xcoord=>null(), ycoord=>null(), zcoord=>null()
+ integer :: npoin=0
+ end type dataXZ_type
+
+
type, EXTENDS (fault_type) :: bc_kinflt_type
private
+ type(dataT_type) :: dataT
+ type(dataXZ_type) :: dataXZ,dataXZ_all
real(kind=CUSTOM_REAL) :: kin_dt
integer :: kin_it
real(kind=CUSTOM_REAL), dimension(:,:), pointer :: v_kin_t1,v_kin_t2
@@ -44,6 +74,9 @@
type(bc_kinflt_type), allocatable, save :: faults(:)
+ !Number of time steps defined by the user : NTOUT
+ integer, save :: NTOUT,NSNAP
+
logical, save :: SIMULATION_TYPE_KIN = .false.
public :: BC_KINFLT_init, BC_KINFLT_set_all, SIMULATION_TYPE_KIN
@@ -138,102 +171,10 @@
NAMELIST / KINPAR / kindt
- read(IIN_BIN) bc%nspec,bc%nglob
- if (.NOT.PARALLEL_FAULT .and. bc%nspec==0) return
- if (bc%nspec>0) then
+ call initialize_simulation(bc,IIN_BIN,Minv,dt)
- 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
-
- endif
-
- if (PARALLEL_FAULT) then
-
- accel=0._CUSTOM_REAL
- 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)
-
- accel=0._CUSTOM_REAL
- if (bc%nspec>0) then
- 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, &
- 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) then
- nx = accel(1,bc%ibulk1)
- ny = accel(2,bc%ibulk1)
- nz = accel(3,bc%ibulk1)
- endif
-
- endif
-
if (bc%nspec>0) then
- 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)
- deallocate(nx,ny,nz)
- ! inverse M 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 ))
- ! WARNING: In non-split nodes at fault edges Minv is assembled across the fault
- ! hence invM1+invM2=2/(M1+M2) instead of 1/M1+1/M2
- ! For instance, in a symmetric mesh (M1=M2) Z is twice its intended value
-
allocate(bc%T(3,bc%nglob))
allocate(bc%D(3,bc%nglob))
allocate(bc%V(3,bc%nglob))
@@ -464,5 +405,232 @@
end subroutine load_vslip_snapshots
!---------------------------------------------------------------
+
+!===============================================================
+! OUTPUTS
+subroutine init_dataT(DataT,coord,nglob,NT,iflt)
+
+ use specfem_par, only : NPROC,myrank
+ ! NT = total number of time steps
+
+ integer, intent(in) :: nglob,NT,iflt
+ real(kind=CUSTOM_REAL), intent(in) :: coord(3,nglob)
+ type (dataT_type), intent(out) :: DataT
+
+ real(kind=CUSTOM_REAL) :: xtarget,ytarget,ztarget,dist,distkeep
+ integer :: i, iglob , IIN, ier, jflt, np, k
+ character(len=70) :: tmpname
+
+ integer :: ipoin, ipoin_local, iproc
+
+ ! 1. read fault output coordinates from user file,
+ ! 2. define iglob: the fault global index of the node nearest to user
+ ! requested coordinate
+
+ IIN = 251
+ open(IIN,file='DATA/FAULT/FAULT_STATIONS',status='old',action='read',iostat=ier)
+ read(IIN,*) np
+ DataT%npoin =0
+ do i=1,np
+ read(IIN,*) xtarget,ytarget,ztarget,tmpname,jflt
+ if (jflt==iflt) DataT%npoin = DataT%npoin +1
+ enddo
+ close(IIN)
+
+ if (DataT%npoin == 0) return
+
+ allocate(DataT%iglob(DataT%npoin))
+ allocate(DataT%name(DataT%npoin))
+ allocate(DataT%dist(DataT%npoin)) !Surendra : for parallel fault
+
+ open(IIN,file='DATA/FAULT/FAULT_STATIONS',status='old',action='read',iostat=ier)
+ if( ier /= 0 ) stop 'error opening FAULT_STATIONS file'
+ read(IIN,*) np
+ k = 0
+ do i=1,np
+ read(IIN,*) xtarget,ytarget,ztarget,tmpname,jflt
+ if (jflt/=iflt) cycle
+ k = k+1
+ DataT%name(k) = tmpname
+ !search nearest node
+ distkeep = huge(distkeep)
+
+ do iglob=1,nglob
+ dist = sqrt((coord(1,iglob)-xtarget)**2 &
+ + (coord(2,iglob)-ytarget)**2 &
+ + (coord(3,iglob)-ztarget)**2)
+ if (dist < distkeep) then
+ distkeep = dist
+ DataT%iglob(k) = iglob
+ endif
+ enddo
+ DataT%dist(k) = distkeep !Surendra : for parallel fault
+ enddo
+
+ !Surendra : for parallel fault
+ if (PARALLEL_FAULT) then
+ allocate(DataT%islice(DataT%npoin))
+ allocate(DataT%iglob_all(DataT%npoin,0:NPROC-1))
+ allocate(DataT%dist_all(DataT%npoin,0:NPROC-1))
+ call gather_all_i(DataT%iglob,DataT%npoin,DataT%iglob_all,DataT%npoin,NPROC)
+ call gather_all_cr(DataT%dist,DataT%npoin,DataT%dist_all,DataT%npoin,NPROC)
+ if(myrank==0) then
+ do ipoin = 1,DataT%npoin
+ iproc = minloc(DataT%dist_all(ipoin,:), 1) - 1
+ DataT%islice(ipoin) = iproc
+ DataT%iglob(ipoin) = DataT%iglob_all(ipoin,iproc)
+ enddo
+ endif
+
+ call bcast_all_i(DataT%islice,DataT%npoin)
+ call bcast_all_i(DataT%iglob,DataT%npoin)
+
+ DataT%npoin_local = 0
+ do ipoin = 1,DataT%npoin
+ if(myrank == DataT%islice(ipoin)) DataT%npoin_local = DataT%npoin_local + 1
+ enddo
+ allocate(DataT%glob_indx(DataT%npoin_local))
+ do ipoin = 1,DataT%npoin
+ if(myrank == DataT%islice(ipoin)) then
+ ipoin_local = ipoin_local + 1
+ DataT%glob_indx(ipoin_local) = ipoin
+ endif
+ enddo
+ else
+ DataT%npoin_local = DataT%npoin
+ endif !Parallel_fault
+
+
+ ! 3. allocate arrays and set to zero
+ allocate(DataT%d1(NT,DataT%npoin))
+ allocate(DataT%v1(NT,DataT%npoin))
+ allocate(DataT%t1(NT,DataT%npoin))
+ allocate(DataT%d2(NT,DataT%npoin))
+ allocate(DataT%v2(NT,DataT%npoin))
+ allocate(DataT%t2(NT,DataT%npoin))
+ allocate(DataT%t3(NT,DataT%npoin))
+ allocate(DataT%theta(NT,DataT%npoin))
+ DataT%d1 = 0e0_CUSTOM_REAL
+ DataT%v1 = 0e0_CUSTOM_REAL
+ DataT%t1 = 0e0_CUSTOM_REAL
+ DataT%d2 = 0e0_CUSTOM_REAL
+ DataT%v2 = 0e0_CUSTOM_REAL
+ DataT%t2 = 0e0_CUSTOM_REAL
+ DataT%t3 = 0e0_CUSTOM_REAL
+ DataT%theta = 0e0_CUSTOM_REAL
+
+ close(IIN)
+
+end subroutine init_dataT
+
+!---------------------------------------------------------------
+subroutine store_dataT(dataT,d,v,t,theta,itime)
+
+ type(dataT_type), intent(inout) :: dataT
+ real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: d,v,t
+ real(kind=CUSTOM_REAL), dimension(:), intent(in) :: theta
+ integer, intent(in) :: itime
+
+ integer :: i,k,ipoin
+
+ do ipoin=1,dataT%npoin_local
+ if (PARALLEL_FAULT) then
+ i = DataT%glob_indx(ipoin)
+ else
+ i = ipoin
+ endif
+ k = dataT%iglob(i)
+ dataT%d1(itime,i) = d(1,k)
+ dataT%d2(itime,i) = d(2,k)
+ dataT%v1(itime,i) = v(1,k)
+ dataT%v2(itime,i) = v(2,k)
+ dataT%t1(itime,i) = t(1,k)
+ dataT%t2(itime,i) = t(2,k)
+ dataT%t3(itime,i) = t(3,k)
+ dataT%theta(itime,i) = theta(k)
+ enddo
+
+end subroutine store_dataT
+
+!-----------------------------------------------------------------
+subroutine write_dataT_all(nt,statelaw)
+
+ integer, intent(in) :: nt
+ integer, intent(in) :: statelaw
+
+ integer :: i
+
+ if (.not.allocated(faults)) return
+ do i = 1,size(faults)
+ call SCEC_write_dataT(faults(i)%dataT,faults(i)%dt,nt,statelaw)
+ enddo
+
+end subroutine write_dataT_all
+
+!------------------------------------------------------------------------
+subroutine SCEC_write_dataT(dataT,DT,NT,statelaw)
+
+ type(dataT_type), intent(in) :: dataT
+ real(kind=CUSTOM_REAL), intent(in) :: DT
+ integer, intent(in) :: NT
+ integer, intent(in) :: statelaw
+
+ integer :: i,k,IOUT,ipoin
+ character :: NTchar*5
+ 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
+
+ IOUT = 121 !WARNING: not very robust. Could instead look for an available ID
+
+ write(NTchar,1) NT
+ NTchar = adjustl(NTchar)
+
+1 format(I5)
+ do ipoin=1,dataT%npoin_local
+ if (PARALLEL_FAULT) then
+ i = DataT%glob_indx(ipoin)
+ else
+ i = ipoin
+ endif
+
+ open(IOUT,file='OUTPUT_FILES/'//trim(dataT%name(i))//'.dat',status='replace')
+ write(IOUT,*) "# problem=Kinmeatic solver"
+ write(IOUT,*) "# author=Surendra Nadh Somala"
+ 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=100 m (*5 GLL nodes)"
+ write(IOUT,*) "# time_step=",DT
+ write(IOUT,*) "# location=",trim(dataT%name(i))
+ write(IOUT,*) "# Column #1 = Time (s)"
+ write(IOUT,*) "# Column #2 = horizontal right-lateral slip (m)"
+ write(IOUT,*) "# Column #3 = horizontal right-lateral slip rate (m/s)"
+ write(IOUT,*) "# Column #4 = horizontal right-lateral shear stress (MPa)"
+ write(IOUT,*) "# Column #5 = vertical up-dip slip (m)"
+ write(IOUT,*) "# Column #6 = vertical up-dip slip rate (m/s)"
+ write(IOUT,*) "# Column #7 = vertical up-dip shear stress (MPa)"
+ 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,*) "#"
+ 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%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 )
+
+
+ enddo
+
+end subroutine SCEC_write_dataT
+
+!-------------------------------------------------------------------------------------------------
+
end module fault_solver_kinematic
-
More information about the CIG-COMMITS
mailing list