elastic.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 elastic
! Elastic properties database
! For each element there is two possible databases:
!   1. homogeneous
!   2. heterogeneous

  use distribution_general, only: distribution_type

  implicit none
  private
  
  type elast_type
    private
    integer :: ndof
    double precision, pointer :: a(:,:,:,:) => null() ! coefs for internal forces
    type(mat_elem_type), pointer :: elem(:) => null() ! element number --> material properties
    type(mat_input_type), pointer :: input(:) => null() ! material index --> input properties
  end type elast_type

  type mat_input_type
    integer :: kind ! 0= homogeneous element
                    ! 1= heterogeneous element, pointwise definition
    type(matpro_type), pointer :: homo => null()  ! elementwise
    type(hete_input_type), pointer :: hete => null() 
  end type mat_input_type

  type mat_elem_type
    type(matpro_type), pointer :: homo => null(), hete(:,:) => null()
  end type mat_elem_type

  type matpro_type  ! material properties
    double precision :: cp=0d0,cs=0d0,rho=0d0,coef(4)=0d0
  end type matpro_type

  type hete_input_type  ! heterogeneous domain, info needed to construct the model
    type(distribution_type) :: cp,cs,rho
  end type hete_input_type

  interface ELAST_inquire
    module procedure ELAST_inquire_node, ELAST_inquire_element !, &
!      ELAST_inquire_domain
  end interface ELAST_inquire

  public :: elast_type,ELAST_read,ELAST_init,ELAST_init_mass,ELAST_inquire &
           ,ELAST_cpminmax,ELAST_csminmax &
           ,ELAST_KD, ELAST_strain_stress, ELAST_inquire_domain

contains


!=======================================================================
!
! BEGIN INPUT BLOCK
!
! NAME   : MATERIAL
! PURPOSE: Define elastic material properties of a tagged domain
! SYNTAX : &MATERIAL tag, mode /
!          Followed by material data, with format depending on the mode (see
!          below)
!
! ARG: tag      [int] [none] The number assigned to a domain during mesh generation
! ARG: mode     [char*5] ['ISOTR'] Type of material and/or spatial distribution.
!               The following modes are implemented and this is their data
!               format:
!
!               'ISOTR' homogeneous isotropic elastic
!               One line, dble(3):
!               density, P-wave-velocity, S-wave-velocity
!
!               'ANISO' homogeneous anisotropic
!               One line, dble(5):
!               density, c11, c13, c33, c44
!
!               'XXXXX' isotropic with any 2D distribution
!               Three $DIST_XXXXX blocks: 
!               density, P-velocity, S-velocity
!
!
! END INPUT BLOCK

! Read properties of a two-dimensional
! isotropic or anisotropic linear elastic element
!
  subroutine ELAST_read(elast,iin)

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

  integer, intent(in) :: iin
  type (elast_type), intent(inout) :: elast

  integer :: i,numat,ntags,tag
  character(10)  :: mode

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

  NAMELIST / MATERIAL / tag,mode

 ! Count material sets 
  numat = 0
  ntags = 0
  rewind(iin)
  do 
    read(iin,MATERIAL,END=50) 
    ntags = max(tag,ntags)
    numat = numat+1
  enddo
  50 if (numat==0)  call IO_abort('MATERIAL parameters not found')
  if (echo_input) write(iout,100) numat

 ! Read properties for each set 
  allocate(elast%input(ntags))
  rewind(iin) 
  do i = 1,numat
    mode = 'ISOTR'
    read(iin,MATERIAL)
    if (echo_input) write(iout,200) tag,mode
    select case(mode)
     
      case('ISOTR') ! Isotropic material: RHO, VP and VS given 
        elast%input(tag)%kind = 0
        allocate(elast%input(tag)%homo)
        call ISOTR_read(elast%input(tag)%homo,iin)

      case('ANISO') ! anisotropic material: RHO,c11, c13, c33 et c44 given (Pa)
        elast%input(tag)%kind = 0
        allocate(elast%input(tag)%homo)
        call ANISO_read(elast%input(tag)%homo,iin)
       
      case default ! isotropic with arbitrary distribution
        elast%input(tag)%kind = 1
        allocate(elast%input(tag)%hete)
        call HETE_read(elast%input(tag)%hete,mode,iin)

!      case default 
!        call IO_abort('Improper value while reading material sets')

    end select

  enddo

  return

