[cig-commits] r20970 - in seismo/3D/FAULT_SOURCE: . branches/new_fault_db/src
ampuero at geodynamics.org
ampuero at geodynamics.org
Mon Oct 29 22:30:39 PDT 2012
Author: ampuero
Date: 2012-10-29 22:30:39 -0700 (Mon, 29 Oct 2012)
New Revision: 20970
Modified:
seismo/3D/FAULT_SOURCE/TO_DO
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_common.f90
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_dynamic.f90
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_kinematic.f90
Log:
refactored dataT; added statelaw to call funcd
Modified: seismo/3D/FAULT_SOURCE/TO_DO
===================================================================
--- seismo/3D/FAULT_SOURCE/TO_DO 2012-10-29 18:56:50 UTC (rev 20969)
+++ seismo/3D/FAULT_SOURCE/TO_DO 2012-10-30 05:30:39 UTC (rev 20970)
@@ -8,8 +8,6 @@
+ fault_solver_dynamic:
- many hard-coded ad hoc features need to be set instead as user options
-+ dataT: move common stuff to fault_solver_common, use inheritance in the dyn/kin modules
-
+ rate-and-state friction:
- make it a user-friendly option (currently a flag in fault_solver)
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-29 18:56:50 UTC (rev 20969)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_common.f90 2012-10-30 05:30:39 UTC (rev 20970)
@@ -1,6 +1,7 @@
! Base module for kinematic and dynamic fault solvers
!
-! Written by: Percy Galvez, Surendra Somala and Jean-Paul Ampuero
+! Authors:
+! Percy Galvez, Surendra Somala, Jean-Paul Ampuero
module fault_solver_common
@@ -9,22 +10,31 @@
include 'constants.h'
private
-
+
type fault_type
integer :: nspec=0, nglob=0
- 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 :: B=>null(),invM1=>null(),invM2=>null(),Z=>null()
- real(kind=CUSTOM_REAL) :: dt
+ real(kind=CUSTOM_REAL), dimension(:,:), pointer :: T=>null(),V=>null(),D=>null(),coord=>null()
+ real(kind=CUSTOM_REAL), dimension(:,:,:), pointer :: R=>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()
- integer, dimension(:), pointer :: poin_offset=>null(),npoin_perproc=>null()
end type fault_type
+ ! outputs on selected fault nodes at every time step:
+ type dataT_type
+ integer :: npoin=0, ndat=0, nt=0
+ real(kind=CUSTOM_REAL) :: dt
+ integer, dimension(:), pointer :: iglob=>null() ! on-fault global index of output nodes
+ real(kind=CUSTOM_REAL), dimension(:,:,:), pointer :: dat=>null()
+ character(len=70), dimension(:), pointer :: name=>null(),longFieldNames=>null()
+ character(len=100) :: shortFieldNames
+ end type dataT_type
+
logical, save :: PARALLEL_FAULT = .false.
public :: fault_type, PARALLEL_FAULT, &
- initialize_fault, get_jump, get_weighted_jump, rotate, add_BT
+ initialize_fault, get_jump, get_weighted_jump, rotate, add_BT, &
+ dataT, init_dataT, store_dataT, SCEC_write_dataT
contains
@@ -179,6 +189,7 @@
end subroutine compute_R
+
!===============================================================
function get_jump (bc,v) result(dv)
@@ -255,4 +266,213 @@
end subroutine add_BT
+
+!===============================================================
+! dataT outputs
+
+subroutine init_dataT(dataT,coord,nglob,NT,DT,ndat,iflt)
+
+ use specfem_par, only : NPROC,myrank
+
+ integer, intent(in) :: nglob,NT,iflt,ndat
+ real(kind=CUSTOM_REAL), intent(in) :: coord(3,nglob),DT
+ class (dataT_type), intent(out) :: dataT
+
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: dist_all
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: dist_loc
+ integer, dimension(:,:), allocatable :: iglob_all
+ integer, dimension(:), allocatable :: iproc,iglob_tmp,glob_indx
+ real(kind=CUSTOM_REAL) :: xtarget,ytarget,ztarget,dist,distkeep
+ integer :: i, iglob , IIN, ier, jflt, np, k
+ character(len=70) :: tmpname
+ character(len=70), dimension(:), allocatable :: name_tmp
+ integer :: ipoin, ipoin_local, npoin_local
+
+ ! 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 ! WARNING: not safe, should check that unit is not aleady opened
+
+ ! count the number of output points on the current fault (#iflt)
+ open(IIN,file='DATA/FAULT_STATIONS',status='old',action='read',iostat=ier)
+ if( ier /= 0 ) stop 'error opening FAULT_STATIONS file'
+ 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(dist_loc(dataT%npoin)) !Surendra : for parallel fault
+
+ open(IIN,file='DATA/FAULT_STATIONS',status='old',action='read')
+ 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
+ dist_loc(k) = distkeep
+
+ enddo
+ close(IIN)
+
+ if (PARALLEL_FAULT) then
+
+ ! For each output point, find the processor that contains the nearest node
+ allocate(iproc(dataT%npoin))
+ allocate(iglob_all(dataT%npoin,0:NPROC-1))
+ allocate(dist_all(dataT%npoin,0:NPROC-1))
+ call gather_all_i(dataT%iglob,dataT%npoin,iglob_all,dataT%npoin,NPROC)
+ call gather_all_cr(dist_loc,dataT%npoin,dist_all,dataT%npoin,NPROC)
+ if (myrank==0) then
+ ! NOTE: output points lying at an interface between procs are assigned to a unique proc
+ iproc = minloc(dist_all,2) - 1
+ dataT%iglob = iglob_all(:,iproc)
+ endif
+ call bcast_all_i(iproc,dataT%npoin)
+ call bcast_all_i(dataT%iglob,dataT%npoin)
+
+ ! Number of output points contained in the current processor
+ npoin_local = count( iproc == myrank )
+
+ if (npoin_local>0) then
+ ! Make a list of output points contained in the current processor
+ allocate(glob_indx(npoin_local))
+ ipoin_local = 0
+ do ipoin = 1,dataT%npoin
+ if (myrank == iproc(ipoin)) then
+ ipoin_local = ipoin_local + 1
+ glob_indx(ipoin_local) = ipoin
+ endif
+ enddo
+ ! Consolidate the output information (remove output points outside current proc)
+ allocate(iglob_tmp(dataT%npoin))
+ allocate(name_tmp(dataT%npoin))
+ iglob_tmp = dataT%iglob
+ name_tmp = dataT%name
+ deallocate(dataT%iglob)
+ deallocate(dataT%name)
+ dataT%npoin = npoin_local
+ allocate(dataT%iglob(dataT%npoin))
+ allocate(dataT%name(dataT%npoin))
+ dataT%iglob = iglob_tmp(glob_indx)
+ dataT%name = name_tmp(glob_indx)
+ deallocate(glob_indx,iglob_tmp,name_tmp)
+
+ else
+ dataT%npoin = 0
+ deallocate(dataT%iglob)
+ deallocate(dataT%name)
+ endif
+
+ deallocate(iproc,iglob_all,dist_all)
+ endif
+
+ ! 3. initialize arrays
+ if (dataT%npoin>0) then
+ dataT%ndat = ndat
+ dataT%nt = NT
+ dataT%dt = DT
+ allocate(dataT%dat(dataT%ndat,dataT%nt,dataT%npoin))
+ dataT%dat = 0e0_CUSTOM_REAL
+ allocate(dataT%longFieldNames(dataT%ndat))
+ dataT%longFieldNames(1) = "horizontal right-lateral slip (m)"
+ dataT%longFieldNames(2) = "horizontal right-lateral slip rate (m/s)"
+ dataT%longFieldNames(3) = "horizontal right-lateral shear stress (MPa)"
+ dataT%longFieldNames(4) = "vertical up-dip slip (m)"
+ dataT%longFieldNames(5) = "vertical up-dip slip rate (m/s)"
+ dataT%longFieldNames(6) = "vertical up-dip shear stress (MPa)"
+ dataT%longFieldNames(7) = "normal stress (MPa)"
+ dataT%shortFieldNames = "h-slip h-slip-rate h-shear-stress v-slip v-slip-rate v-shear-stress n-stress"
+ endif
+
+end subroutine init_dataT
+
+!---------------------------------------------------------------
+subroutine store_dataT(dataT,d,v,t,itime)
+
+ class(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%dat(1,itime,i) = d(1,k)
+ dataT%dat(3,itime,i) = v(1,k)
+ dataT%dat(5,itime,i) = t(1,k)/1.0e6_CUSTOM_REAL
+ dataT%dat(2,itime,i) = -d(2,k)
+ dataT%dat(4,itime,i) = -v(2,k)
+ dataT%dat(6,itime,i) = -t(2,k)/1.0e6_CUSTOM_REAL
+ dataT%dat(7,itime,i) = t(3,k)/1.0e6_CUSTOM_REAL
+ enddo
+
+end subroutine store_dataT
+
+!------------------------------------------------------------------------
+subroutine SCEC_write_dataT(dataT)
+
+ class(my_dataT_type), intent(in) :: dataT
+
+ integer :: i,k,IOUT
+ character(len=10) :: my_fmt
+ 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(my_fmt,'(a,i1,a)') '(',dataT%ndat,'(E15.7))'
+
+ do i=1,dataT%npoin
+ open(IOUT,file='OUTPUT_FILES/'//trim(dataT%name(i))//'.dat',status='replace')
+ write(IOUT,*) "# problem=TPV104" ! WARNING: this should be a user input
+ write(IOUT,*) "# author=Surendra Nadh Somala" ! WARNING: this should be a user input
+ 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)" ! WARNING: this should be a user input
+ write(IOUT,*) "# time_step=",dataT%dt
+ write(IOUT,*) "# location=",trim(dataT%name(i))
+ write(IOUT,*) "# Column #1 = Time (s)"
+ do k=1,dataT%ndat
+ write(IOUT,1100) k+1,dataT%longFieldNames(k)
+ enddo
+ write(IOUT,*) "#"
+ write(IOUT,*) "# The line below lists the names of the data fields:"
+ write(IOUT,*) "# t ", dataT%shortFieldNames
+ write(IOUT,*) "#"
+ do k=1,dataT%nt
+ write(IOUT,my_fmt) k*dataT%dt, dataT%dat(:,k,i)
+ enddo
+ close(IOUT)
+ enddo
+
+1000 format ( ' # Date = ', i2.2, '/', i2.2, '/', i4.4, '; time = ',i2.2, ':', i2.2, ':', i2.2 )
+1100 format ( ' # Column #', i1, ' = ',a )
+
+end subroutine SCEC_write_dataT
+
end module fault_solver_common
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_dynamic.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_dynamic.f90 2012-10-29 18:56:50 UTC (rev 20969)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_dynamic.f90 2012-10-30 05:30:39 UTC (rev 20970)
@@ -1,34 +1,13 @@
-!=====================================================================
+! This module implements dynamic faults: spontaneous rupture with prescribed
+! friction laws (slip-weakening or rate-and-state) and heterogeneous initial conditions
!
-! s p e c f e m 3 d v e r s i o n 1 . 4
-! ---------------------------------------
+! Authors:
+! Percy Galvez, Jean-Paul Ampuero, Tarje Nissen-Meyer, Surendra Somala
!
-! dimitri komatitsch and jeroen tromp
-! seismological laboratory - california institute of technology
-! (c) california institute of technology september 2006
-!
-! this program is free software; you can redistribute it and/or modify
-! it under the terms of the gnu general public license as published by
-! the free software foundation; either version 2 of the license, or
-! (at your option) any later version.
-!
-! this program is distributed in the hope that it will be useful,
-! but without any warranty; without even the implied warranty of
-! merchantability or fitness for a particular purpose. see the
-! gnu general public license for more details.
-!
-! you should have received a copy of the gnu general public license along
-! with this program; if not, write to the free software foundation, inc.,
-! 51 franklin street, fifth floor, boston, ma 02110-1301 usa.
-!
-!===============================================================================================================
+! Surendra Nadh Somala : heterogenous initial stress capabilities (based on TPV16)
+! Surendra Nadh Somala : rate and state friction
+! Somala & Ampuero : fault parallelization
-! This module was written by:
-! Percy Galvez, Jean-Paul Ampuero and Tarje Nissen-Meyer
-! Surendra Nadh Somala : Added Heterogenous initial stress capabilities (based on TPV16)
-! Surendra Nadh Somala : Added Rate and State Friction
-! Somala/Ampuero : fault parallelization
-
module fault_solver_dynamic
use fault_solver_common
@@ -39,17 +18,6 @@
private
- ! outputs on selected fault nodes at every time step:
- ! slip, slip velocity, fault stresses
- type dataT_type
- integer :: npoin=0
- integer, dimension(:), pointer :: iglob=>null() ! on-fault global index of output nodes
- 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)
@@ -82,12 +50,13 @@
type, extends (fault_type) :: bc_dynflt_type
private
real(kind=CUSTOM_REAL), dimension(:,:), pointer :: T0=>null()
- real(kind=CUSTOM_REAL), dimension(:), pointer :: MU=>null(), Fload=>null()
- type(dataT_type) :: dataT
- type(dataXZ_type) :: dataXZ,dataXZ_all
- type(swf_type), pointer :: swf => null()
- type(rsf_type), pointer :: rsf => null()
- logical :: allow_opening = .false. ! default : do not allow opening
+ real(kind=CUSTOM_REAL), dimension(:), pointer :: MU=>null(), Fload=>null()
+ integer, dimension(:), pointer :: npoin_perproc=>null(), npoin_offset=>null()
+ type(dataT_type) :: dataT
+ type(dataXZ_type) :: dataXZ,dataXZ_all
+ type(swf_type), pointer :: swf => null()
+ type(rsf_type), pointer :: rsf => null()
+ logical :: allow_opening = .false. ! default : do not allow opening
end type bc_dynflt_type
type(bc_dynflt_type), allocatable, save :: faults(:)
@@ -267,7 +236,17 @@
endif
- call init_dataT(bc%dataT,bc%coord,bc%nglob,NT,iflt)
+ if (RATE_AND_STATE) then
+ call init_dataT(bc%dataT,bc%coord,bc%nglob,NT,dt,8,iflt)
+ bc%dataT%longFieldNames(8) = "log10 of state variable (log-seconds)"
+ if (bc%rsf%StateLaw==1) then
+ bc%dataT%shortFieldNames = trim(bc%dataT%shortFieldNames)//" log-theta"
+ else
+ bc%dataT%shortFieldNames = trim(bc%dataT%shortFieldNames)//" psi"
+ endif
+ else
+ call init_dataT(bc%dataT,bc%coord,bc%nglob,NT,dt,7,iflt)
+ endif
call init_dataXZ(bc%dataXZ,bc)
!--------------------------------------------------------
@@ -578,11 +557,19 @@
call store_dataXZ(bc%dataXZ, strength, theta_old, theta_new, dc, &
Vf_old, Vf_new, it*bc%dt,bc%dt)
- call store_dataT(bc%dataT,bc%D,bc%V,bc%T,theta_new,it)
+ call store_dataT(bc%dataT,bc%D,bc%V,bc%T,it)
+ if (RATE_AND_STATE) then
+ if (bc%rsf%StateLaw==1) then
+ bc%dataT%dat(8,it,:) = log10(theta_new(bc%dataT%iglob))
+ else
+ bc%dataT%dat(8,it,:) = theta_new(bc%dataT%iglob)
+ endif
+ endif
+
!-- outputs --
! write dataT every NTOUT time step or at the end of simulation
- if ( mod(it,NTOUT) == 0 .or. it==NSTEP) call SCEC_write_dataT(bc%dataT,bc%dt,it,bc%rsf%StateLaw)
+ if ( mod(it,NTOUT) == 0 .or. it==NSTEP) call SCEC_write_dataT(bc%dataT,it)
endif
@@ -892,6 +879,7 @@
endif
fn = tStick - Z*x - a*Seff*asinh(x*arg)
df = -Z - a*Seff/sqrt(ONE + (x*arg)**2)*arg
+
end subroutine funcd
!---------------------------------------------------------------------
@@ -906,8 +894,8 @@
real(kind=CUSTOM_REAL) :: tStick,Seff,Z,f0,V0,a,b,L,theta
integer :: statelaw
- call funcd(dble(x1),fl,df,tStick,Seff,Z,f0,V0,a,b,L,theta)
- call funcd(dble(x2),fh,df,tStick,Seff,Z,f0,V0,a,b,L,theta)
+ call funcd(dble(x1),fl,df,tStick,Seff,Z,f0,V0,a,b,L,theta,statelaw)
+ call funcd(dble(x2),fh,df,tStick,Seff,Z,f0,V0,a,b,L,theta,statelaw)
if( (fl>0 .and. fh>0) .or. (fl<0 .and. fh<0) ) stop 'root must be bracketed in rtsafe'
if(fl==0.) then
rtsafe=x2
@@ -926,7 +914,7 @@
rtsafe=0.5d0*(x1+x2)
dxold=abs(x2-x1)
dx=dxold
- call funcd(rtsafe,f,df,tStick,Seff,Z,f0,V0,a,b,L,theta)
+ call funcd(rtsafe,f,df,tStick,Seff,Z,f0,V0,a,b,L,theta,statelaw)
do j=1,MAXIT
if( ((rtsafe-xh)*df-f)*((rtsafe-xl)*df-f)>0 .or. abs(2.*f)>abs(dxold*df) ) then
dxold=dx
@@ -941,7 +929,7 @@
if(temp==rtsafe) return
endif
if(abs(dx)<xacc) return
- call funcd(rtsafe,f,df,tStick,Seff,Z,f0,V0,a,b,L,theta)
+ call funcd(rtsafe,f,df,tStick,Seff,Z,f0,V0,a,b,L,theta,statelaw)
if(f<0.) then
xl=rtsafe
else
@@ -991,268 +979,7 @@
!===============================================================
! 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), dimension(:,:), allocatable :: dist_all
- real(kind=CUSTOM_REAL), dimension(:), allocatable :: dist_loc
- integer, dimension(:,:), allocatable :: iglob_all
- integer, dimension(:), allocatable :: iproc,iglob_tmp,glob_indx
- real(kind=CUSTOM_REAL) :: xtarget,ytarget,ztarget,dist,distkeep
- integer :: i, iglob , IIN, ier, jflt, np, k
- character(len=70) :: tmpname
- character(len=70), dimension(:), allocatable :: name_tmp
- integer :: ipoin, ipoin_local, npoin_local
-
- ! 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 ! WARNING: not safe, should check that unit is not aleady opened
-
- ! count the number of output points on the current fault (#iflt)
- open(IIN,file='DATA/FAULT_STATIONS',status='old',action='read',iostat=ier)
- if( ier /= 0 ) stop 'error opening FAULT_STATIONS file'
- 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(dist_loc(dataT%npoin)) !Surendra : for parallel fault
-
- open(IIN,file='DATA/FAULT_STATIONS',status='old',action='read')
- 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
- dist_loc(k) = distkeep
-
- enddo
- close(IIN)
-
- if (PARALLEL_FAULT) then
-
- ! For each output point, find the processor that contains the nearest node
- allocate(iproc(dataT%npoin))
- allocate(iglob_all(dataT%npoin,0:NPROC-1))
- allocate(dist_all(dataT%npoin,0:NPROC-1))
- call gather_all_i(dataT%iglob,dataT%npoin,iglob_all,dataT%npoin,NPROC)
- call gather_all_cr(dist_loc,dataT%npoin,dist_all,dataT%npoin,NPROC)
- if (myrank==0) then
- ! NOTE: output points lying at an interface between procs are assigned to a unique proc
- iproc = minloc(dist_all,2) - 1
- dataT%iglob = iglob_all(:,iproc)
- endif
- call bcast_all_i(iproc,dataT%npoin)
- call bcast_all_i(dataT%iglob,dataT%npoin)
-
- ! Number of output points contained in the current processor
- npoin_local = count( iproc == myrank )
-
- if (npoin_local>0) then
- ! Make a list of output points contained in the current processor
- allocate(glob_indx(npoin_local))
- ipoin_local = 0
- do ipoin = 1,dataT%npoin
- if (myrank == iproc(ipoin)) then
- ipoin_local = ipoin_local + 1
- glob_indx(ipoin_local) = ipoin
- endif
- enddo
- ! Consolidate the output information (remove output points outside current proc)
- allocate(iglob_tmp(dataT%npoin))
- allocate(name_tmp(dataT%npoin))
- iglob_tmp = dataT%iglob
- name_tmp = dataT%name
- deallocate(dataT%iglob)
- deallocate(dataT%name)
- dataT%npoin = npoin_local
- allocate(dataT%iglob(dataT%npoin))
- allocate(dataT%name(dataT%npoin))
- dataT%iglob = iglob_tmp(glob_indx)
- dataT%name = name_tmp(glob_indx)
- deallocate(glob_indx,iglob_tmp,name_tmp)
-
- else
- dataT%npoin = 0
- deallocate(dataT%iglob)
- deallocate(dataT%name)
- endif
-
- deallocate(iproc,iglob_all,dist_all)
- endif
-
- ! 3. allocate arrays and set to zero
- if (dataT%npoin>0) then
- 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
- endif
-
-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
-
- 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)
- 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
- 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
-1 format(I5)
- NTchar = adjustl(NTchar)
-
- do i=1,dataT%npoin
-
- 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-29 18:56:50 UTC (rev 20969)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_kinematic.f90 2012-10-30 05:30:39 UTC (rev 20970)
@@ -1,87 +1,47 @@
-!=====================================================================
+! This module implements kinematic faults: prescribed spatio-temporal slip history
!
-! s p e c f e m 3 d v e r s i o n 1 . 4
-! ---------------------------------------
-!
-! dimitri komatitsch and jeroen tromp
-! seismological laboratory - california institute of technology
-! (c) california institute of technology september 2006
-!
-! this program is free software; you can redistribute it and/or modify
-! it under the terms of the gnu general public license as published by
-! the free software foundation; either version 2 of the license, or
-! (at your option) any later version.
-!
-! this program is distributed in the hope that it will be useful,
-! but without any warranty; without even the implied warranty of
-! merchantability or fitness for a particular purpose. see the
-! gnu general public license for more details.
-!
-! you should have received a copy of the gnu general public license along
-! with this program; if not, write to the free software foundation, inc.,
-! 51 franklin street, fifth floor, boston, ma 02110-1301 usa.
-!
-!===============================================================================================================
+! Authors:
+! Percy Galvez, Jean-Paul Ampuero, Javier Ruiz, Surendra Somala
-! This module was written by:
-! Percy Galvez , Jean-Paul Ampuero and Javier Ruiz
-! based on fault_solver.f90
-
module fault_solver_kinematic
- use fault_solver_common
+ use fault_solver_common
- implicit none
+ implicit none
- include 'constants.h'
+ include 'constants.h'
- private
+ 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
+ integer :: npoin=0
+ real(kind=CUSTOM_REAL), dimension(:), pointer :: d1=>null(), d2=>null(), &
+ v1=>null(), v2=>null(), &
+ t1=>null(), t2=>null(), t3=>null(), &
+ xcoord=>null(), ycoord=>null(), zcoord=>null()
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
- end type bc_kinflt_type
-
- 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.
+ type, extends (fault_type) :: bc_kinflt_type
+ private
+ 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
+ end type bc_kinflt_type
- public :: BC_KINFLT_init, BC_KINFLT_set_all, SIMULATION_TYPE_KIN
+ 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
contains
@@ -190,7 +150,7 @@
bc%v_kin_t1 = 0e0_CUSTOM_REAL
bc%v_kin_t2 = 0e0_CUSTOM_REAL
- call init_dataT(bc%dataT,bc%coord,bc%nglob,NT,iflt)
+ call init_dataT(bc%dataT,bc%coord,bc%nglob,NT,dt,7,iflt)
call init_dataXZ(bc%dataXZ,bc%nglob)
endif
@@ -323,7 +283,7 @@
!-- OUTPUTS --
! write dataT every NTOUT time steps or at the end of simulation
- if ( mod(it,NTOUT) == 0 .or. it==NSTEP) call SCEC_write_dataT(bc%dataT,bc%dt,it)
+ if ( mod(it,NTOUT) == 0 .or. it==NSTEP) call SCEC_write_dataT(bc%dataT,it)
! write dataXZ every NSNAP time steps
! if ( mod(it,NSNAP) == 0) call write_dataXZ(bc,it,iflt)
@@ -390,225 +350,4 @@
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_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_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))
- 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 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,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)
- enddo
-
-end subroutine store_dataT
-
-!-----------------------------------------------------------------
-subroutine write_dataT_all(nt)
-
- integer, intent(in) :: nt
-
- integer :: i
-
- if (.not.allocated(faults)) return
- do i = 1,size(faults)
- call SCEC_write_dataT(faults(i)%dataT,faults(i)%dt,nt)
- enddo
-
-end subroutine write_dataT_all
-
-!------------------------------------------------------------------------
-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,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