elastic.f90

Spectral Element Method for wave propagation and rupture dynamics.

   if (OPT_HOMOG) then
     ea = 1
   else
     ea = e
   endif

   if (grid%flat) then

    tmp = elast%a(:,:,1,ea)*dU_dxi
    Uloc = MATMUL( hprime, tmp )

    tmp = elast%a(:,:,2,ea)*dU_deta
    Uloc = Uloc + MATMUL( tmp, hTprime )
    
   else

    tmp = elast%a(:,:,1,e)*dU_dxi + elast%a(:,:,3,e)*dU_deta
    Uloc = MATMUL( hprime, tmp )

    tmp = elast%a(:,:,3,e)*dU_dxi + elast%a(:,:,2,e)*dU_deta
    Uloc = Uloc + MATMUL( tmp, hTprime )

   endif

!-- Assemble internal forces
    do j=1,grid%ngll
    do i=1,grid%ngll
      k = grid%ibool(i,j,e)
      KD(k) = KD(k) + Uloc(i,j)
    enddo
    enddo

  enddo

  end subroutine ELAST_KD1_SH


!=======================================================================
! Version 2: OPT_NGLL declared statically, allows for compiler optimizations
! 
  subroutine ELAST_KD2_PSV(displ,KD,a,H,Ht,ibool &
                          ,nelast,nelem,npoin)

  use constants, only : OPT_NGLL 

  integer, intent(in) :: nelast,nelem,npoin, ibool(OPT_NGLL,OPT_NGLL,nelem)
  double precision, intent(in) :: a(OPT_NGLL,OPT_NGLL,nelast,nelem) &
                  , displ(npoin,2), H(OPT_NGLL,OPT_NGLL), Ht(OPT_NGLL,OPT_NGLL)
  double precision, intent(out):: KD(npoin,2)

  double precision, dimension(OPT_NGLL,OPT_NGLL) :: Uxloc, Uzloc &
                  , tmp, dUx_dxi, dUx_deta, dUz_dxi, dUz_deta
  integer :: i,j,e,k

  KD = 0.d0

  do e=1,nelem

!-- Elementwise storage of DISPL field
    do j=1,OPT_NGLL
    do i=1,OPT_NGLL
      k = ibool(i,j,e)
      Uxloc(i,j) = displ(k,1)
      Uzloc(i,j) = displ(k,2)
    enddo
    enddo

!-- Local gradient
    dUx_dxi  = My_MATMUL( Ht, Uxloc )
    dUz_dxi  = My_MATMUL( Ht, Uzloc )
    dUx_deta = My_MATMUL( Uxloc, H )
    dUz_deta = My_MATMUL( Uzloc, H )

!-- Elementwise forces
    if (nelast==6) then

      tmp = a(:,:,1,e)*dUx_dxi + a(:,:,2,e)*dUz_deta
      Uxloc = My_MATMUL( H, tmp )
  
      tmp = a(:,:,4,e)*( dUx_deta + dUz_dxi )
      Uxloc = Uxloc + My_MATMUL( tmp, Ht )

      tmp = a(:,:,5,e)*dUx_deta + a(:,:,6,e)*dUz_dxi
      Uzloc = My_MATMUL( H, tmp )
  
      tmp = a(:,:,2,e)*dUx_dxi + a(:,:,3,e)*dUz_deta
      Uzloc = Uzloc + My_MATMUL( tmp, Ht )

    elseif (nelast==10) then

      tmp = a(:,:,1,e)*dUx_dxi + a(:,:,7,e)*dUx_deta &
          + a(:,:,8,e)*dUz_dxi + a(:,:,2,e)*dUz_deta
      Uxloc = My_MATMUL( H, tmp )
  
      tmp = a(:,:,7,e)*dUx_dxi + a(:,:,4,e)*dUx_deta &
          + a(:,:,5,e)*dUz_dxi + a(:,:,9,e)*dUz_deta
      Uxloc = Uxloc + My_MATMUL( tmp, Ht )
  
      tmp = a(:,:,8,e)*dUx_dxi + a(:,:,5,e)*dUx_deta &
          + a(:,:,6,e)*dUz_dxi + a(:,:,10,e)*dUz_deta
      Uzloc = My_MATMUL( H, tmp )
  
      tmp = a(:,:,2,e)*dUx_dxi + a(:,:,9,e)*dUx_deta &
          + a(:,:,10,e)*dUz_dxi + a(:,:,3,e)*dUz_deta
      Uzloc = Uzloc + My_MATMUL( tmp, Ht )

    endif

