bc_swff.f90

Spectral Element Method for wave propagation and rupture dynamics.

! SEM2DPACK version 2.2.11 -- A Spectral Element Method for 2D wave propagation and fracture dynamics,
!                             with emphasis on computational seismology and earthquake source dynamics.
! 
! Copyright (C) 2003-2007 Jean-Paul Ampuero
! All Rights Reserved
! 
! Jean-Paul Ampuero
! 
! ETH Zurich (Swiss Federal Institute of Technology)
! Institute of Geophysics
! Seismology and Geodynamics Group
! ETH H鰊ggerberg HPP O 13.1
! CH-8093 Z黵ich
! Switzerland
! 
! [e-mail]
! [phone] (office)
! [phone] (fax)
! 
! http://www.sg.geophys.ethz.ch/geodynamics/ampuero/
! 
! 
! This software is freely available for scientific research purposes. 
! If you use this software in writing scientific papers include proper 
! attributions to its author, Jean-Paul Ampuero.
! 
! 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., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
! 
module bc_swff
! slip weakening friction fault

  use bnd_grid, only: bnd_grid_type
  use distribution_general, only: cd_type

  implicit none
  private

  type bc_swff_input_type
    type(cd_type) :: Dc,MuS,MuD,T,N,Sxx,Sxy,Sxz,Syz,Szz
  end type bc_swff_input_type

  type bc_swff_type
    private
    logical :: UpdateFrictionFirst
    double precision :: CoefA2V,CoefA2D,OutputT1,OutputDt
    integer :: OutputNextTime,OutputNdt,OutputUnit,OutputIFir &
              ,OutputILas,OutputILag,OutputNpoin
    double precision, pointer :: n1(:,:)
    double precision, dimension(:), pointer:: &
      B,invM1=>null(),invM2=>null() ,Z,Tt0,Tn0 &
     ,MU,Tt,Tn,Dc,MuStatic,MuDynamic
    type(bnd_grid_type), pointer :: bc1 => null(), bc2 => null()
    type(bc_swff_input_type) :: input
  end type bc_swff_type

  public :: BC_SWFF_type, BC_SWFF_read, BC_SWFF_init, BC_SWFF_set, BC_SWFF_write

contains

