spec_grid.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 spec_grid

  use stdio, only : IO_abort
  use constants
  use fem_grid
  use bnd_grid

  implicit none
  private


!-----------------------------------------------------------------------
!-- Spectral element grid type
!
!------ Topology
!     nelem        = total number of elements
!     ngll         = number of Gauss-Lobato-Legendre nodes per unit segment
!     npoin        = total number of GLL nodes in the mesh
!     ibool        = local to global numbering for the GLL mesh
!                    (igll,jgll,element) -> bulk node index
!     coord        = coordinates of the GLL nodes
!
!------ Working data
!     shape        = control nodes to GLL nodes shape functions
!     dshape       = derivatives of shape functions
!     hprime       = lagrange derivatives on the unit square
!     hTprime      = transpose of hprime
!     wgll         = GLL integration weights
!     wgll2        = 2D GLL integration weights
!     xgll         = GLL points on the unit segment
!
!------ Other
!     tag          = element -> domain tags
!     flat         = flag for cartesian grid (allows simplifications)
!     fmax         = Highest frequency to be resolved by the grid
!
!
! NOTE: After initialization, the FEM mesh database is only used for plotting
!       purposes, so if no plots are needed you can save them on disk and 
!       clean it from memory.
!
! NOTE: COORD is used in the solver only for computation of the incident
!       wavefield, the boundary coordinates could be stored instead in the 
!       bc structure (boundary_condition).


  type sem_grid_type
    integer :: ngll=0,nelem=0,npoin=0
    type(fem_grid_type) :: fem
    double precision :: fmax=0d0
    double precision, pointer :: coord(:,:)      =>null(), &
                                 hprime(:,:)     =>null(), &
                                 hTprime(:,:)    =>null(), &
                                 wgll(:)         =>null(), &
                                 xgll(:)         =>null(), &
                                 wgll2(:,:)      =>null(), &
                                 shape(:,:,:)    =>null(), &
                                 dshape(:,:,:,:) =>null()
    integer, pointer :: ibool(:,:,:)  =>null(), &
                        tag(:) =>null()
    type (bnd_grid_type), pointer :: bounds(:) =>null()
    logical :: flat=.false.
  end type sem_grid_type

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

  interface SE_Jacobian
    module procedure SE_Jacobian_e, SE_Jacobian_eij
  end interface SE_Jacobian

  interface SE_InverseJacobian
    module procedure SE_InverseJacobian_e, SE_InverseJacobian_eij
  end interface SE_InverseJacobian

  interface SE_VolumeWeights
    module procedure SE_VolumeWeights_e, SE_VolumeWeights_eij
  end interface SE_VolumeWeights

  ! (i,j) GLL for element corner nodes
  !icorner(1) = 1    ; jcorner(1) = 1
  !icorner(2) = ngll ; jcorner(2) = 1
  !icorner(3) = ngll ; jcorner(3) = ngll
  !icorner(4) = 1    ; jcorner(4) = ngll


  public :: sem_grid_type ,&
         SE_init ,&
         SE_init_interpol ,&
         SE_find_nearest_node ,&
         SE_node_belongs_to, &
         SE_find_point ,&
         SE_get_edge_nodes ,&
         SE_inquire ,&
         SE_VolumeWeights ,&
         SE_InverseJacobian ,&
         SE_Jacobian ,&
         BC_inquire ,&
         BC_get_normal_and_weights &
          , edge_D,edge_R,edge_U,edge_L &
          , side_D,side_R,side_U,side_L
 
contains

!=======================================================================
!
  subroutine SE_init(se)

  use stdio, only : IO_new_unit
  type (sem_grid_type), intent(inout) :: se
  integer :: ounit

  call SE_init_gll(se) ! GLL quadrature points, weights 
                            ! and lagrange coefficients,
                            ! shape functions and their derivatives
  call SE_init_numbering(se) ! generate the global node numbering
  call SE_init_coord(se) ! get the coordinates of the GLL nodes,
                              ! and export the grid to a file
  call SE_BcTopoInit(se) ! initialize boundaries generic data

  ! export grid parameters
  ounit = IO_new_unit()
  open(ounit,file='grid_sem2d.hdr',status='replace')
  write(ounit,*) 'NELEM  NPGEO  NGNOD  NPOIN  NGLL'
  write(ounit,*) se%nelem, FE_GetNbNodes(se%fem), FE_GetNodesPerElement(se%fem), &
                 se%npoin, se%ngll
  close(ounit)

  end subroutine SE_init

