elastic.f90

Spectral Element Method for wave propagation and rupture dynamics.

    case(6) ! P-SV flat
      elast%a(:,:,1,e) = Kx * DxiDx*DxiDx
      elast%a(:,:,2,e) = la * DxiDx*DetaDz
      elast%a(:,:,3,e) = Kz * DetaDz*DetaDz
      elast%a(:,:,4,e) = mu * DetaDz*DetaDz
      elast%a(:,:,5,e) = mu * DxiDx*DetaDz
      elast%a(:,:,6,e) = mu * DxiDx*DxiDx
  
    case(10) ! P-SV general
      elast%a(:,:,1,e) = Kx * DxiDx*DxiDx   + mu * DxiDz*DxiDz
      elast%a(:,:,2,e) = la * DxiDx*DetaDz  + mu * DxiDz*DetaDx
      elast%a(:,:,3,e) = Kz * DetaDz*DetaDz + mu * DetaDx*DetaDx
      elast%a(:,:,4,e) = Kx * DetaDx*DetaDx + mu * DetaDz*DetaDz
      elast%a(:,:,5,e) = la * DxiDz*DetaDx  + mu * DxiDx*DetaDz
      elast%a(:,:,6,e) = Kz * DxiDz*DxiDz   + mu * DxiDx*DxiDx
      elast%a(:,:,7,e)  = Kx * DxiDx*DetaDx  + mu * DxiDz*DetaDz
      elast%a(:,:,8,e)  = (la+mu) * DxiDx*DxiDz
      elast%a(:,:,9,e)  = (la+mu) * DetaDx*DetaDz
      elast%a(:,:,10,e) = Kz * DxiDz*DetaDz  + mu * DxiDx*DetaDx

    end select

    weights = SE_VolumeWeights(grid,e)
    do k=1,nelast
      elast%a(:,:,k,e) = - weights*elast%a(:,:,k,e)
    enddo

  enddo


  if (echo_init) write(iout,fmtok)
   
  end subroutine ELAST_init

!=====================================================================
  subroutine HETE_init_elem(matpro,matin,ibool,gcoord)

  use distribution_general, only : DIST_generate

  type(matpro_type), intent(out) :: matpro(:,:)
  type(hete_input_type), intent(in) :: matin
  integer, intent(in) :: ibool(:,:)
  double precision, intent(in) :: gcoord(:,:)

  ! DIST_generate work with vector input/outputs
  double precision :: estuff(size(ibool)),ecoord(2,size(ibool))
  integer :: sha(2)
  
  sha = shape(ibool)
  ecoord = gcoord(:, pack(ibool,mask=.true.) )
  call DIST_generate( estuff, ecoord, matin%cp ) 
  matpro%cp = reshape(estuff,sha)
  call DIST_generate( estuff, ecoord, matin%cs ) 
  matpro%cs = reshape(estuff,sha)
  call DIST_generate( estuff, ecoord, matin%rho ) 
  matpro%rho= reshape(estuff,sha)

  matpro%coef(1) = matpro%rho*(matpro%cp**2-2d0*matpro%cs**2)
  matpro%coef(2) = matpro%rho*matpro%cs**2
  matpro%coef(3) = matpro%rho*matpro%cp**2
  matpro%coef(4) = matpro%coef(3)
      
  end subroutine HETE_init_elem

!=====================================================================
!
! Allocate and compute the mass matrix by summing the contribution of each point
!
  subroutine ELAST_init_mass(mass,grid,elast,ndof)

  use memory_info
  use spec_grid, only : sem_grid_type, SE_VolumeWeights

  double precision, pointer :: mass(:,:)
  type(sem_grid_type), intent(in) :: grid
  type(elast_type), intent(in) :: elast
  integer, intent(in) :: ndof
  
  double precision :: massloc(grid%ngll,grid%ngll)
  integer :: e,i,j,iglob

  allocate(mass(grid%npoin,ndof))
  call storearray('mass',size(mass),idouble)
  mass = 0.d0

  do e = 1,grid%nelem
    call ELAST_inquire(elast,e,rho=massloc)
    massloc = massloc * SE_VolumeWeights(grid,e)
    do j=1,grid%ngll
    do i=1,grid%ngll
      iglob = grid%ibool(i,j,e)
      mass(iglob,1) = mass(iglob,1) + massloc(i,j)
    enddo
    enddo
  enddo

  if (ndof==2) mass(:,2) = mass(:,1)

  end subroutine ELAST_init_mass