!-- Assemble internal forces
    do j=1,OPT_NGLL
    do i=1,OPT_NGLL
      k = ibool(i,j,e)
      KD(k,1) = KD(k,1) + Uxloc(i,j)
      KD(k,2) = KD(k,2) + Uzloc(i,j)
    enddo
    enddo

  enddo

  end subroutine ELAST_KD2_PSV


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

  subroutine ELAST_KD2_SH(displ,KD,a,H,Ht,ibool &
                          ,nelast,nelem,npoin)

  use constants, only : OPT_NGLL 

  integer, intent(in) :: nelast,nelem,npoin, ibool(OPT_NGLL,OPT_NGLL,nelem)
  double precision, intent(in) :: a(OPT_NGLL,OPT_NGLL,nelast,nelem) &
                  , displ(npoin), H(OPT_NGLL,OPT_NGLL), Ht(OPT_NGLL,OPT_NGLL)
  double precision, intent(out):: KD(npoin)

  double precision, dimension(OPT_NGLL,OPT_NGLL) :: Uloc, tmp, dU_dxi, dU_deta
  integer :: i,j,e,k

  KD = 0.d0

  do e=1,nelem

!-- Elementwise storage of DISPL field
    do j=1,OPT_NGLL
    do i=1,OPT_NGLL
      k = ibool(i,j,e)
      Uloc(i,j) = displ(k)
    enddo
    enddo

!-- Local gradient
    dU_dxi  = My_MATMUL( Ht, Uloc )
    dU_deta = My_MATMUL( Uloc, H)

!-- Elementwise forces

   if (nelast==2) then

    tmp = a(:,:,1,e)*dU_dxi
    Uloc = MATMUL( H, tmp )

    tmp = a(:,:,2,e)*dU_deta
    Uloc = Uloc + MATMUL( tmp, Ht )
    
   elseif (nelast==3) then

    tmp = a(:,:,1,e)*dU_dxi + a(:,:,3,e)*dU_deta
    Uloc = MATMUL( H, tmp )

    tmp = a(:,:,3,e)*dU_dxi + a(:,:,2,e)*dU_deta
    Uloc = Uloc + MATMUL( tmp, Ht )

   endif

    
!-- Assemble internal forces
    do j=1,OPT_NGLL
    do i=1,OPT_NGLL
      k = ibool(i,j,e)
      KD(k) = KD(k) + Uloc(i,j)
    enddo
    enddo

  enddo

  end subroutine ELAST_KD2_SH


!----------------------------------
    function My_MATMUL(A,B) result(C)

  use constants, only : OPT_NGLL 

    double precision, dimension(OPT_NGLL,OPT_NGLL), intent(in) :: A,B
    double precision, dimension(OPT_NGLL,OPT_NGLL) :: C

    double precision :: Cij
    integer :: i,j,k

    do j=1,OPT_NGLL
    do i=1,OPT_NGLL
      Cij = 0d0
      do k=1,OPT_NGLL
        Cij = Cij + A(i,k)*B(k,j)
      enddo 
      C(i,j) = Cij
    enddo
    enddo

    end function My_MATMUL


!=======================================================================
!
  subroutine ELAST_strain_stress(elast,grid,displ,store_strain,store_stress,dataout)

  use spec_grid, only : sem_grid_type

  type (elast_type)   , intent(in) :: elast
  type (sem_grid_type), intent(in) :: grid
  logical             , intent(in) :: store_strain,store_stress
  double precision    , intent(in) :: displ(:,:)
  double precision    , intent(out):: dataout(:,:,:,:)
  
  if (.not.(store_strain .or. store_stress)) return
  if ( elast%ndof==1) then
    call ELAST_strain_stress_SH(elast,grid,displ(:,1),store_strain,store_stress,dataout)
  else
    call ELAST_strain_stress_PSV(elast,grid,displ,store_strain,store_stress,dataout)
  endif

  end subroutine ELAST_strain_stress

!-----------------------------------------------------------------------
!
  subroutine ELAST_strain_stress_SH(elast,grid,displ,store_strain,store_stress,dataout)

  use spec_grid, only : sem_grid_type, SE_InverseJacobian

  type (elast_type)   , intent(in) :: elast
  type (sem_grid_type), intent(in) :: grid
  logical             , intent(in) :: store_strain,store_stress
  double precision    , intent(in) :: displ(:)
  double precision    , intent(out):: dataout(:,:,:,:)
  
  double precision, dimension(grid%ngll,grid%ngll) :: &
    Uloc, dU_dxi, dU_deta, e13, e23
  double precision, dimension(:,:), pointer :: dxi_dx, dxi_dz, deta_dx, deta_dz
  double precision, dimension(grid%ngll,grid%ngll,4), target :: coefs
  double precision, dimension(2,2,grid%ngll,grid%ngll), target :: xjaci
  double precision, pointer :: two_mu(:,:)
  integer :: is1,i,j,e

  if (store_strain) then
    is1 = 2
  else
    is1 = 0
  endif

  do e=1,grid%nelem

!-- Elementwise storage of DISPL field
    do j=1,grid%ngll
    do i=1,grid%ngll
      Uloc(i,j) = displ( grid%ibool(i,j,e) )
    enddo
    enddo

