[cig-commits] r20809 - seismo/3D/FAULT_SOURCE/branches/new_fault_db/src

[surendra at geodynamics.org]

Author: surendra
Date: 2012-10-06 11:53:15 -0700 (Sat, 06 Oct 2012)
New Revision: 20809

Modified:
   seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90
   seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/prepare_timerun.f90
Log:
Added bulk velocities for Rate and State (vey specific to tpv102, not robust enough)

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-06 15:27:33 UTC (rev 20808)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/fault_solver.f90	2012-10-06 18:53:15 UTC (rev 20809)
@@ -127,7 +127,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
+       SIMULATION_TYPE_DYN,PARALLEL_FAULT,RATE_AND_STATE
 
 
 contains
@@ -139,13 +139,14 @@
 ! Minv          inverse mass matrix
 ! dt            global time step
 !
-subroutine BC_DYNFLT_init(prname,Minv,DTglobal,nt,accel)
+subroutine BC_DYNFLT_init(prname,Minv,DTglobal,nt,accel,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
   integer :: iflt,ier,dummy_idfault
@@ -198,7 +199,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,accel)
+    call init_one_fault(faults(iflt),IIN_BIN,IIN_PAR,Minv,dt,nt,vel,iflt,accel)
   enddo
   close(IIN_BIN)
   close(IIN_PAR)
@@ -223,11 +224,13 @@
 
 !---------------------------------------------------------------------
 
-subroutine init_one_fault(bc,IIN_BIN,IIN_PAR,Minv,dt_tmp,NT,iflt,accel)
+subroutine init_one_fault(bc,IIN_BIN,IIN_PAR,Minv,dt_tmp,NT,vel,iflt,accel)
 
   use specfem_par
 
   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
@@ -273,9 +276,10 @@
   real(kind=CUSTOM_REAL) :: x,z
   logical :: c1,c2,c3,c4
   real(kind=CUSTOM_REAL) :: b11,b12,b21,b22,B1,B2
-  integer iglob
+  integer iglob,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
@@ -554,6 +558,16 @@
       call init_2d_distribution(bc%rsf%C  ,bc%coord,IIN_PAR,nC)
       call init_2d_distribution(bc%rsf%T  ,bc%coord,IIN_PAR,nForcedRup)
 
+    ! 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
     
       W1=15000._CUSTOM_REAL
       W2=7500._CUSTOM_REAL

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-06 15:27:33 UTC (rev 20808)
+++ seismo/3D/FAULT_SOURCE/branches/new_fault_db/src/prepare_timerun.f90	2012-10-06 18:53:15 UTC (rev 20809)
@@ -114,8 +114,9 @@
   ! sets up mass matrices
   call prepare_timerun_mass_matrices()
 
+veloc(:,:) = 0._CUSTOM_REAL
   ! Loading kinematic and dynamic fault solvers.
-  call BC_DYNFLT_init(prname,rmass,DT,NSTEP,accel)
+  call BC_DYNFLT_init(prname,rmass,DT,NSTEP,accel,veloc)
 
   call BC_KINFLT_init(prname,rmass,DT,NSTEP)
 
@@ -131,7 +132,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