[cig-commits] r8892 - in seismo/2D/SPECFEM2D/trunk: . DATA
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Mon Dec 17 13:56:39 PST 2007
Author: dkomati1
Date: 2007-12-17 13:56:38 -0800 (Mon, 17 Dec 2007)
New Revision: 8892
Modified:
seismo/2D/SPECFEM2D/trunk/DATA/Par_file
seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Abel_Balanche_bathy_source_solid
seismo/2D/SPECFEM2D/trunk/DATA/Par_file_M2_UPPA
seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Ronan_SV_30
seismo/2D/SPECFEM2D/trunk/DATA/Par_file_unstruct
seismo/2D/SPECFEM2D/trunk/compute_Bielak_conditions.f90
seismo/2D/SPECFEM2D/trunk/constants.h
seismo/2D/SPECFEM2D/trunk/create_color_image.f90
seismo/2D/SPECFEM2D/trunk/specfem2D.F90
Log:
cleaned the plane wave code and related Bielak conditions
Modified: seismo/2D/SPECFEM2D/trunk/DATA/Par_file
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file 2007-12-03 13:11:39 UTC (rev 8891)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file 2007-12-17 21:56:38 UTC (rev 8892)
@@ -3,19 +3,6 @@
title = Test for M2 UPPA
interfacesfile = interfaces_M2_UPPA_curved.dat
-# data concerning mesh, when generated using third-party app (more info in README)
-read_external_mesh = .false.
-mesh_file = ./DATA/unstructured_fluide_solide_test/mesh # file containing the mesh
-nodes_coords_file = ./DATA/unstructured_fluide_solide_test/nodes_coords # file containing the nodes coordinates
-materials_file = ./DATA/unstructured_fluide_solide_test/mat # file containing the material number for each element
-free_surface_file = ./DATA/unstructured_fluide_solide_test/surface_free # file containing the free surface
-absorbing_surface_file = ./DATA/unstructured_fluide_solide_test/surface_abs # file containing the absorbing surface
-
-# parameters concerning partitionning
-nproc = 1 # number of processes
-partionning_method = 1 # ascending order = 1, Metis = 2, Scotch = 3
-partitionning_strategy = 01110 #b{strat=m{asc=g,low=g,rat=1.0,type=h,vert=4}} #01110 # options concerning partitionning strategy.
-
# geometry of the model (origin lower-left corner = 0,0) and mesh description
xmin = 0.d0 # abscissa of left side of the model
xmax = 4000.d0 # abscissa of right side of the model
@@ -28,8 +15,7 @@
TURN_ATTENUATION_ON = .false. # turn attenuation on or off for solid medium
# absorbing boundaries parameters
-absorbing_conditions = .true. # absorbing boundary active or not
-absorbbottom = .true.
+absorbbottom = .true. # absorbing boundary active or not
absorbright = .true.
absorbtop = .false.
absorbleft = .true.
@@ -51,15 +37,8 @@
Mxz = 0. # Mxz component (for a moment tensor source only)
factor = 1.d10 # amplification factor
-# constants for attenuation
-N_SLS = 2 # number of standard linear solids for attenuation
-Qp_attenuation = 136.4376068115 # quality factor P for attenuation
-Qs_attenuation = 136.4376068115 # quality factor S for attenuation
-f0_attenuation = 5.196152422706633 # (Hz) relevant only if source is a Dirac or a Heaviside, else it is f0
-
# receiver line parameters for seismograms
seismotype = 1 # record 1=displ 2=veloc 3=accel 4=pressure
-generate_STATIONS = .true. # creates a STATION file in ./DATA
nreceiverlines = 1 # number of receiver lines
anglerec = 0.d0 # angle to rotate components at receivers
@@ -93,7 +72,7 @@
# set vs to zero to make a given model acoustic
# the mesh can contain both acoustic and elastic models simultaneously
1 1 2700.d0 3000.d0 1732.051d0 0 0
-2 1 2500.d0 2700.d0 0.d0 0 0 #1558.89d0 0 0
+2 1 2500.d0 2700.d0 0 0 0 #1558.89d0 0 0
3 1 2200.d0 2500.d0 1443.375d0 0 0
4 1 2200.d0 2200.d0 1343.375d0 0 0
# define the different regions of the model in the (nx,nz) spectral element mesh
Modified: seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Abel_Balanche_bathy_source_solid
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Abel_Balanche_bathy_source_solid 2007-12-03 13:11:39 UTC (rev 8891)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Abel_Balanche_bathy_source_solid 2007-12-17 21:56:38 UTC (rev 8892)
@@ -9,6 +9,7 @@
nx = 134 # number of elements along X
ngnod = 9 # number of control nodes per element (4 or 9)
initialfield = .false. # use a plane wave as source or not
+add_Bielak_conditions = .false. # add Bielak conditions or not if initial plane wave
assign_external_model = .false. # define external earth model or not
TURN_ANISOTROPY_ON = .false. # turn anisotropy on or off for solid medium
TURN_ATTENUATION_ON = .false. # turn attenuation on or off for solid medium
Modified: seismo/2D/SPECFEM2D/trunk/DATA/Par_file_M2_UPPA
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file_M2_UPPA 2007-12-03 13:11:39 UTC (rev 8891)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file_M2_UPPA 2007-12-17 21:56:38 UTC (rev 8892)
@@ -9,6 +9,7 @@
nx = 80 # number of elements along X
ngnod = 9 # number of control nodes per element (4 or 9)
initialfield = .false. # use a plane wave as source or not
+add_Bielak_conditions = .false. # add Bielak conditions or not if initial plane wave
assign_external_model = .false. # define external earth model or not
TURN_ANISOTROPY_ON = .false. # turn anisotropy on or off for solid medium
TURN_ATTENUATION_ON = .false. # turn attenuation on or off for solid medium
Modified: seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Ronan_SV_30
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Ronan_SV_30 2007-12-03 13:11:39 UTC (rev 8891)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Ronan_SV_30 2007-12-17 21:56:38 UTC (rev 8892)
@@ -35,8 +35,8 @@
absorbleft = .true.
# time step parameters
-nt = 20000 # total number of time steps
-deltat = 6.25e-4 # duration of a time step
+nt = 2000 # total number of time steps
+deltat = 1.e-2 # duration of a time step
# source parameters
source_surf = .true. # source inside the medium or at the surface
@@ -72,7 +72,7 @@
enreg_surf = .true. # receivers inside the medium or at the surface
# display parameters
-NTSTEP_BETWEEN_OUTPUT_INFO = 500 # display frequency in time steps
+NTSTEP_BETWEEN_OUTPUT_INFO = 100 # display frequency in time steps
output_postscript_snapshot = .true. # output Postscript snapshot of the results
output_color_image = .true. # output color image of the results
imagetype = 1 # display 1=displ 2=veloc 3=accel 4=pressure
Modified: seismo/2D/SPECFEM2D/trunk/DATA/Par_file_unstruct
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file_unstruct 2007-12-03 13:11:39 UTC (rev 8891)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file_unstruct 2007-12-17 21:56:38 UTC (rev 8892)
@@ -22,6 +22,7 @@
nx = 80 # number of elements along X
ngnod = 4 # number of control nodes per element (4 or 9)
initialfield = .false. # use a plane wave as source or not
+add_Bielak_conditions = .false. # add Bielak conditions or not if initial plane wave
assign_external_model = .false. # define external earth model or not
TURN_ANISOTROPY_ON = .false. # turn anisotropy on or off for solid medium
TURN_ATTENUATION_ON = .false. # turn attenuation on or off for solid medium
@@ -101,4 +102,4 @@
1 80 21 40 4
1 80 41 60 3
60 80 21 40 4
-30 40 50 60 2
\ No newline at end of file
+30 40 50 60 2
Modified: seismo/2D/SPECFEM2D/trunk/compute_Bielak_conditions.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_Bielak_conditions.f90 2007-12-03 13:11:39 UTC (rev 8891)
+++ seismo/2D/SPECFEM2D/trunk/compute_Bielak_conditions.f90 2007-12-17 21:56:38 UTC (rev 8892)
@@ -65,7 +65,6 @@
z = coord(2,iglob)
! times for velocity and traction are staggered i.e. separated by deltat/2.d0
-!! DK DK YYYYYYYYYYYYYYYY UGLY il faudra peut-etre ajouter ou enlever deltat ici
time_veloc = (it-1)*deltat + deltat/2.d0 + time_offset
time_traction = time_veloc + deltat/2.d0
@@ -106,11 +105,11 @@
t = time_veloc
! analytical expression of the two components of the velocity vector
- veloc_horiz = rac3sur2 * ricker_Bielak_veloc(t - x/2.d0 + (9 - z) * rac3sur2) &
- + rac3sur2 * ricker_Bielak_veloc(t - x/2.d0 - (9 - z) * rac3sur2) &
- + rac3 * ricker_Bielak_veloc(t - x/2.d0)
- veloc_vert = - HALF * ricker_Bielak_veloc(t - x/2.d0 + (9 - z) * rac3sur2) &
- + HALF * ricker_Bielak_veloc(t - x/2.d0 - (9 - z) * rac3sur2)
+ veloc_horiz = (sqrt(3.d0)/2.d0) * ricker_Bielak_veloc(t - x/2.d0 + (9 - z) * (sqrt(3.d0)/2.d0)) &
+ + (sqrt(3.d0)/2.d0) * ricker_Bielak_veloc(t - x/2.d0 - (9 - z) * (sqrt(3.d0)/2.d0)) &
+ + sqrt(3.d0) * ricker_Bielak_veloc(t - x/2.d0)
+ veloc_vert = - HALF * ricker_Bielak_veloc(t - x/2.d0 + (9 - z) * (sqrt(3.d0)/2.d0)) &
+ + HALF * ricker_Bielak_veloc(t - x/2.d0 - (9 - z) * (sqrt(3.d0)/2.d0))
end subroutine compute_Bielak_conditions
Modified: seismo/2D/SPECFEM2D/trunk/constants.h
===================================================================
--- seismo/2D/SPECFEM2D/trunk/constants.h 2007-12-03 13:11:39 UTC (rev 8891)
+++ seismo/2D/SPECFEM2D/trunk/constants.h 2007-12-17 21:56:38 UTC (rev 8892)
@@ -87,8 +87,6 @@
double precision, parameter :: time_offset = 3.575d0
double precision, parameter :: f0_ricker_Bielak = 1.d0
double precision, parameter :: a = PI*PI*f0_ricker_Bielak*f0_ricker_Bielak
- double precision, parameter :: rac3 = 1.7320508075688772935d0
- double precision, parameter :: rac3sur2 = rac3 / 2.d0
! non linear display to enhance small amplitudes in color images
double precision, parameter :: POWER_DISPLAY_COLOR = 0.30d0
Modified: seismo/2D/SPECFEM2D/trunk/create_color_image.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/create_color_image.f90 2007-12-03 13:11:39 UTC (rev 8891)
+++ seismo/2D/SPECFEM2D/trunk/create_color_image.f90 2007-12-17 21:56:38 UTC (rev 8892)
@@ -159,10 +159,10 @@
vpmin = min(vpmin,image_color_vp_display(ix,iy))
vpmax = max(vpmax,image_color_vp_display(ix,iy))
endif
-
+
enddo
enddo
-
+
! in the PNM format, the image starts in the upper-left corner
do iy=NY,1,-1
do ix=1,NX
Modified: seismo/2D/SPECFEM2D/trunk/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2007-12-03 13:11:39 UTC (rev 8891)
+++ seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2007-12-17 21:56:38 UTC (rev 8892)
@@ -1471,7 +1471,8 @@
else
! compute analytical initial plane wave field
- print *,'computing analytical initial plane wave field for SV wave at 30 degrees and Poisson = 0.25'
+! the analytical expression below is specific to an SV wave at 30 degrees and Poisson = 0.3333
+ print *,'computing analytical initial plane wave field for SV wave at 30 degrees and Poisson = 0.3333'
do i = 1,npoin
@@ -1482,25 +1483,25 @@
t = 0.d0 + time_offset
! initial analytical displacement
- displ_elastic(1,i) = rac3sur2 * ricker_Bielak_displ(t - x/2.d0 + (9 - z) * rac3sur2) &
- + rac3sur2 * ricker_Bielak_displ(t - x/2.d0 - (9 - z) * rac3sur2) &
- + rac3 * ricker_Bielak_displ(t - x/2.d0)
- displ_elastic(2,i) = - HALF * ricker_Bielak_displ(t - x/2.d0 + (9 - z) * rac3sur2) &
- + HALF * ricker_Bielak_displ(t - x/2.d0 - (9 - z) * rac3sur2)
+ displ_elastic(1,i) = (sqrt(3.d0)/2.d0) * ricker_Bielak_displ(t - x/2.d0 + (9 - z) * (sqrt(3.d0)/2.d0)) &
+ + (sqrt(3.d0)/2.d0) * ricker_Bielak_displ(t - x/2.d0 - (9 - z) * (sqrt(3.d0)/2.d0)) &
+ + sqrt(3.d0) * ricker_Bielak_displ(t - x/2.d0)
+ displ_elastic(2,i) = - HALF * ricker_Bielak_displ(t - x/2.d0 + (9 - z) * (sqrt(3.d0)/2.d0)) &
+ + HALF * ricker_Bielak_displ(t - x/2.d0 - (9 - z) * (sqrt(3.d0)/2.d0))
! initial analytical velocity
- veloc_elastic(1,i) = rac3sur2 * ricker_Bielak_veloc(t - x/2.d0 + (9 - z) * rac3sur2) &
- + rac3sur2 * ricker_Bielak_veloc(t - x/2.d0 - (9 - z) * rac3sur2) &
- + rac3 * ricker_Bielak_veloc(t - x/2.d0)
- veloc_elastic(2,i) = - HALF * ricker_Bielak_veloc(t - x/2.d0 + (9 - z) * rac3sur2) &
- + HALF * ricker_Bielak_veloc(t - x/2.d0 - (9 - z) * rac3sur2)
+ veloc_elastic(1,i) = (sqrt(3.d0)/2.d0) * ricker_Bielak_veloc(t - x/2.d0 + (9 - z) * (sqrt(3.d0)/2.d0)) &
+ + (sqrt(3.d0)/2.d0) * ricker_Bielak_veloc(t - x/2.d0 - (9 - z) * (sqrt(3.d0)/2.d0)) &
+ + sqrt(3.d0) * ricker_Bielak_veloc(t - x/2.d0)
+ veloc_elastic(2,i) = - HALF * ricker_Bielak_veloc(t - x/2.d0 + (9 - z) * (sqrt(3.d0)/2.d0)) &
+ + HALF * ricker_Bielak_veloc(t - x/2.d0 - (9 - z) * (sqrt(3.d0)/2.d0))
! initial analytical acceleration
- accel_elastic(1,i) = rac3sur2 * ricker_Bielak_accel(t - x/2.d0 + (9 - z) * rac3sur2) &
- + rac3sur2 * ricker_Bielak_accel(t - x/2.d0 - (9 - z) * rac3sur2) &
- + rac3 * ricker_Bielak_accel(t - x/2.d0)
- accel_elastic(2,i) = - HALF * ricker_Bielak_accel(t - x/2.d0 + (9 - z) * rac3sur2) &
- + HALF * ricker_Bielak_accel(t - x/2.d0 - (9 - z) * rac3sur2)
+ accel_elastic(1,i) = (sqrt(3.d0)/2.d0) * ricker_Bielak_accel(t - x/2.d0 + (9 - z) * (sqrt(3.d0)/2.d0)) &
+ + (sqrt(3.d0)/2.d0) * ricker_Bielak_accel(t - x/2.d0 - (9 - z) * (sqrt(3.d0)/2.d0)) &
+ + sqrt(3.d0) * ricker_Bielak_accel(t - x/2.d0)
+ accel_elastic(2,i) = - HALF * ricker_Bielak_accel(t - x/2.d0 + (9 - z) * (sqrt(3.d0)/2.d0)) &
+ + HALF * ricker_Bielak_accel(t - x/2.d0 - (9 - z) * (sqrt(3.d0)/2.d0))
enddo
More information about the cig-commits
mailing list