!---- formats
  100   format(//,' M a t e r i a l   s e t s :   2 D  e l a s t i c i t y', &
         /1x,54('='),//5x, &
         'Number of material sets . . . . . . . . . = ',i5)
  200   format(/5x, &
         'Material number . . . . . . . . . . (tag) = ',i5/5x, &
         'Type  . . . . . . . . . . . . . . .(mode) = ',a)

  end subroutine ELAST_read

!=======================================================================
! Homogeneous block, user input
! Isotropic :  lambda, mu, K (= lambda + 2*mu), zero
  subroutine ISOTR_read(homo,iin)
  
  use echo, only : echo_input, iout
  use stdio, only : IO_abort

  type(matpro_type), intent(out) :: homo
  integer, intent(in) :: iin

  double precision :: Kmod,Kvol,young,poiss,denst,cp,cs,mu,lambda,cs2,cp2

  read(iin ,*) denst,cp,cs
  cp2 = cp*cp
  cs2 = cs*cs
  mu   = denst*cs2
  Kmod  = denst*cp2 
  lambda  = denst*(cp2 - 2d0*cs2)

  poiss = 0.5d0*(cp2-2d0*cs2)/(cp2-cs2)
  if (poiss < 0.d0 .or. poiss > 0.5d0) call IO_abort('Poisson''s ratio out of range !')

  if(echo_input) then
    Kvol  = lambda + 2d0*mu/3.d0
    young = 2d0*mu*(1d0+poiss)
    write(iout,200) cp,cs,denst,poiss,lambda,mu,Kvol,young
  endif

  homo%coef(1) = lambda
  homo%coef(2) = mu
  homo%coef(3) = Kmod
  homo%coef(4) = Kmod
  homo%cp   = cp
  homo%cs   = cs
  homo%rho  = denst

  200   format(5x, &
         'P-wave velocity . . . . . . . . . . .(cp) =',EN12.3,/5x, &
         'S-wave velocity . . . . . . . . . . .(cs) =',EN12.3,/5x, &
         'Mass density. . . . . . . . . . . (denst) =',EN12.3,/5x, &
         'Poisson''s ratio . . . . . . . . . (poiss) =',EN12.3,/5x, &
         'First Lame parameter Lambda . . . .(alam) =',EN12.3,/5x, &
         'Second Lame parameter Mu. . . . . . (amu) =',EN12.3,/5x, &
         'Bulk modulus K. . . . . . . . . . .(Kvol) =',EN12.3,/5x, &
         'Young''s modulus E. . . . . . . . (young) =',EN12.3)

  end subroutine ISOTR_read

!=======================================================================
! Homogeneous block, user input
! Transverse anisotropic :  c11, c13, c33, c44
!
  subroutine ANISO_read(homo,iin)

  use echo, only : echo_input, iout

  type(matpro_type), intent(out) :: homo
  integer, intent(in) :: iin

  double precision :: denst,c11,c13,c33,c44,cp,cs

  read(iin ,*) denst, c11, c13, c33, c44

 ! juste une idee des proprietes
  cp  = sqrt(c11/denst)
  cs  = sqrt(c44/denst)

  if(echo_input) write(iout,300) c11,c13,c33,c44,denst, &
  sqrt(c33/denst),sqrt(c11/denst),sqrt(c44/denst),sqrt(c44/denst)

  homo%coef(1) = c13
  homo%coef(2) = c44
  homo%coef(3) = c11
  homo%coef(4) = c33
  homo%cp     = cp
  homo%cs     = cs
  homo%rho    = denst

 300   format(5x, &
         'c11 coefficient (Pascal). . . . . . (c11) =',EN12.3,/5x, &
         'c13 coefficient (Pascal). . . . . . (c13) =',EN12.3,/5x, &
         'c33 coefficient (Pascal). . . . . . (c33) =',EN12.3,/5x, &
         'c44 coefficient (Pascal). . . . . . (c44) =',EN12.3,/5x, &
         'Mass density. . . . . . . . . . . (denst) =',EN12.3,/5x, &
         'Velocity of qP along vertical axis. . . . =',EN12.3,/5x, &
         'Velocity of qP along horizontal axis. . . =',EN12.3,/5x, &
         'Velocity of qSV along vertical axis . . . =',EN12.3,/5x, &
         'Velocity of qSV along horizontal axis . . =',EN12.3)

  end subroutine ANISO_read