!=====================================================================
! BEGIN INPUT BLOCK
!
! NAME   : BC_SWFFLT
! GROUP  : BOUNDARY_CONDITION
! PURPOSE: Slip weakening friction fault
! SYNTAX : &BC_SWFFLT Dc | DcHet, MuS | MuSHet , MuD | MuDHet, 
!                     Tn | TnHet, Tt | TtHet,
!                     Sxx | SxxHet, Sxy | SxyHet, Sxz | SxzHet, 
!		      Syz | SyzHet, Szz | SzzHet
!                     FirstOutput, DtOutput, IxOut /
!             followed eventually by distribution input blocks &DIST_XXX
!             for Dc,MuS,MuD,Tn and/or Tt (the order is important)
!
! NOTE: for better results, use dynamic faults with the leapfrog time scheme
!       and with a layer of damping material (Kelvin-Voigt) near the fault.
!
! Friction law:
! ARG: Dc       [dble] [0.5d0] Critical slip 
! ARG: MuS      [dble] [0.6d0] Static friction coefficient
! ARG: MuD      [dble] [0.5d0] Dynamic friction coefficient
!
! Initial stress, can be a superposition of tractions and background stress:
! ARG: Tn       [dble] [-100d6] Normal traction (positive = tensile)
! ARG: Tt       [dble] [55d6] Tangential traction (positive antiplane: y>0)
! ARG: Sxx      [dble] [0d0] sigma_xx
! ARG: Sxy      [dble] [0d0] sigma_xy
! ARG: Sxz      [dble] [0d0] sigma_xz
! ARG: Syz      [dble] [0d0] sigma_yz
! ARG: Szz      [dble] [0d0] sigma_zz
!
! NOTE: arguments with the suffix "Het" are used to give
!       friction and initial stress parameters non uniform values.
!	For instance, DcHet='GAUSSIAN' followed by a DIST_GAUSSIAN block
!	sets a gaussian distribution of Dc.
!	Several heterogeneous distributions are available, 
!       See DIST_XXX for their syntax.
!
! For outputs in FltXX_sem2d.dat:
! ARG: DtOutput [dble] [0.d0] Time lag between outputs (in seconds)
!               Default resets DtOutput = global timestep
! ARG: FirstOutput [dble] [0.d0] Start output at this time
! ARG: IxOut    [int(3)] [(1,huge,1)] First node, last node and stride
!               Default resets Ixout(2)= last point
!
! NOTE: DtOutput is internally adjusted to the nearest multiple 
!       of the global timestep
!
! END INPUT BLOCK

  subroutine BC_SWFF_read(bc,iin)

  use echo, only : echo_input,iout
  use stdio, only: IO_abort
  use distribution_general, only: DIST_Read_CD

  type(bc_swff_type), intent(out) :: bc
  integer, intent(in) :: iin
  double precision :: Dc,MuS,MuD,Tt,Tn,FirstOutput,DtOutput &
                     ,Sxx,Sxy,Sxz,Syz,Szz
  character(20) :: DcHet,MuSHet,MuDHet,TtHet,TnHet &
                  ,SxxHet,SxyHet,SxzHet,SyzHet,SzzHet
  integer :: i,IxOut(3)

  NAMELIST / BC_SWFFLT /  Dc,MuS,MuD,Tt,Tn &
                         ,DcHet,MuSHet,MuDHet,TtHet,TnHet &
                         ,Sxx,Sxy,Sxz,Syz,Szz &
                         ,SxxHet,SxyHet,SxzHet,SyzHet,SzzHet &
                         ,FirstOutput,DtOutput,IxOut

  Dc = 0.5d0
  MuS = 0.6d0
  MuD = 0.5d0
  Tt = 55d6
  Tn = -100.d6
  Sxx = 0d0
  Sxy = 0d0
  Sxz = 0d0
  Syz = 0d0
  Szz = 0d0

  DcHet = ''
  MuSHet = ''
  MuDHet = ''
  TtHet = ''
  TnHet = ''
  SxxHet = ''
  SxyHet = ''
  SxzHet = ''
  SyzHet = ''
  SzzHet = ''

  FirstOutput = 0.d0
  DtOutput = 0.d0
  IxOut(1) = 1
  IxOut(2) = huge(i)
  IxOut(3) = 1

  read(iin,BC_SWFFLT,END=100)
  
  bc%OutputT1 = FirstOutput
  bc%OutputDt = DtOutput
  bc%OutputIFir = IxOut(1)
  bc%OutputILas = IxOut(2) 
  bc%OutputILag = IxOut(3) 

  call DIST_Read_CD(Dc,DcHet,bc%input%Dc,iin)
  call DIST_Read_CD(MuS,MuSHet,bc%input%MuS,iin)
  call DIST_Read_CD(MuD,MuDHet,bc%input%MuD,iin)
  call DIST_Read_CD(Tn,TnHet,bc%input%N,iin)
  call DIST_Read_CD(Tt,TtHet,bc%input%T,iin)
  call DIST_Read_CD(Sxx,SxxHet,bc%input%Sxx,iin)
  call DIST_Read_CD(Sxy,SxyHet,bc%input%Sxy,iin)
  call DIST_Read_CD(Sxz,SxzHet,bc%input%Sxz,iin)
  call DIST_Read_CD(Syz,SyzHet,bc%input%Syz,iin)
  call DIST_Read_CD(Szz,SzzHet,bc%input%Szz,iin)

  if (echo_input) then
    write(iout,200) DcHet,MuSHet,MuDHet,TnHet,TtHet &
                   ,SxxHet,SxyHet,SxzHet,SyzHet,SzzHet &
                   ,FirstOutput,DtOutput,IxOut
    write(iout,*) 'NOTE: FirstOutput and DtOutput may be modified later'
    if (DtOutput==0d0) &
     write(iout,*) 'NOTE: DtOutput will be reset to the global timestep'
    write(iout,*) '      Final values can be found in FltXX_sem2d.hdr'
    if (IxOut(2)==huge(i)) &
     write(iout,*) 'NOTE: Last output node, IxOut(2), will be reset later'
  endif

  return
  100 call IO_abort('BC_SWFF_read: BC_SWFFLT input block not found')
  200 format(5x,'Type   = Slip weakening fault', &
            /5x,'  Critical slip         = ',A,&
            /5x,'  Static friction       = ',A,&
            /5x,'  Dynamic friction      = ',A,&
            /5x,'  Initial traction T_n  = ',A,&
            /5x,'  Initial traction T_t  = ',A,&
            /5x,'  Initial stress S_xx   = ',A,&
            /5x,'  Initial stress S_xy   = ',A,&
            /5x,'  Initial stress S_xz   = ',A,&
            /5x,'  Initial stress S_yz   = ',A,&
            /5x,'  Initial stress S_zz   = ',A,&
            /5x,'  First output time     = ',EN12.3,&
            /5x,'  Output time stride    = ',EN12.3,&
            /5x,'  First output node     = ',I0,&
            /5x,'  Last output node      = ',I0,&
            /5x,'  Output node stride    = ',I0)

  end subroutine BC_SWFF_read


