[cig-commits] r21110 - seismo/3D/SPECFEM3D/trunk/src/specfem3D
surendra at geodynamics.org
surendra at geodynamics.org
Sat Dec 8 13:33:14 PST 2012
Author: surendra
Date: 2012-12-08 13:33:13 -0800 (Sat, 08 Dec 2012)
New Revision: 21110
Modified:
seismo/3D/SPECFEM3D/trunk/src/specfem3D/fault_solver_common.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/fault_solver_dynamic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/fault_solver_kinematic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.F90
Log:
Fixed some issues with merging
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/fault_solver_common.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/fault_solver_common.f90 2012-12-07 04:28:13 UTC (rev 21109)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/fault_solver_common.f90 2012-12-08 21:33:13 UTC (rev 21110)
@@ -42,12 +42,12 @@
!---------------------------------------------------------------------
-subroutine initialize_fault (bc,IIN_BIN,Minv,dt_tmp)
+subroutine initialize_fault (bc,IIN_BIN,dt_tmp)
use specfem_par
+ use specfem_par_elastic, only : rmassx,rmassy,rmassz
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
@@ -84,7 +84,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(nxyz(3,bc%nglob))
@@ -124,15 +124,16 @@
call normalize_3d_vector(nxyz)
call compute_R(bc%R,bc%nglob,nxyz)
+ !SURENDRA : WARNING! Assuming rmassx=rmassy=rmassz
! Needed in dA_Free = -K2*d2/M2 + K1*d1/M1
- bc%invM1 = Minv(bc%ibulk1)
- bc%invM2 = Minv(bc%ibulk2)
+ bc%invM1 = rmassx(bc%ibulk1)
+ bc%invM2 = rmassx(bc%ibulk2)
! Fault impedance, Z in : Trac=T_Stick-Z*dV
! Z = 1/( B1/M1 + B2/M2 ) / (0.5*dt)
! T_stick = Z*Vfree traction as if the fault was stuck (no displ discontinuity)
! NOTE: same Bi on both sides, see note above
- 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*bc%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
@@ -297,7 +298,7 @@
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)
+ open(IIN,file='../DATA/FAULT_STATIONS',status='old',action='read',iostat=ier)
if (ier /= 0) then
if (myrank==0) write(IMAIN,*) 'Fatal error opening FAULT_STATIONS file. Abort.'
stop
@@ -316,7 +317,7 @@
allocate(dataT%name(dataT%npoin))
allocate(dist_loc(dataT%npoin)) !Surendra : for parallel fault
- open(IIN,file='DATA/FAULT_STATIONS',status='old',action='read')
+ open(IIN,file='../DATA/FAULT_STATIONS',status='old',action='read')
read(IIN,*) np
k = 0
do i=1,np
@@ -418,6 +419,7 @@
!---------------------------------------------------------------
subroutine store_dataT(dataT,d,v,t,itime)
+ use specfem_par, only : myrank
class(dataT_type), intent(inout) :: dataT
real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: d,v,t
integer, intent(in) :: itime
@@ -454,7 +456,7 @@
write(my_fmt,'(a,i1,a)') '(',dataT%ndat+1,'(E15.7))'
do i=1,dataT%npoin
- open(IOUT,file='OUTPUT_FILES/'//trim(dataT%name(i))//'.dat',status='replace')
+ 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
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/fault_solver_dynamic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/fault_solver_dynamic.f90 2012-12-07 04:28:13 UTC (rev 21109)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/fault_solver_dynamic.f90 2012-12-08 21:33:13 UTC (rev 21110)
@@ -91,12 +91,12 @@
! Minv inverse mass matrix
! dt global time step
!
-subroutine BC_DYNFLT_init(prname,Minv,DTglobal,nt,vel,myrank)
+subroutine BC_DYNFLT_init(prname,DTglobal,vel,myrank)
+ use specfem_par, only : nt=>NSTEP
character(len=256), intent(in) :: prname ! 'proc***'
- real(kind=CUSTOM_REAL), intent(in) :: Minv(:)
double precision, intent(in) :: DTglobal
- integer, intent(in) :: nt,myrank
+ integer, intent(in) :: myrank
real(kind=CUSTOM_REAL), intent(inout) :: vel(:,:)
real(kind=CUSTOM_REAL) :: dt
@@ -112,7 +112,7 @@
dummy_idfault = 0
- open(unit=IIN_PAR,file='DATA/Par_file_faults',status='old',iostat=ier)
+ open(unit=IIN_PAR,file='../DATA/Par_file_faults',status='old',iostat=ier)
if( ier /= 0 ) then
if (myrank==0) write(IMAIN,*) 'File DATA/Par_file_faults not found: assume no faults'
close(IIN_PAR)
@@ -158,7 +158,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)
+ call init_one_fault(faults(iflt),IIN_BIN,IIN_PAR,dt,nt,vel,iflt)
enddo
close(IIN_BIN)
close(IIN_PAR)
@@ -186,12 +186,11 @@
!---------------------------------------------------------------------
-subroutine init_one_fault(bc,IIN_BIN,IIN_PAR,Minv,dt,NT,vel,iflt)
+subroutine init_one_fault(bc,IIN_BIN,IIN_PAR,dt,NT,vel,iflt)
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
@@ -200,7 +199,7 @@
NAMELIST / INIT_STRESS / S1,S2,S3,n1,n2,n3
- call initialize_fault(bc,IIN_BIN,Minv,dt)
+ call initialize_fault(bc,IIN_BIN,dt)
if (bc%nspec>0) then
@@ -236,7 +235,7 @@
allocate(bc%MU(bc%nglob))
if (RATE_AND_STATE) then
allocate(bc%rsf)
- call rsf_init(bc%rsf,bc%T0,bc%V,bc%Fload,bc%coord,IIN_PAR)
+ call rsf_init(bc%rsf,bc%T0,bc%V,vel,bc%Fload,bc%coord,IIN_PAR)
else
allocate(bc%swf)
call swf_init(bc%swf,bc%MU,bc%coord,IIN_PAR)
@@ -694,11 +693,12 @@
!=====================================================================
-subroutine rsf_init(f,T0,V,nucFload,coord,IIN_PAR)
+subroutine rsf_init(f,T0,V,vel,nucFload,coord,IIN_PAR)
type(rsf_type), intent(out) :: f
real(kind=CUSTOM_REAL), intent(in) :: T0(:,:)
real(kind=CUSTOM_REAL), intent(inout) :: V(:,:)
+ real(kind=CUSTOM_REAL), intent(inout) :: vel(:,:)
real(kind=CUSTOM_REAL), intent(in) :: coord(:,:)
real(kind=CUSTOM_REAL), pointer :: nucFload(:)
integer, intent(in) :: IIN_PAR
@@ -783,6 +783,7 @@
call init_2d_distribution(f%Vw,coord,IIN_PAR,nVw)
!!$ ! WARNING : Not general enough
+!!$ vel = 0._CUSTOM_REAL
!!$ nglob_bulk = size(vel,2)
!!$ allocate(init_vel(3,nglob_bulk))
!!$ init_vel = 0._CUSTOM_REAL
@@ -1005,7 +1006,7 @@
call itime(now) ! now(1)=hour, (2)=minute, (3)=second
- write(filename,"('OUTPUT_FILES/RuptureTime_Fault',I0)") iflt
+ write(filename,"('../OUTPUT_FILES/RuptureTime_Fault',I0)") iflt
IOUT = 121 !WARNING: not very robust. Could instead look for an available ID
@@ -1172,7 +1173,7 @@
character(len=70) :: filename
- write(filename,"('OUTPUT_FILES/Snapshot',I0,'_F',I0,'.bin')") itime,iflt
+ write(filename,"('../OUTPUT_FILES/Snapshot',I0,'_F',I0,'.bin')") itime,iflt
open(unit=IOUT, file= trim(filename), status='replace', form='unformatted',action='write')
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/fault_solver_kinematic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/fault_solver_kinematic.f90 2012-12-07 04:28:13 UTC (rev 21109)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/fault_solver_kinematic.f90 2012-12-08 21:33:13 UTC (rev 21110)
@@ -53,12 +53,12 @@
! Minv inverse mass matrix
! dt global time step
!
-subroutine BC_KINFLT_init(prname,Minv,DTglobal,nt,myrank)
+subroutine BC_KINFLT_init(prname,DTglobal,myrank)
+ use specfem_par, only : nt=>NSTEP
character(len=256), intent(in) :: prname ! 'proc***'
- real(kind=CUSTOM_REAL), intent(in) :: Minv(:)
double precision, intent(in) :: DTglobal
- integer, intent(in) :: nt,myrank
+ integer, intent(in) :: myrank
real(kind=CUSTOM_REAL) :: dt
integer :: iflt,ier,dummy_idfault
@@ -73,7 +73,7 @@
dummy_idfault = 0
- open(unit=IIN_PAR,file='DATA/Par_file_faults',status='old',iostat=ier)
+ open(unit=IIN_PAR,file='../DATA/Par_file_faults',status='old',iostat=ier)
if( ier /= 0 ) then
if (myrank==0) write(IMAIN,*) 'File DATA/Par_file_faults not found: assume no faults'
close(IIN_PAR)
@@ -111,7 +111,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,iflt)
+ call init_one_fault(faults(iflt),IIN_BIN,IIN_PAR,dt,nt,iflt)
enddo
endif
close(IIN_BIN)
@@ -128,10 +128,9 @@
!---------------------------------------------------------------------
-subroutine init_one_fault(bc,IIN_BIN,IIN_PAR,Minv,dt,NT,iflt)
+subroutine init_one_fault(bc,IIN_BIN,IIN_PAR,dt,NT,iflt)
type(bc_kinflt_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
@@ -139,7 +138,7 @@
NAMELIST / KINPAR / kindt
- call initialize_fault(bc,IIN_BIN,Minv,dt)
+ call initialize_fault(bc,IIN_BIN,dt)
if (bc%nspec>0) then
@@ -346,7 +345,7 @@
IIN_BIN=101
IOUT = 102
- write(filename,"('OUTPUT_FILES/Snapshot',I0,'_F',I0,'.bin')") itime,iflt
+ write(filename,"('../OUTPUT_FILES/Snapshot',I0,'_F',I0,'.bin')") itime,iflt
print*, trim(filename)
open(unit=IIN_BIN, file= trim(filename), status='old', form='formatted',&
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.F90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.F90 2012-12-07 04:28:13 UTC (rev 21109)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.F90 2012-12-08 21:33:13 UTC (rev 21110)
@@ -33,6 +33,8 @@
use specfem_par_elastic
use specfem_par_poroelastic
use specfem_par_movie
+ use fault_solver_dynamic, only : BC_DYNFLT_init, RATE_AND_STATE
+ use fault_solver_kinematic, only : BC_KINFLT_init
implicit none
! local parameters
double precision :: tCPU
@@ -49,6 +51,12 @@
! initializes arrays
call prepare_timerun_init_wavefield()
+ call sync_all()
+ ! Loading kinematic and dynamic fault solvers.
+ call BC_DYNFLT_init(prname,rmass,DT,NSTEP,veloc,myrank)
+
+ call BC_KINFLT_init(prname,rmass,DT,NSTEP,myrank)
+
! sets up time increments
call prepare_timerun_constants()
@@ -326,6 +334,7 @@
use specfem_par_acoustic
use specfem_par_elastic
use specfem_par_poroelastic
+ use fault_solver_dynamic, only : RATE_AND_STATE
implicit none
! initialize acoustic arrays to zero
@@ -340,7 +349,7 @@
! initialize elastic arrays to zero/verysmallvall
if( ELASTIC_SIMULATION ) then
displ(:,:) = 0._CUSTOM_REAL
- veloc(:,:) = 0._CUSTOM_REAL
+ if(.not. RATE_AND_STATE) veloc(:,:) = 0._CUSTOM_REAL
accel(:,:) = 0._CUSTOM_REAL
! put negligible initial value to avoid very slow underflow trapping
if(FIX_UNDERFLOW_PROBLEM) displ(:,:) = VERYSMALLVAL
More information about the CIG-COMMITS
mailing list