!=======================================================================
! Heterogeneous block
  subroutine HETE_read(hete,mode,iin)

  use distribution_general, only : DIST_read

  type(hete_input_type), intent(out) :: hete
  character(*), intent(in) :: mode
  integer, intent(in) :: iin

  call DIST_read(hete%rho,mode,iin)
  call DIST_read(hete%cp,mode,iin)
  call DIST_read(hete%cs,mode,iin)

  end subroutine HETE_read


!=======================================================================
!
!  Define arrays a1 to a10 for the computation of elastic internal forces 
!
! mode = 1 : SH
!      = 2 : P-SV
  subroutine ELAST_init(elast,grid,mode)

  use memory_info
  use spec_grid, only : sem_grid_type, SE_InverseJacobian, SE_VolumeWeights
  use stdio, only : IO_abort
  use echo, only : echo_init,iout,fmt1,fmtok
  use constants, only : OPT_HOMOG

  type(elast_type), intent(inout), target :: elast
  type(sem_grid_type), intent(in) :: grid
  integer, intent(in) :: mode

  double precision, dimension(grid%ngll,grid%ngll,4), target :: coefs
  double precision, dimension(2,2,grid%ngll,grid%ngll), target :: xjaci
  double precision, dimension(grid%ngll,grid%ngll) :: weights
  double precision, dimension(:,:), pointer :: DxiDx,DxiDz,DetaDx,DetaDz
  double precision, dimension(:,:), pointer :: Kx,Kz,la,mu
  integer :: k,e,nelast,tag,nelem_a

  if (echo_init) then
    write(iout,*) 
    write(iout,'(a)') ' M a t e r i a l   p r o p e r t i e s'
    write(iout,'(a)') ' ====================================='
    write(iout,*) 
    write(iout,fmt1,advance='no') 'Translating input velocity model'
  endif

  allocate(elast%elem(grid%nelem))
  do e=1,grid%nelem
    tag = grid%tag(e)
    if ( tag > size(elast%input) .or. tag<1 ) &
      call IO_abort('ELAST_init: tag in mesh does not correspond to a defined material')
    if ( elast%input(tag)%kind == 0 ) then
      elast%elem(e)%homo => elast%input(tag)%homo
    else
      allocate( elast%elem(e)%hete(grid%ngll,grid%ngll) )
      call HETE_init_elem( elast%elem(e)%hete, elast%input(tag)%hete, grid%ibool(:,:,e), grid%coord)
    endif
  enddo
  if (echo_init) write(iout,fmtok)

  if (echo_init) write(iout,fmt1,advance='no') 'Defining elasticity work arrays'
  if (mode==1) then
    nelast = 3
    if (grid%flat) nelast = 2
  elseif (mode==2) then
    nelast = 10
    if (grid%flat) nelast = 6
  else
    call IO_abort('ELAST_init: mode must be 1 or 2 (SH or P-SV)')
  endif
  elast%ndof = mode
  if (OPT_HOMOG) then
    nelem_a = 1
  else
    nelem_a = grid%nelem
  endif
  allocate(elast%a(grid%ngll,grid%ngll,nelast,nelem_a))
  call storearray('elast%a',size(elast%a),idouble)
  
  do e=1,nelem_a 
   
    xjaci = SE_InverseJacobian(grid,e)
    DxiDx  => xjaci(1,1,:,:)
    DxiDz  => xjaci(1,2,:,:)
    DetaDx => xjaci(2,1,:,:)
    DetaDz => xjaci(2,2,:,:)
    
    call ELAST_inquire(elast,e,coefs=coefs)
    la => coefs(:,:,1) ! c13
    mu => coefs(:,:,2) ! c44
    Kx => coefs(:,:,3) ! c11
    Kz => coefs(:,:,4) ! c33

    select case(nelast)

    case(2) ! SH flat
      elast%a(:,:,1,e) = mu * DxiDx*DxiDx
      elast%a(:,:,2,e) = mu * DetaDz*DetaDz
      
    case(3) ! SH general
      elast%a(:,:,1,e) = mu *( DxiDx*DxiDx + DxiDz*DxiDz )
      elast%a(:,:,2,e) = mu *( DetaDx*DetaDx + DetaDz*DetaDz )
      elast%a(:,:,3,e) = mu *( DxiDx*DetaDx + DxiDz*DetaDz )