[cig-commits] r8419 - in seismo/2D/SPECFEM2D/trunk: MAILLE90
SPECFEM90
walter at geodynamics.org
walter at geodynamics.org
Fri Dec 7 15:44:45 PST 2007
Author: walter
Date: 2007-12-07 15:44:44 -0800 (Fri, 07 Dec 2007)
New Revision: 8419
Added:
seismo/2D/SPECFEM2D/trunk/SPECFEM90/define_derivative_matrices.f90
Removed:
seismo/2D/SPECFEM2D/trunk/MAILLE90/Par_rouen_ok
seismo/2D/SPECFEM2D/trunk/SPECFEM90/qmasspec.f90
seismo/2D/SPECFEM2D/trunk/SPECFEM90/qsumspec.f90
Modified:
seismo/2D/SPECFEM2D/trunk/MAILLE90/Par
seismo/2D/SPECFEM2D/trunk/MAILLE90/meshfem2D.f90
seismo/2D/SPECFEM2D/trunk/SPECFEM90/Makefile
seismo/2D/SPECFEM2D/trunk/SPECFEM90/constants.h
seismo/2D/SPECFEM2D/trunk/SPECFEM90/createnum_fast.f90
seismo/2D/SPECFEM2D/trunk/SPECFEM90/createnum_slow.f90
seismo/2D/SPECFEM2D/trunk/SPECFEM90/defarrays.f90
seismo/2D/SPECFEM2D/trunk/SPECFEM90/plotgll.f90
seismo/2D/SPECFEM2D/trunk/SPECFEM90/plotpost.f90
seismo/2D/SPECFEM2D/trunk/SPECFEM90/positsource.f90
seismo/2D/SPECFEM2D/trunk/SPECFEM90/q49shape.f90
seismo/2D/SPECFEM2D/trunk/SPECFEM90/q49spec.f90
seismo/2D/SPECFEM2D/trunk/SPECFEM90/specfem2D.f90
Log:
further cleaned version 5.0 of mesher and solver, added weak formulation based
explicitly on stress tensor, and added full anisotropy
Modified: seismo/2D/SPECFEM2D/trunk/MAILLE90/Par
===================================================================
--- seismo/2D/SPECFEM2D/trunk/MAILLE90/Par 2004-05-16 01:48:11 UTC (rev 8418)
+++ seismo/2D/SPECFEM2D/trunk/MAILLE90/Par 2007-12-07 23:44:44 UTC (rev 8419)
@@ -65,7 +65,7 @@
imeshvect =.true. ! display mesh on vector plots or not
imodelvect =.false. ! display velocity model on vector plots or not
iboundvect =.true. ! display boundary conditions on vector plots
-interpol =.false. ! interpolation of the display or not
+interpol =.true. ! interpolation of the display or not
iptsdisp =6 ! nb of points for interpolation of the display
isubsamp =2 ! subsampling of color snapshots
ignuplot =.false. ! generate a GNUPLOT file for the grid
Deleted: seismo/2D/SPECFEM2D/trunk/MAILLE90/Par_rouen_ok
===================================================================
--- seismo/2D/SPECFEM2D/trunk/MAILLE90/Par_rouen_ok 2004-05-16 01:48:11 UTC (rev 8418)
+++ seismo/2D/SPECFEM2D/trunk/MAILLE90/Par_rouen_ok 2007-12-07 23:44:44 UTC (rev 8419)
@@ -1,81 +0,0 @@
-# ----------------------------------------------------------------
-#
-# This is the parameter file
-# Put variable names first and actual value after 15th column
-#
-# ----------------------------------------------------------------
-<- ->
-#
-# File names and path for different outputs
-#
-title =Test pour INSA Rouen
-topofile =topo_rouen.dat
-interffile =none
-#
-# geometry of the model (origin lower-left corner = 0,0) and mesh description
-#
-xmin =0.0d0 ! abscissa of left side of the model
-xmax =4000.d0 ! abscissa of right side of the model
-nx =80 ! number of elements along X
-nz =60 ! number of elements along Z
-ngnod =4 ! nb noeuds de controle pour blocs generes (4 ou 9)
-ratio =0.967741 ! ratio pour separation en deux zones
-initialfield =.false. ! use a plane wave as source or not
-ireadmodel =.false. ! use a plane wave as source or not
-#
-# absorbing boundaries parameters
-#
-absorbhaut =.false. ! Absorbing boundary active or not
-absorbbas =.true.
-absorbgauche =.true.
-absorbdroite =.true.
-#
-# time step parameters
-#
-nt =1600 ! nb total de pas de temps
-dt =1.d-3 ! valeur du pas de temps
-#
-# source parameters
-#
-isource_surf =.false. ! sources dans le volume ou a la surface
-xs =2000. ! source location x in meters
-zs =1500. ! source location z in meters
-f0 =10.0 ! central source frequency (Hz)
-t0 =0.11 ! time delay of the source in seconds
-source_type =1 ! source type : force=1 or explosion=2
-angle =0. ! angle of the source (for a force only)
-factor =1.d10 ! amplification factor
-#
-# receiver line parameters
-#
-ienreg_surf =.false. ! enregistrement dans le volume ou a la surface
-isismostype =2 ! record 1=displacement 2=velocity 3=acceleration
-nrec =11 ! number of receivers
-xdeb =0. ! first receiver x in meters
-zdeb =2200. ! first receiver z in meters
-xfin =4000. ! last receiver x in meters
-zfin =2200. ! last receiver z in meters
-anglerec =0.d0 ! angle to rotate the components at the receivers
-#
-# display parameters
-#
-itaff =100 ! display frequency in time steps
-ivecttype =2 ! display 1=displacement 2=velocity 3=acceleration
-cutvect =1. ! amplitude min affichee en % pour vector plots
-imeshvect =.true. ! display mesh on vector plots or not
-imodelvect =.false. ! display velocity model on vector plots or not
-iboundvect =.true. ! display boundary conditions on vector plots
-interpol =.false. ! interpolation of the display or not
-iptsdisp =6 ! nb of points for interpolation of the display
-isubsamp =2 ! subsampling of color snapshots
-ignuplot =.false. ! generate a GNUPLOT file for the grid
-ioutputgrid =.false. ! save the grid in a text file or not
-#
-# velocity and density model (nx,nz)
-#
-nbmodels =2 ! nb de modeles differents (rho,vp,vs)
-1 0 2200.d0 2500.d0 1443.375d0 0 0
-2 0 2200.d0 2500.d0 1443.375d0 0 0
-nbzone =2 ! nb of zones and model number for each zone
-1 80 1 60 1
-3 51 3 38 2
Modified: seismo/2D/SPECFEM2D/trunk/MAILLE90/meshfem2D.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/MAILLE90/meshfem2D.f90 2004-05-16 01:48:11 UTC (rev 8418)
+++ seismo/2D/SPECFEM2D/trunk/MAILLE90/meshfem2D.f90 2007-12-07 23:44:44 UTC (rev 8419)
@@ -53,7 +53,6 @@
character(len=15) junk
integer imatnum,inumabs,inumelem,netyp
- integer icodehaut,icodebas,icodegauche,icodedroite
integer nelemabs,npgeo,nspec,ninterf,ntopo
integer k,icol,ili,istepx,istepz,ncut,ix,iz,irec,i,j
integer ixdebzone,ixfinzone,izdebzone,izfinzone,imodnum
@@ -63,6 +62,8 @@
integer ngnod,nt,nx,nz,nxread,nzread
integer icodematread
+ logical codehaut,codebas,codegauche,codedroite
+
double precision ratio
double precision tang1,tangN,vpzone,vszone
double precision cutvect
@@ -81,12 +82,6 @@
integer, external :: num
double precision, external :: bottom,spl,dens
-! --- code des numeros d'aretes pour les bords absorbants
- integer, parameter :: iaretebas = 1
- integer, parameter :: iaretedroite = 2
- integer, parameter :: iaretehaut = 3
- integer, parameter :: iaretegauche = 4
-
! ***
! *** read the parameter file
! ***
@@ -763,18 +758,18 @@
inumabs = 0
do iz=1,nzread
do ix=1,nxread
- icodehaut = 0
- icodebas = 0
- icodegauche = 0
- icodedroite = 0
+ codehaut = .false.
+ codebas = .false.
+ codegauche = .false.
+ codedroite = .false.
inumelem = (iz-1)*nxread + ix
- if(abshaut .and. iz==nzread) icodehaut = iaretehaut
- if(absbas .and. iz== 1) icodebas = iaretebas
- if(absgauche .and. ix== 1) icodegauche = iaretegauche
- if(absdroite .and. ix==nxread) icodedroite = iaretedroite
- if(icodehaut>0 .or. icodebas>0 .or. icodegauche>0 .or. icodedroite>0) then
+ if(abshaut .and. iz==nzread) codehaut = .true.
+ if(absbas .and. iz== 1) codebas = .true.
+ if(absgauche .and. ix== 1) codegauche = .true.
+ if(absdroite .and. ix==nxread) codedroite = .true.
+ if(codehaut .or. codebas .or. codegauche .or. codedroite) then
inumabs = inumabs + 1
- write(15,*) inumabs,inumelem,icodehaut,icodebas,icodegauche,icodedroite
+ write(15,*) inumabs,inumelem,codehaut,codebas,codegauche,codedroite
endif
enddo
enddo
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/Makefile
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/Makefile 2004-05-16 01:48:11 UTC (rev 8418)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/Makefile 2007-12-07 23:44:44 UTC (rev 8419)
@@ -24,10 +24,9 @@
LINK = $(F90)
EXEC = xspecfem2D
OBJS = $O/checkgrid.o $O/datim.o $O/defarrays.o\
- $O/lagrange_poly.o $O/gmat01.o $O/gll_library.o $O/plotgll.o\
+ $O/lagrange_poly.o $O/gmat01.o $O/gll_library.o $O/plotgll.o $O/define_derivative_matrices.o\
$O/plotpost.o $O/positrec.o $O/positsource.o $O/q49spec.o\
- $O/qmasspec.o $O/qsumspec.o $O/specfem2D.o\
- $O/write_seismograms.o $O/createnum_fast.o $O/createnum_slow.o $O/q49shape.o
+ $O/specfem2D.o $O/write_seismograms.o $O/createnum_fast.o $O/createnum_slow.o $O/q49shape.o
default: all
@@ -61,6 +60,9 @@
$O/gll_library.o: gll_library.f90 constants.h
${F90} $(FLAGS) -c -o $O/gll_library.o gll_library.f90
+$O/define_derivative_matrices.o: define_derivative_matrices.f90 constants.h
+ ${F90} $(FLAGS) -c -o $O/define_derivative_matrices.o define_derivative_matrices.f90
+
$O/plotgll.o: plotgll.f90 constants.h
${F90} $(FLAGS) -c -o $O/plotgll.o plotgll.f90
@@ -79,12 +81,6 @@
$O/q49spec.o: q49spec.f90 constants.h
${F90} $(FLAGS) -c -o $O/q49spec.o q49spec.f90
-$O/qmasspec.o: qmasspec.f90 constants.h
- ${F90} $(FLAGS) -c -o $O/qmasspec.o qmasspec.f90
-
-$O/qsumspec.o: qsumspec.f90 constants.h
- ${F90} $(FLAGS) -c -o $O/qsumspec.o qsumspec.f90
-
$O/specfem2D.o: specfem2D.f90 constants.h
${F90} $(FLAGS) -c -o $O/specfem2D.o specfem2D.f90
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/constants.h
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/constants.h 2004-05-16 01:48:11 UTC (rev 8418)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/constants.h 2007-12-07 23:44:44 UTC (rev 8419)
@@ -1,7 +1,7 @@
! polynomial degree
integer, parameter :: NGLLX = 5
- integer, parameter :: NGLLY = NGLLX
+ integer, parameter :: NGLLZ = NGLLX
! select fast (Paul Fischer) or slow (topology only) global numbering algorithm
logical, parameter :: FAST_NUMBERING = .true.
@@ -18,20 +18,18 @@
! integer, parameter :: IOUT = 41
! flags for absorbing boundaries
- integer, parameter :: IHAUT = 1
- integer, parameter :: IBAS = 2
- integer, parameter :: IGAUCHE = 3
- integer, parameter :: IDROITE = 4
+ integer, parameter :: ITOP = 1
+ integer, parameter :: IBOTTOM = 2
+ integer, parameter :: ILEFT = 3
+ integer, parameter :: IRIGHT = 4
- integer, parameter :: IARETEBAS = 1
- integer, parameter :: IARETEDROITE = 2
- integer, parameter :: IARETEHAUT = 3
- integer, parameter :: IARETEGAUCHE = 4
-
! a few useful constants
double precision, parameter :: ZERO = 0.d0,ONE = 1.d0
- double precision, parameter :: HALF = 0.5d0,TWO = 2.0d0
+ double precision, parameter :: HALF = 0.5d0,TWO = 2.0d0,QUART = 0.25d0
+! pi
double precision, parameter :: PI = 3.141592653589793d0
+! 4/3
+ double precision, parameter :: FOUR_THIRDS = 4.d0/3.d0
! parameters to define the Gauss-Lobatto-Legendre points
double precision, parameter :: GAUSSALPHA = ZERO,GAUSSBETA = ZERO
@@ -75,3 +73,40 @@
! taille de la fenetre de display Postscript en pourcentage de la feuille
double precision, parameter :: RPERCENTX = 70.0d0,RPERCENTZ = 77.0d0
+!-----------------------------------------------------------------------
+
+!! DK DK
+!! DK DK anisotropic copper crystal (cubic symmetry)
+!! DK DK
+
+! switch anisotropy on or off
+ logical, parameter :: TURN_ANISOTROPY_ON = .false.
+
+!! DK DK regular c_ijkl with no rotation
+ double precision, parameter :: c11val = 169.d9
+ double precision, parameter :: c12val = 122.d9
+ double precision, parameter :: c13val = c12val
+ double precision, parameter :: c14val = 0.d0
+ double precision, parameter :: c15val = 0.d0
+ double precision, parameter :: c16val = 0.d0
+
+ double precision, parameter :: c22val = c11val
+ double precision, parameter :: c23val = c12val
+ double precision, parameter :: c24val = 0.d0
+ double precision, parameter :: c25val = 0.d0
+ double precision, parameter :: c26val = 0.d0
+
+ double precision, parameter :: c33val = c11val
+ double precision, parameter :: c34val = 0.d0
+ double precision, parameter :: c35val = 0.d0
+ double precision, parameter :: c36val = 0.d0
+
+ double precision, parameter :: c44val = 75.3d9
+ double precision, parameter :: c45val = 0.d0
+ double precision, parameter :: c46val = 0.d0
+
+ double precision, parameter :: c55val = c44val
+ double precision, parameter :: c56val = 0.d0
+
+ double precision, parameter :: c66val = c44val
+
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/createnum_fast.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/createnum_fast.f90 2004-05-16 01:48:11 UTC (rev 8418)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/createnum_fast.f90 2007-12-07 23:44:44 UTC (rev 8419)
@@ -20,7 +20,7 @@
include "constants.h"
integer npoin,npgeo,nspec,ngnod
- integer knods(ngnod,nspec),ibool(NGLLX,NGLLY,nspec)
+ integer knods(ngnod,nspec),ibool(NGLLX,NGLLZ,nspec)
double precision shape(ngnod,NGLLX,NGLLX)
double precision coorg(NDIME,npgeo)
@@ -43,7 +43,7 @@
print *,'Generating global mesh numbering (fast version)...'
print *
- nxyz = NGLLX*NGLLY
+ nxyz = NGLLX*NGLLZ
ntot = nxyz*nspec
allocate(loc(ntot))
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/createnum_slow.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/createnum_slow.f90 2004-05-16 01:48:11 UTC (rev 8418)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/createnum_slow.f90 2007-12-07 23:44:44 UTC (rev 8419)
@@ -21,7 +21,7 @@
integer npoin,nspec,ngnod
- integer knods(ngnod,nspec),ibool(NGLLX,NGLLY,nspec)
+ integer knods(ngnod,nspec),ibool(NGLLX,NGLLZ,nspec)
integer i,j,num2,i2,j2,ipos,ipos2,iloc,jloc,kloc
integer ngnodloc,ngnodother,nedgeloc,nedgeother,npedge,numelem,npcorn
@@ -62,7 +62,7 @@
do numelem = 1,nspec
do i=1,NGLLX
- do j=1,NGLLY
+ do j=1,NGLLZ
! verifier que le point n'a pas deja ete genere
@@ -71,7 +71,7 @@
!
!---- point interieur a un element, donc forcement unique
!
- if(i /= 1 .and. i /= NGLLX .and. j /= 1 .and. j /= NGLLY) then
+ if(i /= 1 .and. i /= NGLLX .and. j /= 1 .and. j /= NGLLZ) then
npoin = npoin + 1
ibool(i,j,numelem) = npoin
@@ -79,17 +79,17 @@
!
!---- point au coin d'un element, rechercher les coins des autres elements
!
- else if((i == 1 .and. j == 1) .or. (i == 1 .and. j == NGLLY) .or. &
- (i == NGLLX .and. j == 1) .or. (i == NGLLX .and. j == NGLLY)) then
+ else if((i == 1 .and. j == 1) .or. (i == 1 .and. j == NGLLZ) .or. &
+ (i == NGLLX .and. j == 1) .or. (i == NGLLX .and. j == NGLLZ)) then
! trouver numero local du coin
if(i == 1 .and. j == 1) then
ngnodloc = 1
else if(i == NGLLX .and. j == 1) then
ngnodloc = 2
- else if(i == NGLLX .and. j == NGLLY) then
+ else if(i == NGLLX .and. j == NGLLZ) then
ngnodloc = 3
- else if(i == 1 .and. j == NGLLY) then
+ else if(i == 1 .and. j == NGLLZ) then
ngnodloc = 4
endif
@@ -117,10 +117,10 @@
j2 = 1
else if(ngnodother == 3) then
i2 = NGLLX
- j2 = NGLLY
+ j2 = NGLLZ
else if(ngnodother == 4) then
i2 = 1
- j2 = NGLLY
+ j2 = NGLLZ
else
stop 'bad corner'
endif
@@ -156,7 +156,7 @@
nedgeloc = 1
else if(i == NGLLX) then
nedgeloc = 2
- else if(j == NGLLY) then
+ else if(j == NGLLZ) then
nedgeloc = 3
else if(i == 1) then
nedgeloc = 4
@@ -189,7 +189,7 @@
alreadyexist = .true.
! obtenir la numerotation dans l'autre element
-! maillage conforme donc on doit supposer que NGLLX == NGLLY
+! maillage conforme donc on doit supposer que NGLLX == NGLLZ
! generer toute l'arete pour eviter des recherches superflues
do kloc = 2,NGLLX-1
@@ -205,12 +205,12 @@
ipos = jloc
else if(nedgeloc == 3) then
iloc = kloc
- jloc = NGLLY
+ jloc = NGLLZ
ipos = NGLLX - iloc + 1
else if(nedgeloc == 4) then
iloc = 1
jloc = kloc
- ipos = NGLLY - jloc + 1
+ ipos = NGLLZ - jloc + 1
else
stop 'bad nedgeloc'
endif
@@ -229,10 +229,10 @@
j2 = ipos2
else if(nedgeother == 3) then
i2 = NGLLX - ipos2 + 1
- j2 = NGLLY
+ j2 = NGLLZ
else if(nedgeother == 4) then
i2 = 1
- j2 = NGLLY - ipos2 + 1
+ j2 = NGLLZ - ipos2 + 1
else
stop 'bad nedgeother'
endif
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/defarrays.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/defarrays.f90 2004-05-16 01:48:11 UTC (rev 8418)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/defarrays.f90 2007-12-07 23:44:44 UTC (rev 8419)
@@ -12,74 +12,63 @@
!========================================================================
subroutine defarrays(vpext,vsext,rhoext,density,elastcoef, &
- xigll,yigll,wxgll,wygll,hprime,hTprime, &
- a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13x,a13z, &
- ibool,kmato,dvolu,xjaci,coord,gltfu, &
- numabs,codeabs,anyabs,npoin,rsizemin,rsizemax, &
- cpoverdxmin,cpoverdxmax,rlamdaSmin,rlamdaSmax, &
- rlamdaPmin,rlamdaPmax,vpmin,vpmax,ireadmodel,nelemabs,nspec,numat)
+ xigll,zigll,xix,xiz,gammax,gammaz,a11,a12, &
+ ibool,kmato,coord,gltfu,npoin,rsizemin,rsizemax, &
+ cpoverdxmin,cpoverdxmax,rlambdaSmin,rlambdaSmax, &
+ rlambdaPmin,rlambdaPmax,vpmin,vpmax,ireadmodel,nspec,numat)
-
! define all the arrays for the variational formulation
implicit none
include "constants.h"
- integer i,j,ispec,material,ipointnum,npoin,nelemabs,nspec,numat
- integer isourx,isourz,ielems,ir,is,ip
+ integer i,j,ispec,material,ipointnum,npoin,nspec,numat
+ integer isourx,isourz,ielems,ir,is
integer kmato(nspec),ibool(NGLLX,NGLLX,nspec)
- double precision density(numat),elastcoef(4,numat), &
- xigll(NGLLX),yigll(NGLLY),wxgll(NGLLX),wygll(NGLLX), &
- dvolu(nspec,NGLLX,NGLLX), &
- xjaci(nspec,NDIME,NDIME,NGLLX,NGLLX), &
- hprime(NGLLX,NGLLX), hTprime(NGLLX,NGLLX)
+ double precision xix(NGLLX,NGLLZ,nspec)
+ double precision xiz(NGLLX,NGLLZ,nspec)
+ double precision gammax(NGLLX,NGLLZ,nspec)
+ double precision gammaz(NGLLX,NGLLZ,nspec)
+ double precision density(numat),elastcoef(4,numat)
+
double precision coord(NDIME,npoin)
- double precision a1(NGLLX,NGLLX,nspec),a2(NGLLX,NGLLX,nspec), &
- a3(NGLLX,NGLLX,nspec),a4(NGLLX,NGLLX,nspec), &
- a5(NGLLX,NGLLX,nspec),a6(NGLLX,NGLLX,nspec), &
- a7(NGLLX,NGLLX,nspec),a8(NGLLX,NGLLX,nspec), &
- a9(NGLLX,NGLLX,nspec),a10(NGLLX,NGLLX,nspec)
- double precision a13x(NGLLX,NGLLX,nelemabs),a13z(NGLLX,NGLLX,nelemabs)
+
double precision a11(NGLLX,NGLLX),a12(NGLLX,NGLLX)
+ double precision xigll(NGLLX),zigll(NGLLZ)
+
double precision gltfu(20)
double precision vpext(npoin)
double precision vsext(npoin)
double precision rhoext(npoin)
double precision vsmin,vsmax,densmin,densmax
- double precision rKmod,rlamda,rmu,xix,xiz,etax,etaz,denst,rjacob
- double precision rKvol,cploc,csloc,xxi,zeta,rwgll,x0,z0
- double precision c11,c13,c33,c44
+ double precision rKmod,rlambda,rmu,denst
+ double precision rKvol,cploc,csloc,x0,z0
double precision x1,z1,x2,z2,rdist1,rdist2,rapportmin,rapportmax
- double precision rlambmin,rlambmax,coefintegr
+ double precision rlambmin,rlambmax
double precision flagxprime,flagzprime,sig0
double precision rsizemin,rsizemax,cpoverdxmin,cpoverdxmax, &
- rlamdaSmin,rlamdaSmax,rlamdaPmin,rlamdaPmax,vpmin,vpmax
+ rlambdaSmin,rlambdaSmax,rlambdaPmin,rlambdaPmax,vpmin,vpmax
- logical anyabs,anisotrope,ireadmodel
+ logical ireadmodel
- integer numabs(nelemabs),codeabs(4,nelemabs)
-
double precision, external :: lagrange_deriv_GLL
!
!-----------------------------------------------------------------------
!
-!---- compute parameters a1 to a13 for the spectral elements
+!---- compute parameters for the spectral elements
a11 = zero
a12 = zero
- a13x = zero
- a13z = zero
-
vpmin = HUGEVAL
vsmin = HUGEVAL
vpmax = -HUGEVAL
@@ -93,50 +82,27 @@
cpoverdxmin = HUGEVAL
cpoverdxmax = -HUGEVAL
- rlamdaPmin = HUGEVAL
- rlamdaSmin = HUGEVAL
- rlamdaPmax = -HUGEVAL
- rlamdaSmax = -HUGEVAL
+ rlambdaPmin = HUGEVAL
+ rlambdaSmin = HUGEVAL
+ rlambdaPmax = -HUGEVAL
+ rlambdaSmax = -HUGEVAL
do ispec=1,nspec
-
material = kmato(ispec)
- rlamda = elastcoef(1,material)
+ rlambda = elastcoef(1,material)
rmu = elastcoef(2,material)
rKmod = elastcoef(3,material)
denst = density(material)
- rKvol = rlamda + 2.d0*rmu/3.d0
+ rKvol = rlambda + 2.d0*rmu/3.d0
cploc = dsqrt((rKvol + 4.d0*rmu/3.d0)/denst)
csloc = dsqrt(rmu/denst)
-! determiner si le materiau est anisotrope ou non
- if(elastcoef(4,material) > 0.00001d0) then
- anisotrope = .true.
- c11 = elastcoef(1,material)
- c13 = elastcoef(2,material)
- c33 = elastcoef(3,material)
- c44 = elastcoef(4,material)
- else
- anisotrope = .false.
- endif
-
- do j=1,NGLLY
+ do j=1,NGLLZ
do i=1,NGLLX
- xix = xjaci(ispec,1,1,i,j)
- xiz = xjaci(ispec,1,2,i,j)
- etax = xjaci(ispec,2,1,i,j)
- etaz = xjaci(ispec,2,2,i,j)
- rjacob = dvolu(ispec,i,j)
-
- xxi = etaz * rjacob
- zeta = xix * rjacob
-
- rwgll = - wxgll(i)*wygll(j)
-
!--- si formulation heterogene pour un modele de vitesse externe
if(ireadmodel) then
ipointnum = ibool(i,j,ispec)
@@ -144,16 +110,10 @@
csloc = vsext(ipointnum)
denst = rhoext(ipointnum)
rmu = denst*csloc*csloc
- rlamda = denst*cploc*cploc - 2.d0*rmu
- rKmod = rlamda + 2.d0*rmu
+ rlambda = denst*cploc*cploc - 2.d0*rmu
+ rKmod = rlambda + 2.d0*rmu
endif
-!--- si materiau transverse isotrope, donner une idee des proprietes
- if(anisotrope) then
- cploc = sqrt(c11/denst)
- csloc = sqrt(c44/denst)
- endif
-
!--- calculer min et max du modele de vitesse et densite
vpmin = dmin1(vpmin,cploc)
vpmax = dmax1(vpmax,cploc)
@@ -165,7 +125,7 @@
densmax = dmax1(densmax,denst)
!--- stocker parametres pour verifs diverses
- if(i < NGLLX .and. j < NGLLY) then
+ if(i < NGLLX .and. j < NGLLZ) then
x0 = coord(1,ibool(i,j,ispec))
z0 = coord(2,ibool(i,j,ispec))
@@ -191,51 +151,24 @@
z0 = coord(2,ibool(1,1,ispec))
x1 = coord(1,ibool(NGLLX,1,ispec))
z1 = coord(2,ibool(NGLLX,1,ispec))
- x2 = coord(1,ibool(1,NGLLY,ispec))
- z2 = coord(2,ibool(1,NGLLY,ispec))
+ x2 = coord(1,ibool(1,NGLLZ,ispec))
+ z2 = coord(2,ibool(1,NGLLZ,ispec))
rdist1 = dsqrt((x1-x0)**2 + (z1-z0)**2)
rdist2 = dsqrt((x2-x0)**2 + (z2-z0)**2)
rlambmin = cploc/dmax1(rdist1,rdist2)
rlambmax = cploc/dmin1(rdist1,rdist2)
- rlamdaPmin = dmin1(rlamdaPmin,rlambmin)
- rlamdaPmax = dmax1(rlamdaPmax,rlambmax)
+ rlambdaPmin = dmin1(rlambdaPmin,rlambmin)
+ rlambdaPmax = dmax1(rlambdaPmax,rlambmax)
rlambmin = csloc/dmax1(rdist1,rdist2)
rlambmax = csloc/dmin1(rdist1,rdist2)
- rlamdaSmin = dmin1(rlamdaSmin,rlambmin)
- rlamdaSmax = dmax1(rlamdaSmax,rlambmax)
+ rlambdaSmin = dmin1(rlambdaSmin,rlambmin)
+ rlambdaSmax = dmax1(rlambdaSmax,rlambmax)
endif
-!--- definir tableaux
- if(.not. anisotrope) then
- a1(i,j,ispec) = rwgll*(rKmod*xix*xix + rmu*xiz*xiz)*rjacob
- a2(i,j,ispec) = rwgll*(rKmod*etax*xix + rmu*etaz*xiz)*rjacob
- a3(i,j,ispec) = rwgll*(rlamda+rmu)*xiz*xix*rjacob
- a4(i,j,ispec) = rwgll*(rlamda*etaz*xix + rmu*etax*xiz)*rjacob
- a5(i,j,ispec) = rwgll*(rKmod*etaz*etaz + rmu*etax*etax)*rjacob
- a6(i,j,ispec) = rwgll*(rKmod*etax*etax + rmu*etaz*etaz)*rjacob
- a7(i,j,ispec) = rwgll*(rlamda*etax*xiz + rmu*etaz*xix)*rjacob
- a8(i,j,ispec) = rwgll*(rlamda+rmu)*etax*etaz*rjacob
- a9(i,j,ispec) = rwgll*(rKmod*xiz*xiz + rmu*xix*xix)*rjacob
- a10(i,j,ispec) = rwgll*(rKmod*etaz*xiz + rmu*etax*xix)*rjacob
- else
- a3(i,j,ispec) = rwgll*(c13+c44)*xiz*xix*rjacob
- a4(i,j,ispec) = rwgll*(c13*etaz*xix + c44*etax*xiz)*rjacob
- a7(i,j,ispec) = rwgll*(c13*etax*xiz + c44*etaz*xix)*rjacob
- a8(i,j,ispec) = rwgll*(c13+c44)*etax*etaz*rjacob
-
- a1(i,j,ispec) = rwgll*(c11*xix*xix + c44*xiz*xiz)*rjacob
- a2(i,j,ispec) = rwgll*(c11*etax*xix + c44*etaz*xiz)*rjacob
- a6(i,j,ispec) = rwgll*(c11*etax*etax + c44*etaz*etaz)*rjacob
-
- a5(i,j,ispec) = rwgll*(c33*etaz*etaz + c44*etax*etax)*rjacob
- a9(i,j,ispec) = rwgll*(c33*xiz*xiz + c44*xix*xix)*rjacob
- a10(i,j,ispec) = rwgll*(c33*etaz*xiz + c44*etax*xix)*rjacob
- endif
-
enddo
enddo
enddo
@@ -248,87 +181,6 @@
print *,'********'
print *
-!
-!--- definition coefficients pour bords absorbants
-!
-
- if(anyabs) then
-
- do ispec=1,nelemabs
-
- material = kmato(numabs(ispec))
-
- rlamda = elastcoef(1,material)
- rmu = elastcoef(2,material)
- rKmod = elastcoef(3,material)
- denst = density(material)
-
- rKvol = rlamda + 2.d0*rmu/3.d0
- cploc = dsqrt((rKvol + 4.d0*rmu/3.d0)/denst)
- csloc = dsqrt(rmu/denst)
-
- do i=1,NGLLX
- do j=1,NGLLY
-
-!--- si formulation heterogene pour un modele de vitesse externe
- if(ireadmodel) then
- ipointnum = ibool(i,j,numabs(ispec))
- cploc = vpext(ipointnum)
- csloc = vsext(ipointnum)
- denst = rhoext(ipointnum)
- rmu = denst*csloc*csloc
- rlamda = denst*cploc*cploc - 2.d0*rmu
- rKmod = rlamda + 2.d0*rmu
- endif
-
- xix = xjaci(numabs(ispec),1,1,i,j)
- xiz = xjaci(numabs(ispec),1,2,i,j)
- etax = xjaci(numabs(ispec),2,1,i,j)
- etaz = xjaci(numabs(ispec),2,2,i,j)
- rjacob = dvolu(numabs(ispec),i,j)
-
- xxi = etaz * rjacob
- zeta = xix * rjacob
-
- rwgll = - wxgll(i)*wygll(j)
-
-!---- sommer les contributions dans les coins pour l'ancienne formulation
-!---- ne pas sommer les contributions dans les coins pour la nouvelle
-
-! bord absorbant du bas
- if(codeabs(ibas,ispec) /= 1 .and. j == 1) then
- coefintegr = wxgll(i)*xxi
- a13x(i,j,ispec) = denst*csloc*coefintegr
- a13z(i,j,ispec) = denst*cploc*coefintegr
- endif
-
-! bord absorbant du haut (signe moins)
- if(codeabs(ihaut,ispec) /= 1 .and. j == NGLLY) then
- coefintegr = wxgll(i)*xxi
- a13x(i,j,ispec) = denst*csloc*coefintegr
- a13z(i,j,ispec) = denst*cploc*coefintegr
- endif
-
-! bord absorbant de gauche
- if(codeabs(igauche,ispec) /= 1 .and. i == 1) then
- coefintegr = wygll(j)*zeta
- a13x(i,j,ispec) = denst*cploc*coefintegr
- a13z(i,j,ispec) = denst*csloc*coefintegr
- endif
-
-! bord absorbant de droite
- if(codeabs(idroite,ispec) /= 1 .and. i == NGLLX) then
- coefintegr = wygll(j)*zeta
- a13x(i,j,ispec) = denst*cploc*coefintegr
- a13z(i,j,ispec) = denst*csloc*coefintegr
- endif
-
- enddo
- enddo
- enddo
-
- endif
-
! seulement si source explosive
if(nint(gltfu(2)) == 2) then
@@ -336,40 +188,26 @@
isourz = nint(gltfu(11))
ielems = nint(gltfu(12))
- if(isourx == 1.or.isourx == NGLLX.or.isourz == 1 .or.isourz == NGLLX) &
+ if(isourx == 1 .or. isourx == NGLLX .or. isourz == 1 .or. isourz == NGLLX) &
stop 'Explosive source on element edge'
!---- definir a11 et a12 - dirac (schema en croix)
- xix = xjaci(ielems,1,1,isourx,isourz)
- xiz = xjaci(ielems,1,2,isourx,isourz)
- etax = xjaci(ielems,2,1,isourx,isourz)
- etaz = xjaci(ielems,2,2,isourx,isourz)
-
sig0 = one
do ir=1,NGLLX
flagxprime = lagrange_deriv_GLL(ir-1,isourx-1,xigll,NGLLX)
- a11(ir,isourz) = a11(ir,isourz) + sig0*xix*flagxprime
- a12(ir,isourz) = a12(ir,isourz) + sig0*xiz*flagxprime
+ a11(ir,isourz) = a11(ir,isourz) + sig0*xix(isourx,isourz,ielems)*flagxprime
+ a12(ir,isourz) = a12(ir,isourz) + sig0*xiz(isourx,isourz,ielems)*flagxprime
enddo
- do is=1,NGLLY
- flagzprime = lagrange_deriv_GLL(is-1,isourz-1,yigll,NGLLY)
- a11(isourx,is) = a11(isourx,is) + sig0*etax*flagzprime
- a12(isourx,is) = a12(isourx,is) + sig0*etaz*flagzprime
+ do is=1,NGLLZ
+ flagzprime = lagrange_deriv_GLL(is-1,isourz-1,zigll,NGLLZ)
+ a11(isourx,is) = a11(isourx,is) + sig0*gammax(isourx,isourz,ielems)*flagzprime
+ a12(isourx,is) = a12(isourx,is) + sig0*gammaz(isourx,isourz,ielems)*flagzprime
enddo
endif
-!---- compute hprime coefficients (derivatives of Lagrange polynomials)
-!---- (works only if NGLLX = NGLLY)
- do ip=1,NGLLX
- do i=1,NGLLX
- hprime(ip,i) = lagrange_deriv_GLL(ip-1,i-1,xigll,NGLLX)
- hTprime(i,ip) = hprime(ip,i)
- enddo
- enddo
-
end subroutine defarrays
Added: seismo/2D/SPECFEM2D/trunk/SPECFEM90/define_derivative_matrices.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/define_derivative_matrices.f90 2004-05-16 01:48:11 UTC (rev 8418)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/define_derivative_matrices.f90 2007-12-07 23:44:44 UTC (rev 8419)
@@ -0,0 +1,62 @@
+
+!========================================================================
+!
+! S P E C F E M 2 D Version 5.0
+! ------------------------------
+!
+! Dimitri Komatitsch
+! Universite de Pau et des Pays de l'Adour, France
+!
+! (c) May 2004
+!
+!========================================================================
+
+ subroutine define_derivative_matrices(xigll,zigll,wxgll,wzgll,hprime_xx,hprime_zz)
+
+ implicit none
+
+ include "constants.h"
+
+! Gauss-Lobatto-Legendre points of integration
+ double precision, dimension(NGLLX) :: xigll
+ double precision, dimension(NGLLZ) :: zigll
+
+! weights
+ double precision, dimension(NGLLX) :: wxgll
+ double precision, dimension(NGLLZ) :: wzgll
+
+! array with derivatives of Lagrange polynomials
+ double precision, dimension(NGLLX,NGLLX) :: hprime_xx
+ double precision, dimension(NGLLZ,NGLLZ) :: hprime_zz
+
+! function for calculating derivatives of Lagrange polynomials
+ double precision, external :: lagrange_deriv_GLL
+
+ integer i1,i2,k1,k2
+
+! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
+! set up coordinates of the Gauss-Lobatto-Legendre points
+ call zwgljd(xigll,wxgll,NGLLX,GAUSSALPHA,GAUSSBETA)
+ call zwgljd(zigll,wzgll,NGLLZ,GAUSSALPHA,GAUSSBETA)
+
+! if number of points is odd, the middle abscissa is exactly zero
+ if(mod(NGLLX,2) /= 0) xigll((NGLLX-1)/2+1) = ZERO
+ if(mod(NGLLZ,2) /= 0) zigll((NGLLZ-1)/2+1) = ZERO
+
+! 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,NGLLX
+ do i2=1,NGLLX
+ hprime_xx(i1,i2) = lagrange_deriv_GLL(i1-1,i2-1,xigll,NGLLX)
+ enddo
+ enddo
+
+ do k1=1,NGLLZ
+ do k2=1,NGLLZ
+ hprime_zz(k1,k2) = lagrange_deriv_GLL(k1-1,k2-1,zigll,NGLLZ)
+ enddo
+ enddo
+
+ end subroutine define_derivative_matrices
+
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/plotgll.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/plotgll.f90 2004-05-16 01:48:11 UTC (rev 8418)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/plotgll.f90 2007-12-07 23:44:44 UTC (rev 8419)
@@ -60,7 +60,7 @@
!
!---- plot the lines in xi-direction
!
- do iy = 1,NGLLY
+ do iy = 1,NGLLZ
do ix = 1,NGLLX-1
!
!---- get the global point number
@@ -75,7 +75,7 @@
write(20,15) coord(1,iglobnum2),coord(2,iglobnum2)
write(20,10)
- if(iy == 1 .or. iy == NGLLY) then
+ if(iy == 1 .or. iy == NGLLZ) then
write(21,15) coord(1,iglobnum),coord(2,iglobnum)
write(21,15) coord(1,iglobnum2),coord(2,iglobnum2)
write(21,10)
@@ -88,7 +88,7 @@
!---- plot the lines in eta-direction
!
do ix = 1,NGLLX
- do iy = 1,NGLLY-1
+ do iy = 1,NGLLZ-1
!
!---- get the global point number
!
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/plotpost.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/plotpost.f90 2004-05-16 01:48:11 UTC (rev 8418)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/plotpost.f90 2007-12-07 23:44:44 UTC (rev 8419)
@@ -34,7 +34,7 @@
integer i,iglobrec,iglobsource,npoin,npgeo,ngnod
integer kmato(nspec),knods(ngnod,nspec)
- integer ibool(NGLLX,NGLLY,nspec)
+ integer ibool(NGLLX,NGLLZ,nspec)
double precision xinterp(iptsdisp,iptsdisp),zinterp(iptsdisp,iptsdisp)
double precision shapeint(ngnod,iptsdisp,iptsdisp)
@@ -181,8 +181,7 @@
zmin=0.d0
! rapport taille page/taille domaine physique
- rapp_page = min(rpercentz*sizez/(zmax-zmin),rpercentx*sizex/(xmax-xmin)) &
- / 100.d0
+ rapp_page = min(rpercentz*sizez/(zmax-zmin),rpercentx*sizex/(xmax-xmin)) / 100.d0
! recherche de la valeur maximum de la norme du deplacement
dispmax = maxval(sqrt(displ(1,:)**2 + displ(2,:)**2))
@@ -284,17 +283,17 @@
write(24,*) '/Times-Roman findfont'
write(24,*) '.5 CM SCSF'
- if (legendes) then
+ if(legendes) then
write(24,*) '24. CM 1.2 CM MV'
write(24,610) usoffset,it
write(24,*) '%'
write(24,*) '24. CM 1.95 CM MV'
timeval = it*dt
- if(timeval >= 1.d-3) then
- write(24,600) usoffset,timeval
+ if(timeval >= 1.d-3) then
+ write(24,600) usoffset,timeval
else
- write(24,601) usoffset,timeval
+ write(24,601) usoffset,timeval
endif
write(24,*) '%'
write(24,*) '24. CM 2.7 CM MV'
@@ -306,7 +305,7 @@
write(24,*) '%'
write(24,*) '/Times-Roman findfont'
write(24,*) '.6 CM SCSF'
- if (icolor == 1) write(24,*) '.4 .9 .9 setrgbcolor'
+ if(icolor == 1) write(24,*) '.4 .9 .9 setrgbcolor'
write(24,*) '11 CM 1.1 CM MV'
write(24,*) '(X axis) show'
write(24,*) '%'
@@ -317,15 +316,15 @@
write(24,*) '%'
write(24,*) '/Times-Roman findfont'
write(24,*) '.7 CM SCSF'
- if (icolor == 1) write(24,*) '.8 0 .8 setrgbcolor'
+ if (icolor == 1) write(24,*) '.8 0 .8 setrgbcolor'
write(24,*) '24.35 CM 18.9 CM MV'
write(24,*) usoffset,' CM 2 div neg 0 MR'
write(24,*) 'currentpoint gsave translate -90 rotate 0 0 moveto'
- if (ivecttype == 1) then
+ if (ivecttype == 1) then
write(24,*) '(Displacement vector field) show'
- else if (ivecttype == 2) then
+ else if (ivecttype == 2) then
write(24,*) '(Velocity vector field) show'
- else if (ivecttype == 3) then
+ else if (ivecttype == 3) then
write(24,*) '(Acceleration vector field) show'
else
stop 'Bad field code in PostScript display'
@@ -365,7 +364,7 @@
do i=1,NGLLX-isubsamp,isubsamp
do j=1,NGLLX-isubsamp,isubsamp
- if((vpmax-vpmin)/vpmin > 0.02d0) then
+ if((vpmax-vpmin)/vpmin > 0.02d0) then
if(ireadmodel) then
x1 = (vpext(ibool(i,j,ispec))-vpmin)/ (vpmax-vpmin)
else
@@ -383,7 +382,7 @@
! rescaler pour eviter gris trop sombre
x1 = x1*0.7 + 0.2
- if (x1 > 1.d0) x1=1.d0
+ if (x1 > 1.d0) x1=1.d0
! inverser echelle : blanc = vpmin, gris = vpmax
x1 = 1.d0 - x1
@@ -459,7 +458,7 @@
write(24,*) 'MK'
write(24,681) x1,z1
- if (ngnod == 4) then
+ if (ngnod == 4) then
! tracer des droites si elements 4 noeuds
@@ -536,7 +535,7 @@
write(24,*) 'CO'
- if (icolor == 1) then
+ if (icolor == 1) then
! For the moment 20 different colors max
nbcols = 20
@@ -557,7 +556,7 @@
! write the element number, the group number and the
! material number inside the element
- if (inumber == 1) then
+ if (inumber == 1) then
xw = (coorg(1,knods(1,ispec)) + coorg(1,knods(2,ispec)) + &
coorg(1,knods(3,ispec)) + coorg(1,knods(4,ispec))) / 4.d0
@@ -567,7 +566,7 @@
zw = (zw-zmin)*rapp_page + orig_z
xw = xw * centim
zw = zw * centim
- if (icolor == 1) write(24,*) '1 setgray'
+ if (icolor == 1) write(24,*) '1 setgray'
write(24,500) xw,zw
@@ -607,21 +606,21 @@
do ibord = 1,4
- if(codeabs(ibord,ispecabs) /= 0) then
+ if(codeabs(ibord,ispecabs) /= 0) then
- if(ibord == ihaut) then
+ if(ibord == ITOP) then
write(24,*) '1. .85 0. RG'
ideb = 3
ifin = 4
- else if(ibord == ibas) then
+ else if(ibord == IBOTTOM) then
write(24,*) '.4 1. .4 RG'
ideb = 1
ifin = 2
- else if(ibord == igauche) then
+ else if(ibord == ILEFT) then
write(24,*) '1. .43 1. RG'
ideb = 4
ifin = 1
- else if(ibord == idroite) then
+ else if(ibord == IRIGHT) then
write(24,*) '.25 1. 1. RG'
ideb = 2
ifin = 3
@@ -656,7 +655,7 @@
!
! return if the maximum displacement equals zero (no source)
- if (dispmax == 0.d0) then
+ if (dispmax == 0.d0) then
print *,' null displacement : returning !'
return
endif
@@ -679,8 +678,7 @@
do ispec=1,nspec
! interpolation sur grille reguliere
- if(mod(ispec,100) == 0) &
- write(*,*) 'Interpolation uniform grid element ',ispec
+ if(mod(ispec,1000) == 0) write(*,*) 'Interpolation uniform grid element ',ispec
do i=1,iptsdisp
do j=1,iptsdisp
@@ -716,17 +714,17 @@
d = dsqrt(x2**2 + z2**2)
! ignorer si vecteur trop petit
- if (d > cutvect*sizemax) then
+ if (d > cutvect*sizemax) then
d1 = d * rapport
d2 = d1 * dcos(angle*convert)
dummy = x2/d
- if (dummy > 0.9999d0) dummy = 0.9999d0
- if (dummy < -0.9999d0) dummy = -0.9999d0
+ if (dummy > 0.9999d0) dummy = 0.9999d0
+ if (dummy < -0.9999d0) dummy = -0.9999d0
theta = dacos(dummy)
- if(z2 < 0.d0) theta = 360.d0*convert - theta
+ if(z2 < 0.d0) theta = 360.d0*convert - theta
thetaup = theta - angle*convert
thetadown = theta + angle*convert
@@ -750,13 +748,13 @@
indice = 1
first = .false.
do ii=1,longueur-1
- if(ch1(ii) /= ' '.or.first) then
- if(ch1(ii) /= ' '.or.ch1(ii+1) /= ' ') then
- ch2(indice) = ch1(ii)
- indice = indice + 1
- first = .true.
+ if(ch1(ii) /= ' ' .or. first) then
+ if(ch1(ii) /= ' ' .or. ch1(ii+1) /= ' ') then
+ ch2(indice) = ch1(ii)
+ indice = indice + 1
+ first = .true.
+ endif
endif
- endif
enddo
ch2(indice) = ch1(longueur)
write(24,200) (ch2(ii),ii=1,indice)
@@ -781,17 +779,17 @@
d = dsqrt(x2**2 + z2**2)
! ignorer si vecteur trop petit
- if (d > cutvect*sizemax) then
+ if (d > cutvect*sizemax) then
d1 = d * rapport
d2 = d1 * dcos(angle*convert)
dummy = x2/d
- if (dummy > 0.9999d0) dummy = 0.9999d0
- if (dummy < -0.9999d0) dummy = -0.9999d0
+ if (dummy > 0.9999d0) dummy = 0.9999d0
+ if (dummy < -0.9999d0) dummy = -0.9999d0
theta = dacos(dummy)
- if(z2 < 0.d0) theta = 360.d0*convert - theta
+ if(z2 < 0.d0) theta = 360.d0*convert - theta
thetaup = theta - angle*convert
thetadown = theta + angle*convert
@@ -815,13 +813,13 @@
indice = 1
first = .false.
do ii=1,longueur-1
- if(ch1(ii) /= ' '.or.first) then
- if(ch1(ii) /= ' '.or.ch1(ii+1) /= ' ') then
- ch2(indice) = ch1(ii)
- indice = indice + 1
- first = .true.
+ if(ch1(ii) /= ' ' .or. first) then
+ if(ch1(ii) /= ' ' .or. ch1(ii+1) /= ' ') then
+ ch2(indice) = ch1(ii)
+ indice = indice + 1
+ first = .true.
+ endif
endif
- endif
enddo
ch2(indice) = ch1(longueur)
write(24,200) (ch2(ii),ii=1,indice)
@@ -863,8 +861,8 @@
!---- write position of the receivers
!
do i=1,nrec
- if(i == 1) write(24,*) '% debut ligne recepteurs'
- if(i == nrec) write(24,*) '% fin ligne recepteurs'
+ if(i == 1) write(24,*) '% debut ligne recepteurs'
+ if(i == nrec) write(24,*) '% fin ligne recepteurs'
iglobrec = nint(posrec(1,i))
xw = coord(1,iglobrec)
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/positsource.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/positsource.f90 2004-05-16 01:48:11 UTC (rev 8418)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/positsource.f90 2007-12-07 23:44:44 UTC (rev 8419)
@@ -24,7 +24,7 @@
integer npoin,nspec
double precision coord(NDIME,npoin)
double precision gltfu(20)
- integer ibool(NGLLX,NGLLY,nspec)
+ integer ibool(NGLLX,NGLLZ,nspec)
double precision dminmax,dmin,xs,zs,xp,zp,dist
integer ip,ipoint,ix,iy,numelem,ilowx,ilowy,ihighx,ihighy
@@ -42,14 +42,14 @@
ilowx = 1
ilowy = 1
ihighx = NGLLX
- ihighy = NGLLY
+ ihighy = NGLLZ
! on ne fait la recherche que sur l'interieur de l'element si source explosive
if(nint(gltfu(2)) == 2) then
ilowx = 2
ilowy = 2
ihighx = NGLLX-1
- ihighy = NGLLY-1
+ ihighy = NGLLZ-1
endif
! recherche du point de grille le plus proche
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/q49shape.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/q49shape.f90 2004-05-16 01:48:11 UTC (rev 8418)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/q49shape.f90 2007-12-07 23:44:44 UTC (rev 8419)
@@ -11,11 +11,11 @@
!
!========================================================================
- subroutine q49shape(shape,dershape,xi,yi,ngnod,NGLLX,NGLLY,NDIME)
+ subroutine q49shape(shape,dershape,xi,yi,ngnod)
!=======================================================================
!
-! "q 4 9 s h a p e" : set up the shape functions and their derivatives
+! set up the shape functions and their derivatives
! for the isoparametric transformation of the
! spectral macroblocs.
! The routine is able to deal with
@@ -38,20 +38,17 @@
implicit none
- integer ngnod,NGLLX,NGLLY,NDIME
+ include "constants.h"
+ integer ngnod
+
double precision shape(ngnod,NGLLX,NGLLX)
double precision dershape(NDIME,ngnod,NGLLX,NGLLX)
- double precision xi(NGLLX),yi(NGLLY)
+ double precision xi(NGLLX),yi(NGLLZ)
- double precision, parameter :: &
- zero=0.d0,one=1.d0,two=2.d0,half=0.5d0,quart=0.25d0
-
integer l1,l2
double precision s,sp,sm,t,tp,tm,s2,t2,ss,tt,st
- double precision, external :: hgll
-
!
!---- set up the shape functions and their local derivatives
!
@@ -59,7 +56,7 @@
!
!---- 4-noded rectangular element
!
- do l2 = 1,NGLLY
+ do l2 = 1,NGLLZ
t = yi(l2)
@@ -97,7 +94,7 @@
!
!---- 9-noded rectangular element
!
- do l2 = 1,NGLLY
+ do l2 = 1,NGLLZ
t = yi(l2)
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/q49spec.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/q49spec.f90 2004-05-16 01:48:11 UTC (rev 8418)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/q49spec.f90 2007-12-07 23:44:44 UTC (rev 8419)
@@ -11,17 +11,16 @@
!
!========================================================================
- subroutine q49spec(shapeint,dershape,dvolu,xjaci,xi, &
- coorg,knods,ngnod,NGLLX,NGLLY,NDIME,nspec,npgeo, &
+ subroutine q49spec(shapeint,dershape,xix,xiz,gammax,gammaz,jacobian,xi, &
+ coorg,knods,ngnod,nspec,npgeo, &
xirec,etarec,flagrange,iptsdisp)
!=======================================================================
!
-! "q 4 9 s p e c" : set up the jacobian matrix
+! set up the jacobian matrix
! for the isoparametric transformation of the
! spectral macroblocs.
-! The routine is able to deal with
-! 4 or 9 control nodes for each bloc.
+! The routine can handle 4 or 9 control nodes.
! The control nodes are defined as follows:
!
! 4 . . . . 7 . . . . 3
@@ -40,20 +39,23 @@
implicit none
- integer ngnod,NGLLX,NGLLY,NDIME,nspec,npgeo,iptsdisp
+ include "constants.h"
+ integer ngnod,nspec,npgeo,iptsdisp
+
integer knods(ngnod,nspec)
double precision shapeint(ngnod,iptsdisp,iptsdisp)
double precision dershape(NDIME,ngnod,NGLLX,NGLLX)
- double precision dvolu(nspec,NGLLX,NGLLX)
- double precision xjaci(nspec,NDIME,NDIME,NGLLX,NGLLX)
double precision coorg(NDIME,npgeo)
double precision xi(NGLLX)
double precision xirec(iptsdisp),etarec(iptsdisp)
double precision flagrange(NGLLX,iptsdisp)
- double precision, parameter :: &
- zero=0.d0,one=1.d0,two=2.d0,half=0.5d0,quart=0.25d0
+ double precision xix(NGLLX,NGLLZ,nspec)
+ double precision xiz(NGLLX,NGLLZ,nspec)
+ double precision gammax(NGLLX,NGLLZ,nspec)
+ double precision gammaz(NGLLX,NGLLZ,nspec)
+ double precision jacobian(NGLLX,NGLLZ,nspec)
integer l1,l2,ispec,in,nnum,ip1,ip2,i,k
double precision s,sp,sm,t,tp,tm,s2,t2,ss,tt,st
@@ -68,36 +70,36 @@
do ispec = 1,nspec
do ip1 = 1,NGLLX
- do ip2 = 1,NGLLY
+ do ip2 = 1,NGLLZ
- xjac2_11 = zero
- xjac2_21 = zero
- xjac2_12 = zero
- xjac2_22 = zero
+ xjac2_11 = ZERO
+ xjac2_21 = ZERO
+ xjac2_12 = ZERO
+ xjac2_22 = ZERO
- do in = 1,ngnod
+ do in = 1,ngnod
- nnum = knods(in,ispec)
+ nnum = knods(in,ispec)
- xjac2_11 = xjac2_11 + dershape(1,in,ip1,ip2)*coorg(1,nnum)
- xjac2_21 = xjac2_21 + dershape(1,in,ip1,ip2)*coorg(2,nnum)
- xjac2_12 = xjac2_12 + dershape(2,in,ip1,ip2)*coorg(1,nnum)
- xjac2_22 = xjac2_22 + dershape(2,in,ip1,ip2)*coorg(2,nnum)
+ xjac2_11 = xjac2_11 + dershape(1,in,ip1,ip2)*coorg(1,nnum)
+ xjac2_21 = xjac2_21 + dershape(1,in,ip1,ip2)*coorg(2,nnum)
+ xjac2_12 = xjac2_12 + dershape(2,in,ip1,ip2)*coorg(1,nnum)
+ xjac2_22 = xjac2_22 + dershape(2,in,ip1,ip2)*coorg(2,nnum)
- enddo
+ enddo
!
!---- invert the jacobian matrix
!
- dvolu(ispec,ip1,ip2) = xjac2_11*xjac2_22 - xjac2_12*xjac2_21
+ jacobian(ip1,ip2,ispec) = xjac2_11*xjac2_22 - xjac2_12*xjac2_21
- if (dvolu(ispec,ip1,ip2) <= zero) stop 'Error : Jacobian undefined !!'
+ if(jacobian(ip1,ip2,ispec) <= zero) stop 'Error: Jacobian undefined!'
- xjaci(ispec,1,1,ip1,ip2) = xjac2_22 / dvolu(ispec,ip1,ip2)
- xjaci(ispec,2,1,ip1,ip2) = - xjac2_21 / dvolu(ispec,ip1,ip2)
- xjaci(ispec,1,2,ip1,ip2) = - xjac2_12 / dvolu(ispec,ip1,ip2)
- xjaci(ispec,2,2,ip1,ip2) = xjac2_11 / dvolu(ispec,ip1,ip2)
+ xix(ip1,ip2,ispec) = xjac2_22 / jacobian(ip1,ip2,ispec)
+ gammax(ip1,ip2,ispec) = - xjac2_21 / jacobian(ip1,ip2,ispec)
+ xiz(ip1,ip2,ispec) = - xjac2_12 / jacobian(ip1,ip2,ispec)
+ gammaz(ip1,ip2,ispec) = xjac2_11 / jacobian(ip1,ip2,ispec)
enddo
enddo
@@ -112,7 +114,7 @@
etarec(i) = 2.d0*dble(i-1)/dble(iptsdisp-1) - 1.d0
enddo
-!---- calcul des interpolateurs de Lagrange (suppose NGLLX = NGLLY)
+!---- calcul des interpolateurs de Lagrange (suppose NGLLX = NGLLZ)
do i=1,NGLLX
do k=1,iptsdisp
flagrange(i,k) = hgll(i-1,xirec(k),xi,NGLLX)
Deleted: seismo/2D/SPECFEM2D/trunk/SPECFEM90/qmasspec.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/qmasspec.f90 2004-05-16 01:48:11 UTC (rev 8418)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/qmasspec.f90 2007-12-07 23:44:44 UTC (rev 8419)
@@ -1,60 +0,0 @@
-
-!========================================================================
-!
-! S P E C F E M 2 D Version 5.0
-! ------------------------------
-!
-! Dimitri Komatitsch
-! Universite de Pau et des Pays de l'Adour, France
-!
-! (c) May 2004
-!
-!========================================================================
-
- subroutine qmasspec(rhoext,wxgll,wygll,ibool,dvolu,rmass,density,kmato,npoin,ireadmodel,nspec,numat)
-
-! build the mass matrix
-
- implicit none
-
- include "constants.h"
-
- integer npoin,nspec,numat
-
- double precision wxgll(NGLLX),wygll(NGLLY),rmass(npoin),dvolu(nspec,NGLLX,NGLLX),density(numat)
- double precision rhoext(npoin)
-
- integer kmato(nspec),ibool(NGLLX,NGLLX,nspec)
-
- integer numelem,material,i,j,iglobnum
- logical ireadmodel
- double precision denst
-
-!
-!---- compute the mass matrix by summing the contribution of each point
-!
-
- rmass = zero
-
- do numelem = 1,nspec
-
- material = kmato(numelem)
- denst = density(material)
-
- do i=1,NGLLX
- do j=1,NGLLY
-
- iglobnum = ibool(i,j,numelem)
-
-!--- si formulation heterogene pour un modele de densite externe
- if(ireadmodel) denst = rhoext(iglobnum)
-
- rmass(iglobnum) = rmass(iglobnum) + denst * wxgll(i) * wygll(j) * dvolu(numelem,i,j)
-
- enddo
- enddo
-
- enddo
-
- end subroutine qmasspec
-
Deleted: seismo/2D/SPECFEM2D/trunk/SPECFEM90/qsumspec.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/qsumspec.f90 2004-05-16 01:48:11 UTC (rev 8418)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/qsumspec.f90 2007-12-07 23:44:44 UTC (rev 8419)
@@ -1,198 +0,0 @@
-
-!========================================================================
-!
-! S P E C F E M 2 D Version 5.0
-! ------------------------------
-!
-! Dimitri Komatitsch
-! Universite de Pau et des Pays de l'Adour, France
-!
-! (c) May 2004
-!
-!========================================================================
-
- subroutine qsumspec(hprime,hTprime, &
- a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13x,a13z, &
- ibool,displ,veloc,accel,Uxnewloc,Uznewloc, &
- rmass,npoin,nspec,gltfu,initialfield, &
- is_bordabs,nelemabs,anyabs,time)
-
- implicit none
-
- include "constants.h"
-
- integer npoin,nspec,nelemabs
- logical anyabs
-
- double precision hprime(NGLLX,NGLLX),hTprime(NGLLX,NGLLX)
- double precision a1(NGLLX,NGLLX,nspec),a2(NGLLX,NGLLX,nspec), &
- a3(NGLLX,NGLLX,nspec),a4(NGLLX,NGLLX,nspec),a5(NGLLX,NGLLX,nspec), &
- a6(NGLLX,NGLLX,nspec),a7(NGLLX,NGLLX,nspec), &
- a8(NGLLX,NGLLX,nspec),a9(NGLLX,NGLLX,nspec),a10(NGLLX,NGLLX,nspec)
- double precision a11(NGLLX,NGLLX),a12(NGLLX,NGLLX)
- double precision a13x(NGLLX,NGLLX,nelemabs)
- double precision a13z(NGLLX,NGLLX,nelemabs)
- double precision Uxnewloc(NGLLX,NGLLX,nspec)
- double precision Uznewloc(NGLLX,NGLLX,nspec)
-
- integer is_bordabs(nspec)
-
-! local arrays
- double precision Uxoldloc(NGLLX,NGLLX)
- double precision Uzoldloc(NGLLX,NGLLX)
- double precision t1(NGLLX,NGLLX)
- double precision t2(NGLLX,NGLLX)
- double precision t3(NGLLX,NGLLX)
- double precision t4(NGLLX,NGLLX)
-
- double precision dUx_dxi,dUz_dxi,dUx_deta,dUz_deta
- double precision hprimex,hTprimex,hprimez,hTprimez
-
- integer ibool(NGLLX,NGLLX,nspec)
-
- double precision rmass(npoin)
- double precision displ(NDIME,npoin),veloc(NDIME,npoin),accel(NDIME,npoin)
-
- double precision gltfu(20)
-
- integer i,j,k,l,ielems,iglobsource,iglobnum,numer_abs
- double precision ricker,time
- logical initialfield
-
- double precision f0,t0,factor,a,angleforce
-
-! main loop on all the spectral elements
- do k=1,nspec
-
-! map the global displacement field to the local mesh
- do j=1,NGLLX
- do i=1,NGLLX
- iglobnum = ibool(i,j,k)
- Uxoldloc(i,j) = displ(1,iglobnum)
- Uzoldloc(i,j) = displ(2,iglobnum)
- enddo
- enddo
-
- do j=1,NGLLX
- do i=1,NGLLX
-
-! compute the gradient of the displacement field (matrix products)
- dUx_dxi = zero
- dUz_dxi = zero
- dUx_deta = zero
- dUz_deta = zero
-
- do l=1,NGLLX
-
- hTprimex = hTprime(i,l)
- hprimez = hprime(l,j)
-
- dUx_dxi = dUx_dxi + hTprimex*Uxoldloc(l,j)
- dUz_dxi = dUz_dxi + hTprimex*Uzoldloc(l,j)
- dUx_deta = dUx_deta + Uxoldloc(i,l)*hprimez
- dUz_deta = dUz_deta + Uzoldloc(i,l)*hprimez
-
- enddo
-
-! compute the local arrays using the components of the stiffness matrix
- t1(i,j) = a1(i,j,k)*dUx_dxi + a2(i,j,k)*dUx_deta + &
- a3(i,j,k)*dUz_dxi + a4(i,j,k)*dUz_deta
- t2(i,j) = a2(i,j,k)*dUx_dxi + a6(i,j,k)*dUx_deta + &
- a7(i,j,k)*dUz_dxi + a8(i,j,k)*dUz_deta
- t3(i,j)= a3(i,j,k)*dUx_dxi + a7(i,j,k)*dUx_deta + &
- a9(i,j,k)*dUz_dxi + a10(i,j,k)*dUz_deta
- t4(i,j)= a4(i,j,k)*dUx_dxi + a8(i,j,k)*dUx_deta + &
- a10(i,j,k)*dUz_dxi + a5(i,j,k)*dUz_deta
-
- enddo
- enddo
-
-! compute the local forces (sum of two matrix products)
- do j=1,NGLLX
- do i=1,NGLLX
-
- Uxnewloc(i,j,k) = zero
- Uznewloc(i,j,k) = zero
-
- do l=1,NGLLX
- hprimex = hprime(i,l)
- hTprimez = hTprime(l,j)
-
- Uxnewloc(i,j,k) = Uxnewloc(i,j,k) + hprimex*t1(l,j) + t2(i,l)*hTprimez
- Uznewloc(i,j,k) = Uznewloc(i,j,k) + hprimex*t3(l,j) + t4(i,l)*hTprimez
-
- enddo
-
- enddo
- enddo
-
-! conditions absorbantes nouvelle formulation
-! pas de dependance par l'adressage indirect
-! car chaque element absorbant est mappe sur un element spectral different
- if(anyabs) then
- numer_abs = is_bordabs(k)
- if(numer_abs .gt. 0) then
- do j=1,NGLLX
- do i=1,NGLLX
- if(a13x(i,j,numer_abs) .ne. zero) then
- iglobnum = ibool(i,j,k)
- Uxnewloc(i,j,k) = Uxnewloc(i,j,k) - a13x(i,j,numer_abs)*veloc(1,iglobnum)
- Uznewloc(i,j,k) = Uznewloc(i,j,k) - a13z(i,j,numer_abs)*veloc(2,iglobnum)
- endif
- enddo
- enddo
- endif
- endif
-
-! assemblage des contributions des differents elements
- do j=1,NGLLX
- do i=1,NGLLX
- iglobnum = ibool(i,j,k)
- accel(1,iglobnum) = accel(1,iglobnum) + Uxnewloc(i,j,k)
- accel(2,iglobnum) = accel(2,iglobnum) + Uznewloc(i,j,k)
- enddo
- enddo
-
- enddo
-
-! --- add the source
- if(.not. initialfield) then
-
- f0 = gltfu(5)
- t0 = gltfu(6)
- factor = gltfu(7)
- angleforce = gltfu(8)
-
-! Ricker wavelet for the source time function
- a = pi*pi*f0*f0
- ricker = - factor * (1.d0-2.d0*a*(time-t0)**2)*exp(-a*(time-t0)**2)
-
-! --- collocated force
- if(nint(gltfu(2)) == 1) then
- iglobsource = nint(gltfu(9))
- accel(1,iglobsource) = accel(1,iglobsource) - dsin(angleforce)*ricker
- accel(2,iglobsource) = accel(2,iglobsource) + dcos(angleforce)*ricker
-
-!---- explosion
- else if(nint(gltfu(2)) == 2) then
-! recuperer numero d'element de la source
- ielems = nint(gltfu(12))
- do i=1,NGLLX
- do j=1,NGLLX
- iglobnum = ibool(i,j,ielems)
- accel(1,iglobnum) = accel(1,iglobnum) + a11(i,j)*ricker
- accel(2,iglobnum) = accel(2,iglobnum) + a12(i,j)*ricker
- enddo
- enddo
- endif
-
- else
- stop 'wrong source type'
- endif
-
-! --- multiplier par l'inverse de la matrice de masse
- accel(1,:) = accel(1,:)*rmass(:)
- accel(2,:) = accel(2,:)*rmass(:)
-
- end subroutine qsumspec
-
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/specfem2D.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/specfem2D.f90 2004-05-16 01:48:11 UTC (rev 8418)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/specfem2D.f90 2007-12-07 23:44:44 UTC (rev 8419)
@@ -17,12 +17,14 @@
!
!========================================================================
-! cleaned version 5.0 is based on SPECFEM2D version 4.2
+! version 5.0 : - got rid of useless routines, suppressed commons etc.
+! - weak formulation based explicitly on stress tensor
+! - implementation of full anisotropy
+
+! version 5.0 is based on SPECFEM2D version 4.2
! (c) June 1998 by Dimitri Komatitsch, Harvard University, USA
! and Jean-Pierre Vilotte, Institut de Physique du Globe de Paris, France
-! version 5.0 : got rid of useless routines, suppressed commons etc.
-
program specfem2D
implicit none
@@ -42,43 +44,69 @@
logical anyabs
integer i,j,it,irec,iglobrec,ipoin,ip,id,ip1,ip2,in,nnum
- integer nbpoin,inump,n,npoinext,netyp,ispec,npoin,npgeo
+ integer nbpoin,inump,n,npoinext,netyp,ispec,npoin,npgeo,iglob
double precision valux,valuz,rhoextread,vpextread,vsextread
+ double precision cpl,csl,rhol
double precision dcosrot,dsinrot,xcor,zcor
! coefficients of the explicit Newmark time scheme
integer NSTEP
- double precision deltatover2,deltatsqover2,time,deltat
+ double precision deltatover2,deltatsquareover2,time,deltat
- double precision, dimension(:), allocatable :: xigll,yigll,wxgll,wygll,xirec,etarec
+! Gauss-Lobatto-Legendre points and weights
+ double precision, dimension(NGLLX) :: xigll,wxgll
+ double precision, dimension(NGLLZ) :: yigll,wzgll
+! array with derivatives of Lagrange polynomials
+ double precision, dimension(NGLLX,NGLLX) :: hprime_xx
+ double precision, dimension(NGLLZ,NGLLZ) :: hprime_zz
+
+! space derivatives
+ double precision tempx1l,tempx2l,tempz1l,tempz2l
+ double precision fac1,fac2,hp1,hp2
+ double precision duxdxl,duzdxl,duxdzl,duzdzl
+ double precision sigma_xx,sigma_xz,sigma_zx,sigma_zz
+
+ double precision, dimension(NGLLX,NGLLZ) :: tempx1,tempx2,tempz1,tempz2
+
+! for anisotropy
+ double precision duydyl,duydzl,duzdyl,duxdyl,duydxl
+ double precision duxdxl_plus_duydyl,duxdxl_plus_duzdzl,duydyl_plus_duzdzl
+ double precision duxdyl_plus_duydxl,duzdxl_plus_duxdzl,duzdyl_plus_duydzl
+
+! Jacobian matrix and determinant
+ double precision xixl,xizl,gammaxl,gammazl,jacobianl
+
+! material properties of the elastic medium
+ double precision mul,lambdal,lambdalplus2mul
+
+ double precision, dimension(:), allocatable :: xirec,etarec
+
double precision, dimension(:,:), allocatable :: coord,accel,veloc,displ, &
- hprime,hTprime,flagrange,xinterp,zinterp,Uxinterp,Uzinterp,elastcoef
+ flagrange,xinterp,zinterp,Uxinterp,Uzinterp,elastcoef
double precision, dimension(:), allocatable :: rmass, &
fglobx,fglobz,density,vpext,vsext,rhoext,displread,velocread,accelread
- double precision, dimension(:,:,:), allocatable :: shapeint,shape,dvolu, &
- a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a13x,a13z,Uxnewloc,Uznewloc
+ double precision, dimension(:,:,:), allocatable :: shapeint,shape, &
+ xix,xiz,gammax,gammaz,jacobian,a13x,a13z
double precision, dimension(:,:), allocatable :: a11,a12
double precision, dimension(:,:,:,:), allocatable :: dershape
- double precision, dimension(:,:,:,:,:), allocatable :: xjaci
-
integer, dimension(:,:,:), allocatable :: ibool
- integer, dimension(:,:), allocatable :: knods,codeabs
- integer, dimension(:), allocatable :: kmato,numabs,is_bordabs
+ integer, dimension(:,:), allocatable :: knods
+ integer, dimension(:), allocatable :: kmato,numabs
- integer ie,k
- integer inum,itourne,ntourne,idummy,numabsread
- integer isource,iexplo
- integer codeabsread(4)
+ integer ie,k,isource,iexplo
+ integer ielems,iglobsource
+ double precision f0,t0,factor,a,angleforce,ricker
+
double precision rsizemin,rsizemax,cpoverdxmin,cpoverdxmax, &
- rlamdaSmin,rlamdaSmax,rlamdaPmin,rlamdaPmax,vpmin,vpmax
+ lambdalSmin,lambdalSmax,lambdalPmin,lambdalPmax,vpmin,vpmax
integer icolor,inumber,isubsamp,ivecttype,itaff,nrec,isismostype
integer numat,ngnod,nspec,iptsdisp,nelemabs
@@ -88,10 +116,16 @@
double precision cutvect,anglerec
+! for absorbing conditions
+ integer ispecabs,inum,numabsread
+ logical codeabsread(4)
+ double precision nx,nz,vx,vz,vn,rho_vp,rho_vs,tx,tz,weight,xxi,zeta,rKvol
+
+ logical, dimension(:,:), allocatable :: codeabs
+
! title of the plot
character(len=60) stitle
-!
!***********************************************************************
!
! i n i t i a l i z a t i o n p h a s e
@@ -203,7 +237,7 @@
!
isource = nint(gltfu(1))
iexplo = nint(gltfu(2))
- if(isource >= 4.and.isource <= 6.and.iexplo == 1) gltfu(8) = gltfu(8) * pi / 180.d0
+ if(isource >= 4 .and. isource <= 6 .and. iexplo == 1) gltfu(8) = gltfu(8) * pi / 180.d0
!
!---- read receiver locations
@@ -240,45 +274,29 @@
!
!---- allocate arrays
!
+ allocate(shape(ngnod,NGLLX,NGLLZ))
+ allocate(shapeint(ngnod,iptsdisp,iptsdisp))
+ allocate(dershape(NDIME,ngnod,NGLLX,NGLLZ))
+ allocate(xix(NGLLX,NGLLZ,nspec))
+ allocate(xiz(NGLLX,NGLLZ,nspec))
+ allocate(gammax(NGLLX,NGLLZ,nspec))
+ allocate(gammaz(NGLLX,NGLLZ,nspec))
+ allocate(jacobian(NGLLX,NGLLZ,nspec))
+ allocate(a11(NGLLX,NGLLZ))
+ allocate(a12(NGLLX,NGLLZ))
+ allocate(xirec(iptsdisp))
+ allocate(etarec(iptsdisp))
+ allocate(flagrange(NGLLX,iptsdisp))
+ allocate(xinterp(iptsdisp,iptsdisp))
+ allocate(zinterp(iptsdisp,iptsdisp))
+ allocate(Uxinterp(iptsdisp,iptsdisp))
+ allocate(Uzinterp(iptsdisp,iptsdisp))
+ allocate(density(numat))
+ allocate(elastcoef(4,numat))
+ allocate(kmato(nspec))
+ allocate(knods(ngnod,nspec))
+ allocate(ibool(NGLLX,NGLLZ,nspec))
-allocate(shape(ngnod,NGLLX,NGLLY))
-allocate(shapeint(ngnod,iptsdisp,iptsdisp))
-allocate(dershape(NDIME,ngnod,NGLLX,NGLLY))
-allocate(dvolu(nspec,NGLLX,NGLLY))
-allocate(xjaci(nspec,NDIME,NDIME,NGLLX,NGLLY))
-allocate(hprime(NGLLX,NGLLY))
-allocate(hTprime(NGLLX,NGLLY))
-allocate(a1(NGLLX,NGLLY,nspec))
-allocate(a2(NGLLX,NGLLY,nspec))
-allocate(a3(NGLLX,NGLLY,nspec))
-allocate(a4(NGLLX,NGLLY,nspec))
-allocate(a5(NGLLX,NGLLY,nspec))
-allocate(a6(NGLLX,NGLLY,nspec))
-allocate(a7(NGLLX,NGLLY,nspec))
-allocate(a8(NGLLX,NGLLY,nspec))
-allocate(a9(NGLLX,NGLLY,nspec))
-allocate(a10(NGLLX,NGLLY,nspec))
-allocate(a11(NGLLX,NGLLY))
-allocate(a12(NGLLX,NGLLY))
-allocate(xigll(NGLLX))
-allocate(yigll(NGLLY))
-allocate(wxgll(NGLLX))
-allocate(wygll(NGLLY))
-allocate(Uxnewloc(NGLLX,NGLLY,nspec))
-allocate(Uznewloc(NGLLX,NGLLY,nspec))
-allocate(xirec(iptsdisp))
-allocate(etarec(iptsdisp))
-allocate(flagrange(NGLLX,iptsdisp))
-allocate(xinterp(iptsdisp,iptsdisp))
-allocate(zinterp(iptsdisp,iptsdisp))
-allocate(Uxinterp(iptsdisp,iptsdisp))
-allocate(Uzinterp(iptsdisp,iptsdisp))
-allocate(density(numat))
-allocate(elastcoef(4,numat))
-allocate(kmato(nspec))
-allocate(knods(ngnod,nspec))
-allocate(ibool(NGLLX,NGLLY,nspec))
-
! --- allocate arrays for absorbing boundary conditions
if(nelemabs <= 0) then
nelemabs = 1
@@ -286,7 +304,6 @@
else
anyabs = .true.
endif
- allocate(is_bordabs(nspec))
allocate(numabs(nelemabs))
allocate(codeabs(4,nelemabs))
@@ -294,21 +311,12 @@
!---- print element group main parameters
!
write(IOUT,107)
- write(IOUT,207) nspec,ngnod,NGLLX,NGLLY,NGLLX*NGLLY,iptsdisp,numat,nelemabs
+ write(IOUT,207) nspec,ngnod,NGLLX,NGLLZ,NGLLX*NGLLZ,iptsdisp,numat,nelemabs
-!
-!---- set up coordinates of the Gauss-Lobatto-Legendre points
-!
- call zwgljd(xigll,wxgll,NGLLX,GAUSSALPHA,GAUSSBETA)
- call zwgljd(yigll,wygll,NGLLY,GAUSSALPHA,GAUSSBETA)
+! set up Gauss-Lobatto-Legendre derivation matrices
+ call define_derivative_matrices(xigll,yigll,wxgll,wzgll,hprime_xx,hprime_zz)
!
-!---- if nb of points is odd, the middle abscissa is exactly zero
-!
- if(mod(NGLLX,2) /= 0) xigll((NGLLX-1)/2+1) = ZERO
- if(mod(NGLLY,2) /= 0) yigll((NGLLY-1)/2+1) = ZERO
-
-!
!---- read the material properties
!
call gmat01(density,elastcoef,numat)
@@ -326,79 +334,24 @@
!---- read absorbing boundary data
!
if(anyabs) then
- read(IIN ,40) datlin
- do n=1,nelemabs
- read(IIN ,*) inum,numabsread,codeabsread(1), &
- codeabsread(2),codeabsread(3),codeabsread(4)
- if(inum < 1 .or. inum > nelemabs) stop 'Wrong absorbing element number'
- numabs(inum) = numabsread
- codeabs(ihaut,inum) = codeabsread(1)
- codeabs(ibas,inum) = codeabsread(2)
- codeabs(igauche,inum) = codeabsread(3)
- codeabs(idroite,inum) = codeabsread(4)
-
-!---- eventuellement tourner element counterclockwise si condition absorbante
-
- if(codeabs(ibas,inum) == iaretebas .or. &
- codeabs(ihaut,inum) == iaretehaut .or. &
- codeabs(igauche,inum) == iaretegauche .or. &
- codeabs(idroite,inum) == iaretedroite) then
- ntourne = 0
-
- else if(codeabs(ibas,inum) == iaretegauche .or. &
- codeabs(ihaut,inum) == iaretedroite .or. &
- codeabs(igauche,inum) == iaretehaut .or. &
- codeabs(idroite,inum) == iaretebas) then
- ntourne = 3
-
- else if(codeabs(ibas,inum) == iaretehaut .or. &
- codeabs(ihaut,inum) == iaretebas .or. &
- codeabs(igauche,inum) == iaretedroite .or. &
- codeabs(idroite,inum) == iaretegauche) then
- ntourne = 2
-
- else if(codeabs(ibas,inum) == iaretedroite .or. &
- codeabs(ihaut,inum) == iaretegauche .or. &
- codeabs(igauche,inum) == iaretebas .or. &
- codeabs(idroite,inum) == iaretehaut) then
- ntourne = 1
- else
- stop 'Error in absorbing conditions numbering'
+ read(IIN ,40) datlin
+ do n=1,nelemabs
+ read(IIN ,*) inum,numabsread,codeabsread(1),codeabsread(2),codeabsread(3),codeabsread(4)
+ if(inum < 1 .or. inum > nelemabs) stop 'Wrong absorbing element number'
+ numabs(inum) = numabsread
+ codeabs(ITOP,inum) = codeabsread(1)
+ codeabs(IBOTTOM,inum) = codeabsread(2)
+ codeabs(ILEFT,inum) = codeabsread(3)
+ codeabs(IRIGHT,inum) = codeabsread(4)
+ enddo
+ write(*,*)
+ write(*,*) 'Number of absorbing elements: ',nelemabs
endif
-!---- rotate element counterclockwise
- if(ntourne /= 0) then
-
- do itourne = 1,ntourne
-
- idummy = knods(1,numabs(inum))
- knods(1,numabs(inum)) = knods(2,numabs(inum))
- knods(2,numabs(inum)) = knods(3,numabs(inum))
- knods(3,numabs(inum)) = knods(4,numabs(inum))
- knods(4,numabs(inum)) = idummy
-
- if(ngnod == 9) then
- idummy = knods(5,numabs(inum))
- knods(5,numabs(inum)) = knods(6,numabs(inum))
- knods(6,numabs(inum)) = knods(7,numabs(inum))
- knods(7,numabs(inum)) = knods(8,numabs(inum))
- knods(8,numabs(inum)) = idummy
- endif
-
- enddo
-
- endif
-
- enddo
- write(*,*)
- write(*,*) 'Number of absorbing elements: ',nelemabs
- endif
-
-
!
!---- compute the spectral element shape functions and their local derivatives
!
- call q49shape(shape,dershape,xigll,yigll,ngnod,NGLLX,NGLLY,NDIME)
+ call q49shape(shape,dershape,xigll,yigll,ngnod)
!
!---- generate the global numbering
@@ -415,8 +368,8 @@
!---- compute the spectral element jacobian matrix
!
- call q49spec(shapeint,dershape,dvolu,xjaci,xigll,coorg,knods,ngnod, &
- NGLLX,NGLLY,NDIME,nspec,npgeo,xirec,etarec,flagrange,iptsdisp)
+ call q49spec(shapeint,dershape,xix,xiz,gammax,gammaz,jacobian,xigll, &
+ coorg,knods,ngnod,nspec,npgeo,xirec,etarec,flagrange,iptsdisp)
!
!---- close input file
@@ -428,10 +381,13 @@
!
allocate(coord(NDIME,npoin))
+
allocate(accel(NDIME,npoin))
allocate(displ(NDIME,npoin))
allocate(veloc(NDIME,npoin))
+
allocate(rmass(npoin))
+
allocate(fglobx(npoin))
allocate(fglobz(npoin))
@@ -444,15 +400,15 @@
allocate(vsext(npoinext))
allocate(rhoext(npoinext))
- allocate(a13x(NGLLX,NGLLY,nelemabs))
- allocate(a13z(NGLLX,NGLLY,nelemabs))
+ allocate(a13x(NGLLX,NGLLZ,nelemabs))
+ allocate(a13z(NGLLX,NGLLZ,nelemabs))
!
!---- set the coordinates of the points of the global grid
!
do ispec = 1,nspec
do ip1 = 1,NGLLX
- do ip2 = 1,NGLLY
+ do ip2 = 1,NGLLZ
xcor = zero
zcor = zero
@@ -492,8 +448,8 @@
!
!---- define coefficients of the Newmark time scheme
!
- deltatover2 = 0.5d0*deltat
- deltatsqover2 = deltat*deltat/2.d0
+ deltatover2 = HALF*deltat
+ deltatsquareover2 = HALF*deltat*deltat
!
!---- definir la position reelle des points source et recepteurs
@@ -522,144 +478,502 @@
endif
!
-!---- build the mass matrix for spectral elements
+!---- define all arrays
!
- call qmasspec(rhoext,wxgll,wygll,ibool,dvolu,rmass,density,kmato,npoin,ireadmodel,nspec,numat)
-
-!
-!---- definir les tableaux a1 a a13
-!
call defarrays(vpext,vsext,rhoext,density,elastcoef, &
- xigll,yigll,wxgll,wygll,hprime,hTprime, &
- a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13x,a13z, &
- ibool,kmato,dvolu,xjaci,coord,gltfu, &
- numabs,codeabs,anyabs,npoin,rsizemin,rsizemax, &
- cpoverdxmin,cpoverdxmax,rlamdaSmin,rlamdaSmax, &
- rlamdaPmin,rlamdaPmax,vpmin,vpmax,ireadmodel,nelemabs,nspec,numat)
+ xigll,yigll,xix,xiz,gammax,gammaz,a11,a12, &
+ ibool,kmato,coord,gltfu,npoin,rsizemin,rsizemax, &
+ cpoverdxmin,cpoverdxmax,lambdalSmin,lambdalSmax, &
+ lambdalPmin,lambdalPmax,vpmin,vpmax,ireadmodel,nspec,numat)
-! initialiser les tableaux a zero
- accel = zero
- veloc = zero
- displ = zero
-
-!
-!--- precalculer l'inverse de la matrice de masse pour efficacite
-!
- rmass(:) = one / rmass(:)
-
-! calculer la numerotation inverse pour les bords absorbants
- is_bordabs(:) = 0
- if(anyabs) then
- do ispec = 1,nelemabs
- is_bordabs(numabs(ispec)) = ispec
+! build the global mass matrix once and for all
+ rmass(:) = ZERO
+ do ispec = 1,nspec
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
+ iglob = ibool(i,j,ispec)
+!--- if external density model
+ if(ireadmodel) then
+ rhol = rhoext(iglob)
+ else
+ rhol = density(kmato(ispec))
+ endif
+ rmass(iglob) = rmass(iglob) + wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec)
+ enddo
enddo
- endif
+ enddo
! convertir angle recepteurs en radians
anglerec = anglerec * pi / 180.d0
!
-!---- eventuellement lecture des champs initiaux dans un fichier
-!
- if(initialfield) then
- print *
- print *,'Reading initial fields from external file...'
- print *
- open(unit=55,file='wavefields.txt',status='unknown')
- read(55,*) nbpoin
- if(nbpoin /= npoin) stop 'Wrong number of points in input file'
- allocate(displread(NDIME))
- allocate(velocread(NDIME))
- allocate(accelread(NDIME))
- do n = 1,npoin
- read(55,*) inump, (displread(i), i=1,NDIME), &
- (velocread(i), i=1,NDIME), (accelread(i), i=1,NDIME)
- if(inump<1 .or. inump>npoin) stop 'Wrong point number'
- displ(:,inump) = displread
- veloc(:,inump) = velocread
- accel(:,inump) = accelread
- enddo
- deallocate(displread)
- deallocate(velocread)
- deallocate(accelread)
- close(55)
- endif
-
-!
-!---- afficher le max du deplacement initial
-!
- print *,'Max norme U initial = ',maxval(sqrt(displ(1,:)**2 + displ(2,:)**2))
-
-!
!---- verifier le maillage, la stabilite et le nb de points par lambda
!
call checkgrid(deltat,gltfu,initialfield,rsizemin,rsizemax, &
- cpoverdxmin,cpoverdxmax,rlamdaSmin,rlamdaSmax,rlamdaPmin,rlamdaPmax)
+ cpoverdxmin,cpoverdxmax,lambdalSmin,lambdalSmax,lambdalPmin,lambdalPmax)
!
!---- initialiser sismogrammes
!
- sisux = zero
- sisuz = zero
+ sisux = ZERO
+ sisuz = ZERO
dcosrot = dcos(anglerec)
dsinrot = dsin(anglerec)
+! initialiser les tableaux a zero
+ accel = ZERO
+ veloc = ZERO
+ displ = ZERO
+
!
+!---- eventuellement lecture des champs initiaux dans un fichier
+!
+ if(initialfield) then
+ print *
+ print *,'Reading initial fields from external file...'
+ print *
+ open(unit=55,file='wavefields.txt',status='unknown')
+ read(55,*) nbpoin
+ if(nbpoin /= npoin) stop 'Wrong number of points in input file'
+ allocate(displread(NDIME))
+ allocate(velocread(NDIME))
+ allocate(accelread(NDIME))
+ do n = 1,npoin
+ read(55,*) inump, (displread(i), i=1,NDIME), &
+ (velocread(i), i=1,NDIME), (accelread(i), i=1,NDIME)
+ if(inump<1 .or. inump>npoin) stop 'Wrong point number'
+ displ(:,inump) = displread
+ veloc(:,inump) = velocread
+ accel(:,inump) = accelread
+ enddo
+ deallocate(displread)
+ deallocate(velocread)
+ deallocate(accelread)
+ close(55)
+ print *,'Max norm of initial displacement = ',maxval(sqrt(displ(1,:)**2 + displ(2,:)**2))
+ endif
+
+!
!---- s t a r t t i m e i t e r a t i o n s
!
write(IOUT,400)
! boucle principale d'evolution en temps
- do it=1,NSTEP
+ do it = 1,NSTEP
! compute current time
- time = (it-1)*deltat
+ time = (it-1)*deltat
- if(mod(it-1,itaff) == 0) then
- if(time >= 1.d-3) then
- write(IOUT,100) it,time
- else
- write(IOUT,101) it,time
+ if(mod(it,itaff) == 0) then
+ write(IOUT,*)
+ if(time >= 1.d-3) then
+ write(IOUT,100) it,time
+ else
+ write(IOUT,101) it,time
+ endif
endif
+
+! `predictor' update displacement using finite-difference time scheme (Newmark)
+ displ(:,:) = displ(:,:) + deltat*veloc(:,:) + deltatsquareover2*accel(:,:)
+ veloc(:,:) = veloc(:,:) + deltatover2*accel(:,:)
+ accel(:,:) = ZERO
+
+! integration over spectral elements
+ do ispec = 1,NSPEC
+
+! get elastic parameters of current spectral element
+ lambdal = elastcoef(1,kmato(ispec))
+ mul = elastcoef(2,kmato(ispec))
+ lambdalplus2mul = lambdal + 2.d0*mul
+
+! first double loop over GLL to compute and store gradients
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
+
+!--- if external medium, get elastic parameters of current grid point
+ if(ireadmodel) then
+ iglob = ibool(i,j,ispec)
+ cpl = vpext(iglob)
+ csl = vsext(iglob)
+ rhol = rhoext(iglob)
+ mul = rhol*csl*csl
+ lambdal = rhol*cpl*cpl - 2.d0*mul
+ lambdalplus2mul = lambdal + 2.d0*mul
+ endif
+
+! derivative along x
+ tempx1l = ZERO
+ tempz1l = ZERO
+ do k = 1,NGLLX
+ hp1 = hprime_xx(k,i)
+ iglob = ibool(k,j,ispec)
+ tempx1l = tempx1l + displ(1,iglob)*hp1
+ tempz1l = tempz1l + displ(2,iglob)*hp1
+ enddo
+
+! derivative along z
+ tempx2l = ZERO
+ tempz2l = ZERO
+ do k = 1,NGLLZ
+ hp2 = hprime_zz(k,j)
+ iglob = ibool(i,k,ispec)
+ tempx2l = tempx2l + displ(1,iglob)*hp2
+ tempz2l = tempz2l + displ(2,iglob)*hp2
+ enddo
+
+ xixl = xix(i,j,ispec)
+ xizl = xiz(i,j,ispec)
+ gammaxl = gammax(i,j,ispec)
+ gammazl = gammaz(i,j,ispec)
+
+! derivatives of displacement
+ duxdxl = tempx1l*xixl + tempx2l*gammaxl
+ duxdzl = tempx1l*xizl + tempx2l*gammazl
+
+ duzdxl = tempz1l*xixl + tempz2l*gammaxl
+ duzdzl = tempz1l*xizl + tempz2l*gammazl
+
+! compute stress tensor
+ sigma_xx = lambdalplus2mul*duxdxl + lambdal*duzdzl
+ sigma_xz = mul*(duzdxl + duxdzl)
+ sigma_zz = lambdalplus2mul*duzdzl + lambdal*duxdxl
+
+! full anisotropy
+ if(TURN_ANISOTROPY_ON) then
+
+! implement anisotropy in 2D
+ duydyl = ZERO
+ duydzl = ZERO
+ duzdyl = ZERO
+ duxdyl = ZERO
+ duydxl = ZERO
+
+! precompute some sums
+ duxdxl_plus_duydyl = duxdxl + duydyl
+ duxdxl_plus_duzdzl = duxdxl + duzdzl
+ duydyl_plus_duzdzl = duydyl + duzdzl
+ duxdyl_plus_duydxl = duxdyl + duydxl
+ duzdxl_plus_duxdzl = duzdxl + duxdzl
+ duzdyl_plus_duydzl = duzdyl + duydzl
+
+ sigma_xx = c11val*duxdxl + c16val*duxdyl_plus_duydxl + c12val*duydyl + &
+ c15val*duzdxl_plus_duxdzl + c14val*duzdyl_plus_duydzl + c13val*duzdzl
+
+! sigma_yy = c12val*duxdxl + c26val*duxdyl_plus_duydxl + c22val*duydyl + &
+! c25val*duzdxl_plus_duxdzl + c24val*duzdyl_plus_duydzl + c23val*duzdzl
+
+ sigma_zz = c13val*duxdxl + c36val*duxdyl_plus_duydxl + c23val*duydyl + &
+ c35val*duzdxl_plus_duxdzl + c34val*duzdyl_plus_duydzl + c33val*duzdzl
+
+! sigma_xy = c16val*duxdxl + c66val*duxdyl_plus_duydxl + c26val*duydyl + &
+! c56val*duzdxl_plus_duxdzl + c46val*duzdyl_plus_duydzl + c36val*duzdzl
+
+ sigma_xz = c15val*duxdxl + c56val*duxdyl_plus_duydxl + c25val*duydyl + &
+ c55val*duzdxl_plus_duxdzl + c45val*duzdyl_plus_duydzl + c35val*duzdzl
+
+! sigma_yz = c14val*duxdxl + c46val*duxdyl_plus_duydxl + c24val*duydyl + &
+! c45val*duzdxl_plus_duxdzl + c44val*duzdyl_plus_duydzl + c34val*duzdzl
+
endif
-! calculer le predictor
- displ(:,:) = displ(:,:) + deltat*veloc(:,:) + deltatsqover2*accel(:,:)
- veloc(:,:) = veloc(:,:) + deltatover2*accel(:,:)
- accel(:,:) = zero
+! stress tensor is symmetric
+ sigma_zx = sigma_xz
+ jacobianl = jacobian(i,j,ispec)
+
+! weak formulation term based on stress tensor (non-symmetric form)
+ tempx1(i,j) = jacobianl*(sigma_xx*xixl+sigma_zx*xizl)
+ tempz1(i,j) = jacobianl*(sigma_xz*xixl+sigma_zz*xizl)
+
+ tempx2(i,j) = jacobianl*(sigma_xx*gammaxl+sigma_zx*gammazl)
+ tempz2(i,j) = jacobianl*(sigma_xz*gammaxl+sigma_zz*gammazl)
+
+ enddo
+ enddo
+
!
-!---- calcul du residu d'acceleration pour le corrector
-!---- retourne dans accel le terme Fext - M*A(i,n+1) - K*D(i,n+1)
+! second double-loop over GLL to compute all terms
!
- call qsumspec(hprime,hTprime,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13x,a13z, &
- ibool,displ,veloc,accel,Uxnewloc,Uznewloc,rmass,npoin, &
- nspec,gltfu,initialfield,is_bordabs,nelemabs,anyabs,time)
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
+! along x direction
+ tempx1l = ZERO
+ tempz1l = ZERO
+ do k = 1,NGLLX
+ fac1 = wxgll(k)*hprime_xx(i,k)
+ tempx1l = tempx1l + tempx1(k,j)*fac1
+ tempz1l = tempz1l + tempz1(k,j)*fac1
+ enddo
+
+! along z direction
+ tempx2l = ZERO
+ tempz2l = ZERO
+ do k = 1,NGLLZ
+ fac2 = wzgll(k)*hprime_zz(j,k)
+ tempx2l = tempx2l + tempx2(i,k)*fac2
+ tempz2l = tempz2l + tempz2(i,k)*fac2
+ enddo
+
+ fac1 = wzgll(j)
+ fac2 = wxgll(i)
+
+ iglob = ibool(i,j,ispec)
+ accel(1,iglob) = accel(1,iglob) - (fac1*tempx1l + fac2*tempx2l)
+ accel(2,iglob) = accel(2,iglob) - (fac1*tempz1l + fac2*tempz2l)
+
+ enddo ! second loop over the GLL points
+ enddo
+
+ enddo ! end of loop over all spectral elements
+
!
-!---- mise a jour globale du deplacement par corrector
+!--- absorbing boundaries
!
+ if(anyabs) then
+
+ do ispecabs=1,nelemabs
+
+ ispec = numabs(ispecabs)
+
+! get elastic parameters of current spectral element
+ lambdal = elastcoef(1,kmato(ispec))
+ mul = elastcoef(2,kmato(ispec))
+ rhol = density(kmato(ispec))
+ rKvol = lambdal + 2.d0*mul/3.d0
+ cpl = dsqrt((rKvol + 4.d0*mul/3.d0)/rhol)
+ csl = dsqrt(mul/rhol)
+
+
+!--- left absorbing boundary
+ if(codeabs(ILEFT,ispecabs)) then
+
+ i = 1
+
+ do j=1,NGLLZ
+
+ iglob = ibool(i,j,ispec)
+
+ zeta = xix(i,j,ispec) * jacobian(i,j,ispec)
+
+! external velocity model
+ if(ireadmodel) then
+ cpl = vpext(iglob)
+ csl = vsext(iglob)
+ rhol = rhoext(iglob)
+ endif
+
+ rho_vp = rhol*cpl
+ rho_vs = rhol*csl
+
+ nx = -1.d0
+ nz = 0.d0
+
+ vx = veloc(1,iglob)
+ vz = veloc(2,iglob)
+
+ vn = nx*vx+nz*vz
+
+ tx = rho_vp*vn*nx+rho_vs*(vx-vn*nx)
+ tz = rho_vp*vn*nz+rho_vs*(vz-vn*nz)
+
+ weight = zeta*wzgll(j)
+
+ accel(1,iglob) = accel(1,iglob) - tx*weight
+ accel(2,iglob) = accel(2,iglob) - tz*weight
+
+ enddo
+
+ endif ! end of left absorbing boundary
+
+!--- right absorbing boundary
+ if(codeabs(IRIGHT,ispecabs)) then
+
+ i = NGLLX
+
+ do j=1,NGLLZ
+
+ iglob = ibool(i,j,ispec)
+
+ zeta = xix(i,j,ispec) * jacobian(i,j,ispec)
+
+! external velocity model
+ if(ireadmodel) then
+ cpl = vpext(iglob)
+ csl = vsext(iglob)
+ rhol = rhoext(iglob)
+ endif
+
+ rho_vp = rhol*cpl
+ rho_vs = rhol*csl
+
+ nx = 1.d0
+ nz = 0.d0
+
+ vx = veloc(1,iglob)
+ vz = veloc(2,iglob)
+
+ vn = nx*vx+nz*vz
+
+ tx = rho_vp*vn*nx+rho_vs*(vx-vn*nx)
+ tz = rho_vp*vn*nz+rho_vs*(vz-vn*nz)
+
+ weight = zeta*wzgll(j)
+
+ accel(1,iglob) = accel(1,iglob) - tx*weight
+ accel(2,iglob) = accel(2,iglob) - tz*weight
+
+ enddo
+
+ endif ! end of right absorbing boundary
+
+!--- bottom absorbing boundary
+ if(codeabs(IBOTTOM,ispecabs)) then
+
+ j = 1
+
+ do i=1,NGLLX
+
+ iglob = ibool(i,j,ispec)
+
+ xxi = gammaz(i,j,ispec) * jacobian(i,j,ispec)
+
+! external velocity model
+ if(ireadmodel) then
+ cpl = vpext(iglob)
+ csl = vsext(iglob)
+ rhol = rhoext(iglob)
+ endif
+
+ rho_vp = rhol*cpl
+ rho_vs = rhol*csl
+
+ nx = 0.d0
+ nz = -1.d0
+
+ vx = veloc(1,iglob)
+ vz = veloc(2,iglob)
+
+ vn = nx*vx+nz*vz
+
+ tx = rho_vp*vn*nx+rho_vs*(vx-vn*nx)
+ tz = rho_vp*vn*nz+rho_vs*(vz-vn*nz)
+
+ weight = xxi*wxgll(i)
+
+ accel(1,iglob) = accel(1,iglob) - tx*weight
+ accel(2,iglob) = accel(2,iglob) - tz*weight
+
+ enddo
+
+ endif ! end of bottom absorbing boundary
+
+!--- top absorbing boundary
+ if(codeabs(ITOP,ispecabs)) then
+
+ j = NGLLZ
+
+ do i=1,NGLLX
+
+ iglob = ibool(i,j,ispec)
+
+ xxi = gammaz(i,j,ispec) * jacobian(i,j,ispec)
+
+! external velocity model
+ if(ireadmodel) then
+ cpl = vpext(iglob)
+ csl = vsext(iglob)
+ rhol = rhoext(iglob)
+ endif
+
+ rho_vp = rhol*cpl
+ rho_vs = rhol*csl
+
+ nx = 0.d0
+ nz = 1.d0
+
+ vx = veloc(1,iglob)
+ vz = veloc(2,iglob)
+
+ vn = nx*vx+nz*vz
+
+ tx = rho_vp*vn*nx+rho_vs*(vx-vn*nx)
+ tz = rho_vp*vn*nz+rho_vs*(vz-vn*nz)
+
+ weight = xxi*wxgll(i)
+
+ accel(1,iglob) = accel(1,iglob) - tx*weight
+ accel(2,iglob) = accel(2,iglob) - tz*weight
+
+ enddo
+
+ endif ! end of top absorbing boundary
+
+ enddo
+
+ endif ! end of absorbing boundaries
+
+
+! --- add the source
+ if(.not. initialfield) then
+
+ f0 = gltfu(5)
+ t0 = gltfu(6)
+ factor = gltfu(7)
+ angleforce = gltfu(8)
+
+! Ricker wavelet for the source time function
+ a = pi*pi*f0*f0
+ ricker = - factor * (1.d0-2.d0*a*(time-t0)**2)*exp(-a*(time-t0)**2)
+
+! --- collocated force
+ if(nint(gltfu(2)) == 1) then
+ iglobsource = nint(gltfu(9))
+ accel(1,iglobsource) = accel(1,iglobsource) - dsin(angleforce)*ricker
+ accel(2,iglobsource) = accel(2,iglobsource) + dcos(angleforce)*ricker
+
+!---- explosion
+ else if(nint(gltfu(2)) == 2) then
+! recuperer numero d'element de la source
+ ielems = nint(gltfu(12))
+ do i=1,NGLLX
+ do j=1,NGLLX
+ iglob = ibool(i,j,ielems)
+ accel(1,iglob) = accel(1,iglob) + a11(i,j)*ricker
+ accel(2,iglob) = accel(2,iglob) + a12(i,j)*ricker
+ enddo
+ enddo
+ endif
+
+ else
+ stop 'wrong source type'
+ endif
+
+! divide by the mass matrix
+ accel(1,:) = accel(1,:) / rmass(:)
+ accel(2,:) = accel(2,:) / rmass(:)
+
+! `corrector' update velocity
veloc(:,:) = veloc(:,:) + deltatover2*accel(:,:)
!
!---- display max of norm of displacement
!
- if(mod(it-1,itaff) == 0) &
- print *,'Max norme U = ',maxval(sqrt(displ(1,:)**2 + displ(2,:)**2))
+ if(mod(it,itaff) == 0) print *,'Max norm of displacement = ',maxval(sqrt(displ(1,:)**2 + displ(2,:)**2))
! store the seismograms
do irec=1,nrec
iglobrec = nint(posrec(1,irec))
- if(isismostype == 1) then
+ if(isismostype == 1) then
valux = displ(1,iglobrec)
valuz = displ(2,iglobrec)
- else if(isismostype == 2) then
+ else if(isismostype == 2) then
valux = veloc(1,iglobrec)
valuz = veloc(2,iglobrec)
- else if(isismostype == 3) then
+ else if(isismostype == 3) then
valux = accel(1,iglobrec)
valuz = accel(2,iglobrec)
else
@@ -678,7 +992,7 @@
if(mod(it,itaff) == 0 .or. it == 5 .or. it == NSTEP) then
write(IOUT,*)
- if(time >= 1.d-3) then
+ if(time >= 1.d-3) then
write(IOUT,110) time
else
write(IOUT,111) time
@@ -689,7 +1003,7 @@
!---- affichage postscript
!
write(IOUT,*) 'Dump PostScript'
- if(ivecttype == 1) then
+ if(ivecttype == 1) then
write(IOUT,*) 'drawing displacement field...'
call plotpost(displ,coord,vpext,gltfu,posrec, &
it,deltat,coorg,xinterp,zinterp,shapeint, &
@@ -697,7 +1011,7 @@
numabs,codeabs,anyabs,stitle,npoin,npgeo,vpmin,vpmax,nrec, &
icolor,inumber,isubsamp,ivecttype,interpol,imeshvect,imodelvect, &
iboundvect,ireadmodel,cutvect,nelemabs,numat,iptsdisp,nspec,ngnod)
- else if(ivecttype == 2) then
+ else if(ivecttype == 2) then
write(IOUT,*) 'drawing velocity field...'
call plotpost(veloc,coord,vpext,gltfu,posrec, &
it,deltat,coorg,xinterp,zinterp,shapeint, &
@@ -705,7 +1019,7 @@
numabs,codeabs,anyabs,stitle,npoin,npgeo,vpmin,vpmax,nrec, &
icolor,inumber,isubsamp,ivecttype,interpol,imeshvect,imodelvect, &
iboundvect,ireadmodel,cutvect,nelemabs,numat,iptsdisp,nspec,ngnod)
- else if(ivecttype == 3) then
+ else if(ivecttype == 3) then
write(IOUT,*) 'drawing acceleration field...'
call plotpost(accel,coord,vpext,gltfu,posrec, &
it,deltat,coorg,xinterp,zinterp,shapeint, &
@@ -781,8 +1095,8 @@
'Number of spectral elements . . . . . (nspec) =',i7,/5x, &
'Number of control nodes per element . (ngnod) =',i7,/5x, &
'Number of points in X-direction . . . (NGLLX) =',i7,/5x, &
- 'Number of points in Y-direction . . . (NGLLY) =',i7,/5x, &
- 'Number of points per element. . .(NGLLX*NGLLY) =',i7,/5x, &
+ 'Number of points in Y-direction . . . (NGLLZ) =',i7,/5x, &
+ 'Number of points per element. . .(NGLLX*NGLLZ) =',i7,/5x, &
'Number of points for display . . . .(iptsdisp) =',i7,/5x, &
'Number of element material sets . . . (numat) =',i7,/5x, &
'Number of absorbing elements . . . .(nelemabs) =',i7)
More information about the cig-commits
mailing list