!=====================================================================
!
  subroutine BC_SWFF_init(bc,tags,grid,M,time,v,perio)
  
  use spec_grid, only : sem_grid_type,BC_inquire,BC_get_normal_and_weights
  use stdio, only: IO_abort,IO_new_unit
  use time_evol, only: timescheme_type
  use distribution_general, only: DIST_Init_CD
  use constants, only : TINY_XABS
  use bc_periodic, only : bc_periodic_type,BC_PERIO_intersects

  type(bc_swff_type)  , intent(inout) :: bc
  type(sem_grid_type), intent(in) :: grid
  double precision, intent(in) :: M(:)
  integer, intent(in) :: tags(2)
  type(timescheme_type), intent(in) :: time
  double precision, intent(inout) :: v(:,:)
  type(bc_periodic_type), pointer :: perio

  double precision, allocatable :: FltCoord(:,:)
  double precision, dimension(:), pointer :: Sxx,Sxy,Sxz,Syz,Szz,Tx,Ty,Tz,nx,nz
  integer :: i,HdrUnit,npoin,NSAMP,NDAT
  character(15) :: OutputFileName,HdrFileName
  logical :: two_sides
  
! if fault is on a domain boundary and symmetry is required, 
! bc2 and invM2 are not defined
  two_sides = tags(2)>0

! bc1 --> grid%bounds(i) boundary corresponding to TAG
  call BC_inquire( grid%bounds, tag = tags(1), bc_topo_ptr = bc%bc1 )
  if ( two_sides ) then
    call BC_inquire( grid%bounds, tag = tags(2), bc_topo_ptr = bc%bc2 )

    if (bc%bc1%nelem/=bc%bc2%nelem) &
     call IO_abort('bc_swff_init: number of boundary elements do not match')

    if (bc%bc1%npoin/=bc%bc2%npoin) &
     call IO_abort('bc_swff_init: number of nodes on boundaries do not match')
  endif
  npoin = bc%bc1%npoin

  allocate(FltCoord(2,npoin))
  FltCoord = grid%coord(:,bc%bc1%node)
  if ( two_sides ) then
    if ( any(abs(FltCoord(1,:)-grid%coord(1,bc%bc2%node))>TINY_XABS) &
      .OR.any(abs(FltCoord(2,:)-grid%coord(2,bc%bc2%node))>TINY_XABS) )&
      call IO_abort('bc_swff_init: coordinates on boundaries do not match properly')
  endif