!=======================================================================
!
! get coordinates and weights of the Gauss-Lobatto-Legendre points
! and derivatives of Lagrange polynomials  H_ij = h'_i(xgll(j))
  subroutine SE_init_gll(se)

  use gll, only : get_GLL_info

  type (sem_grid_type), intent(inout) :: se
  double precision :: xi,eta 
  integer :: ngll,i,j, ngnod

  ngll = se%ngll

  allocate(se%xgll(ngll))
  allocate(se%wgll(ngll))
  allocate(se%hprime(ngll,ngll))
  allocate(se%hTprime(ngll,ngll))
  allocate(se%wgll2(ngll,ngll))

  call get_GLL_info(ngll,se%xgll,se%wgll,se%hprime)
  se%hTprime = transpose(se%hprime)
  
 ! wgll2(i,j) = wgll(i) * wgll(j)
  do j = 1,ngll
    se%wgll2(:,j) = se%wgll * se%wgll(j)
  enddo

!-----------------------------------------------------------------------
!-- set up the shape functions and their local derivatives
  ngnod = FE_GetNodesPerElement(se%fem)

  allocate(se%shape(ngnod,ngll,ngll))
  allocate(se%dshape(ngnod,NDIME,ngll,ngll))

  do j = 1,ngll
    eta = se%xgll(j)
    do i = 1,ngll
      xi = se%xgll(i)
      se%shape(:,i,j) = FE_getshape(xi,eta,ngnod)
      se%dshape(:,:,i,j) = FE_getdershape(xi,eta,ngnod)
    enddo
  enddo

  end subroutine SE_init_gll



!=======================================================================
!
!-- generate the global numbering
!
  subroutine SE_init_numbering(se)

  use memory_info
  use echo, only : echo_init,iout,fmt1,fmtok
  use stdio, only : IO_new_unit

  type(sem_grid_type), intent(inout) :: se

  integer, pointer :: ibool(:,:,:)
  integer, dimension(se%ngll,4) :: iedg,jedg,iedgR,jedgR
  integer, dimension(4) :: ivtx,jvtx
  integer :: i,j,k,e,n,ii,jj,ee,nn,ngll,npoin,iol,ounit
  integer, pointer :: ees(:),nns(:)

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

  if (echo_init) then
    write(iout,*) 
    write(iout,'(a)') ' S p e c t r a l   e l e m e n t s   g r i d'
    write(iout,'(a)') ' ==========================================='
    write(iout,*) 
  endif

  ngll  = se%ngll
  se%nelem = FE_GetNbElements(se%fem)
  se%tag => se%fem%tag

! global numbering table
  allocate(se%ibool(ngll,ngll,se%nelem))
  call storearray('ibool',size(se%ibool),iinteg)
  ibool => se%ibool
  ibool = 0

!---- start numbering
  if (echo_init) write(iout,fmt1,advance='no') 'Numbering GLL points'
  npoin = 0

 ! GLL index (i,j) of edge nodes, counterclockwise
  do n = 1,4
    call SE_inquire(se,edge=n,itab=iedg(:,n),jtab=jedg(:,n))
  enddo
 ! reverse order for matching edges
  iedgR = iedg(ngll:1:-1,:)
  jedgR = jedg(ngll:1:-1,:)

 ! GLL index (i,j) for vertex nodes
  ivtx(1) = 1;    jvtx(1) = 1
  ivtx(2) = ngll; jvtx(2) = 1
  ivtx(3) = ngll; jvtx(3) = ngll
  ivtx(4) = 1;    jvtx(4) = ngll

  do e = 1,se%nelem

   !-- interior nodes are unique
    do j=2,ngll-1
    do i=2,ngll-1
      npoin = npoin + 1
      ibool(i,j,e) = npoin
    enddo
    enddo
  
   !-- edge nodes
    do n = 1,4
      if ( ibool(iedg(2,n),jedg(2,n),e)>0 ) cycle ! skip if already processed
      call FE_GetEdgeConn(se%fem,e,n,ee,nn)
      do k = 2,ngll-1
        npoin = npoin+1
        ibool(iedg(k,n),jedg(k,n),e) = npoin
        if (ee>0) ibool(iedgR(k,nn),jedgR(k,nn),ee) = npoin
      enddo
    enddo

   !-- vertex nodes
    do n = 1,4
      i = ivtx(n)
      j = jvtx(n)
      if ( ibool(i,j,e) > 0 ) cycle
      npoin = npoin +1
     ! scan Vertex_Conn for shared elements
      call FE_GetVertexConn(se%fem,e,n,ees,nns)
      do k= 1,size(ees)
        ii = ivtx(nns(k))
        jj = jvtx(nns(k))
        ibool(ii,jj,ees(k)) = npoin
      enddo
    enddo 

  enddo

  call FE_UnsetConnectivity(se%fem)

  se%npoin = npoin

  if (echo_init) then
    write(iout,fmtok)
    write(iout,100) 'Total number of GLL points. . . . . . . = ',npoin
    write(iout,*)
    write(iout,fmt1,advance='no') 'Saving element/node table in binary file ibool_sem2d.dat'
  endif

  inquire( IOLENGTH=iol ) se%ibool
  ounit = IO_new_unit()
  open(ounit,file='ibool_sem2d.dat',status='replace',access='direct',recl=iol)
  write(ounit,rec=1) se%ibool
  close(ounit)
  if (echo_init) write(iout,fmtok)

  return

