[cig-commits] r20900 - seismo/3D/FAULT_SOURCE/branches/new_fault_db/src
surendra at geodynamics.org
surendra at geodynamics.org
Wed Oct 24 06:46:55 PDT 2012
Author: surendra
Date: 2012-10-24 06:46:55 -0700 (Wed, 24 Oct 2012)
New Revision: 20900
Modified:
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_kinematic.f90
Log:
Started working on re-factoring common fault routines of dyn & kin types
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 13:37:51 UTC (rev 20899)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90 2012-10-24 13:46:55 UTC (rev 20900)
@@ -37,33 +37,7 @@
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 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
@@ -87,22 +61,12 @@
real(kind=CUSTOM_REAL), dimension(:), pointer :: Fload=>null()
end type asperity_type
- type bc_dynflt_type
+ type, EXTENDS (fault_type) :: bc_dynflt_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 :: 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) :: dt
- integer, dimension(:), pointer :: ibulk1=>null(), ibulk2=>null()
type(swf_type), pointer :: swf => null()
type(rsf_type), pointer :: rsf => null()
type(asperity_type), pointer :: asp => null()
logical :: allow_opening = .false. ! default : do not allow opening
- type(dataT_type) :: dataT
- type(dataXZ_type) :: dataXZ,dataXZ_all
- integer, dimension(:), pointer :: poin_offset=>null(),npoin_perproc=>null()
end type bc_dynflt_type
type(bc_dynflt_type), allocatable, save :: faults(:)
@@ -114,13 +78,8 @@
!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 :: PARALLEL_FAULT = .false.
-
logical, save :: TPV16 = .false.
logical, save :: RATE_AND_STATE = .false.
@@ -128,7 +87,7 @@
real(kind=CUSTOM_REAL), allocatable, save :: Kelvin_Voigt_eta(:)
public :: BC_DYNFLT_init, BC_DYNFLT_set3d_all, Kelvin_Voigt_eta, &
- SIMULATION_TYPE_DYN,PARALLEL_FAULT,RATE_AND_STATE
+ SIMULATION_TYPE_DYN,,RATE_AND_STATE
contains
@@ -669,128 +628,6 @@
end subroutine init_one_fault
-!---------------------------------------------------------------------
-subroutine compute_R(R,nglob,nx,ny,nz)
-
- integer :: nglob
- real(kind=CUSTOM_REAL), intent(out) :: R(3,3,nglob)
- real(kind=CUSTOM_REAL), dimension(nglob), intent(in) :: nx,ny,nz
-
- real(kind=CUSTOM_REAL), dimension(nglob) :: sx,sy,sz,dx,dy,dz,norm
-
- ! Percy , defining fault directions (in concordance with SCEC conventions) .
- ! fault coordinates (s,d,n) = (1,2,3)
- ! s = strike , d = dip , n = n.
- ! 1 = strike , 2 = dip , 3 = n.
- norm = sqrt(nx*nx+ny*ny)
- sx = ny/norm
- sy = -nx/norm
- sz = 0.e0_CUSTOM_REAL
-
- norm = sqrt(sy*sy*nz*nz+sx*sx*nz*nz+(sy*nx-ny*sx)*(nx*sy-ny*sx))
- dx = -sy*nz/norm
- dy = sx*nz/norm
- dz = (sy*nx-ny*sx)/norm
- !Percy, dz is always dipwards = -1/norm , because (nx*sy-ny*sx)= - 1
-
- R(1,1,:)=sx
- R(1,2,:)=sy
- R(1,3,:)=sz
- R(2,1,:)=dx
- R(2,2,:)=dy
- R(2,3,:)=dz
- R(3,1,:)=nx
- R(3,2,:)=ny
- R(3,3,:)=nz
-
-end subroutine compute_R
-
-!---------------------------------------------------------------------
-! 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.
!
@@ -1030,33 +867,7 @@
end function get_weighted_jump
-!----------------------------------------------------------------------
-function rotate(bc,v,fb) result(vr)
- type(bc_dynflt_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)
-
- ! Percy, tangential direction Vt, equation 7 of Pablo's notes in agreement with SPECFEM3D
-
- ! forward rotation
- if (fb==1) then
- vr(1,:) = v(1,:)*bc%R(1,1,:)+v(2,:)*bc%R(1,2,:)+v(3,:)*bc%R(1,3,:) ! vs
- vr(2,:) = v(1,:)*bc%R(2,1,:)+v(2,:)*bc%R(2,2,:)+v(3,:)*bc%R(2,3,:) ! vd
- vr(3,:) = v(1,:)*bc%R(3,1,:)+v(2,:)*bc%R(3,2,:)+v(3,:)*bc%R(3,3,:) ! vn
-
- ! backward rotation
- else
- vr(1,:) = v(1,:)*bc%R(1,1,:)+v(2,:)*bc%R(2,1,:)+v(3,:)*bc%R(3,1,:) !vx
- vr(2,:) = v(1,:)*bc%R(1,2,:)+v(2,:)*bc%R(2,2,:)+v(3,:)*bc%R(3,2,:) !vy
- vr(3,:) = v(1,:)*bc%R(1,3,:)+v(2,:)*bc%R(2,3,:)+v(3,:)*bc%R(3,3,:) !vz
-
- endif
-
-end function rotate
-
-
!=====================================================================
subroutine swf_update_state(dold,dnew,vold,f)
@@ -1209,253 +1020,6 @@
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_kinematic.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_kinematic.f90 2012-10-24 13:37:51 UTC (rev 20899)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_kinematic.f90 2012-10-24 13:46:55 UTC (rev 20900)
@@ -35,33 +35,8 @@
private
-! outputs on selected fault nodes at every time step:
-! slip, slip velocity, fault stresses
- type dataT_type
- integer :: npoin
- integer, dimension(:), pointer :: iglob
- real(kind=CUSTOM_REAL), dimension(:,:), pointer :: d1,v1,t1,d2,v2,t2,t3
- character(len=70), dimension(:), pointer :: name
- end type dataT_type
-
-! DATAXZ_type used to read snapshots (temporal)
- type dataXZ_type
- real(kind=CUSTOM_REAL), dimension(:), pointer :: d1, d2, v1, v2, & !Slip and Slip rate.
- t1, t2, t3 !Tractions.
- real(kind=CUSTOM_REAL), dimension(:), pointer :: xcoord,ycoord,zcoord
- integer :: npoin
- end type dataXZ_type
-
- type bc_kinflt_type
+ type, EXTENDS (fault_type) :: bc_kinflt_type
private
- integer :: nspec,nglob
- real(kind=CUSTOM_REAL) :: dt
- real(kind=CUSTOM_REAL), dimension(:), pointer :: B,invM1,invM2,Z
- real(kind=CUSTOM_REAL), dimension(:,:), pointer :: T,D,V,coord
- real(kind=CUSTOM_REAL), dimension(:,:,:), pointer :: R
- integer, dimension(:), pointer :: ibulk1, ibulk2
- type(dataT_type) :: dataT
- type(dataXZ_type) :: dataXZ
real(kind=CUSTOM_REAL) :: kin_dt
integer :: kin_it
real(kind=CUSTOM_REAL), dimension(:,:), pointer :: v_kin_t1,v_kin_t2
@@ -69,9 +44,6 @@
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
@@ -427,140 +399,7 @@
end subroutine BC_KINFLT_set_single
!===============================================================
-function get_jump(bc,v) result(dv)
- type(bc_kinflt_type), intent(in) :: bc
- real(kind=CUSTOM_REAL), intent(in) :: v(:,:)
- real(kind=CUSTOM_REAL) :: dv(3,bc%nglob)
-
- ! diference between side 2 and side 1 of fault nodes. dv
- dv(1,:) = v(1,bc%ibulk2)-v(1,bc%ibulk1)
- dv(2,:) = v(2,bc%ibulk2)-v(2,bc%ibulk1)
- dv(3,:) = v(3,bc%ibulk2)-v(3,bc%ibulk1)
-
-end function get_jump
-
-!---------------------------------------------------------------------
-function get_weighted_jump(bc,f) result(da)
-
- type(bc_kinflt_type), intent(in) :: bc
- real(kind=CUSTOM_REAL), intent(in) :: f(:,:)
- real(kind=CUSTOM_REAL) :: da(3,bc%nglob)
-
- ! diference between side 2 and side 1 of fault nodes. M-1 * F
- da(1,:) = bc%invM2*f(1,bc%ibulk2)-bc%invM1*f(1,bc%ibulk1)
- da(2,:) = bc%invM2*f(2,bc%ibulk2)-bc%invM1*f(2,bc%ibulk1)
- da(3,:) = bc%invM2*f(3,bc%ibulk2)-bc%invM1*f(3,bc%ibulk1)
-
- ! NOTE: In non-split nodes at fault edges da=0
- ! (invM and f are assembled across the fault, hence f1=f2 and invM1=invM2)
-
-end function get_weighted_jump
-
-!----------------------------------------------------------------------
-function rotate(bc,v,fb) result(vr)
-
- type(bc_kinflt_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)
-
- ! Percy, tangential direction Vt, equation 7 of Pablo's notes in agreement with SPECFEM3D
-
- ! forward rotation
- if (fb==1) then
- vr(1,:) = v(1,:)*bc%R(1,1,:)+v(2,:)*bc%R(1,2,:)+v(3,:)*bc%R(1,3,:) ! vs
- vr(2,:) = v(1,:)*bc%R(2,1,:)+v(2,:)*bc%R(2,2,:)+v(3,:)*bc%R(2,3,:) ! vd
- vr(3,:) = v(1,:)*bc%R(3,1,:)+v(2,:)*bc%R(3,2,:)+v(3,:)*bc%R(3,3,:) ! vn
-
- ! backward rotation
- else
- vr(1,:) = v(1,:)*bc%R(1,1,:)+v(2,:)*bc%R(2,1,:)+v(3,:)*bc%R(3,1,:) !vx
- vr(2,:) = v(1,:)*bc%R(1,2,:)+v(2,:)*bc%R(2,2,:)+v(3,:)*bc%R(3,2,:) !vy
- vr(3,:) = v(1,:)*bc%R(1,3,:)+v(2,:)*bc%R(2,3,:)+v(3,:)*bc%R(3,3,:) !vz
-
- endif
-
-end function rotate
-
-
-!===============================================================
-! OUTPUTS
-
-subroutine init_dataT(DataT,coord,nglob,NT,iflt)
- ! 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
-
- ! 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))
-
- 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
- enddo
-
- ! 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))
- 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
-
- close(IIN)
-
-end subroutine init_dataT
-
-
!---------------------------------------------------------------
subroutine init_dataXZ(dataXZ,nglob)
@@ -582,78 +421,8 @@
end subroutine init_dataXZ
-!---------------------------------------------------------------
-subroutine store_dataT(dataT,d,v,t,itime)
- type(dataT_type), intent(inout) :: dataT
- real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: d,v,t
- integer, intent(in) :: itime
- integer :: i,k
-
- do i=1,dataT%npoin
- 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)
- enddo
-
-end subroutine store_dataT
-
-
-!-----------------------------------------------------------------
-
-subroutine SCEC_write_dataT(dataT,DT,NT)
-
- type(dataT_type), intent(in) :: dataT
- real(kind=CUSTOM_REAL), intent(in) :: DT
- integer, intent(in) :: NT
-
- integer :: i,k,IOUT
-
- IOUT = 121 !WARNING: not very robust. Could instead look for an available ID
-
- do i=1,dataT%npoin
-
- open(IOUT,file='OUTPUT_FILES/'//trim(dataT%name(i))//'.dat',status='replace')
- write(IOUT,*) "% problem=TPV5"
- write(IOUT,*) "% author=Galvez, Ampuero, Nissen-Meyer"
- write(IOUT,*) "% date=2010/xx/xx"
- write(IOUT,*) "% code=SPECFEM3D_FAULT "
- write(IOUT,*) "% code_version=1.1"
- write(IOUT,*) "% element_size=100 m (*4 GLL nodes)"
- write(IOUT,*) "% time_step=",DT
- write(IOUT,*) "% num_time_steps=",NT
- write(IOUT,*) "% location=",trim(dataT%name(i))
- write(IOUT,*) "% Time series in 8 column of E15.7"
- 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,*) "%"
- write(IOUT,*) "% Here is the time-series data."
- 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)
- enddo
-
-end subroutine SCEC_write_dataT
-
-
!---------------------------------------------------------------
!LOAD_VSLIP_SNAPSHOTS(v,dataXZ,itime,coord,npoin,nglob,iflt)
!Loading slip velocity from snapshots.
@@ -688,22 +457,6 @@
read(IIN_BIN,"(5F24.15)") dataXZ%xcoord,dataXZ%ycoord,dataXZ%zcoord,dataXZ%v1,dataXZ%v2
-! read(IOUT) dataXZ%xcoord
-! read(IOUT) dataXZ%ycoord
-! read(IOUT) dataXZ%zcoord
-! write(IOUT) dataXZ%d1
-! write(IOUT) dataXZ%d2
-! read(IOUT) dataXZ%v1
-! read(IOUT) dataXZ%v2
-! write(IOUT) dataXZ%t1
-! write(IOUT) dataXZ%t2
-! write(IOUT) dataXZ%t3
-! write(IOUT) dataXZ%sta
-! write(IOUT) dataXZ%stg
-! write(IOUT) dataXZ%tRUP
-! write(IOUT) dataXZ%tPZ
- close(IOUT)
-
close(IIN_BIN)
end subroutine load_vslip_snapshots
More information about the CIG-COMMITS
mailing list