[cig-commits] r8614 - seismo/1D/SPECFEM1D/trunk
walter at geodynamics.org
walter at geodynamics.org
Fri Dec 7 16:24:13 PST 2007
Author: walter
Date: 2007-12-07 16:24:12 -0800 (Fri, 07 Dec 2007)
New Revision: 8614
Added:
seismo/1D/SPECFEM1D/trunk/Makefile
seismo/1D/SPECFEM1D/trunk/constants.h
seismo/1D/SPECFEM1D/trunk/define_derivative_matrix.f90
seismo/1D/SPECFEM1D/trunk/diffusion.f90
seismo/1D/SPECFEM1D/trunk/gll_library.f90
seismo/1D/SPECFEM1D/trunk/jeroen_teaching_software.txt
seismo/1D/SPECFEM1D/trunk/lagrange_poly.f90
seismo/1D/SPECFEM1D/trunk/source_time_function.f90
seismo/1D/SPECFEM1D/trunk/wave.f90
Log:
Initial revision
Added: seismo/1D/SPECFEM1D/trunk/Makefile
===================================================================
--- seismo/1D/SPECFEM1D/trunk/Makefile 2007-12-08 00:21:23 UTC (rev 8613)
+++ seismo/1D/SPECFEM1D/trunk/Makefile 2007-12-08 00:24:12 UTC (rev 8614)
@@ -0,0 +1,63 @@
+O = obj
+
+# Portland compiler
+#F90 = pgf90
+#FLAGS = -fast -Mbounds -Mneginfo -Mdclchk -Mstandard -Knoieee
+#FLAGS = -fast -Mnobounds -Mneginfo -Mdclchk -Munroll=c:6 -Mstandard -Knoieee
+
+# Absoft compiler
+#F90 = f90
+#FLAGS = -O2 -W132 -YEXT_NAMES=LCS -s -B108 -YCFRL=1
+
+# Intel Linux compiler
+F90 = ifort
+FLAGS = -O3 -e95 -implicitnone
+
+wave: constants.h \
+ $O/wave.o \
+ $O/gll_library.o \
+ $O/lagrange_poly.o \
+ $O/source_time_function.o \
+ $O/define_derivative_matrix.o
+ ${F90} $(FLAGS) -o xwave \
+ $O/wave.o \
+ $O/gll_library.o \
+ $O/lagrange_poly.o \
+ $O/source_time_function.o \
+ $O/define_derivative_matrix.o
+
+diffusion: constants.h \
+ $O/diffusion.o \
+ $O/gll_library.o \
+ $O/lagrange_poly.o \
+ $O/define_derivative_matrix.o
+ ${F90} $(FLAGS) -o xdiffusion \
+ $O/diffusion.o \
+ $O/gll_library.o \
+ $O/lagrange_poly.o \
+ $O/define_derivative_matrix.o
+
+bak:
+ cp *f90 *h Makefile bak
+
+clean:
+ rm -f $O/*.o *.o snapshot* seismogram xwave xdiffusion
+
+$O/wave.o: constants.h wave.f90
+ ${F90} $(FLAGS) -c -o $O/wave.o wave.f90
+
+$O/diffusion.o: constants.h diffusion.f90
+ ${F90} $(FLAGS) -c -o $O/diffusion.o diffusion.f90
+
+$O/define_derivative_matrix.o: constants.h define_derivative_matrix.f90
+ ${F90} $(FLAGS) -c -o $O/define_derivative_matrix.o define_derivative_matrix.f90
+
+$O/gll_library.o: gll_library.f90
+ ${F90} $(FLAGS) -c -o $O/gll_library.o gll_library.f90
+
+$O/lagrange_poly.o: lagrange_poly.f90
+ ${F90} $(FLAGS) -c -o $O/lagrange_poly.o lagrange_poly.f90
+
+$O/source_time_function.o: source_time_function.f90
+ ${F90} $(FLAGS) -c -o $O/source_time_function.o source_time_function.f90
+
Added: seismo/1D/SPECFEM1D/trunk/constants.h
===================================================================
--- seismo/1D/SPECFEM1D/trunk/constants.h 2007-12-08 00:21:23 UTC (rev 8613)
+++ seismo/1D/SPECFEM1D/trunk/constants.h 2007-12-08 00:24:12 UTC (rev 8614)
@@ -0,0 +1,11 @@
+! number of spectral elements
+ integer, parameter :: NSPEC = 11
+
+! number of GLL points (polynomial degree plus one)
+ integer, parameter :: NGLL = 7
+
+! number of global points
+ integer, parameter :: NGLOB = (NGLL-1)*NSPEC + 1
+
+! for the Gauss-Lobatto-Legendre points and weights
+ double precision, parameter :: GAUSSALPHA = 0.d0,GAUSSBETA = 0.d0
Added: seismo/1D/SPECFEM1D/trunk/define_derivative_matrix.f90
===================================================================
--- seismo/1D/SPECFEM1D/trunk/define_derivative_matrix.f90 2007-12-08 00:21:23 UTC (rev 8613)
+++ seismo/1D/SPECFEM1D/trunk/define_derivative_matrix.f90 2007-12-08 00:24:12 UTC (rev 8614)
@@ -0,0 +1,39 @@
+ subroutine define_derivative_matrix(xigll,wgll,hprime)
+
+ implicit none
+
+ include "constants.h"
+
+! Gauss-Lobatto-Legendre points of integration
+ double precision, dimension(NGLL) :: xigll
+
+! weights
+ double precision, dimension(NGLL) :: wgll
+
+! array with derivatives of Lagrange polynomials
+ double precision, dimension(NGLL,NGLL) :: hprime
+
+! function for calculating derivatives of Lagrange polynomials
+ double precision, external :: lagrange_deriv_GLL
+
+ integer i1,i2
+
+! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
+
+! set up coordinates of the Gauss-Lobatto-Legendre points
+ call zwgljd(xigll,wgll,NGLL,GAUSSALPHA,GAUSSBETA)
+
+! if number of points is odd, the middle abscissa is exactly zero
+ if(mod(NGLL,2) /= 0) xigll((NGLL-1)/2+1) = 0.d0
+
+! calculate derivatives of the Lagrange polynomials
+! and precalculate some products in double precision
+! hprime(i,j) = h'_i(xigll_j) by definition of the derivative matrix
+ do i1=1,NGLL
+ do i2=1,NGLL
+ hprime(i1,i2) = lagrange_deriv_GLL(i1-1,i2-1,xigll,NGLL)
+ enddo
+ enddo
+
+ end subroutine define_derivative_matrix
+
Added: seismo/1D/SPECFEM1D/trunk/diffusion.f90
===================================================================
--- seismo/1D/SPECFEM1D/trunk/diffusion.f90 2007-12-08 00:21:23 UTC (rev 8613)
+++ seismo/1D/SPECFEM1D/trunk/diffusion.f90 2007-12-08 00:24:12 UTC (rev 8614)
@@ -0,0 +1,227 @@
+ program diffusion
+
+ implicit none
+
+ include "constants.h"
+
+! number of timesteps
+ integer, parameter :: NSTEP = 30000
+
+! time step in seconds
+ double precision, parameter :: DT = 100000000. ! s
+
+! fixed boundary conditions
+ logical, parameter :: FIXED_BC = .true.
+
+! model parameters (SI)
+ double precision, parameter :: LENGTH = 3.0d+03 ! m
+ double precision, parameter :: DENSITY = 2.5d+03 ! kg/m^3
+ double precision, parameter :: THERMALCONDUCTIVITY = 10.0d-01 ! cal/m/s/K
+ double precision, parameter :: HEATCAPACITY = 0.3d+03 ! cal/kg/K
+
+ integer ispec,i,j,iglob,itime
+
+! Gauss-Lobatto-Legendre points of integration
+ double precision, dimension(NGLL) :: xigll
+
+! weights
+ double precision, dimension(NGLL) :: wgll
+
+! array with derivatives of Lagrange polynomials
+ double precision, dimension(NGLL,NGLL) :: hprime
+
+! anchors
+ double precision, dimension(NSPEC) :: x1,x2
+
+! global grid points
+ double precision, dimension(NGLOB) :: x
+
+! material properties
+ double precision, dimension(NGLL,NSPEC) :: rho,heat_capacity,thermal_conductivity
+
+! Jacobian `matrix' and Jacobian
+ double precision, dimension(NGLL,NSPEC) :: dxidx,jacobian
+
+! local mass matrix
+ double precision mass_local
+
+! global mass matrix
+ double precision, dimension(NGLOB) :: mass_global
+
+! temperature and temperature gradient
+ double precision, dimension(NGLOB) :: temperature,grad_temperature
+
+! local to global numbering
+ integer, dimension(NGLL,NSPEC) :: ibool
+
+! time marching
+ double precision deltat,deltatover2,deltatsqover2
+ double precision dh,diffusivity,time_step
+
+! end fluxes
+ double precision flux_1,flux_NGLOB
+
+! end temperatures
+ double precision temperature_1,temperature_NGLOB
+
+! derivatives
+ double precision dtdx,flux,templ,temp(NGLL)
+
+! movie
+ character(len=50) moviefile
+
+!++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ call define_derivative_matrix(xigll,wgll,hprime)
+
+! evenly spaced achors between 0 and 1
+ do ispec = 1,NSPEC
+ x1(ispec) = LENGTH*dble(ispec-1)/dble(NSPEC)
+ x2(ispec) = LENGTH*dble(ispec)/dble(NSPEC)
+ enddo
+
+! set up the mesh properties
+ do ispec = 1,NSPEC
+ do i = 1,NGLL
+ rho(i,ispec) = DENSITY
+ thermal_conductivity(i,ispec) = THERMALCONDUCTIVITY
+ heat_capacity(i,ispec) = HEATCAPACITY
+ dxidx(i,ispec) = 2. / (x2(ispec)-x1(ispec))
+ jacobian(i,ispec) = (x2(ispec)-x1(ispec)) / 2.
+ enddo
+ enddo
+
+! set up local to global numbering
+ iglob = 1
+ do ispec = 1,NSPEC
+ do i = 1,NGLL
+ if(i > 1) iglob = iglob+1
+ ibool(i,ispec) = iglob
+ enddo
+ enddo
+
+! get the global grid points
+ do ispec = 1,NSPEC
+ do i = 1,NGLL
+ iglob = ibool(i,ispec)
+ x(iglob) = 0.5*(1.-xigll(i))*x1(ispec)+0.5*(1.+xigll(i))*x2(ispec)
+ enddo
+ enddo
+
+! calculate the global mass matrix
+ mass_global(:) = 0.
+ do ispec = 1,NSPEC
+ do i = 1,NGLL
+ mass_local = wgll(i)*rho(i,ispec)*heat_capacity(i,ispec)*jacobian(i,ispec)
+ iglob = ibool(i,ispec)
+ mass_global(iglob) = mass_global(iglob) + mass_local
+ enddo
+ enddo
+
+! estimate the time step
+ dh = LENGTH/dble(NGLOB-1)
+ diffusivity = THERMALCONDUCTIVITY/(HEATCAPACITY*DENSITY)
+ time_step = 0.5*dh*dh/diffusivity
+! print *,'time step estimate: ',time_step,' seconds'
+
+ if(FIXED_BC) then
+! set up the temperatures at the ends
+ temperature_1 = 10.
+ temperature_NGLOB = 0.
+ else
+! set up the fluxes at the ends
+ flux_1 = 0.
+ flux_NGLOB = 0.
+ endif
+
+! time marching parameters
+ deltat = DT
+ deltatover2 = deltat/2.
+ deltatsqover2 = deltat*deltat/2.
+
+! initialize
+ temperature(:) = 0.
+ grad_temperature(:) = 0.
+
+ do itime = 1,NSTEP
+
+! update temperature
+ temperature(:) = temperature(:) + deltatover2*grad_temperature(:)
+ grad_temperature(:) = 0.
+
+ do ispec = 1,NSPEC
+
+ do i = 1,NGLL
+! get dtdx
+ templ = 0.
+ do j = 1,NGLL
+ iglob = ibool(j,ispec)
+ templ = templ + temperature(iglob)*hprime(j,i)
+ enddo
+ dtdx = templ*dxidx(i,ispec)
+
+! heat flux
+ flux = -thermal_conductivity(i,ispec)*dtdx
+
+ temp(i) = jacobian(i,ispec)*flux*dxidx(i,ispec)
+ enddo ! first loop over the GLL points
+
+ do i = 1,NGLL
+ templ = 0.
+ do j = 1,NGLL
+ templ = templ + temp(j)*hprime(i,j)*wgll(j)
+ enddo
+
+! update the gradient of temperature
+ iglob = ibool(i,ispec)
+ grad_temperature(iglob) = grad_temperature(iglob) + templ
+ enddo ! second loop over the GLL points
+
+ enddo ! end loop over all spectral elements
+
+! fixed boundary conditions at global level
+ if(FIXED_BC) then
+! left side
+ temperature(1) = temperature_1
+ templ = 0.
+ do i = 1,NGLL
+ iglob = ibool(i,1)
+ templ = templ + temperature(iglob)*hprime(i,1)
+ enddo
+ dtdx = templ*dxidx(1,1)
+ flux_1 = -thermal_conductivity(1,1)*dtdx
+! right side
+ temperature(NGLOB) = temperature_NGLOB
+ templ = 0.
+ do i = 1,NGLL
+ iglob = ibool(i,NSPEC)
+ templ = templ + temperature(iglob)*hprime(i,NGLL)
+ enddo
+ dtdx = templ*dxidx(NGLL,NSPEC)
+ flux_NGLOB = -thermal_conductivity(NGLL,NSPEC)*dtdx
+ endif
+
+! add in the end fluxes
+ grad_temperature(1) = grad_temperature(1) + flux_1
+ grad_temperature(NGLOB) = grad_temperature(NGLOB) - flux_NGLOB
+
+! divide by the mass matrix
+ grad_temperature(:) = grad_temperature(:)/mass_global(:)
+
+! update temperature
+ temperature(:) = temperature(:) + deltatover2*grad_temperature(:)
+
+! write out snapshots
+ if(mod(itime-1,1000) == 0) then
+ write(moviefile,10) itime
+10 format('snapshot',i5.5)
+ open(unit=10,file=moviefile,status='unknown')
+ do iglob = 1,NGLOB
+ write(10,*) sngl(x(iglob)),sngl(temperature(iglob))
+ enddo
+ close(10)
+ endif
+
+ enddo ! end time loop
+
+ end program diffusion
Added: seismo/1D/SPECFEM1D/trunk/gll_library.f90
===================================================================
--- seismo/1D/SPECFEM1D/trunk/gll_library.f90 2007-12-08 00:21:23 UTC (rev 8613)
+++ seismo/1D/SPECFEM1D/trunk/gll_library.f90 2007-12-08 00:24:12 UTC (rev 8614)
@@ -0,0 +1,528 @@
+!=======================================================================
+!
+! Library to compute the Gauss-Lobatto-Legendre points and weights
+! Based on Gauss-Lobatto routines from M.I.T.
+! Department of Mechanical Engineering
+!
+!=======================================================================
+
+ double precision function endw1(n,alpha,beta)
+
+ implicit none
+
+ 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
+ double precision, external :: gammaf
+ 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
+
+ end function endw1
+
+!
+!=======================================================================
+!
+
+ double precision function endw2(n,alpha,beta)
+
+ implicit none
+
+ 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
+ double precision, external :: gammaf
+ 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
+
+ end function endw2
+
+!
+!=======================================================================
+!
+
+ double precision function gammaf (x)
+
+ implicit none
+
+ double precision, parameter :: pi = 3.141592653589793d0
+
+ double precision x
+
+ double precision, parameter :: half=0.5d0,one=1.d0,two=2.d0
+
+ gammaf = one
+
+ if (x == -half) gammaf = -two*dsqrt(pi)
+ if (x == half) gammaf = dsqrt(pi)
+ if (x == one ) gammaf = one
+ if (x == two ) gammaf = one
+ if (x == 1.5d0) gammaf = dsqrt(pi)/2.d0
+ if (x == 2.5d0) gammaf = 1.5d0*dsqrt(pi)/2.d0
+ if (x == 3.5d0) gammaf = 2.5d0*1.5d0*dsqrt(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
+
+ end function gammaf
+
+!
+!=====================================================================
+!
+
+ subroutine jacg (xjac,np,alpha,beta)
+
+!=======================================================================
+!
+! computes np Gauss points, which are the zeros of the
+! Jacobi polynomial with parameters alpha and beta
+!
+! .alpha = beta = 0.0 -> Legendre points
+! .alpha = beta = -0.5 -> Chebyshev points
+!
+!=======================================================================
+
+ implicit none
+
+ 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 :: K_MAX_ITER = 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 j=1,np
+ if(j == 1) then
+ x = dcos((2.d0*(dble(j)-1.d0)+1.d0)*dth)
+ else
+ x1 = dcos((2.d0*(dble(j)-1.d0)+1.d0)*dth)
+ x2 = xlast
+ x = (x1+x2)/2.d0
+ endif
+ do k=1,K_MAX_ITER
+ call jacobf (p,pd,pm1,pdm1,pm2,pdm2,np,alpha,beta,x)
+ recsum = 0.d0
+ jm = j-1
+ do i=1,jm
+ recsum = recsum+1.d0/(x-xjac(np-i+1))
+ enddo
+ delx = -p/(pd-recsum*p)
+ x = x+delx
+ if(abs(delx) < eps) goto 31
+ enddo
+ 31 continue
+ xjac(np-j+1) = x
+ xlast = x
+ enddo
+ do i=1,np
+ xmin = 2.d0
+ do j=i,np
+ if(xjac(j) < xmin) then
+ xmin = xjac(j)
+ jmin = j
+ endif
+ enddo
+ if(jmin /= i) then
+ swap = xjac(i)
+ xjac(i) = xjac(jmin)
+ xjac(jmin) = swap
+ endif
+ enddo
+
+ end subroutine jacg
+
+!
+!=====================================================================
+!
+
+ subroutine jacobf (poly,pder,polym1,pderm1,polym2,pderm2,n,alp,bet,x)
+
+!=======================================================================
+!
+! Computes the Jacobi polynomial of degree n and its derivative at x
+!
+!=======================================================================
+
+ implicit none
+
+ double precision poly,pder,polym1,pderm1,polym2,pderm2,alp,bet,x
+ integer n
+
+ double precision apb,polyl,pderl,dk,a1,a2,b3,a3,a4,polyn,pdern,psave,pdsave
+ integer k
+
+ apb = alp+bet
+ poly = 1.d0
+ pder = 0.d0
+ psave = 0.d0
+ pdsave = 0.d0
+
+ if (n == 0) return
+
+ polyl = poly
+ pderl = pder
+ poly = (alp-bet+(apb+2.d0)*x)/2.d0
+ pder = (apb+2.d0)/2.d0
+ if (n == 1) return
+
+ do k=2,n
+ dk = dble(k)
+ a1 = 2.d0*dk*(dk+apb)*(2.d0*dk+apb-2.d0)
+ a2 = (2.d0*dk+apb-1.d0)*(alp**2-bet**2)
+ b3 = (2.d0*dk+apb-2.d0)
+ a3 = b3*(b3+1.d0)*(b3+2.d0)
+ a4 = 2.d0*(dk+alp-1.d0)*(dk+bet-1.d0)*(2.d0*dk+apb)
+ polyn = ((a2+a3*x)*poly-a4*polyl)/a1
+ pdern = ((a2+a3*x)*pder-a4*pderl+a3*poly)/a1
+ psave = polyl
+ pdsave = pderl
+ polyl = poly
+ poly = polyn
+ pderl = pder
+ pder = pdern
+ enddo
+
+ polym1 = polyl
+ pderm1 = pderl
+ polym2 = psave
+ pderm2 = pdsave
+
+ end subroutine jacobf
+
+!
+!------------------------------------------------------------------------
+!
+
+ double precision FUNCTION PNDLEG (Z,N)
+
+!------------------------------------------------------------------------
+!
+! Compute the derivative of the Nth order Legendre polynomial at Z.
+! Based on the recursion formula for the Legendre polynomials.
+!
+!------------------------------------------------------------------------
+ implicit none
+
+ double precision z
+ integer n
+
+ double precision P1,P2,P1D,P2D,P3D,FK,P3
+ integer k
+
+ P1 = 1.d0
+ P2 = Z
+ P1D = 0.d0
+ P2D = 1.d0
+ P3D = 1.d0
+
+ do K = 1, N-1
+ FK = dble(K)
+ P3 = ((2.d0*FK+1.d0)*Z*P2 - FK*P1)/(FK+1.d0)
+ P3D = ((2.d0*FK+1.d0)*P2 + (2.d0*FK+1.d0)*Z*P2D - FK*P1D) / (FK+1.d0)
+ P1 = P2
+ P2 = P3
+ P1D = P2D
+ P2D = P3D
+ enddo
+
+ PNDLEG = P3D
+
+ end function pndleg
+
+!
+!------------------------------------------------------------------------
+!
+
+ double precision FUNCTION PNLEG (Z,N)
+
+!------------------------------------------------------------------------
+!
+! Compute the value of the Nth order Legendre polynomial at Z.
+! Based on the recursion formula for the Legendre polynomials.
+!
+!------------------------------------------------------------------------
+ implicit none
+
+ double precision z
+ integer n
+
+ double precision P1,P2,P3,FK
+ integer k
+
+ P1 = 1.d0
+ P2 = Z
+ P3 = P2
+
+ do K = 1, N-1
+ FK = dble(K)
+ P3 = ((2.d0*FK+1.d0)*Z*P2 - FK*P1)/(FK+1.d0)
+ P1 = P2
+ P2 = P3
+ enddo
+
+ PNLEG = P3
+
+ end function pnleg
+
+!
+!------------------------------------------------------------------------
+!
+
+ double precision function pnormj (n,alpha,beta)
+
+ implicit none
+
+ double precision alpha,beta
+ integer n
+
+ double precision one,two,dn,const,prod,dindx,frac
+ double precision, external :: gammaf
+ integer i
+
+ one = 1.d0
+ two = 2.d0
+ dn = dble(n)
+ const = alpha+beta+one
+
+ if (n <= 1) then
+ prod = gammaf(dn+alpha)*gammaf(dn+beta)
+ prod = prod/(gammaf(dn)*gammaf(dn+alpha+beta))
+ pnormj = prod * two**const/(two*dn+const)
+ return
+ endif
+
+ prod = gammaf(alpha+one)*gammaf(beta+one)
+ prod = prod/(two*(one+const)*gammaf(const+one))
+ prod = prod*(one+alpha)*(two+alpha)
+ prod = prod*(one+beta)*(two+beta)
+
+ do i=3,n
+ dindx = dble(i)
+ frac = (dindx+alpha)*(dindx+beta)/(dindx*(dindx+alpha+beta))
+ prod = prod*frac
+ enddo
+
+ pnormj = prod * two**const/(two*dn+const)
+
+ end function pnormj
+
+!
+!------------------------------------------------------------------------
+!
+
+ subroutine zwgjd(z,w,np,alpha,beta)
+
+!=======================================================================
+!
+! Z w g j d : Generate np Gauss-Jacobi points and weights
+! associated with Jacobi polynomial of degree n = np-1
+!
+! Note : Coefficients alpha and beta must be greater than -1.
+! ----
+!=======================================================================
+
+ implicit none
+
+ double precision, parameter :: zero=0.d0,one=1.d0,two=2.d0
+
+ integer np
+ double precision z(np),w(np)
+ double precision alpha,beta
+
+ integer n,np1,np2,i
+ double precision p,pd,pm1,pdm1,pm2,pdm2
+ double precision apb,dnp1,dnp2,fac1,fac2,fac3,fnorm,rcoef
+ double precision, external :: gammaf,pnormj
+
+ pd = zero
+ pm1 = zero
+ pm2 = zero
+ pdm1 = zero
+ pdm2 = zero
+
+ n = np-1
+ apb = alpha+beta
+ p = zero
+ pdm1 = zero
+
+ if (np <= 0) stop 'minimum number of Gauss points is 1'
+
+ if ((alpha <= -one) .or. (beta <= -one)) stop 'alpha and beta must be greater than -1'
+
+ if (np == 1) then
+ z(1) = (beta-alpha)/(apb+two)
+ w(1) = gammaf(alpha+one)*gammaf(beta+one)/gammaf(apb+two) * two**(apb+one)
+ return
+ endif
+
+ call jacg(z,np,alpha,beta)
+
+ np1 = n+1
+ np2 = n+2
+ dnp1 = dble(np1)
+ dnp2 = dble(np2)
+ fac1 = dnp1+alpha+beta+one
+ fac2 = fac1+dnp1
+ fac3 = fac2+one
+ fnorm = pnormj(np1,alpha,beta)
+ rcoef = (fnorm*fac2*fac3)/(two*fac1*dnp2)
+ do i=1,np
+ call jacobf(p,pd,pm1,pdm1,pm2,pdm2,np2,alpha,beta,z(i))
+ w(i) = -rcoef/(p*pdm1)
+ enddo
+
+ end subroutine zwgjd
+
+!
+!------------------------------------------------------------------------
+!
+
+ subroutine zwgljd(z,w,np,alpha,beta)
+
+!=======================================================================
+!
+! Z w g l j d : Generate np Gauss-Lobatto-Jacobi points and the
+! ----------- weights associated with Jacobi polynomials of degree
+! n = np-1.
+!
+! Note : alpha and beta coefficients must be greater than -1.
+! Legendre polynomials are special case of Jacobi polynomials
+! just by setting alpha and beta to 0.
+!
+!=======================================================================
+
+ implicit none
+
+ double precision, parameter :: zero=0.d0,one=1.d0,two=2.d0
+
+ integer np
+ double precision alpha,beta
+ double precision z(np), w(np)
+
+ integer n,nm1,i
+ double precision p,pd,pm1,pdm1,pm2,pdm2
+ double precision alpg,betg
+ double precision, external :: endw1,endw2
+
+ p = zero
+ pm1 = zero
+ pm2 = zero
+ pdm1 = zero
+ pdm2 = zero
+
+ n = np-1
+ nm1 = n-1
+ pd = zero
+
+ if (np <= 1) stop 'minimum number of Gauss-Lobatto points is 2'
+
+! with spectral elements, use at least 3 points
+ if (np <= 2) stop 'minimum number of Gauss-Lobatto points for the SEM is 3'
+
+ if ((alpha <= -one) .or. (beta <= -one)) stop 'alpha and beta must be greater than -1'
+
+ if (nm1 > 0) then
+ alpg = alpha+one
+ betg = beta+one
+ call zwgjd(z(2),w(2),nm1,alpg,betg)
+ endif
+
+ z(1) = - one
+ z(np) = one
+
+ do i=2,np-1
+ w(i) = w(i)/(one-z(i)**2)
+ enddo
+
+ call jacobf(p,pd,pm1,pdm1,pm2,pdm2,n,alpha,beta,z(1))
+ w(1) = endw1(n,alpha,beta)/(two*pd)
+ call jacobf(p,pd,pm1,pdm1,pm2,pdm2,n,alpha,beta,z(np))
+ w(np) = endw2(n,alpha,beta)/(two*pd)
+
+ end subroutine zwgljd
+
Added: seismo/1D/SPECFEM1D/trunk/jeroen_teaching_software.txt
===================================================================
--- seismo/1D/SPECFEM1D/trunk/jeroen_teaching_software.txt 2007-12-08 00:21:23 UTC (rev 8613)
+++ seismo/1D/SPECFEM1D/trunk/jeroen_teaching_software.txt 2007-12-08 00:24:12 UTC (rev 8614)
@@ -0,0 +1,24 @@
+
+From: "Jeroen Tromp"
+To: komatitsch
+Subject: course software
+Date: Wed, 8 Jan 2003 14:04:45 -0800
+
+Hi Dimitri:
+
+I am teaching part of a course in computational geophysics this
+quarter. For
+this reason I put together the basic software for students to do 1D
+simulations of the wave equation and the diffusion equation based upon
+the
+spectral-element method. Attached the package. What I have in mind is to
+give
+them the basic subroutines lagrange_poly.f90, gll_library.f90, and
+define_derivative_matrix.f90, and to ask them to write the small
+programs
+wave.f90 and diffusion.f90 based upon class notes. This way they can
+get a
+pretty good feel for how the method works.
+
+Jeroen
+
Added: seismo/1D/SPECFEM1D/trunk/lagrange_poly.f90
===================================================================
--- seismo/1D/SPECFEM1D/trunk/lagrange_poly.f90 2007-12-08 00:21:23 UTC (rev 8613)
+++ seismo/1D/SPECFEM1D/trunk/lagrange_poly.f90 2007-12-08 00:24:12 UTC (rev 8614)
@@ -0,0 +1,83 @@
+ subroutine lagrange_any(xi,NGLL,xigll,h,hprime)
+
+! subroutine to compute the Lagrange interpolants based upon the GLL points
+! and their first derivatives at any point xi in [-1,1]
+
+ implicit none
+
+ integer NGLL
+ double precision xi,xigll(NGLL),h(NGLL),hprime(NGLL)
+
+ integer dgr,i,j
+ double precision prod1,prod2
+
+ do dgr=1,NGLL
+
+ prod1 = 1.0d0
+ prod2 = 1.0d0
+ do i=1,NGLL
+ if(i /= dgr) then
+ prod1 = prod1*(xi-xigll(i))
+ prod2 = prod2*(xigll(dgr)-xigll(i))
+ endif
+ enddo
+ h(dgr)=prod1/prod2
+
+ hprime(dgr)=0.0d0
+ do i=1,NGLL
+ if(i /= dgr) then
+ prod1=1.0d0
+ do j=1,NGLL
+ if(j /= dgr .and. j /= i) prod1 = prod1*(xi-xigll(j))
+ enddo
+ hprime(dgr) = hprime(dgr)+prod1
+ endif
+ enddo
+ hprime(dgr) = hprime(dgr)/prod2
+
+ enddo
+
+ end subroutine lagrange_any
+
+!
+!=====================================================================
+!
+
+! subroutine to compute the derivative of the Lagrange interpolants
+! at the GLL points at any given GLL point
+
+ double precision function lagrange_deriv_GLL(I,j,ZGLL,NZ)
+
+!------------------------------------------------------------------------
+!
+! Compute the value of the derivative of the I-th
+! Lagrange interpolant through the
+! NZ Gauss-Lobatto Legendre points ZGLL at point ZGLL(j)
+!
+!------------------------------------------------------------------------
+
+ implicit none
+
+ integer i,j,nz
+ double precision zgll(0:nz-1)
+
+ integer degpoly
+
+ double precision, external :: pnleg,pndleg
+
+ degpoly = nz - 1
+ if (i == 0 .and. j == 0) then
+ lagrange_deriv_GLL = - dble(degpoly)*(dble(degpoly)+1.d0) / 4.d0
+ else if (i == degpoly .and. j == degpoly) then
+ lagrange_deriv_GLL = dble(degpoly)*(dble(degpoly)+1.d0) / 4.d0
+ else if (i == j) then
+ lagrange_deriv_GLL = 0.d0
+ else
+ lagrange_deriv_GLL = pnleg(zgll(j),degpoly) / &
+ (pnleg(zgll(i),degpoly)*(zgll(j)-zgll(i))) &
+ + (1.d0-zgll(j)*zgll(j))*pndleg(zgll(j),degpoly) / (dble(degpoly)* &
+ (dble(degpoly)+1.d0)*pnleg(zgll(i),degpoly)*(zgll(j)-zgll(i))*(zgll(j)-zgll(i)))
+ endif
+
+ end function lagrange_deriv_GLL
+
Added: seismo/1D/SPECFEM1D/trunk/source_time_function.f90
===================================================================
--- seismo/1D/SPECFEM1D/trunk/source_time_function.f90 2007-12-08 00:21:23 UTC (rev 8613)
+++ seismo/1D/SPECFEM1D/trunk/source_time_function.f90 2007-12-08 00:24:12 UTC (rev 8614)
@@ -0,0 +1,21 @@
+ double precision function source_time_function(t,hdur)
+
+ implicit none
+
+! source decay rate (also change in source spectrum if needed)
+ double precision, parameter :: decay_rate = 2.628d0
+ double precision, parameter :: PI = 3.141592653589793d0
+
+ double precision t,hdur
+
+ double precision alpha
+
+ alpha = decay_rate/hdur
+
+! Gaussian
+ source_time_function = alpha*dexp(-alpha*alpha*t*t)/dsqrt(PI)
+
+! Ricker wavelet
+! source_time_function = -2.*(alpha**3)*t*dexp(-alpha*alpha*t*t)/dsqrt(PI)
+
+ end function source_time_function
Added: seismo/1D/SPECFEM1D/trunk/wave.f90
===================================================================
--- seismo/1D/SPECFEM1D/trunk/wave.f90 2007-12-08 00:21:23 UTC (rev 8613)
+++ seismo/1D/SPECFEM1D/trunk/wave.f90 2007-12-08 00:24:12 UTC (rev 8614)
@@ -0,0 +1,231 @@
+ program wave
+
+ implicit none
+
+ include "constants.h"
+
+! number of timesteps
+ integer, parameter :: NSTEP = 10000
+
+! time step in seconds
+ double precision, parameter :: DT = 0.0002 ! s
+
+! fixed boundary conditions
+ logical, parameter :: FIXED_BC = .true.
+
+! model parameters (SI)
+ double precision, parameter :: LENGTH = 3.0d+03 ! m
+ double precision, parameter :: DENSITY = 2.5d+03 ! kg/m^3
+ double precision, parameter :: RIGIDITY = 3.0d+10 ! Pa
+
+! pi
+ double precision, parameter :: PI = 3.141592653589793d0
+
+ integer ispec,i,j,iglob,itime
+
+! Gauss-Lobatto-Legendre points of integration
+ double precision, dimension(NGLL) :: xigll
+
+! weights
+ double precision, dimension(NGLL) :: wgll
+
+! array with derivatives of Lagrange polynomials
+ double precision, dimension(NGLL,NGLL) :: hprime
+
+! anchors
+ double precision, dimension(NSPEC) :: x1,x2
+
+! global grid points
+ double precision, dimension(NGLOB) :: x
+
+! material properties
+ double precision, dimension(NGLL,NSPEC) :: rho,mu
+
+! Jacobian `matrix' and Jacobian
+ double precision, dimension(NGLL,NSPEC) :: dxidx,jacobian
+
+! local mass matrix
+ double precision mass_local
+
+! global mass matrix
+ double precision, dimension(NGLOB) :: mass_global
+
+! displacement, velocity and acceleration
+ double precision, dimension(NGLOB) :: displ,veloc,accel
+
+! local to global numbering
+ integer, dimension(NGLL,NSPEC) :: ibool
+
+! time marching
+ double precision dh,v,courant,time_step
+ double precision deltat,deltatover2,deltatsqover2
+
+! source
+ integer ispec_source,i_source,iglob_source
+ double precision stf,hdur,source_amp
+ double precision, external :: source_time_function
+
+! receiver
+ integer ireceiver
+ double precision seismogram(NSTEP)
+
+! derivatives
+ double precision dudx,sigma,templ,temp(NGLL)
+
+! movie
+ character(len=50) moviefile
+
+!++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ call define_derivative_matrix(xigll,wgll,hprime)
+
+! evenly spaced achors between 0 and 1
+ do ispec = 1,NSPEC
+ x1(ispec) = LENGTH*dble(ispec-1)/dble(NSPEC)
+ x2(ispec) = LENGTH*dble(ispec)/dble(NSPEC)
+ enddo
+
+! set up the mesh properties
+ do ispec = 1,NSPEC
+ do i = 1,NGLL
+ rho(i,ispec) = DENSITY
+ mu(i,ispec) = RIGIDITY
+ dxidx(i,ispec) = 2. / (x2(ispec)-x1(ispec))
+ jacobian(i,ispec) = (x2(ispec)-x1(ispec)) / 2.
+ enddo
+ enddo
+
+! set up local to global numbering
+ iglob = 1
+ do ispec = 1,NSPEC
+ do i = 1,NGLL
+ if(i > 1) iglob = iglob+1
+ ibool(i,ispec) = iglob
+ enddo
+ enddo
+
+! get the global grid points
+ do ispec = 1,NSPEC
+ do i = 1,NGLL
+ iglob = ibool(i,ispec)
+ x(iglob) = 0.5*(1.-xigll(i))*x1(ispec)+0.5*(1.+xigll(i))*x2(ispec)
+ enddo
+ enddo
+
+! calculate the global mass matrix
+ mass_global(:) = 0.
+ do ispec = 1,NSPEC
+ do i = 1,NGLL
+ mass_local = wgll(i)*rho(i,ispec)*jacobian(i,ispec)
+ iglob = ibool(i,ispec)
+ mass_global(iglob) = mass_global(iglob) + mass_local
+ enddo
+ enddo
+
+! estimate the time step
+ dh = LENGTH/dble(NGLOB-1)
+ v = dsqrt(RIGIDITY/DENSITY)
+ courant = 0.2
+ time_step = courant*dh/v
+! print *,'time step estimate: ',time_step,' seconds'
+
+! set the source
+ ispec_source = (NSPEC+1)/2
+ i_source = (NGLL+1)/2
+ hdur = 50.*DT
+ source_amp = 0.001
+
+! set the receiver
+ ireceiver = 2*NGLL-2
+
+! time marching parameters
+ deltat = DT
+ deltatover2 = deltat/2.
+ deltatsqover2 = deltat*deltat/2.
+
+! initialize
+ displ(:) = 0.
+ veloc(:) = 0.
+ accel(:) = 0.
+
+ do iglob = 1,NGLOB
+ displ(iglob) = dsin(PI*x(iglob)/LENGTH)
+ enddo
+
+ do itime = 1,NSTEP
+
+! `predictor' update displacement using finite-difference time scheme (Newmark)
+ displ(:) = displ(:) + deltat*veloc(:) + deltatsqover2*accel(:)
+ veloc(:) = veloc(:) + deltatover2*accel(:)
+ accel(:) = 0.
+
+ do ispec = 1,NSPEC
+
+ do i = 1,NGLL
+! get dudx
+ templ = 0.
+ do j = 1,NGLL
+ iglob = ibool(j,ispec)
+ templ = templ + displ(iglob)*hprime(j,i)
+ enddo
+ dudx = templ*dxidx(i,ispec)
+
+! stress
+ sigma = mu(i,ispec)*dudx
+
+ temp(i) = jacobian(i,ispec)*sigma*dxidx(i,ispec)
+ enddo ! first loop over the GLL points
+
+ do i = 1,NGLL
+ templ = 0.
+ do j = 1,NGLL
+ templ = templ + temp(j)*hprime(i,j)*wgll(j)
+ enddo
+
+! `corrector' update acceleration
+ iglob = ibool(i,ispec)
+ accel(iglob) = accel(iglob) - templ
+ enddo ! second loop over the GLL points
+
+ enddo ! end loop over all spectral elements
+
+! add source at global level
+! iglob_source = ibool(ispec_source,i_source)
+! stf = source_time_function(dble(itime-1)*DT-hdur,hdur)
+! accel(iglob_source) = accel(iglob_source) + stf*source_amp
+
+! fixed boundary conditions at global level
+ if(FIXED_BC) then
+ accel(1) = 0.
+ accel(NGLOB) = 0.
+ endif
+
+! divide by the mass matrix
+ accel(:) = accel(:)/mass_global(:)
+
+! `corrector' update velocity
+ veloc(:) = veloc(:) + deltatover2*accel(:)
+
+! write out snapshots
+ if(mod(itime-1,1000) == 0) then
+ write(moviefile,10) itime
+10 format('snapshot',i5.5)
+ open(unit=10,file=moviefile,status='unknown')
+ do iglob = 1,NGLOB
+ write(10,*) sngl(x(iglob)),sngl(displ(iglob))
+ enddo
+ close(10)
+ endif
+
+ seismogram(itime) = displ(ireceiver)
+
+ enddo ! end time loop
+
+ open(unit=12,file='seismogram',status='unknown')
+ do itime = 1,NSTEP
+ write(12,*) sngl(dble(itime-1)*DT-hdur),sngl(seismogram(itime))
+ enddo
+ close(12)
+
+ end program wave
+
More information about the cig-commits
mailing list