100 format(5X,A,I0)
  
  end subroutine SE_init_numbering

  
!=======================================================================
!
!  set the global nodal coordinates
!
  subroutine SE_init_coord(se)

  use echo, only : echo_check,echo_init,iout,fmt1,fmtok
  use memory_info
  use stdio, only : IO_new_unit

  type (sem_grid_type), intent(inout) :: se

  integer :: i,j,e,ounit,iol
  double precision, pointer :: coorg(:,:)

  allocate(se%coord(NDIME,se%npoin))
  call storearray('coord',size(se%coord),idouble)

!---- Coordinates of the global points 
  if (echo_init) write(iout,fmt1,advance='no') 'Defining nodes coordinates'
  do e = 1,se%nelem
    coorg => FE_GetElementCoord(se%fem,e)
    do j = 1,se%ngll
    do i = 1,se%ngll
      se%coord(:,se%ibool(i,j,e)) = matmul( coorg, se%shape(:,i,j) )
    enddo
    enddo
    deallocate(coorg)
  enddo
  if (echo_init) write(iout,fmtok)

  if (echo_check) then

!----  Save the grid in a text file
    write(iout,*) 
    write(iout,fmt1,advance='no') 'Saving the grid coordinates (coord) in a text file'
    ounit = IO_new_unit()
    open(ounit,file='coord_sem2d.tab',status='unknown')
    write(ounit,*) se%npoin
    do i = 1,se%npoin
      write(ounit,*) i, se%coord(:,i)
    enddo
    close(ounit)
    write(iout,fmtok)
  endif
    
!----  Always save the grid in a binary file
  if (echo_init) write(iout,fmt1,advance='no') 'Saving the grid coordinates (coord) in a binary file'
  inquire( IOLENGTH=iol ) se%coord
  iol=iol/2
  ounit = IO_new_unit()
  open(ounit,file='coord_sem2d.dat',status='replace',access='direct',recl=iol)
 ! NOTE: crashes here with segmentation fault if se%coord exceeds stack size
 !       FIX: (bash) ulimit -s unlimited
  write(ounit,rec=1) real(se%coord)
  close(ounit)
  if (echo_init) write(iout,fmtok)

  end subroutine SE_init_coord



!=======================================================================
!
! lagrange interpolation functions at (xi,eta)
! interp(ngll*ngll) is then used as
!   interpolated_value = matmul(interp,nodal_values)
!   interpolated_vector(1:ndof) = matmul(interp,nodal_vectors(:,1:ndof))
!
  subroutine SE_init_interpol(xi,eta,interp,grid)

  use gll, only : hgll

  double precision, intent(in) :: xi,eta
  type (sem_grid_type), intent(in) :: grid
  double precision, intent(out) :: interp(grid%ngll*grid%ngll)

  double precision :: fi,fj
  integer :: i,j,k

  k=0
  do j=1,grid%ngll
    fj = hgll(j-1,eta,grid%xgll,grid%ngll)
    do i=1,grid%ngll
      k=k+1
      fi = hgll(i-1,xi,grid%xgll,grid%ngll)
      interp(k) = fi*fj
    enddo
  enddo

  end subroutine SE_init_interpol


!=====================================================================
!
  subroutine SE_find_nearest_node(coord_in,grid,iglob,coord,distance)

  type(sem_grid_type), intent(in)  :: grid
  double precision   , intent(in)  :: coord_in(:)
  double precision, optional, intent(out) :: coord(:),distance
  integer, intent(out) :: iglob

  double precision :: d2min,d2
  integer :: ip

  d2min = huge(d2min)
  do ip = 1,grid%npoin
    d2 = (coord_in(1)-grid%coord(1,ip))**2 + (coord_in(2)-grid%coord(2,ip))**2
    if (d2 <= d2min) then
      d2min  = d2
      iglob = ip
    endif
  enddo

  if (present(coord)) coord(:) = grid%coord(:,iglob)
  if (present(distance)) distance = sqrt(d2min)

  end subroutine SE_find_nearest_node

!=====================================================================
!
  subroutine SE_node_belongs_to(iglob,etab,itab,jtab,grid)

  integer, optional, intent(in) :: iglob
  integer, pointer, dimension(:) :: itab,jtab,etab
  type(sem_grid_type), intent(in)  :: grid