bc_swff.f90

Spectral Element Method for wave propagation and rupture dynamics.

!    v(bc%bc1%node,1) = -sign( max( abs(bc%Tt0)+bc%MU*bc%Tn0, 0d0), bc%Tt0 ) &
!                         / (2670d0*3464d0)
!    if (two_sides) v(bc%bc2%node,1) = -v(bc%bc1%node,1)
!  else ! P-SV
!! WARNING: not implemented yet
!  endif


! Output arrays
  bc%OutputIFir = max(bc%OutputIFir, 1)
  bc%OutputILas = min(bc%OutputILas, npoin)
  bc%OutputNpoin = (bc%OutputILas-bc%OutputIFir)/bc%OutputILag +1
  allocate(bc%Tt(bc%OutputNpoin))
  allocate(bc%Tn(bc%OutputNpoin))
  bc%Tt=bc%Tt0(bc%OutputIFir:bc%OutputILas:bc%OutputILag)
  bc%Tn=bc%Tn0(bc%OutputIFir:bc%OutputILas:bc%OutputILag)
  if (size(bc%Tn0(bc%OutputIFir:bc%OutputILas:bc%OutputILag))/=bc%OutputNpoin) &
    stop 'couille avec npoin'

!-- Open output files -----
! Set output parameters

! NOTE: file names limited to tags(1)<100
  write(OutputFileName,'("Flt",I2.2,"_sem2d.dat")') tags(1)
  bc%OutputUnit = IO_new_unit()
  open(bc%OutputUnit,file=OutputFileName,status='replace',form='unformatted')

 ! adjust output timestep to the nearest multiple of dt:
  if (bc%OutputDt==0d0) then
    bc%OutputNdt=1
  else
    bc%OutputNdt = nint(bc%OutputDt/time%dt)
  endif
  bc%OutputDt = time%dt * bc%OutputNdt
  bc%OutputNextTime = nint(bc%OutputT1/time%dt)

 ! HEADER FILE FORMAT:
 !      Name            FltXX_sem2d.hdr where XX=tags(1) of the BC_SWFFLT
 !                      input block, the tag of the first side of the fault.
 !      Line 1          Name of size/other variables
 !      Line 2          Value of size/other variables
 !      Line 3          Name of data fields, separated by ":"
 !      Line 4          Name of coordinate axis
 !      Line 5:end      Two column table of nodes coordinates
 !
 ! OUTPUT FILE FORMAT:
 !      At each time (lag DELT, total NSAMP), NDAT data lines with
 !      NPTS columns, one per fault node, are written.
 !      Stress fields are relative to initial values and
 !      the first NDAT lines contain the initial conditions.
 !
 ! NOTE: when the binary file is written one line at a time
 !      a record tag is also written at the beginning and end of line,
 !      => number of columns in output file = nb of data columns +2
 !
  write(HdrFileName,'("Flt",I2.2,"_sem2d")') tags(1)
  NDAT = 5
  NSAMP = (time%nt -bc%OutputNextTime)/bc%OutputNdt +1
  HdrUnit = IO_new_unit()

  open(HdrUnit,file=trim(HdrFileName)//'.hdr',status='replace')
  write(HdrUnit,*) 'NPTS NDAT NSAMP DELT'
  write(HdrUnit,*) bc%OutputNpoin,NDAT,NSAMP,bc%OutputDt
                  
  write(HdrUnit,*) 'Slip:Slip Rate:Shear Stress:Normal Stress:Friction'
  write(HdrUnit,*) 'XPTS ZPTS'
  do i=bc%OutputIFir,bc%OutputILas,bc%OutputILag
    write(HdrUnit,*) FltCoord(:,i)
  enddo
  close(HdrUnit)
  deallocate(FltCoord)

! Sample SU script subset|b2a
 ! NOTE: how to convert to a human-readable ASCII table
 !      the Slip at time 10*DELT (assuming first output at t=0), 
 !      using Seismic Unix's subset and b2a :
 !              subset <Flt05_sem2d.hdr n1=(NPTS+2) n2=NDAT n3=NSAMP \
 !              if1s=1 n1s=NPTS ix2s=0 ix3s=10 \
 !              | b2a n1=NPTS > slip_t10.tab
 !
  open(HdrUnit,file=trim(HdrFileName)//'.csh',status='replace')
  write(HdrUnit,*) '#!/bin/csh -f'
  write(HdrUnit,*) '#USAGE: '//trim(HdrFileName)//'.csh node field'
  write(HdrUnit,*) '@ ifs = $2 - 1'
  write(HdrUnit,'(A,I0,A,I0,A,I0,A)') 'subset <'//trim(HdrFileName)//'.dat n1=',bc%OutputNpoin+2, &
    ' n2=',NDAT,' n3=',NSAMP,' ix1s=$1 ix2s=$ifs | b2a n1=1 > '//trim(HdrFileName)//'.$1.$ifs.tab'
  close(HdrUnit)

  end subroutine BC_SWFF_init


!=====================================================================
! Computes the boundary term B*T with T satisfying a friction law
!
! Convention:   T = sigma*n1
!               D = d2-d1
! => T and D have same sign
!
! NOTE: this is an EXPLICIT scheme FOR FRICTION, 
!       we use the displacement at the PREVIOUS timestep 
!       to compute the friction coefficient
!
! NOTE: possible periodicity does not need to be enforced at this level
!       because it is assumed that MxA and D are already periodic
!
  subroutine BC_SWFF_set(bc,MxA,V,D)

  type(bc_swff_type) :: bc
  double precision, intent(in) :: V(:,:),D(:,:)
  double precision, intent(inout) :: MxA(:,:)

  if (size(V,2) == 1) then
    call BC_SWFF_set_SH(bc,MxA(:,1),V(:,1),D(:,1))
  else
    call BC_SWFF_set_PSV(bc,MxA,V,D)
  endif

  end subroutine BC_SWFF_set

!---------------------------------------------------------------------

  subroutine BC_SWFF_set_SH(bc,MxA,V,D)

  type(bc_swff_type) :: bc
  double precision, intent(in) :: V(:),D(:)
  double precision, intent(inout) :: MxA(:)

  double precision, dimension(bc%bc1%npoin) :: T,dD
  integer, pointer :: i1(:),i2(:)
  logical :: two_sides

  two_sides = associated(bc%bc2)
  i1 => bc%bc1%node

! A test traction is computed as if the fault was sticking at n+1 (no slip rate)
!           = Z * SlipRateFree 
!	    = Z * (SlipRate_predicted + coef*SlipAccelFree)
! where: 
!   Z = fault impedance matrix (diagonal)
!   SlipRateFree, SlipAccelFree = slip velocity and acceleration 
!     				  as if the fault was stress free
!   coef = depends on the time integration scheme
!
  if (two_sides) then
    i2 => bc%bc2%node
    T = bc%Z* ( V(i2)-V(i1) + bc%CoefA2V*( bc%invM2*MxA(i2)-bc%invM1*MxA(i1) ))
    dD = D(i2) - D(i1)
  else
    T  = -2d0*bc%Z* ( V(i1) + bc%CoefA2V * bc%invM1*MxA(i1) )
    dD = -2d0*D(i1)
  endif

! Update slip weakening (only for leapfrog scheme)
  if (bc%UpdateFrictionFirst) bc%MU = friction(dD,bc%MuStatic,bc%MuDynamic,bc%Dc)

! Coulomb friction, applied to the absolute stress
  T = T + bc%Tt0
  T = sign( min(abs(T), -bc%Tn0*bc%MU), T)
  T = T - bc%Tt0

! Add boundary term B*T to M*a
  MxA(i1) = MxA(i1) + bc%B*T
  if (two_sides) MxA(i2) = MxA(i2) - bc%B*T

! Store relative tractions for output
  bc%Tt = T(bc%OutputIFir:bc%OutputILas:bc%OutputILag)
  bc%Tn = 0d0

! Newmark scheme: update friction for next time step
! Introduces a delay of dt, but keeps the friction solver fully explicit
  if (.not. bc%UpdateFrictionFirst) then
    if (two_sides) then
      dD = dD + bc%CoefA2D *( bc%invM2*MxA(i2)-bc%invM1*MxA(i1) )
    else
      dD = dD -2d0*bc%CoefA2D * bc%invM1*MxA(i1)
    endif
    bc%MU = friction(dD, bc%MuStatic, bc%MuDynamic, bc%Dc)
  endif

  end subroutine BC_SWFF_set_SH


!---------------------------------------------------------------------
! WARNING: the opening conditions are incomplete
!          there can be a problem if the fault opens and closes again

  subroutine BC_SWFF_set_PSV(bc,MxA,V,D)

  type(bc_swff_type) :: bc
  double precision, intent(in) :: V(:,:),D(:,:)
  double precision, intent(inout) :: MxA(:,:)

  double precision, dimension(bc%bc1%npoin) :: Slip,Tt,Tn
  double precision, dimension(bc%bc1%npoin,2) :: BxT,T_Stick,dD,dA

  integer, pointer :: i1(:),i2(:)
  logical :: two_sides

  two_sides = associated(bc%bc2)

! Test traction
  i1 => bc%bc1%node
  if (two_sides) then
    i2 => bc%bc2%node
    dD = D(i2,:) - D(i1,:)
    dA(:,1) = bc%invM2*MxA(i2,1) - bc%invM1*MxA(i1,1)
    dA(:,2) = bc%invM2*MxA(i2,2) - bc%invM1*MxA(i1,2)
    T_Stick = V(i2,:)-V(i1,:) + bc%CoefA2V * dA
  else
    dD = - 2d0*D(i1,:) ! used only to compute tangential discontinuity (slip)
    dA(:,1) = - 2d0*bc%invM1*MxA(i1,1)
    dA(:,2) = - 2d0*bc%invM1*MxA(i1,2)
    T_Stick = -2d0*V(i1,:) + bc%CoefA2V * dA
  endif
  T_Stick(:,1) = bc%Z * T_Stick(:,1)
  T_Stick(:,2) = bc%Z * T_Stick(:,2)

 ! Update slip weakening:
!WARNING: there should be no update if there is opening, 
!         if (Norm>0) return
!         should be implemented here if required (unlikely)
!
  Slip = bc%n1(:,2)*dD(:,1) - bc%n1(:,1)*dD(:,2)
  bc%MU = friction(Slip, bc%MuStatic, bc%MuDynamic, bc%Dc)

! Rotate tractions to fault frame (Tt,Tn)
  Tt = bc%n1(:,2)*T_Stick(:,1) - bc%n1(:,1)*T_Stick(:,2)
  if (two_sides) then
    Tn = bc%n1(:,1)*T_Stick(:,1) + bc%n1(:,2)*T_Stick(:,2)
  else
    Tn = 0d0    ! if symmetry: costant normal stress
  endif

! Add initial stress
  Tt = Tt +bc%Tt0
  Tn = Tn +bc%Tn0

 !Opening-->free stress
  Tn = min(Tn,0.d0) ! negative normal stress is compressive
 !Coulomb friction
  Tt = sign( min(abs(Tt),-Tn*bc%MU), Tt)

! Subtract initial stress
  Tt = Tt -bc%Tt0
  Tn = Tn -bc%Tn0

! Rotate tractions back to (x,z) frame and apply boundary weights
  BxT(:,1) = bc%B*(  bc%n1(:,2)*Tt + bc%n1(:,1)*Tn )
  BxT(:,2) = bc%B*( -bc%n1(:,1)*Tt + bc%n1(:,2)*Tn )

! Add boundary term B*T to M*a
  MxA(i1,:) = MxA(i1,:) + BxT
  if (two_sides) MxA(i2,:) = MxA(i2,:) - BxT

! Store tractions for output
  bc%Tt = Tt(bc%OutputIFir:bc%OutputILas:bc%OutputILag)
  bc%Tn = Tn(bc%OutputIFir:bc%OutputILas:bc%OutputILag)

  end subroutine BC_SWFF_set_PSV



!=====================================================================
! OUTPUT FORMAT: at each time, 5 data lines, one column per fault node
! 1: Slip
! 2: Slip rate
! 3: Shear stress 
! 4: Normal stress 
! 5: Friction coefficient
!
  subroutine BC_SWFF_write(bc,displ,veloc,itime)

  type(bc_swff_type) :: bc
  double precision, intent(in) :: Veloc(:,:),Displ(:,:)
  integer, intent(in) :: itime

  double precision, pointer :: n1(:,:)
  integer, pointer :: i1(:),i2(:)

  if (itime<bc%OutputNextTime) return
  bc%OutputNextTime = bc%OutputNextTime + bc%OutputNdt

  if (associated(bc%bc2)) then

    i1 => bc%bc1%node(bc%OutputIFir:bc%OutputILas:bc%OutputILag)
    i2 => bc%bc2%node(bc%OutputIFir:bc%OutputILas:bc%OutputILag)

    if ( size(displ,2)==1 ) then
      write(bc%OutputUnit) real(displ(i2,1)-displ(i1,1))
      write(bc%OutputUnit) real(veloc(i2,1)-veloc(i1,1))
    else
      n1 => bc%n1(bc%OutputIFir:bc%OutputILas:bc%OutputILag,:)
      write(bc%OutputUnit) real( n1(:,2)*(displ(i2,1)-displ(i1,1)) &
                               - n1(:,1)*(displ(i2,2)-displ(i1,2)) )
      write(bc%OutputUnit) real( n1(:,2)*(veloc(i2,1)-veloc(i1,1)) &
                               - n1(:,1)*(veloc(i2,2)-veloc(i1,2)) )
    endif

  else

    i1 => bc%bc1%node(bc%OutputIFir:bc%OutputILas:bc%OutputILag)
    if ( size(displ,2)==1 ) then
      write(bc%OutputUnit) -real(2d0*displ(i1,1))
      write(bc%OutputUnit) -real(2d0*veloc(i1,1))
    else
      n1 => bc%n1(bc%OutputIFir:bc%OutputILas:bc%OutputILag,:)
      write(bc%OutputUnit) -2.0*real( n1(:,2)*displ(i1,1) &
                                - n1(:,1)*displ(i1,2) )
      write(bc%OutputUnit) -2.0*real( n1(:,2)*veloc(i1,1) &
                                - n1(:,1)*veloc(i1,2) )
    endif
  endif

  write(bc%OutputUnit) real( bc%Tt )
  write(bc%OutputUnit) real( bc%Tn )
  write(bc%OutputUnit) real( bc%MU(bc%OutputIFir:bc%OutputILas:bc%OutputILag) )

  end subroutine BC_SWFF_write



!=====================================================================
 !-- linear slip weakening:
  function friction(s,mus,mud,dc) result(mu)

  double precision, dimension(:), intent(in) :: s,mus,mud,dc
  double precision :: mu(size(s))

  mu = mus -(mus-mud)/dc *abs(s)
  mu = max( mu, mud)

  end function friction

!---------------------------------------------------------------------
 !-- exponential slip weakening:
  function friction_exponential(s,mus,mud,dc) result(mu)

  double precision, dimension(:), intent(in) :: s,mus,mud,dc
  double precision :: mu(size(s))

  mu = mus -(mus-mud)*exp(-abs(s)/dc)

  end function friction_exponential

  end module bc_swff