[cig-commits] r8422 - in seismo/2D/SPECFEM2D/trunk: MAILLE90
SPECFEM90
walter at geodynamics.org
walter at geodynamics.org
Fri Dec 7 15:45:00 PST 2007
Author: walter
Date: 2007-12-07 15:44:59 -0800 (Fri, 07 Dec 2007)
New Revision: 8422
Added:
seismo/2D/SPECFEM2D/trunk/MAILLE90/Par_poisson_minus_one
Modified:
seismo/2D/SPECFEM2D/trunk/SPECFEM90/gmat01.f90
seismo/2D/SPECFEM2D/trunk/SPECFEM90/plotpost.f90
seismo/2D/SPECFEM2D/trunk/SPECFEM90/specfem2D.f90
Log:
cleaned the FD and SEM source codes following visit to INSA Rouen.
added ability to handle negative Poisson's ratio in SEM.
added vertical force source in FD code (no_pml).
Added: seismo/2D/SPECFEM2D/trunk/MAILLE90/Par_poisson_minus_one
===================================================================
--- seismo/2D/SPECFEM2D/trunk/MAILLE90/Par_poisson_minus_one 2004-05-21 18:05:12 UTC (rev 8421)
+++ seismo/2D/SPECFEM2D/trunk/MAILLE90/Par_poisson_minus_one 2007-12-07 23:44:59 UTC (rev 8422)
@@ -0,0 +1,79 @@
+# ----------------------------------------------------------------
+#
+# 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 =2900. ! 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 =.true. ! 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 =.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
+ioutputgrid =.false. ! save the grid in a text file or not
+#
+# velocity and density model (nx,nz)
+#
+nbmodels =1 ! nb de modeles differents (rho,vp,vs)
+1 0 2200.d0 2500.d0 2165.d0 0 0
+nbzone =1 ! nb of zones and model number for each zone
+1 80 1 60 1
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/gmat01.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/gmat01.f90 2004-05-21 18:05:12 UTC (rev 8421)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/gmat01.f90 2007-12-07 23:44:59 UTC (rev 8422)
@@ -26,7 +26,7 @@
double precision density(numat),elastcoef(4,numat)
integer in,n,indic
- double precision young,poiss,denst,cp,cs,amu,a2mu,alam
+ double precision young,poisson,denst,cp,cs,amu,a2mu,alam
double precision val1,val2,val3,val4
double precision c11,c13,c33,c44
@@ -55,17 +55,19 @@
Kmod = alam + a2mu
Kvol = alam + a2mu/3.d0
young = 9.d0*Kvol*amu/(3.d0*Kvol + amu)
- poiss = half*(3.d0*Kvol-a2mu)/(3.d0*Kvol+amu)
- if (poiss < 0.d0 .or. poiss >= 0.50001d0) stop 'Poisson''s ratio out of range !'
+ poisson = half*(3.d0*Kvol-a2mu)/(3.d0*Kvol+amu)
+! Poisson's ratio must be between -1 and +1/2
+ if (poisson < -1.d0 .or. poisson > 0.5d0) stop 'Poisson''s ratio out of range'
!---- materiau isotrope, module de Young et coefficient de Poisson donnes
else if(indic == 1) then
young = val1
- poiss = val2
- if (poiss < 0.d0 .or. poiss >= 0.50001d0) stop 'Poisson''s ratio out of range !'
- a2mu = young/(one+poiss)
+ poisson = val2
+! Poisson's ratio must be between -1 and +1/2
+ if (poisson < -1.d0 .or. poisson > 0.5d0) stop 'Poisson''s ratio out of range'
+ a2mu = young/(one+poisson)
amu = half*a2mu
- alam = a2mu*poiss/(one-two*poiss)
+ alam = a2mu*poisson/(one-two*poisson)
Kmod = alam + a2mu
Kvol = alam + a2mu/3.d0
cp = dsqrt((Kvol + 4.d0*amu/3.d0)/denst)
@@ -106,7 +108,7 @@
!---- check the input
!
if(indic == 0 .or. indic == 1) then
- write(iout,200) n,cp,cs,denst,poiss,alam,amu,Kvol,young
+ write(iout,200) n,cp,cs,denst,poisson,alam,amu,Kvol,young
else
write(iout,300) n,c11,c13,c33,c44,denst, &
sqrt(c33/denst),sqrt(c11/denst),sqrt(c44/denst),sqrt(c44/denst)
@@ -128,7 +130,7 @@
'P-wave velocity. . . . . . . . . . . (cp) =',1pe15.8,/5x, &
'S-wave velocity. . . . . . . . . . . (cs) =',1pe15.8,/5x, &
'Mass density. . . . . . . . . . . (denst) =',1pe15.8,/5x, &
- 'Poisson''s ratio . . . . . . . . . (poiss) =',1pe15.8,/5x, &
+ 'Poisson''s ratio . . . . . . . .(poisson) =',1pe15.8,/5x, &
'First Lame parameter Lambda. . . . (alam) =',1pe15.8,/5x, &
'Second Lame parameter Mu. . . . . . (amu) =',1pe15.8,/5x, &
'Bulk modulus K . . . . . . . . . . (Kvol) =',1pe15.8,/5x, &
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/plotpost.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/plotpost.f90 2004-05-21 18:05:12 UTC (rev 8421)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/plotpost.f90 2007-12-07 23:44:59 UTC (rev 8422)
@@ -82,7 +82,7 @@
sizex = 27.94d0
sizez = 21.59d0
else
- usoffset = 0.
+ usoffset = 0.d0
sizex = 29.7d0
sizez = 21.d0
endif
@@ -316,18 +316,18 @@
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
- write(24,*) '(Displacement vector field) show'
- else if (ivecttype == 2) then
- write(24,*) '(Velocity vector field) show'
- else if (ivecttype == 3) then
- write(24,*) '(Acceleration vector field) show'
+ if(ivecttype == 1) then
+ write(24,*) '(Displacement vector field) show'
+ else if(ivecttype == 2) then
+ write(24,*) '(Velocity vector field) show'
+ else if(ivecttype == 3) then
+ write(24,*) '(Acceleration vector field) show'
else
- stop 'Bad field code in PostScript display'
+ stop 'Bad field code in PostScript display'
endif
write(24,*) 'grestore'
write(24,*) '25.35 CM 18.9 CM MV'
@@ -368,13 +368,13 @@
if(ireadmodel) then
x1 = (vpext(ibool(i,j,ispec))-vpmin)/ (vpmax-vpmin)
else
- material = kmato(ispec)
- rlamda = elastcoef(1,material)
- rmu = elastcoef(2,material)
- denst = density(material)
- rKvol = rlamda + 2.d0*rmu/3.d0
- cploc = dsqrt((rKvol + 4.d0*rmu/3.d0)/denst)
- x1 = (cploc-vpmin)/(vpmax-vpmin)
+ material = kmato(ispec)
+ rlamda = elastcoef(1,material)
+ rmu = elastcoef(2,material)
+ denst = density(material)
+ rKvol = rlamda + 2.d0*rmu/3.d0
+ cploc = dsqrt((rKvol + 4.d0*rmu/3.d0)/denst)
+ x1 = (cploc-vpmin)/(vpmax-vpmin)
endif
else
x1 = 0.5d0
@@ -382,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
@@ -393,29 +393,32 @@
zw = (zw-zmin)*rapp_page + orig_z
xw = xw * centim
zw = zw * centim
- write(24,500) xw,zw
+ write(24,500) xw,zw
+
xw = coord(1,ibool(i+isubsamp,j,ispec))
zw = coord(2,ibool(i+isubsamp,j,ispec))
xw = (xw-xmin)*rapp_page + orig_x
zw = (zw-zmin)*rapp_page + orig_z
xw = xw * centim
zw = zw * centim
- write(24,499) xw,zw
+ write(24,499) xw,zw
+
xw = coord(1,ibool(i+isubsamp,j+isubsamp,ispec))
zw = coord(2,ibool(i+isubsamp,j+isubsamp,ispec))
xw = (xw-xmin)*rapp_page + orig_x
zw = (zw-zmin)*rapp_page + orig_z
xw = xw * centim
zw = zw * centim
- write(24,499) xw,zw
+ write(24,499) xw,zw
+
xw = coord(1,ibool(i,j+isubsamp,ispec))
zw = coord(2,ibool(i,j+isubsamp,ispec))
xw = (xw-xmin)*rapp_page + orig_x
zw = (zw-zmin)*rapp_page + orig_z
xw = xw * centim
zw = zw * centim
- write(24,499) xw,zw
- write(24,604) x1
+ write(24,499) xw,zw
+ write(24,604) x1
enddo
enddo
@@ -427,7 +430,7 @@
!---- draw spectral element mesh
!
- if (imeshvect) then
+ if(imeshvect) then
write(24,*) '%'
write(24,*) '% spectral element mesh'
@@ -458,7 +461,7 @@
write(24,*) 'MK'
write(24,681) x1,z1
- if (ngnod == 4) then
+ if(ngnod == 4) then
! tracer des droites si elements 4 noeuds
@@ -497,45 +500,45 @@
! tracer des courbes si elements 9 noeuds
do ir=2,iptsdisp
- x2 = (xinterp(ir,is)-xmin)*rapp_page + orig_x
- z2 = (zinterp(ir,is)-zmin)*rapp_page + orig_z
- x2 = x2 * centim
- z2 = z2 * centim
- write(24,681) x2,z2
+ x2 = (xinterp(ir,is)-xmin)*rapp_page + orig_x
+ z2 = (zinterp(ir,is)-zmin)*rapp_page + orig_z
+ x2 = x2 * centim
+ z2 = z2 * centim
+ write(24,681) x2,z2
enddo
ir=iptsdisp
do is=2,iptsdisp
- x2 = (xinterp(ir,is)-xmin)*rapp_page + orig_x
- z2 = (zinterp(ir,is)-zmin)*rapp_page + orig_z
- x2 = x2 * centim
- z2 = z2 * centim
- write(24,681) x2,z2
+ x2 = (xinterp(ir,is)-xmin)*rapp_page + orig_x
+ z2 = (zinterp(ir,is)-zmin)*rapp_page + orig_z
+ x2 = x2 * centim
+ z2 = z2 * centim
+ write(24,681) x2,z2
enddo
is=iptsdisp
do ir=iptsdisp-1,1,-1
- x2 = (xinterp(ir,is)-xmin)*rapp_page + orig_x
- z2 = (zinterp(ir,is)-zmin)*rapp_page + orig_z
- x2 = x2 * centim
- z2 = z2 * centim
- write(24,681) x2,z2
+ x2 = (xinterp(ir,is)-xmin)*rapp_page + orig_x
+ z2 = (zinterp(ir,is)-zmin)*rapp_page + orig_z
+ x2 = x2 * centim
+ z2 = z2 * centim
+ write(24,681) x2,z2
enddo
ir=1
do is=iptsdisp-1,2,-1
- x2 = (xinterp(ir,is)-xmin)*rapp_page + orig_x
- z2 = (zinterp(ir,is)-zmin)*rapp_page + orig_z
- x2 = x2 * centim
- z2 = z2 * centim
- write(24,681) x2,z2
+ x2 = (xinterp(ir,is)-xmin)*rapp_page + orig_x
+ z2 = (zinterp(ir,is)-zmin)*rapp_page + orig_z
+ x2 = x2 * centim
+ z2 = z2 * centim
+ write(24,681) x2,z2
enddo
endif
write(24,*) 'CO'
- if (icolor == 1) then
+ if(icolor == 1) then
! For the moment 20 different colors max
nbcols = 20
@@ -556,7 +559,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
@@ -566,7 +569,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
@@ -655,7 +658,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
@@ -671,7 +674,7 @@
write(24,*) '0 setgray'
endif
- if (interpol) then
+ if(interpol) then
print *,'Interpolating the vector field...'
@@ -714,14 +717,14 @@
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
@@ -779,14 +782,14 @@
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
@@ -851,7 +854,7 @@
xw = xw * centim
zw = zw * centim
write(24,510) xw,zw
- if (isymbols) then
+ if(isymbols) then
write(24,*) 'Cross'
else
write(24,*) '(S) show'
@@ -873,7 +876,7 @@
xw = xw * centim
zw = zw * centim
write(24,510) xw,zw
- if (isymbols) then
+ if(isymbols) then
if(nrec > ndots .and. i /= 1 .and. i /= nrec) then
write(24,*) 'VDot'
else
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/specfem2D.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/specfem2D.f90 2004-05-21 18:05:12 UTC (rev 8421)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/specfem2D.f90 2007-12-07 23:44:59 UTC (rev 8422)
@@ -117,7 +117,7 @@
double precision cutvect,anglerec
! for absorbing conditions
- integer ispecabs,inum,numabsread
+ integer ispecabs,inum,numabsread,i1abs,i2abs
logical codeabsread(4)
double precision nx,nz,vx,vz,vn,rho_vp,rho_vs,tx,tz,weight,xxi,zeta,rKvol
@@ -834,8 +834,14 @@
j = 1
- do i=1,NGLLX
+! exclude corners to make sure there is no contradiction on the normal
+ i1abs = 1
+ i2abs = NGLLX
+ if(codeabs(ILEFT,ispecabs)) i1abs = 2
+ if(codeabs(IRIGHT,ispecabs)) i2abs = NGLLX-1
+ do i=i1abs,i2abs
+
iglob = ibool(i,j,ispec)
xxi = gammaz(i,j,ispec) * jacobian(i,j,ispec)
@@ -875,8 +881,14 @@
j = NGLLZ
- do i=1,NGLLX
+! exclude corners to make sure there is no contradiction on the normal
+ i1abs = 1
+ i2abs = NGLLX
+ if(codeabs(ILEFT,ispecabs)) i1abs = 2
+ if(codeabs(IRIGHT,ispecabs)) i2abs = NGLLX-1
+ do i=i1abs,i2abs
+
iglob = ibool(i,j,ispec)
xxi = gammaz(i,j,ispec) * jacobian(i,j,ispec)
More information about the cig-commits
mailing list