!=======================================================================
!
  subroutine ELAST_inquire_node(elast,i,j,e,rho,coefs,cp,cs,mu,lambda)

    type(elast_type), intent(in) :: elast
    integer, intent(in) :: e,i,j
    double precision, optional, intent(out) :: coefs(4),rho,cp,cs,mu,lambda
    
    type (matpro_type), pointer :: mato

    if ( .not.associated(elast%elem(e)%hete) ) then
      mato => elast%elem(e)%homo
    else
      mato => elast%elem(e)%hete(i,j)
    endif  

    if (present(rho))   rho   = mato%rho
    if (present(coefs)) coefs = mato%coef
    if (present(lambda)) lambda = mato%coef(1)
    if (present(mu))    mu    = mato%coef(2)
    if (present(cp))    cp    = mato%cp
    if (present(cs))    cs    = mato%cs

  end subroutine ELAST_inquire_node

!=======================================================================
!
  subroutine ELAST_inquire_element(elast,e,rho,coefs,cp,cs,mu,lambda)

    type(elast_type), intent(in) :: elast
    integer, intent(in) :: e
    double precision, optional, dimension(:,:,:), intent(out) :: coefs
    double precision, optional, dimension(:,:), intent(out) :: &
      rho,cp,cs,mu,lambda
    
    type (matpro_type), pointer :: mato1,mato2(:,:)

    if ( .not.associated(elast%elem(e)%hete) ) then
      mato1 => elast%elem(e)%homo
      if (present(rho))   rho   = mato1%rho
      if (present(coefs)) then
        coefs(:,:,1) = mato1%coef(1)
        coefs(:,:,2) = mato1%coef(2)
        coefs(:,:,3) = mato1%coef(3)
        coefs(:,:,4) = mato1%coef(4)
      endif
      if (present(lambda)) lambda = mato1%coef(1)
      if (present(mu))    mu    = mato1%coef(2)
      if (present(cp))    cp    = mato1%cp
      if (present(cs))    cs    = mato1%cs
    else
      mato2 => elast%elem(e)%hete(:,:)
      if (present(rho))   rho   = mato2%rho
      if (present(coefs)) then
        coefs(:,:,1) = mato2%coef(1)
        coefs(:,:,2) = mato2%coef(2)
        coefs(:,:,3) = mato2%coef(3)
        coefs(:,:,4) = mato2%coef(4)
      endif
      if (present(lambda)) lambda = mato2%coef(1)
      if (present(mu))    mu    = mato2%coef(2)
      if (present(cp))    cp    = mato2%cp
      if (present(cs))    cs    = mato2%cs
    endif  


  end subroutine ELAST_inquire_element


!=======================================================================
!
  subroutine ELAST_inquire_domain(elast,ndof)

    type(elast_type), intent(in) :: elast
    integer, intent(out) :: ndof

    ndof = elast%ndof

  end subroutine ELAST_inquire_domain

!=======================================================================
  subroutine ELAST_cpminmax(elast,cpmin,cpmax)

    type(elast_type), intent(in) :: elast
    double precision, intent(out) :: cpmin,cpmax

    integer :: i,e

    cpmax = -huge(cpmax)
    cpmin = huge(cpmin)

    do i=1,size(elast%input)
    if (elast%input(i)%kind == 0) then
      cpmax = max(cpmax, elast%input(i)%homo%cp )
      cpmin = min(cpmin, elast%input(i)%homo%cp )
    endif  
    enddo

    do e=1,size(elast%elem)
    if ( associated(elast%elem(e)%hete) ) then
      cpmax = max(cpmax,maxval(elast%elem(e)%hete%cp))
      cpmin = min(cpmin,minval(elast%elem(e)%hete%cp))
    endif      
    enddo

  end subroutine ELAST_cpminmax

!=======================================================================
  subroutine ELAST_csminmax(elast,csmin,csmax)

    type(elast_type), intent(in) :: elast
    double precision, intent(out) :: csmin,csmax

    integer :: i,e

    csmax = -huge(csmax)
    csmin = huge(csmin)

    do i=1,size(elast%input)
    if (elast%input(i)%kind == 0) then
      csmax = max(csmax, elast%input(i)%homo%cs )
      csmin = min(csmin, elast%input(i)%homo%cs )
    endif  
    enddo

    do e=1,size(elast%elem)
    if ( associated(elast%elem(e)%hete) ) then
      csmax = max(csmax,maxval(elast%elem(e)%hete%cs))
      csmin = min(csmin,minval(elast%elem(e)%hete%cs))
    endif      
    enddo

  end subroutine ELAST_csminmax


