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

  implicit none
  private

  public :: get_GLL_info,hgll,hdgll,zwgljd

contains

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

  subroutine get_GLL_info(n,x,w,H)

  integer :: i,ip,n
  double precision :: w(n),x(n),H(n,n)

  call zwgljd(x,w,n,0.d0,0.d0)
  !if n is odd make sure the middle point is exactly zero:
  if(mod(n,2) /= 0) x((n-1)/2+1) = 0.d0

  do i=1,n
  do ip=1,n
    H(ip,i)  = hdgll(ip-1,i-1,x,n) ! in hdggl, indices start at 0
  enddo
  enddo

  end subroutine get_GLL_info
!
!=======================================================================
      
  double precision function endw1 (n,alpha,beta)
!
!=======================================================================
!
!     E n d w 1 :
!     ---------
!
!=======================================================================
!
  

  integer n
  double precision alpha,beta

  double precision, parameter :: zero=0.d0,one=1.d0,two=2.d0, &
        three=3.d0,four=4.d0

  double precision apb,f1,fint1,fint2,f2,di,abn,abnn,a1,a2,a3,f3

  integer i
!
!-----------------------------------------------------------------------
!
  f3 = zero
  apb   = alpha+beta
  if (n == 0) then
   endw1 = zero
   return
  endif
  f1   = gammaf(alpha+two)*gammaf(beta+one)/gammaf(apb+three)
  f1   = f1*(apb+two)*two**(apb+two)/two
  if (n == 1) then
   endw1 = f1
   return
  endif
  fint1 = gammaf(alpha+two)*gammaf(beta+one)/gammaf(apb+three)
  fint1 = fint1*two**(apb+two)
  fint2 = gammaf(alpha+two)*gammaf(beta+two)/gammaf(apb+four)
  fint2 = fint2*two**(apb+three)
  f2    = (-two*(beta+two)*fint1 + (apb+four)*fint2) * (apb+three)/four
  if (n == 2) then
   endw1 = f2
   return
  endif
  do i=3,n
   di   = dble(i-1)
   abn  = alpha+beta+di
   abnn = abn+di
   a1   = -(two*(di+alpha)*(di+beta))/(abn*abnn*(abnn+one))
   a2   =  (two*(alpha-beta))/(abnn*(abnn+two))
   a3   =  (two*(abn+one))/((abnn+two)*(abnn+one))
   f3   =  -(a2*f2+a1*f1)/a3
   f1   = f2
   f2   = f3
  enddo
  endw1  = f3

  return
  end function endw1

  double precision function endw2 (n,alpha,beta)
!
!=======================================================================
!
!     E n d w 2 :
!     ---------
!
!=======================================================================
!
  

  integer n
  double precision alpha,beta

  double precision, parameter :: zero=0.d0,one=1.d0,two=2.d0, &
      three=3.d0,four=4.d0

  double precision apb,f1,fint1,fint2,f2,di,abn,abnn,a1,a2,a3,f3

  integer i

!
!-----------------------------------------------------------------------
!
  apb   = alpha+beta
  f3 = zero
  if (n == 0) then
   endw2 = zero
   return
  endif
  f1   = gammaf(alpha+one)*gammaf(beta+two)/gammaf(apb+three)
  f1   = f1*(apb+two)*two**(apb+two)/two
  if (n == 1) then
   endw2 = f1
   return
  endif
  fint1 = gammaf(alpha+one)*gammaf(beta+two)/gammaf(apb+three)
  fint1 = fint1*two**(apb+two)
  fint2 = gammaf(alpha+two)*gammaf(beta+two)/gammaf(apb+four)
  fint2 = fint2*two**(apb+three)
  f2    = (two*(alpha+two)*fint1 - (apb+four)*fint2) * (apb+three)/four
  if (n == 2) then
   endw2 = f2
   return
  endif
  do i=3,n
   di   = dble(i-1)
   abn  = alpha+beta+di
   abnn = abn+di
   a1   =  -(two*(di+alpha)*(di+beta))/(abn*abnn*(abnn+one))
   a2   =  (two*(alpha-beta))/(abnn*(abnn+two))
   a3   =  (two*(abn+one))/((abnn+two)*(abnn+one))
   f3   =  -(a2*f2+a1*f1)/a3
   f1   = f2
   f2   = f3
  enddo
  endw2  = f3
  return
  end function endw2

  double precision function gammaf (x)
