[cig-commits] r20929 - in seismo/3D/FAULT_SOURCE: . branches/new_fault_db/EXAMPLES/splay_faults branches/new_fault_db/EXAMPLES/tohoku branches/new_fault_db/EXAMPLES/tpv102 branches/new_fault_db/EXAMPLES/tpv103 branches/new_fault_db/EXAMPLES/tpv16.crack branches/new_fault_db/EXAMPLES/tpv5.crack branches/new_fault_db/src
ampuero at geodynamics.org
ampuero at geodynamics.org
Thu Oct 25 11:26:38 PDT 2012
Author: ampuero
Date: 2012-10-25 11:26:37 -0700 (Thu, 25 Oct 2012)
New Revision: 20929
Modified:
seismo/3D/FAULT_SOURCE/TO_DO
seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/splay_faults/splay_faults.py
seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/tohoku/mesh_japan.py
seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/tpv102/TPV102.py
seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/tpv103/TPV103.py
seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/tpv16.crack/TPV16.py
seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/tpv5.crack/tpv5.py
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
seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/prepare_timerun.f90
Log:
factored add_BT; clean up fault solver init
Modified: seismo/3D/FAULT_SOURCE/TO_DO
===================================================================
--- seismo/3D/FAULT_SOURCE/TO_DO 2012-10-25 17:27:11 UTC (rev 20928)
+++ seismo/3D/FAULT_SOURCE/TO_DO 2012-10-25 18:26:37 UTC (rev 20929)
@@ -6,15 +6,17 @@
- in fault_solver.f90: parallelize dataT outputs
- set parallel version as default
++ improve documentation of mesh generation
+
++ fault_solver:
+ - many hard-coded ad hoc features need to be set instead as user options
+
+ rate-and-state friction:
- make it a user-friendly option (currently a flag in fault_solver)
+ cubit interface:
- - merge the common operations into a single function call (and modfiy the examples accordingly)
- consolidate save_fault_nodes_elements.py and functions.py
-+ factor common subroutines of dynamic and kinematic fault solvers into a lower level module
-
+ verify consistency of fault edge nodes in kinematic solver:
currently, these nodes are not split but they are included in the fault database
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/splay_faults/splay_faults.py
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/splay_faults/splay_faults.py 2012-10-25 17:27:11 UTC (rev 20928)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/splay_faults/splay_faults.py 2012-10-25 18:26:37 UTC (rev 20929)
@@ -163,7 +163,6 @@
#SAVING FAULT NODES AND ELEMENTS.
os.system('mkdir -p MESH')
########## FAULT A ##############################################################
-## NOTE : FACE FAULT REFERENCE : face_down.
Au = [8,15] #face_up
Ad = [6,7] #face_down
faultA = fault_input(1,Au,Ad)
@@ -205,6 +204,3 @@
cubit2specfem3d.export2SESAME('MESH')
# all files needed by SCOTCH are now in directory MESH
-
-
-
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/tohoku/mesh_japan.py
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/tohoku/mesh_japan.py 2012-10-25 17:27:11 UTC (rev 20928)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/tohoku/mesh_japan.py 2012-10-25 18:26:37 UTC (rev 20929)
@@ -21,23 +21,22 @@
cubit.cmd(open_mesh)
cubit.cmd('compress ids hex face edge node')
- # THIS IS RELLY IMPORTANT TO AVOID GAPS IN NODES indexing AND HAVING
- # Indices bigger that the total grid points and
- # avoid fault segmentation fault errors
+ # THIS IS REALLY IMPORTANT to avoid gaps in node indices
+ # and indices bigger than the total number of grid points
+ # and avoid segmentation fault errors
########### Slab elements and nodes ###############
######## SLAB ################################################
os.system('mkdir -p MESH')
-Au = [42] # A_up
-Ad = [36] # A_down
-
####### Fault opening #############################################
cubit.cmd('set node constraint off')
cubit.cmd('node in surf 42 move X 0 Y 0 Z 0.001') # delta = 2e-3 km. in Z direction.
cubit.cmd('node in surf 36 move X 0 Y 0 Z -0.001') # In general the shift should be normal to the fault.
####################################################################
+Au = [42] # A_up
+Ad = [36] # A_down
faultA = fault_input(1,Au,Ad)
## FOR THE BULK (Seismic wave propagation part for SESAME)
@@ -210,10 +209,7 @@
cubit2specfem3d.export2SESAME('MESH')
m2km()
-### this mesh is in km . It is more convinient to work in meters
-### changing nodes_coords_file to meters .
-
-
+### This mesh is in km. It is more convenient to work in meters.
+### Changing nodes_coords_file to meters.
-
# all files needed by SCOTCH are now in directory MESH
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/tpv102/TPV102.py
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/tpv102/TPV102.py 2012-10-25 17:27:11 UTC (rev 20928)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/tpv102/TPV102.py 2012-10-25 18:26:37 UTC (rev 20929)
@@ -45,20 +45,6 @@
cubit2specfem3d.export2SESAME('MESH')
-
-
-
-
Au = [8] # A_up
Ad = [3] # A_down
-
faultA = fault_input(1,Au,Ad)
-
-
-
-
-
-
-
-
-
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/tpv103/TPV103.py
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/tpv103/TPV103.py 2012-10-25 17:27:11 UTC (rev 20928)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/tpv103/TPV103.py 2012-10-25 18:26:37 UTC (rev 20929)
@@ -14,7 +14,6 @@
os.system('mkdir -p MESH')
-
xmin = [9,16]
xmax = [11,13]
ymin = [3]
@@ -44,19 +43,6 @@
cubit2specfem3d.export2SESAME('MESH')
-
-
-
-
Au = [8] # A_up
Ad = [3] # A_down
-
faultA = fault_input(1,Au,Ad)
-
-
-
-
-
-
-
-
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/tpv16.crack/TPV16.py
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/tpv16.crack/TPV16.py 2012-10-25 17:27:11 UTC (rev 20928)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/tpv16.crack/TPV16.py 2012-10-25 18:26:37 UTC (rev 20929)
@@ -18,14 +18,10 @@
os.system('mkdir -p MESH')
-
-
Au = [8] # A_up
Ad = [13] # A_down
-
faultA = fault_input(1,Au,Ad)
-
xmin = [4]
xmax = [6]
ymin = [17,10]
@@ -54,9 +50,3 @@
#### Export to SESAME format using cubit2specfem3d.py of GEOCUBIT
cubit2specfem3d.export2SESAME('MESH')
-
-
-
-
-
-
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/tpv5.crack/tpv5.py
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/tpv5.crack/tpv5.py 2012-10-25 17:27:11 UTC (rev 20928)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/EXAMPLES/tpv5.crack/tpv5.py 2012-10-25 18:26:37 UTC (rev 20929)
@@ -90,11 +90,6 @@
cubit.cmd("mesh volume 1")
#cubit.cmd("unmerge surface 2 3")
-########### SIDESETS (NOT USED ) ###############
-#fault_A_elements_up="sideset 200 surface 2"
-#fault_A_elements_down="sideset 201 surface 3"
-#cubit.cmd(fault_A_elements_up)
-#cubit.cmd(fault_A_elements_down)
########### Fault elements and nodes ###############
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-25 17:27:11 UTC (rev 20928)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90 2012-10-25 18:26:37 UTC (rev 20929)
@@ -84,20 +84,13 @@
fw=>null(), Vw=>null()
end type rsf_type
- type asperity_type
+ type, extends (fault_type) :: bc_dynflt_type
private
- real(kind=CUSTOM_REAL), dimension(:), pointer :: Fload=>null()
- end type asperity_type
-
- type, EXTENDS (fault_type) :: bc_dynflt_type
- private
- real(kind=CUSTOM_REAL), dimension(:,:), pointer :: T0=>null()
- real(kind=CUSTOM_REAL), dimension(:), pointer :: MU=>null()
+ real(kind=CUSTOM_REAL), dimension(:,:), pointer :: T0=>null(), 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()
- type(asperity_type), pointer :: asp => null()
logical :: allow_opening = .false. ! default : do not allow opening
end type bc_dynflt_type
@@ -134,13 +127,12 @@
! Minv inverse mass matrix
! dt global time step
!
-subroutine BC_DYNFLT_init(prname,Minv,DTglobal,nt,accel,vel)
+subroutine BC_DYNFLT_init(prname,Minv,DTglobal,nt,vel)
character(len=256), intent(in) :: prname ! 'proc***'
real(kind=CUSTOM_REAL), intent(in) :: Minv(:)
double precision, intent(in) :: DTglobal
integer, intent(in) :: nt
- real(kind=CUSTOM_REAL) :: accel(:,:)
real(kind=CUSTOM_REAL), intent(inout) :: vel(:,:)
real(kind=CUSTOM_REAL) :: dt
@@ -195,7 +187,7 @@
dt = real(DTglobal)
do iflt=1,nbfaults
read(IIN_PAR,nml=BEGIN_FAULT,end=100)
- call init_one_fault(faults(iflt),IIN_BIN,IIN_PAR,Minv,dt,nt,vel,iflt,accel)
+ call init_one_fault(faults(iflt),IIN_BIN,IIN_PAR,Minv,dt,nt,vel,iflt)
enddo
close(IIN_BIN)
close(IIN_PAR)
@@ -220,64 +212,21 @@
!---------------------------------------------------------------------
-subroutine init_one_fault(bc,IIN_BIN,IIN_PAR,Minv,dt_tmp,NT,vel,iflt,accel)
+subroutine init_one_fault(bc,IIN_BIN,IIN_PAR,Minv,dt,NT,vel,iflt)
- real(kind=CUSTOM_REAL) :: accel(:,:)
real(kind=CUSTOM_REAL), intent(inout) :: vel(:,:)
type(bc_dynflt_type), intent(inout) :: bc
real(kind=CUSTOM_REAL), intent(in) :: Minv(:)
integer, intent(in) :: IIN_BIN,IIN_PAR,NT,iflt
- real(kind=CUSTOM_REAL), intent(in) :: dt_tmp
+ real(kind=CUSTOM_REAL), intent(in) :: dt
- real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: jacobian2Dw
- real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: normal
- integer, dimension(:,:), allocatable :: ibool1
- real(kind=CUSTOM_REAL) :: norm
real(kind=CUSTOM_REAL) :: S1,S2,S3
integer :: n1,n2,n3
- real(kind=CUSTOM_REAL) :: mus,mud,dc
- integer :: nmus,nmud,ndc,ij,k,e
- real(kind=CUSTOM_REAL), dimension(:), allocatable :: nx,ny,nz
- real(kind=CUSTOM_REAL) :: V0,f0,a,b,L,theta,theta_init,V_init,fw,Vw
- integer :: nV0,nf0,na,nb,nL,nV_init,ntheta_init,nfw,nVw
- real(kind=CUSTOM_REAL) :: AspTx,Fload
- integer :: nAsp,nFload
- real(kind=CUSTOM_REAL) :: C,T
- integer :: nC,nForcedRup
- integer, parameter :: IIN_NUC =270
- integer :: ipar
- integer :: ier
-
- integer :: relz_num,sub_relz_num
- integer :: num_cell_str,num_cell_dip
- real(kind=CUSTOM_REAL) :: siz_str,siz_dip
- integer :: hypo_cell_str,hypo_cell_dip
- real(kind=CUSTOM_REAL) :: hypo_loc_str,hypo_loc_dip,rad_T_str,rad_T_dip
-
- integer, dimension(:), allocatable :: inp_nx,inp_nz
- real(kind=CUSTOM_REAL), dimension(:), allocatable :: loc_str,loc_dip,sigma0,tau0_str,tau0_dip,Rstress_str,Rstress_dip,static_fc, &
- dyn_fc,swcd,cohes,tim_forcedRup
-
- real(kind=CUSTOM_REAL) :: minX
-
- real(kind=CUSTOM_REAL) :: W1,W2,w,hypo_z
- real(kind=CUSTOM_REAL) :: x,z
- logical :: c1,c2,c3,c4
- real(kind=CUSTOM_REAL) :: b11,b12,b21,b22,B1,B2
- integer :: i,nglob_bulk
-
- double precision, dimension(:), allocatable :: mu1,mu2,mu3
- real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: init_vel
-
-
NAMELIST / INIT_STRESS / S1,S2,S3,n1,n2,n3
- NAMELIST / SWF / mus,mud,dc,nmus,nmud,ndc,C,T,nC,nForcedRup
- 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
- call initialize_fault(bc,IIN_BIN,Minv,dt_tmp)
+ call initialize_fault(bc,IIN_BIN,Minv,dt)
if (bc%nspec>0) then
@@ -300,7 +249,6 @@
bc%T0(1,:) = S1
bc%T0(2,:) = S2
bc%T0(3,:) = S3
-
call init_2d_distribution(bc%T0(1,:),bc%coord,IIN_PAR,n1)
call init_2d_distribution(bc%T0(2,:),bc%coord,IIN_PAR,n2)
call init_2d_distribution(bc%T0(3,:),bc%coord,IIN_PAR,n3)
@@ -310,263 +258,75 @@
! if (abs(bc%zcoord(k)-0.e0_CUSTOM_REAL) <= SMALLVAL) bc%T0(2,k) = 0
! end do
- if (TPV16) then
-
- allocate(inp_nx(bc%nglob))
- allocate(inp_nz(bc%nglob))
- allocate(loc_str(bc%nglob))
- allocate(loc_dip(bc%nglob))
- allocate(sigma0(bc%nglob))
- allocate(tau0_str(bc%nglob))
- allocate(tau0_dip(bc%nglob))
- allocate(Rstress_str(bc%nglob))
- allocate(Rstress_dip(bc%nglob))
- allocate(static_fc(bc%nglob))
- allocate(dyn_fc(bc%nglob))
- allocate(swcd(bc%nglob))
- allocate(cohes(bc%nglob))
- allocate(tim_forcedRup(bc%nglob))
-
- open(unit=IIN_NUC,file='DATA/FAULT/input_file.txt',status='old',iostat=ier)
- read(IIN_NUC,*) relz_num,sub_relz_num
- read(IIN_NUC,*) num_cell_str,num_cell_dip,siz_str,siz_dip
- read(IIN_NUC,*) hypo_cell_str,hypo_cell_dip,hypo_loc_str,hypo_loc_dip,rad_T_str,rad_T_dip
- do ipar=1,bc%nglob
- read(IIN_NUC,*) inp_nx(ipar),inp_nz(ipar),loc_str(ipar),loc_dip(ipar),sigma0(ipar),tau0_str(ipar),tau0_dip(ipar), &
- Rstress_str(ipar),Rstress_dip(ipar),static_fc(ipar),dyn_fc(ipar),swcd(ipar),cohes(ipar),tim_forcedRup(ipar)
- enddo
- close(IIN_NUC)
-
- allocate( bc%swf )
- allocate( bc%swf%mus(bc%nglob) )
- allocate( bc%swf%mud(bc%nglob) )
- allocate( bc%swf%Dc(bc%nglob) )
- allocate( bc%swf%theta(bc%nglob) )
- allocate( bc%swf%T(bc%nglob) )
- allocate( bc%swf%C(bc%nglob) )
-
- minX = minval(bc%coord(1,:))
-
- do i=1,bc%nglob
-
- !loc_dip is negative of Z-coord
- ipar = minloc( (minX+loc_str(:)-bc%coord(1,i))**2 + (-loc_dip(:)-bc%coord(3,i))**2 , 1)
- bc%T0(3,i) = -sigma0(ipar)
- bc%T0(1,i) = tau0_str(ipar)
- bc%T0(2,i) = tau0_dip(ipar)
-
- bc%swf%mus(i) = static_fc(ipar)
- bc%swf%mud(i) = dyn_fc(ipar)
- bc%swf%Dc(i) = swcd(ipar)
- bc%swf%C(i) = cohes(ipar)
- bc%swf%T(i) = tim_forcedRup(ipar)
- enddo
-
- endif
-
-
! Set friction parameters and initialize friction variables
- ! Slip weakening friction
- if(.not. RATE_AND_STATE) then
- allocate( bc%swf )
- allocate( bc%swf%mus(bc%nglob) )
- allocate( bc%swf%mud(bc%nglob) )
- allocate( bc%swf%Dc(bc%nglob) )
- allocate( bc%swf%theta(bc%nglob) )
- allocate( bc%swf%C(bc%nglob) )
- allocate( bc%swf%T(bc%nglob) )
- ! WARNING: if V_HEALING is negative we turn off healing
- bc%swf%healing = (V_HEALING > 0e0_CUSTOM_REAL)
-
- mus = 0.6e0_CUSTOM_REAL
- mud = 0.1e0_CUSTOM_REAL
- dc = 1e0_CUSTOM_REAL
- nmus = 0
- nmud = 0
- ndc = 0
-
- C = 0._CUSTOM_REAL
- T = HUGEVAL
- nC = 0
- nForcedRup = 0
-
- read(IIN_PAR, nml=SWF)
- bc%swf%mus = mus
- bc%swf%mud = mud
- bc%swf%Dc = dc
- bc%swf%C = C
- bc%swf%T = T
- call init_2d_distribution(bc%swf%mus,bc%coord,IIN_PAR,nmus)
- call init_2d_distribution(bc%swf%mud,bc%coord,IIN_PAR,nmud)
- call init_2d_distribution(bc%swf%Dc ,bc%coord,IIN_PAR,ndc)
- call init_2d_distribution(bc%swf%C ,bc%coord,IIN_PAR,nC)
- call init_2d_distribution(bc%swf%T ,bc%coord,IIN_PAR,nForcedRup)
-
- bc%swf%theta = 0e0_CUSTOM_REAL
- allocate(bc%MU(bc%nglob))
- bc%MU = swf_mu(bc%swf)
-
- ! Rate and state friction
+ allocate(bc%MU(bc%nglob))
+ if (RATE_AND_STATE) then
+ allocate(bc%rsf)
+ call rsf_init(bc%rsf,bc%Fload,bc%coord,IIN_PAR)
+ ! WARNING: the line below is only valid for pure strike-slip faulting
+ bc%V(1,:) = f%V_init
else
- allocate( bc%rsf )
- allocate( bc%rsf%V0(bc%nglob) )
- allocate( bc%rsf%f0(bc%nglob) )
- allocate( bc%rsf%a(bc%nglob) )
- allocate( bc%rsf%b(bc%nglob) )
- allocate( bc%rsf%L(bc%nglob) )
- allocate( bc%rsf%V_init(bc%nglob) )
- allocate( bc%rsf%theta(bc%nglob) )
- allocate( bc%rsf%C(bc%nglob) )
- allocate( bc%rsf%T(bc%nglob) )
- allocate( bc%rsf%fw(bc%nglob) )
- allocate( bc%rsf%Vw(bc%nglob) )
+ allocate(bc%swf)
+ call swf_init(bc%swf,bc%MU,bc%coord,IIN_PAR)
+ if (TPV16) call TPV16_init() !WARNING: ad hoc, initializes T0 and swf
+ endif
- V0 =1.e-6_CUSTOM_REAL
- f0 =0.6_CUSTOM_REAL
- a =0.0080_CUSTOM_REAL !0.0080_CUSTOM_REAL
- b =0.0040_CUSTOM_REAL !0.0120_CUSTOM_REAL
- L =0.0135_CUSTOM_REAL
- V_init =1.e-12_CUSTOM_REAL
- theta_init =1.084207680000000e+09_CUSTOM_REAL
+ endif
- nV0 =0
- nf0 =0
- na =0
- nb =0
- nL =0
- nV_init =0
- ntheta_init =0
+! call init_dataT(bc%dataT,bc%coord,bc%nglob,NT,iflt)
+ call init_dataXZ(bc%dataXZ,bc)
- C = 0._CUSTOM_REAL
- T = HUGEVAL
- nC = 0
- nForcedRup = 0
+!--------------------------------------------------------
+contains
- fw = 0.2_CUSTOM_REAL
- Vw = 0.1_CUSTOM_REAL
- nfw =0
- nVw =0
+subroutine TPV16_init
- read(IIN_PAR, nml=RSF)
- bc%rsf%V0 = V0
- bc%rsf%f0 = f0
- bc%rsf%a = a
- bc%rsf%b = b
- bc%rsf%L = L
- bc%rsf%V_init = V_init
- bc%rsf%theta = theta_init
- bc%rsf%C = C
- bc%rsf%T = T
- bc%rsf%fw = fw
- bc%rsf%Vw = Vw
- call init_2d_distribution(bc%rsf%V0,bc%coord,IIN_PAR,nV0)
- call init_2d_distribution(bc%rsf%f0,bc%coord,IIN_PAR,nf0)
- call init_2d_distribution(bc%rsf%a,bc%coord,IIN_PAR,na)
- call init_2d_distribution(bc%rsf%b,bc%coord,IIN_PAR,nb)
- call init_2d_distribution(bc%rsf%L,bc%coord,IIN_PAR,nL)
- call init_2d_distribution(bc%rsf%V_init,bc%coord,IIN_PAR,nV_init)
- call init_2d_distribution(bc%rsf%theta,bc%coord,IIN_PAR,ntheta_init)
- call init_2d_distribution(bc%rsf%C ,bc%coord,IIN_PAR,nC)
- call init_2d_distribution(bc%rsf%T ,bc%coord,IIN_PAR,nForcedRup)
- call init_2d_distribution(bc%rsf%fw ,bc%coord,IIN_PAR,nfw)
- call init_2d_distribution(bc%rsf%Vw ,bc%coord,IIN_PAR,nVw)
+ integer :: ier, ipar
+ integer, parameter :: IIN_NUC =270 ! WARNING: not safe, should look for an available unit
+ real(kind=CUSTOM_REAL), dimension(bc%nglob) :: loc_str,loc_dip,sigma0,tau0_str,tau0_dip,Rstress_str,Rstress_dip,static_fc, &
+ dyn_fc,swcd,cohes,tim_forcedRup
+ integer, dimension(bc%nglob) :: inp_nx,inp_nz
+ real(kind=CUSTOM_REAL) :: minX, siz_str,siz_dip, hypo_loc_str,hypo_loc_dip,rad_T_str,rad_T_dip
+ integer :: relz_num,sub_relz_num, num_cell_str,num_cell_dip, hypo_cell_str,hypo_cell_dip
-!!$ ! WARNING : Not general enough
-!!$ nglob_bulk = size(vel,2)
-!!$ allocate(init_vel(3,nglob_bulk))
-!!$ init_vel = 0._CUSTOM_REAL
-!!$ init_vel(1,bc%ibulk1) = -bc%rsf%V_init/2._CUSTOM_REAL
-!!$ init_vel(1,bc%ibulk2) = bc%rsf%V_init/2._CUSTOM_REAL
-!!$ where(ystore > 0) init_vel(1,:) = -V_init/2._CUSTOM_REAL
-!!$ where(ystore < 0) init_vel(1,:) = V_init/2._CUSTOM_REAL
-!!$ !init_vel = rotate(bc,init_vel,-1) ! directly assigned in global coordinates here as it is a simplified case
-!!$ vel = vel + init_vel
-
- !WARNING: below is an ad-hoc setting of a(x,z) for some SCEC benchmark
- ! This should be instead an option in init_2d_distribution
- W1=15000._CUSTOM_REAL
- W2=7500._CUSTOM_REAL
- w=3000._CUSTOM_REAL
- hypo_z = -7500._CUSTOM_REAL
- do i=1,bc%nglob
- x=bc%coord(1,i)
- z=bc%coord(3,i)
- c1=abs(x)<W1+w
- c2=abs(x)>W1
- c3=abs(z-hypo_z)<W2+w
- c4=abs(z-hypo_z)>W2
- if( (c1 .and. c2 .and. c3) .or. (c3 .and. c4 .and. c1) ) then
+ open(unit=IIN_NUC,file='DATA/FAULT/input_file.txt',status='old',iostat=ier)
+ read(IIN_NUC,*) relz_num,sub_relz_num
+ read(IIN_NUC,*) num_cell_str,num_cell_dip,siz_str,siz_dip
+ read(IIN_NUC,*) hypo_cell_str,hypo_cell_dip,hypo_loc_str,hypo_loc_dip,rad_T_str,rad_T_dip
+ do ipar=1,bc%nglob
+ read(IIN_NUC,*) inp_nx(ipar),inp_nz(ipar),loc_str(ipar),loc_dip(ipar),sigma0(ipar),tau0_str(ipar),tau0_dip(ipar), &
+ Rstress_str(ipar),Rstress_dip(ipar),static_fc(ipar),dyn_fc(ipar),swcd(ipar),cohes(ipar),tim_forcedRup(ipar)
+ enddo
+ close(IIN_NUC)
+
+ minX = minval(bc%coord(1,:))
- if (c1 .and. c2) then
- b11 = w/(abs(x)-W1-w)
- b12 = w/(abs(x)-W1)
- B1 = HALF * (ONE + tanh(b11 + b12))
- elseif(abs(x)<=W1) then
- B1 = 1._CUSTOM_REAL
- else
- B1 = 0._CUSTOM_REAL
- endif
-
- if (c3 .and. c4) then
- b21 = w/(abs(z-hypo_z)-W2-w)
- b22 = w/(abs(z-hypo_z)-W2)
- B2 = HALF * (ONE + tanh(b21 + b22))
- elseif(abs(z-hypo_z)<=W2) then
- B2 = 1._CUSTOM_REAL
- else
- B2 = 0._CUSTOM_REAL
- endif
-
- bc%rsf%a(i) = 0.008 + 0.008 * (ONE - B1*B2)
- bc%rsf%Vw(i) = 0.1 + 0.9 * (ONE - B1*B2)
+ do i=1,bc%nglob
+
+ ! WARNING: nearest neighbor interpolation
+ ipar = minloc( (minX+loc_str(:)-bc%coord(1,i))**2 + (-loc_dip(:)-bc%coord(3,i))**2 , 1)
+ !loc_dip is negative of Z-coord
- elseif( abs(x)<=W1 .and. abs(z-hypo_z)<=W2 ) then
- bc%rsf%a(i) = 0.008
- bc%rsf%Vw(i) = 0.1_CUSTOM_REAL
- else
- bc%rsf%a(i) = 0.016
- bc%rsf%Vw(i) = 1.0_CUSTOM_REAL
- endif
+ bc%T0(3,i) = -sigma0(ipar)
+ bc%T0(1,i) = tau0_str(ipar)
+ bc%T0(2,i) = tau0_dip(ipar)
- enddo
+ bc%swf%mus(i) = static_fc(ipar)
+ bc%swf%mud(i) = dyn_fc(ipar)
+ bc%swf%Dc(i) = swcd(ipar)
+ bc%swf%C(i) = cohes(ipar)
+ bc%swf%T(i) = tim_forcedRup(ipar)
- ! WARNING: The line below scratches an earlier initialization of theta through theta_init
- ! We should implement it as an option for the user
- if(bc%rsf%stateLaw == 1) then
- bc%rsf%theta = bc%rsf%L/bc%rsf%V0 &
- * exp( ( bc%rsf%a * log(TWO*sinh(-sqrt(bc%T0(1,:)**2+bc%T0(2,:)**2)/bc%T0(3,:)/bc%rsf%a)) &
- - bc%rsf%f0 - bc%rsf%a*log(bc%rsf%V_init/bc%rsf%V0) ) &
- / bc%rsf%b )
- else
- bc%rsf%theta = bc%rsf%a * log(TWO*V0/V_init * sinh(-sqrt(bc%T0(1,:)**2+bc%T0(2,:)**2)/bc%T0(3,:)/bc%rsf%a))
- endif
+ enddo
- allocate(bc%MU(bc%nglob))
+subroutine TPV16_init
- ! WARNING: the line below is only valid for pure strike-slip faulting
- bc%V(1,:) = bc%rsf%V_init
-
- allocate( bc%asp )
- allocate( bc%asp%Fload(bc%nglob) )
- Fload =0.e0_CUSTOM_REAL
- nFload =0
- read(IIN_PAR, nml=ASP)
- bc%asp%Fload =Fload
- call init_2d_distribution(bc%asp%Fload,bc%coord,IIN_PAR,nFload)
-
- endif
- endif
-
-
-! call init_dataT(bc%dataT,bc%coord,bc%nglob,NT,iflt)
- call init_dataXZ(bc%dataXZ,bc)
-
end subroutine init_one_fault
!---------------------------------------------------------------------
-! adds a value to a fault parameter inside an area with prescribed shape
+! REPLACES the value of 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
+!JPA refactor: background value should be an argument
real(kind=CUSTOM_REAL), intent(inout) :: a(:)
real(kind=CUSTOM_REAL), intent(in) :: coord(:,:)
@@ -630,6 +390,7 @@
stop 'bc_dynflt_3d::init_2d_distribution:: unknown shape'
end select
+ ! REPLACE the value inside the prescribed area
where (b /= 0e0_CUSTOM_REAL) a = b
enddo
@@ -672,7 +433,7 @@
!---------------------------------------------------------------------
subroutine BC_DYNFLT_set3d(bc,MxA,V,D,iflt)
- use specfem_par, only:it,NSTEP,myrank
+ use specfem_par, only: it,NSTEP,myrank
real(kind=CUSTOM_REAL), intent(inout) :: MxA(:,:)
type(bc_dynflt_type), intent(inout) :: bc
@@ -682,7 +443,7 @@
real(kind=CUSTOM_REAL), dimension(3,bc%nglob) :: T,dD,dV,dA
real(kind=CUSTOM_REAL), dimension(bc%nglob) :: strength,tStick,tnew, &
theta_old, theta_new, dc, &
- ta,Vf_old,Vf_new,tau1,TxExt
+ ta,Vf_old,Vf_new,TxExt
real(kind=CUSTOM_REAL) :: half_dt,TLoad,DTau0,GLoad,time
integer :: i
@@ -732,7 +493,7 @@
else
GLoad = 1.0_CUSTOM_REAL
endif
- TxExt = DTau0 * bc%asp%Fload * GLoad
+ TxExt = DTau0 * bc%Fload * GLoad
T(1,:) = T(1,:) + TxExt
endif
@@ -808,14 +569,8 @@
T = rotate(bc,T,-1)
! Add boundary term B*T to M*a
- MxA(1,bc%ibulk1) = MxA(1,bc%ibulk1) + bc%B*T(1,:)
- MxA(2,bc%ibulk1) = MxA(2,bc%ibulk1) + bc%B*T(2,:)
- MxA(3,bc%ibulk1) = MxA(3,bc%ibulk1) + bc%B*T(3,:)
+ call add_BT(bc,MxA,T)
- MxA(1,bc%ibulk2) = MxA(1,bc%ibulk2) - bc%B*T(1,:)
- MxA(2,bc%ibulk2) = MxA(2,bc%ibulk2) - bc%B*T(2,:)
- MxA(3,bc%ibulk2) = MxA(3,bc%ibulk2) - bc%B*T(3,:)
-
!-- intermediate storage of outputs --
Vf_new = sqrt(bc%V(1,:)*bc%V(1,:)+bc%V(2,:)*bc%V(2,:))
if(.not. RATE_AND_STATE) then
@@ -858,6 +613,62 @@
!===============================================================
+subroutine swf_init(f,mu,coord,IIN_PAR)
+
+ type(swf_type), intent(out) :: f
+ real(kind=CUSTOM_REAL), intent(out) :: mu(:)
+ real(kind=CUSTOM_REAL), intent(in) :: coord(:,:)
+ integer, intent(in) :: IIN_PAR
+
+ integer :: nglob
+ real(kind=CUSTOM_REAL) :: mus,mud,dc,C,T
+ integer :: nmus,nmud,ndc,nC,nForcedRup
+
+ NAMELIST / SWF / mus,mud,dc,nmus,nmud,ndc,C,T,nC,nForcedRup
+
+ nglob = size(mu)
+ allocate( f%mus(nglob) )
+ allocate( f%mud(nglob) )
+ allocate( f%Dc(nglob) )
+ allocate( f%theta(nglob) )
+ allocate( f%C(nglob) )
+ allocate( f%T(nglob) )
+
+ ! WARNING: if V_HEALING is negative we turn off healing
+ f%healing = (V_HEALING > 0e0_CUSTOM_REAL)
+
+ mus = 0.6e0_CUSTOM_REAL
+ mud = 0.1e0_CUSTOM_REAL
+ dc = 1e0_CUSTOM_REAL
+ C = 0._CUSTOM_REAL
+ T = HUGEVAL
+
+ nmus = 0
+ nmud = 0
+ ndc = 0
+ nC = 0
+ nForcedRup = 0
+
+ read(IIN_PAR, nml=SWF)
+
+ f%mus = mus
+ f%mud = mud
+ f%Dc = dc
+ f%C = C
+ f%T = T
+ call init_2d_distribution(f%mus,coord,IIN_PAR,nmus)
+ call init_2d_distribution(f%mud,coord,IIN_PAR,nmud)
+ call init_2d_distribution(f%Dc ,coord,IIN_PAR,ndc)
+ call init_2d_distribution(f%C ,coord,IIN_PAR,nC)
+ call init_2d_distribution(f%T ,coord,IIN_PAR,nForcedRup)
+
+ f%theta = 0e0_CUSTOM_REAL
+
+ mu = swf_mu(f)
+
+subroutine swf_init
+
+!---------------------------------------------------------------------
subroutine swf_update_state(dold,dnew,vold,f)
real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: vold,dold,dnew
@@ -889,6 +700,171 @@
end function swf_mu
!=====================================================================
+
+subroutine rsf_init(f,nucFload,coord,IIN_PAR)
+
+ type(rsf_type), intent(out) :: f
+ real(kind=CUSTOM_REAL), intent(in) :: coord(:,:)
+ real(kind=CUSTOM_REAL), pointer :: nucFload(:,:)
+ integer, intent(in) :: IIN_PAR
+
+ real(kind=CUSTOM_REAL) :: V0,f0,a,b,L,theta,theta_init,V_init,fw,Vw, C,T
+ integer :: nV0,nf0,na,nb,nL,nV_init,ntheta_init,nfw,nVw, nC,nForcedRup
+ real(kind=CUSTOM_REAL) :: W1,W2,w,hypo_z
+ real(kind=CUSTOM_REAL) :: x,z
+ logical :: c1,c2,c3,c4
+ real(kind=CUSTOM_REAL) :: b11,b12,b21,b22,B1,B2
+ integer :: i !,nglob_bulk
+ real(kind=CUSTOM_REAL) :: Fload
+ integer :: nFload
+! real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: init_vel
+
+ 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
+
+ allocate( f%V0(bc%nglob) )
+ allocate( f%f0(bc%nglob) )
+ allocate( f%a(bc%nglob) )
+ allocate( f%b(bc%nglob) )
+ allocate( f%L(bc%nglob) )
+ allocate( f%V_init(bc%nglob) )
+ allocate( f%theta(bc%nglob) )
+ allocate( f%C(bc%nglob) )
+ allocate( f%T(bc%nglob) )
+ allocate( f%fw(bc%nglob) )
+ allocate( f%Vw(bc%nglob) )
+
+ V0 =1.e-6_CUSTOM_REAL
+ f0 =0.6_CUSTOM_REAL
+ a =0.0080_CUSTOM_REAL !0.0080_CUSTOM_REAL
+ b =0.0040_CUSTOM_REAL !0.0120_CUSTOM_REAL
+ L =0.0135_CUSTOM_REAL
+ V_init =1.e-12_CUSTOM_REAL
+ theta_init =1.084207680000000e+09_CUSTOM_REAL
+ C = 0._CUSTOM_REAL
+ T = HUGEVAL
+ fw = 0.2_CUSTOM_REAL
+ Vw = 0.1_CUSTOM_REAL
+
+ nV0 =0
+ nf0 =0
+ na =0
+ nb =0
+ nL =0
+ nV_init =0
+ ntheta_init =0
+ nC = 0
+ nForcedRup = 0
+ nfw =0
+ nVw =0
+
+ read(IIN_PAR, nml=RSF)
+
+ f%V0 = V0
+ f%f0 = f0
+ f%a = a
+ f%b = b
+ f%L = L
+ f%V_init = V_init
+ f%theta = theta_init
+ f%C = C
+ f%T = T
+ f%fw = fw
+ f%Vw = Vw
+
+ call init_2d_distribution(f%V0,coord,IIN_PAR,nV0)
+ call init_2d_distribution(f%f0,coord,IIN_PAR,nf0)
+ call init_2d_distribution(f%a,coord,IIN_PAR,na)
+ call init_2d_distribution(f%b,coord,IIN_PAR,nb)
+ call init_2d_distribution(f%L,coord,IIN_PAR,nL)
+ call init_2d_distribution(f%V_init,coord,IIN_PAR,nV_init)
+ call init_2d_distribution(f%theta,coord,IIN_PAR,ntheta_init)
+ call init_2d_distribution(f%C,coord,IIN_PAR,nC)
+ call init_2d_distribution(f%T,coord,IIN_PAR,nForcedRup)
+ call init_2d_distribution(f%fw,coord,IIN_PAR,nfw)
+ call init_2d_distribution(f%Vw,coord,IIN_PAR,nVw)
+
+!!$ ! WARNING : Not general enough
+!!$ nglob_bulk = size(vel,2)
+!!$ allocate(init_vel(3,nglob_bulk))
+!!$ init_vel = 0._CUSTOM_REAL
+!!$ init_vel(1,bc%ibulk1) = -f%V_init/2._CUSTOM_REAL
+!!$ init_vel(1,bc%ibulk2) = f%V_init/2._CUSTOM_REAL
+!!$ where(ystore > 0) init_vel(1,:) = -V_init/2._CUSTOM_REAL
+!!$ where(ystore < 0) init_vel(1,:) = V_init/2._CUSTOM_REAL
+!!$ !init_vel = rotate(bc,init_vel,-1) ! directly assigned in global coordinates here as it is a simplified case
+!!$ vel = vel + init_vel
+
+ ! WARNING: below is an ad-hoc setting of a(x,z) for some SCEC benchmark
+ ! This should be instead an option in init_2d_distribution
+ W1=15000._CUSTOM_REAL
+ W2=7500._CUSTOM_REAL
+ w=3000._CUSTOM_REAL
+ hypo_z = -7500._CUSTOM_REAL
+ do i=1,bc%nglob
+ x=coord(1,i)
+ z=coord(3,i)
+ c1=abs(x)<W1+w
+ c2=abs(x)>W1
+ c3=abs(z-hypo_z)<W2+w
+ c4=abs(z-hypo_z)>W2
+ if( (c1 .and. c2 .and. c3) .or. (c3 .and. c4 .and. c1) ) then
+
+ if (c1 .and. c2) then
+ b11 = w/(abs(x)-W1-w)
+ b12 = w/(abs(x)-W1)
+ B1 = HALF * (ONE + tanh(b11 + b12))
+ elseif(abs(x)<=W1) then
+ B1 = 1._CUSTOM_REAL
+ else
+ B1 = 0._CUSTOM_REAL
+ endif
+
+ if (c3 .and. c4) then
+ b21 = w/(abs(z-hypo_z)-W2-w)
+ b22 = w/(abs(z-hypo_z)-W2)
+ B2 = HALF * (ONE + tanh(b21 + b22))
+ elseif(abs(z-hypo_z)<=W2) then
+ B2 = 1._CUSTOM_REAL
+ else
+ B2 = 0._CUSTOM_REAL
+ endif
+
+ f%a(i) = 0.008 + 0.008 * (ONE - B1*B2)
+ f%Vw(i) = 0.1 + 0.9 * (ONE - B1*B2)
+
+ elseif( abs(x)<=W1 .and. abs(z-hypo_z)<=W2 ) then
+ f%a(i) = 0.008
+ f%Vw(i) = 0.1_CUSTOM_REAL
+ else
+ f%a(i) = 0.016
+ f%Vw(i) = 1.0_CUSTOM_REAL
+ endif
+
+ enddo
+
+ ! WARNING: The line below scratches an earlier initialization of theta through theta_init
+ ! We should implement it as an option for the user
+ if(f%stateLaw == 1) then
+ f%theta = f%L/f%V0 &
+ * exp( ( f%a * log(TWO*sinh(-sqrt(bc%T0(1,:)**2+bc%T0(2,:)**2)/bc%T0(3,:)/f%a)) &
+ - f%f0 - f%a*log(f%V_init/f%V0) ) &
+ / f%b )
+ else
+ f%theta = f%a * log(TWO*V0/V_init * sinh(-sqrt(bc%T0(1,:)**2+bc%T0(2,:)**2)/bc%T0(3,:)/f%a))
+ endif
+
+ ! WARNING : ad hoc for SCEC benchmark TPV10x
+ allocate( nucFload(bc%nglob) )
+ Fload = 0.e0_CUSTOM_REAL
+ nFload = 0
+ read(IIN_PAR, nml=ASP)
+ nucFload = Fload
+ call init_2d_distribution(nucFload,coord,IIN_PAR,nFload)
+
+end subroutine rsf_init
+
+!---------------------------------------------------------------------
! Rate and state friction coefficient
function rsf_mu(f,V) result(mu)
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-25 17:27:11 UTC (rev 20928)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_common.f90 2012-10-25 18:26:37 UTC (rev 20929)
@@ -21,24 +21,23 @@
integer, dimension(:), pointer :: poin_offset=>null(),npoin_perproc=>null()
end type fault_type
-
-
logical, save :: PARALLEL_FAULT = .false.
- public :: initialize_fault,compute_R,get_jump,get_weighted_jump,rotate,PARALLEL_FAULT,fault_type
+ public :: fault_type, PARALLEL_FAULT, &
+ initialize_fault, get_jump, get_weighted_jump, rotate, add_BT
contains
!---------------------------------------------------------------------
-subroutine initialize_fault (bc,IIN_BIN,Minv,dt_tmp)
+subroutine initialize_fault (bc,IIN_BIN,Minv,dt)
use specfem_par
class(fault_type), intent(inout) :: bc
real(kind=CUSTOM_REAL), intent(in) :: Minv(:)
integer, intent(in) :: IIN_BIN
- real(kind=CUSTOM_REAL), intent(in) :: dt_tmp
+ real(kind=CUSTOM_REAL), intent(in) :: dt
real(kind=CUSTOM_REAL) :: tmp_vec(3,NGLOB_AB)
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: jacobian2Dw
@@ -75,7 +74,7 @@
read(IIN_BIN) bc%coord(2,:)
read(IIN_BIN) bc%coord(3,:)
- bc%dt = dt_tmp
+ bc%dt = dt
bc%B = 0e0_CUSTOM_REAL
allocate(nx(bc%nglob),ny(bc%nglob),nz(bc%nglob))
@@ -141,7 +140,7 @@
! 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
- bc%Z = 1.e0_CUSTOM_REAL/(0.5e0_CUSTOM_REAL*dt_tmp * bc%B *( bc%invM1 + bc%invM2 ))
+ 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
@@ -246,5 +245,20 @@
!----------------------------------------------------------------------
+subroutine add_BT(bc,MxA,T)
+ class(fault_type), intent(in) :: bc
+ real(kind=CUSTOM_REAL), intent(inout) :: MxA(:,:)
+ real(kind=CUSTOM_REAL), dimension(3,bc%nglob), intent(in) :: T
+
+ MxA(1,bc%ibulk1) = MxA(1,bc%ibulk1) + bc%B*T(1,:)
+ MxA(2,bc%ibulk1) = MxA(2,bc%ibulk1) + bc%B*T(2,:)
+ MxA(3,bc%ibulk1) = MxA(3,bc%ibulk1) + bc%B*T(3,:)
+
+ MxA(1,bc%ibulk2) = MxA(1,bc%ibulk2) - bc%B*T(1,:)
+ MxA(2,bc%ibulk2) = MxA(2,bc%ibulk2) - bc%B*T(2,:)
+ MxA(3,bc%ibulk2) = MxA(3,bc%ibulk2) - bc%B*T(3,:)
+
+end subroutine add_BT
+
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-25 17:27:11 UTC (rev 20928)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver_kinematic.f90 2012-10-25 18:26:37 UTC (rev 20929)
@@ -65,12 +65,12 @@
end type dataXZ_type
- type, EXTENDS (fault_type) :: bc_kinflt_type
+ type, extends (fault_type) :: bc_kinflt_type
private
- type(dataT_type) :: dataT
- type(dataXZ_type) :: dataXZ,dataXZ_all
+ type(dataT_type) :: dataT
+ type(dataXZ_type) :: dataXZ,dataXZ_all
real(kind=CUSTOM_REAL) :: kin_dt
- integer :: kin_it
+ integer :: kin_it
real(kind=CUSTOM_REAL), dimension(:,:), pointer :: v_kin_t1,v_kin_t2
end type bc_kinflt_type
@@ -163,12 +163,6 @@
integer, intent(in) :: IIN_BIN,IIN_PAR,NT,iflt
real(kind=CUSTOM_REAL), intent(in) :: dt
- real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: jacobian2Dw
- real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: normal
- integer, dimension(:,:), allocatable :: ibool1
- real(kind=CUSTOM_REAL) :: norm
- integer :: ij,k,e
- real(kind=CUSTOM_REAL), dimension(:), allocatable :: nx,ny,nz
real(kind=CUSTOM_REAL) :: kindt
NAMELIST / KINPAR / kindt
@@ -322,14 +316,8 @@
T = rotate(bc,T,-1)
! Add boundary term B*T to M*a
- MxA(1,bc%ibulk1) = MxA(1,bc%ibulk1) + bc%B*T(1,:)
- MxA(2,bc%ibulk1) = MxA(2,bc%ibulk1) + bc%B*T(2,:)
- MxA(3,bc%ibulk1) = MxA(3,bc%ibulk1) + bc%B*T(3,:)
+ call add_BT(bc,MxA,T)
- MxA(1,bc%ibulk2) = MxA(1,bc%ibulk2) - bc%B*T(1,:)
- MxA(2,bc%ibulk2) = MxA(2,bc%ibulk2) - bc%B*T(2,:)
- MxA(3,bc%ibulk2) = MxA(3,bc%ibulk2) - bc%B*T(3,:)
-
!-- intermediate storage of outputs --
call store_dataT(bc%dataT,bc%D,bc%V,bc%T,it)
@@ -345,7 +333,6 @@
!===============================================================
- !---------------------------------------------------------------
subroutine init_dataXZ(dataXZ,nglob)
type(dataXZ_type), intent(inout) :: dataXZ
@@ -365,9 +352,6 @@
end subroutine init_dataXZ
-
-
-
!---------------------------------------------------------------
!LOAD_VSLIP_SNAPSHOTS(v,dataXZ,itime,coord,npoin,nglob,iflt)
!Loading slip velocity from snapshots.
@@ -405,7 +389,6 @@
close(IIN_BIN)
end subroutine load_vslip_snapshots
-!---------------------------------------------------------------
!===============================================================
@@ -502,7 +485,6 @@
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))
Modified: seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/prepare_timerun.f90
===================================================================
--- seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/prepare_timerun.f90 2012-10-25 17:27:11 UTC (rev 20928)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/prepare_timerun.f90 2012-10-25 18:26:37 UTC (rev 20929)
@@ -116,7 +116,7 @@
veloc(:,:) = 0._CUSTOM_REAL
! Loading kinematic and dynamic fault solvers.
- call BC_DYNFLT_init(prname,rmass,DT,NSTEP,accel,veloc)
+ call BC_DYNFLT_init(prname,rmass,DT,NSTEP,veloc)
call BC_KINFLT_init(prname,rmass,DT,NSTEP)
More information about the CIG-COMMITS
mailing list