!=======================================================================
!
! Computes the elastic internal forces term = -K.displ 
! in a SEM grid using the coefficients in elast.
! On output the result is stored in the field KD (scratched)
!
! Number of multiplications per GLL node ( = total / ngll^2*nelem) :
!                       SH              P-SV
!       flat            4*ngll +2       8*ngll +8
!       general         4*ngll +4       8*ngll +16
!
  subroutine ELAST_KD(elast,grid,displ,KD)

  use spec_grid, only : sem_grid_type
  use constants, only : OPT_NGLL 

  type (elast_type)   , intent(in) :: elast
  type (sem_grid_type), intent(in) :: grid
  double precision    , intent(in) :: displ(:,:)
  double precision    , intent(out):: KD(:,:)

  if ( elast%ndof==1 ) then
    if (OPT_NGLL==grid%ngll) then
      call ELAST_KD2_SH(displ(:,1),KD(:,1),elast%a,grid%Hprime,grid%Htprime,grid%ibool &
                       ,size(elast%a,3),grid%nelem,grid%npoin)
    else
      call ELAST_KD1_SH(elast,grid,displ(:,1),KD(:,1))
    endif
  else
    if (OPT_NGLL==grid%ngll) then
      call ELAST_KD2_PSV(displ,KD,elast%a,grid%Hprime,grid%Htprime,grid%ibool &
                       ,size(elast%a,3),grid%nelem,grid%npoin)
    else
      call ELAST_KD1_PSV(elast,grid,displ,KD)
    endif
  endif
 
  end subroutine ELAST_KD

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

  subroutine ELAST_KD1_PSV(elast,grid,displ,KD)

  use spec_grid, only : sem_grid_type

  type (elast_type)   , intent(in) :: elast
  type (sem_grid_type), intent(in) :: grid
  double precision    , intent(in) :: displ(:,:)
  double precision    , intent(out):: KD(:,:)

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

  KD = 0.d0
  do e=1,grid%nelem

!-- Elementwise storage of DISPL field
    do j=1,grid%ngll
    do i=1,grid%ngll
      k = grid%ibool(i,j,e)
      Uxloc(i,j) = displ(k,1)
      Uzloc(i,j) = displ(k,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 )

!-- Elementwise forces

   if (grid%flat) then

    tmp = elast%a(:,:,1,e)*dUx_dxi + elast%a(:,:,2,e)*dUz_deta
    Uxloc = MATMUL( grid%hprime, tmp )

    tmp = elast%a(:,:,4,e)*dUx_deta + elast%a(:,:,5,e)*dUz_dxi
    Uxloc = Uxloc + MATMUL( tmp, grid%hTprime )

    tmp = elast%a(:,:,5,e)*dUx_deta + elast%a(:,:,6,e)*dUz_dxi
    Uzloc = MATMUL( grid%hprime, tmp )

    tmp = elast%a(:,:,2,e)*dUx_dxi + elast%a(:,:,3,e)*dUz_deta
    Uzloc = Uzloc + MATMUL( tmp, grid%hTprime )

   else

    tmp = elast%a(:,:,1,e)*dUx_dxi + elast%a(:,:,7,e)*dUx_deta &
        + elast%a(:,:,8,e)*dUz_dxi + elast%a(:,:,2,e)*dUz_deta
    Uxloc = MATMUL( grid%hprime, tmp )

    tmp = elast%a(:,:,7,e)*dUx_dxi + elast%a(:,:,4,e)*dUx_deta &
        + elast%a(:,:,5,e)*dUz_dxi + elast%a(:,:,9,e)*dUz_deta
    Uxloc = Uxloc + MATMUL( tmp, grid%hTprime )

    tmp = elast%a(:,:,8,e)*dUx_dxi + elast%a(:,:,5,e)*dUx_deta &
        + elast%a(:,:,6,e)*dUz_dxi + elast%a(:,:,10,e)*dUz_deta
    Uzloc = MATMUL( grid%hprime, tmp )

    tmp = elast%a(:,:,2,e)*dUx_dxi + elast%a(:,:,9,e)*dUx_deta &
        + elast%a(:,:,10,e)*dUz_dxi + elast%a(:,:,3,e)*dUz_deta
    Uzloc = Uzloc + MATMUL( tmp, grid%hTprime )

   endif

!-- Assemble internal forces
    do j=1,grid%ngll
    do i=1,grid%ngll
      k = grid%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_KD1_PSV

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

  subroutine ELAST_KD1_SH(elast,grid,displ,KD)

  use spec_grid, only : sem_grid_type
  use constants, only : OPT_HOMOG

  type (elast_type)   , intent(in) :: elast
  type (sem_grid_type), intent(in) :: grid
  double precision    , intent(in) :: displ(:)
  double precision    , intent(out):: KD(:)

  double precision, dimension(grid%ngll,grid%ngll) :: &
    Uloc, hprime,htprime, tmp, dU_dxi, dU_deta
  integer :: i,j,e,ea,k

  hprime = grid%hprime
  htprime = grid%hTprime

  KD = 0.d0
  do e=1,grid%nelem

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

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

!-- Elementwise forces