!-- Local gradient
    dU_dxi  = MATMUL( grid%hTprime, Uloc )
    dU_deta = MATMUL( Uloc, grid%hprime )

!-- Jacobian matrix
    xjaci = SE_InverseJacobian(grid,e)
    dxi_dx  => xjaci(1,1,:,:)
    dxi_dz  => xjaci(1,2,:,:)
    deta_dx => xjaci(2,1,:,:)
    deta_dz => xjaci(2,2,:,:)

!-- Strain 
    e13 = 0.5d0*( dU_dxi*dxi_dx + dU_deta*deta_dx ) ! = dUy_dx
    e23 = 0.5d0*( dU_dxi*dxi_dz + dU_deta*deta_dz ) ! = dUy_dz

    if (store_strain) then
      dataout(:,:,e,1) = e13
      dataout(:,:,e,2) = e23
    endif

    if (store_stress) then

!-- Shear modulus
      call ELAST_inquire(elast,e,coefs=coefs)
      two_mu => coefs(:,:,2)
      two_mu = 2d0*two_mu

!-- Stress, isotropic medium, antiplane strain
      dataout(:,:,e,is1+1) = two_mu*e13     ! sigma_13
      dataout(:,:,e,is1+2) = two_mu*e23     ! sigma_23

    endif

  enddo

  end subroutine ELAST_strain_stress_SH

!-----------------------------------------------------------------------
!
  subroutine ELAST_strain_stress_PSV(elast,grid,displ,store_strain,store_stress,dataout)

  use spec_grid, only : sem_grid_type, SE_InverseJacobian

  type (elast_type)   , intent(in) :: elast
  type (sem_grid_type), intent(in) :: grid
  logical             , intent(in) :: store_strain,store_stress
  double precision    , intent(in) :: displ(:,:)
  double precision    , intent(out):: dataout(:,:,:,:)
  
  double precision, dimension(grid%ngll,grid%ngll) :: &
    Uxloc, Uzloc, dUx_dxi, dUx_deta, dUz_dxi, dUz_deta, &
    e11, e22, e12, s33
  double precision, dimension(:,:), pointer :: dxi_dx, dxi_dz, deta_dx, deta_dz
  double precision, dimension(grid%ngll,grid%ngll,4), target :: coefs
  double precision, dimension(2,2,grid%ngll,grid%ngll), target :: xjaci
  double precision, pointer :: two_mu(:,:), lambda(:,:)
  integer :: is1,iglob,i,j,e

  if (store_strain) then
    is1 = 3
  else
    is1 = 0
  endif

  do e=1,grid%nelem

!-- Elementwise storage of DISPL field
    do j=1,grid%ngll
    do i=1,grid%ngll
      iglob = grid%ibool(i,j,e)
      Uxloc(i,j) = displ(iglob,1)
      Uzloc(i,j) = displ(iglob,2)
    enddo
    enddo

!-- Local gradient
    dUx_dxi  = MATMUL( grid%hTprime, Uxloc )
    dUz_dxi  = MATMUL( grid%hTprime, Uzloc )
    dUx_deta = MATMUL( Uxloc, grid%hprime )
    dUz_deta = MATMUL( Uzloc, grid%hprime )

!-- Jacobian matrix
    xjaci = SE_InverseJacobian(grid,e)
    dxi_dx  => xjaci(1,1,:,:)
    dxi_dz  => xjaci(1,2,:,:)
    deta_dx => xjaci(2,1,:,:)
    deta_dz => xjaci(2,2,:,:)

!-- Strain 
    e11 = dUx_dxi*dxi_dx + dUx_deta*deta_dx
    e22 = dUz_dxi*dxi_dz + dUz_deta*deta_dz
    e12 = 0.5d0*( dUx_dxi*dxi_dz + dUx_deta*deta_dz  &
                       + dUz_dxi*dxi_dx + dUz_deta*deta_dx  )

    if (store_strain) then
      dataout(:,:,e,1) = e11
      dataout(:,:,e,2) = e22
      dataout(:,:,e,3) = e12
    endif

    if (store_stress) then

!-- Elastic moduli
      call ELAST_inquire(elast,e,coefs=coefs)
      lambda => coefs(:,:,1)
      two_mu => coefs(:,:,2)
      two_mu = 2d0*two_mu

!-- Stress, isotropic medium, plane strain
      s33 = lambda*(e11+e22)
      dataout(:,:,e,is1+1) = s33 +two_mu*e11     ! sigma_11
      dataout(:,:,e,is1+2) = s33 +two_mu*e22     ! sigma_22
      dataout(:,:,e,is1+3) = s33                 ! sigma_33
      dataout(:,:,e,is1+4) = two_mu*e12          ! sigma_12

    endif

  enddo

  end subroutine ELAST_strain_stress_PSV


end module elastic