! NOTE: the mesh being conformal, the weights B=GLL_weights*jac1D are equal on both
!       sides of the fault. 
  allocate( bc%n1(npoin,2) )
  allocate( bc%B(npoin) ) ! assembled[ GLL_weights * jac1D ]
  call BC_get_normal_and_weights(bc%bc1,grid,bc%n1,bc%B, BC_PERIO_intersects(bc%bc1,perio) )

! Update coefficients of Newmark scheme
! vnew = vpredictor +gamma*dt*anew   See solver.f90 
! dnew = dpredictor +beta*dt^2*anew   See solver.f90
  bc%CoefA2V = time%gamma * time%dt
  bc%CoefA2D = time%beta * time%dt**2
  bc%UpdateFrictionFirst = time%kind == 'leapfrog'
!  bc%UpdateFrictionFirst = .true. !WARNING: testing

! Needed in dA_Free = -K2*d2/M2 + K1*d1/M1
  allocate(bc%invM1(npoin))
  bc%invM1 = 1d0/M(bc%bc1%node)
  if ( two_sides ) then
    allocate(bc%invM2(npoin))
    bc%invM2 = 1d0/M(bc%bc2%node)
  endif

! Fault impedance, Z in :  Trac=T_Stick-Z*dV
!   Z = 1/( B1/M1 + B2/M2 )
! T_Stick = Z*Vfree traction as if the fault was stuck (no displ discontinuity) 
! NOTE: same Bi on both sides, see note above
  allocate(bc%Z(npoin))
  if ( two_sides ) then
    bc%Z = 1.d0/(bc%CoefA2V * bc%B *( bc%invM1 + bc%invM2 ))
  else
    bc%Z = 0.5d0/(bc%CoefA2V * bc%B * bc%invM1 )
  endif

! Friction parameters
  call DIST_Init_CD(bc%input%Dc,FltCoord,bc%Dc)
  call DIST_Init_CD(bc%input%MuS,FltCoord,bc%MuStatic)
  call DIST_Init_CD(bc%input%MuD,FltCoord,bc%MuDynamic)

! Initial stress
  call DIST_Init_CD(bc%input%T,FltCoord,bc%Tt0)
  call DIST_Init_CD(bc%input%N,FltCoord,bc%Tn0)
  call DIST_Init_CD(bc%input%Sxx,FltCoord,Sxx)
  call DIST_Init_CD(bc%input%Sxy,FltCoord,Sxy)
  call DIST_Init_CD(bc%input%Sxz,FltCoord,Sxz)
  call DIST_Init_CD(bc%input%Syz,FltCoord,Syz)
  call DIST_Init_CD(bc%input%Szz,FltCoord,Szz)
  nx => bc%n1(:,1)
  nz => bc%n1(:,2)
  allocate(Tx(npoin))
  allocate(Ty(npoin))
  allocate(Tz(npoin))
  Tx = Sxx*nx + Sxz*nz
  Ty = Sxy*nx + Syz*nz
  Tz = Sxz*nx + Szz*nz
  bc%Tn0 = bc%Tn0 + Tx*nx + Tz*nz
  if (size(v,2)==1) then ! SH
    bc%Tt0 = bc%Tt0 + Ty
  else ! P-SV
    bc%Tt0 = bc%Tt0 + Tx*nz - Tz*nx
  endif 
  deallocate(Sxx,Sxy,Sxz,Syz,Szz,Tx,Ty,Tz)

! Initial friction
  allocate(bc%MU(npoin))
  bc%MU = bc%MuStatic
 
! BETA :
! Initial velocities
! To guarantee second order when initial stress drop is abrupt, 
! the initial velocity must be set analytically:
!   impedance*slip_rate = max( initial_stress - initial_strength, 0)
! where impedance = 0.5*rho*cs
!       slip rate = v_2 - v_1
!                 = -2*v_1 if symmetry
!  if (size(v,2)==1) then ! SH