!
!=======================================================================
!
!     G a m m a f :
!     -----------
!
!=======================================================================
!

  double precision x
  double precision, parameter :: pi = 3.141592653589793d0
  double precision, parameter :: half=0.5d0,one=1.d0,two=2.d0

  gammaf = one

  if (x == -half) gammaf = -two*sqrt(pi)
  if (x ==  half) gammaf =  sqrt(pi)
  if (x ==  one ) gammaf =  one
  if (x ==  two ) gammaf =  one
  if (x ==  1.5d0) gammaf =  sqrt(pi)/2.d0
  if (x ==  2.5d0) gammaf =  1.5d0*sqrt(pi)/2.d0
  if (x ==  3.5d0) gammaf =  2.5d0*1.5d0*sqrt(pi)/2.d0
  if (x ==  3.d0 ) gammaf =  2.d0
  if (x ==  4.d0 ) gammaf = 6.d0
  if (x ==  5.d0 ) gammaf = 24.d0
  if (x ==  6.d0 ) gammaf = 120.d0

  return
  end function gammaf

  double precision FUNCTION HDGLL (I,j,ZGLL,NZ)
!-------------------------------------------------------------
!
!     Compute the value of the derivative of the I-th
!     Lagrangian interpolant through the
!     NZ Gauss-Lobatto Legendre points ZGLL at point ZGLL(j).
!
!-------------------------------------------------------------

  

  integer i,j,nz
  double precision zgll(0:nz-1)

  integer idegpoly
  double precision rlegendre1,rlegendre2,rlegendre3

    idegpoly = nz - 1
  if ((i == 0).and.(j == 0)) then
          hdgll = - dble(idegpoly)*(dble(idegpoly)+1.d0)/4.d0
  else if ((i == idegpoly).and.(j == idegpoly)) then
          hdgll = dble(idegpoly)*(dble(idegpoly)+1.d0)/4.d0
  else if (i == j) then
          hdgll = 0.d0
  else
         rlegendre1 = pnleg(zgll(j),idegpoly)
         rlegendre2 = pndleg(zgll(j),idegpoly)
         rlegendre3 = pnleg(zgll(i),idegpoly)
  hdgll = rlegendre1 / (rlegendre3*(zgll(j)-zgll(i))) &
    + (1.d0-zgll(j)*zgll(j))*rlegendre2/(dble(idegpoly)* &
    (dble(idegpoly)+1.d0)*rlegendre3* &
    (zgll(j)-zgll(i))*(zgll(j)-zgll(i)))
  endif

  return
  end FUNCTION hdgll

!=====================================================================

  double precision FUNCTION HGLL (I,Z,ZGLL,NZ)
!-------------------------------------------------------------
!
!     Compute the value of the Lagrangian interpolant L through
!     the NZ Gauss-Lobatto Legendre points ZGLL at the point Z.
!
!-------------------------------------------------------------

  

  integer i,nz
  double precision z
  double precision ZGLL(0:nz-1)

  integer n
  double precision EPS,DZ,ALFAN

  EPS = 1.d-5
  DZ = Z - ZGLL(I)
  IF (ABS(DZ)  <  EPS) THEN
   HGLL = 1.d0
   RETURN
  ENDIF
  N = NZ - 1
  ALFAN = dble(N)*(dble(N)+1.d0)
  HGLL = - (1.d0-Z*Z)*PNDLEG(Z,N)/ (ALFAN*PNLEG(ZGLL(I),N)*(Z-ZGLL(I)))
  RETURN
  end function hgll
!
!=====================================================================

  subroutine jacg (xjac,np,alpha,beta)
!
!=======================================================================
!
!     J a c g : Compute np Gauss points, which are the zeros of the
!               Jacobi polynomial with parameter alpha and beta.
!
!=======================================================================
!
!     Note :
!     ----
!          .Alpha and Beta determines the specific type of gauss points.
!                  .alpha = beta =  0.0  ->  Legendre points
!                  .alpha = beta = -0.5  ->  Chebyshev points
!
!=======================================================================
!
  

  integer np
  double precision alpha,beta
  double precision xjac(np)

  integer k,j,i,jmin,jm,n
  double precision xlast,dth,x,x1,x2,recsum,delx,xmin,swap
  double precision p,pd,pm1,pdm1,pm2,pdm2

  integer, parameter :: kstop = 10
  double precision, parameter :: zero = 0.d0, eps = 1.0d-12

!
!-----------------------------------------------------------------------
!
  pm1 = zero
  pm2 = zero
  pdm1 = zero
  pdm2 = zero

  xlast = 0.d0
  n   = np-1
  dth = 4.d0*datan(1.d0)/(2.d0*dble(n)+2.d0)
  p = 0.d0
  pd = 0.d0
  jmin = 0
  do 40 j=1,np
   if (j == 1) then
      x = cos((2.d0*(dble(j)-1.d0)+1.d0)*dth)