[cig-commits] r18869 - in seismo/2D/SPECFEM2D/trunk: EXAMPLES/canyon src/meshfem2D src/specfem2D
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Thu Sep 1 06:42:51 PDT 2011
Author: dkomati1
Date: 2011-09-01 06:42:50 -0700 (Thu, 01 Sep 2011)
New Revision: 18869
Modified:
seismo/2D/SPECFEM2D/trunk/EXAMPLES/canyon/Par_file_canyon
seismo/2D/SPECFEM2D/trunk/EXAMPLES/canyon/SOURCE_canyon
seismo/2D/SPECFEM2D/trunk/src/meshfem2D/meshfem2D.F90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/noise_tomography.f90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/prepare_initialfield.F90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
Log:
fixed the plane wave code and made sure the "EXAMPLES/canyon" example now works fine;
also removed useless white spaces in source files.
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/canyon/Par_file_canyon
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/canyon/Par_file_canyon 2011-09-01 09:29:08 UTC (rev 18868)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/canyon/Par_file_canyon 2011-09-01 13:42:50 UTC (rev 18869)
@@ -1,14 +1,15 @@
# title of job
-title = Calcul Mexique Alejandro
+title = Canyon Mexico Paco P-SV
# forward or adjoint simulation
SIMULATION_TYPE = 1 # 1 = forward, 2 = adjoint + kernels
-SAVE_FORWARD = .true. # save the last frame, needed for adjoint simulation
+NOISE_TOMOGRAPHY = 0 # 0 = earthquake simulation, 1/2/3 = noise simulation
+SAVE_FORWARD = .false. # save the last frame, needed for adjoint simulation
# parameters concerning partitioning
nproc = 1 # number of processes
partitioning_method = 3 # SCOTCH = 3, ascending order (very bad idea) = 1
-
+
ngnod = 4 # number of control nodes per element (4 or 9)
initialfield = .true. # use a plane wave as source or not
add_Bielak_conditions = .true. # add Bielak conditions or not if initial plane wave
@@ -40,17 +41,17 @@
rec_normal_to_surface = .false. # base anglerec normal to surface (external mesh and curve file needed)
# first receiver line
-nrec = 20 # number of receivers
-xdeb = 1. # first receiver x in meters
+nrec = 5 # number of receivers
+xdeb = 14.5 # first receiver x in meters
zdeb = 9. # first receiver z in meters
-xfin = 18. # last receiver x in meters (ignored if onlyone receiver)
+xfin = 17.5 # last receiver x in meters (ignored if onlyone receiver)
zfin = 9. # last receiver z in meters (ignored if onlyone receiver)
enreg_surf_same_vertical = .false. # receivers inside the medium or at the surface
# display parameters
NTSTEP_BETWEEN_OUTPUT_INFO = 250 # 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
+output_color_image = .false. # output color image of the results
imagetype = 1 # display 1=displ 2=veloc 3=accel 4=pressure
cutsnaps = 1. # minimum amplitude in % for snapshots
meshvect = .true. # display mesh on vector plots or not
@@ -63,17 +64,17 @@
gnuplot = .false. # generate a GNUPLOT file for the grid
output_grid = .false. # save the grid in a text file or not
output_energy = .false. # compute and output acoustic and elastic energy (slows down the code significantly)
-output_wavefield_snapshot = .false. # output Ux,Uy,Uz text file for each output time (big files)
+output_wavefield_snapshot = .false. # output Ux,Uy,Uz text file for each output time (big files)
# velocity and density models
nbmodels = 1 # nb of different models
# define models as
-# I: (model_number 1 rho Vp Vs 0 0 QKappa Qmu 0 0 0 0 0 0) or
+# I: (model_number 1 rho Vp Vs 0 0 Qp Qs 0 0 0 0 0 0) or
# II: (model_number 2 rho c11 c13 c15 c33 c35 c55 0 0 0 0 0 0) or
-# III: (model_number 3 rhos rhof phi c kxx kxz kzz Ks Kf Kfr etaf mufr Qmu).
+# III: (model_number 3 rhos rhof phi c kxx kxz kzz Ks Kf Kfr etaf mufr Qs).
# For istropic elastic/acoustic material use I and set Vs to zero to make a given model acoustic, for anisotropic elastic use II,
# and for isotropic poroelastic material use III. The mesh can contain acoustic, elastic, and poroelastic models simultaneously.
-1 1 1.d0 2.d0 1.d0 0 0 10.d0 10.d0 0 0 0 0 0 0
+1 1 1.d0 2.d0 1.d0 0 0 9999.d0 9999.d0 0 0 0 0 0 0
# external mesh or not
read_external_mesh = .true.
@@ -86,11 +87,11 @@
# data concerning mesh, when generated using third-party app (more info in README)
# (see also absorbing_conditions above)
-mesh_file = ./mesh/canyon_mesh_file # file containing the mesh
-nodes_coords_file = ./mesh/canyon_nodes_coords_file # file containing the nodes coordinates
-materials_file = ./mesh/canyon_materials_file # file containing the material number for each element
-free_surface_file = ./mesh/canyon_free_surface_file # file containing the free surface
-absorbing_surface_file = ./mesh/canyon_absorbing_surface_file # file containing the absorbing surface
+mesh_file = ./DATA/mesh/canyon_mesh_file # file containing the mesh
+nodes_coords_file = ./DATA/mesh/canyon_nodes_coords_file # file containing the nodes coordinates
+materials_file = ./DATA/mesh/canyon_materials_file # file containing the material number for each element
+free_surface_file = ./DATA/mesh/canyon_free_surface_file # file containing the free surface
+absorbing_surface_file = ./DATA/mesh/canyon_absorbing_surface_file # file containing the absorbing surface
tangential_detection_curve_file = ./DATA/courbe_eros_nodes # file containing the curve delimiting the velocity model
#-----------------------------------------------------------------------------
@@ -105,7 +106,7 @@
nx = 95 # number of elements along X
# absorbing boundary parameters (see absorbing_conditions above)
-absorbbottom = .true.
+absorbbottom = .true.
absorbright = .true.
absorbtop = .false.
absorbleft = .true.
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/canyon/SOURCE_canyon
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/canyon/SOURCE_canyon 2011-09-01 09:29:08 UTC (rev 18868)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/canyon/SOURCE_canyon 2011-09-01 13:42:50 UTC (rev 18869)
@@ -1,8 +1,8 @@
# source parameters
source_surf = .false. # source inside the medium or at the surface
-xs = 2. # source location x in meters
-zs = 3. # source location z in meters
-source_type = 1 # 1 for plane P waves, 2 for plane SV waves, 3 for Rayleigh wave
+xs = 8. # source location x in meters
+zs = 9. # source location z in meters
+source_type = 2 # 1 for plane P waves, 2 for plane SV waves, 3 for Rayleigh wave
time_function_type = 1 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4, Heaviside = 5
f0 = 1.0 # dominant source frequency (Hz) if not Dirac or Heaviside
tshift = 0.0 # time shift when multi sources (if one source, must be zero)
Modified: seismo/2D/SPECFEM2D/trunk/src/meshfem2D/meshfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/meshfem2D/meshfem2D.F90 2011-09-01 09:29:08 UTC (rev 18868)
+++ seismo/2D/SPECFEM2D/trunk/src/meshfem2D/meshfem2D.F90 2011-09-01 13:42:50 UTC (rev 18869)
@@ -945,7 +945,23 @@
!--- compute position of the receivers and write the STATIONS file
- if (generate_STATIONS .and. .not. read_external_mesh) then
+ if (generate_STATIONS) then
+
+!! DK DK for now we cannot use both enreg_surf_same_vertical and read_external_mesh
+!! DK DK because we need to know splines to define the shape of the surface of the model
+ if(any(enreg_surf_same_vertical) .and. read_external_mesh) &
+ stop 'for now we cannot use both enreg_surf_same_vertical and read_external_mesh'
+
+!! DK DK if we read an external mesh, the splines are not defined for the shape of the surface and of the interfaces
+!! DK DK therefore let us allocate dummy arrays just to be able to call the "save_stations_file" subroutine
+ if (read_external_mesh) then
+ npoints_interface_top = 1
+ max_npoints_interface = 1
+ allocate(xinterface_top(1))
+ allocate(zinterface_top(1))
+ allocate(coefs_interface_top(1))
+ endif
+
call save_stations_file(nreceiverlines,nrec,xdeb,zdeb,xfin,zfin,enreg_surf_same_vertical, &
xinterface_top,zinterface_top,coefs_interface_top, &
npoints_interface_top,max_npoints_interface)
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/noise_tomography.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/noise_tomography.f90 2011-09-01 09:29:08 UTC (rev 18868)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/noise_tomography.f90 2011-09-01 13:42:50 UTC (rev 18869)
@@ -163,7 +163,7 @@
! moment tensor elements must be zero!
do i=1,NSOURCES
- if ( (Mxx(i) /= 0.d0) .or. (Mxz(i) /= 0.d0) .or. (Mzz(i) /= 0.d0) .or. &
+ if ( (Mxx(i) /= 0.d0) .or. (Mxz(i) /= 0.d0) .or. (Mzz(i) /= 0.d0) .or. &
(factor(i) /= 0.d0)) then
call exit_mpi('For noise simulations, all moment tensor elements must be zero. Exiting.')
endif
@@ -174,7 +174,7 @@
! =============================================================================================================
-! read in time series based on noise spectrum and construct noise "source" array
+! read in time series based on noise spectrum and construct noise "source" array
subroutine compute_source_array_noise(p_sv,NSTEP,deltat,nglob,ibool,ispec_noise, &
xi_noise,gamma_noise,xigll,zigll, &
time_function_noise,source_array_noise)
@@ -203,7 +203,7 @@
integer, parameter :: time_function_type = 1
-! ---------------------------------------------------------------------------------
+! ---------------------------------------------------------------------------------
! A NOTE ABOUT TIME FUNCTIONS FOR NOISE SIMULATIONS
!
! In noise forward modeling and inversion, "time function" is used to refer
@@ -215,7 +215,7 @@
! particular geographic region, you would have to use a time function created
! by inverse Fourier transforming a model of the noise spectrum for that
! region.
-!
+!
! IN CASES SUCH AS THIS--WHERE THE USER REQUIRES A REALISTIC MODEL OF THE
! SEISMIC NOISE FIELD--THE VARIABE "time_function_type" SHOULD BE SET TO 0
! AND A TIME FUNCTION ENCODING THE DESIRED NOISE SPECTRUM SHOULD BE
@@ -230,7 +230,7 @@
! ----------------------------------------------------------------------------------
time_function_noise(:) = 0._CUSTOM_REAL
- t0 = ((NSTEP-1)/2.)*deltat
+ t0 = ((NSTEP-1)/2.)*deltat
f0 = 15.0d0
aval = PI*PI*f0*f0
factor_noise = 1.d3
@@ -267,7 +267,7 @@
t = it*deltat
time_function_noise(it) = factor_noise * &
4.*aval**2. * (3. - 12.*aval*(t-t0)**2. + 4.*aval**2.*(t-t0)**4.) * &
- exp(-aval*(t-t0)**2.)
+ exp(-aval*(t-t0)**2.)
enddo
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/prepare_initialfield.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/prepare_initialfield.F90 2011-09-01 09:29:08 UTC (rev 18868)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/prepare_initialfield.F90 2011-09-01 13:42:50 UTC (rev 18869)
@@ -129,7 +129,7 @@
call exit_MPI('Unrecognized source_type: should be 1 for plane P waves, 2 for plane SV waves, 3 for Rayleigh wave')
endif
- ! allow negative angleforce(1): incidence from the right side of the domain
+ ! allow negative angleforce(1): incidence from the right side of the domain
angleforce_abs=abs(angleforce(1))
if (angleforce_abs > pi/2.d0 .and. source_type(1) /= 3) then
call exit_MPI("incorrect angleforce: must have 0 <= angleforce < 90")
@@ -330,7 +330,7 @@
+ C_plane(2) * ricker_Bielak_accel(t - sin(angleforce_refl)*x/c_refl - cos(angleforce_refl)*z/c_refl,f0(1))
enddo
-
+
endif
end subroutine prepare_initialfield
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2011-09-01 09:29:08 UTC (rev 18868)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2011-09-01 13:42:50 UTC (rev 18869)
@@ -692,7 +692,7 @@
! to compute analytical initial plane wave field
double precision :: angleforce_refl, c_inc, c_refl, cploc, csloc
double precision, dimension(2) :: A_plane, B_plane, C_plane
- double precision :: z0_source, x0_source, time_offset
+ double precision :: time_offset
! beyond critical angle
integer , dimension(:), allocatable :: left_bound,right_bound,bot_bound
@@ -792,7 +792,7 @@
integer :: iglob_target_to_replace, ispec3, i3, j3
!<NOISE_TOMOGRAPHY
- ! NOISE_TOMOGRAPHY = 0 - turn noise tomography subroutines off, resulting in
+ ! NOISE_TOMOGRAPHY = 0 - turn noise tomography subroutines off, resulting in
! an earthquake simulation rather than a noise simulation
! NOISE_TOMOGRAPHY = 1 - compute "generating wavefield" and store the result
@@ -800,7 +800,7 @@
! NOISE_TOMOGRAPHY = 2 - compute "ensemble forward wavefield"; if an adjoint
! simulation is planned, users need to manually extract cross-correlograms
- ! NOISE_TOMOGRAPHY = 3 - carry out adjoint simulation; for noise tomography
+ ! NOISE_TOMOGRAPHY = 3 - carry out adjoint simulation; for noise tomography
! applications, users need to supply adjoint source(s) based on cross-
! -correlograms from previous simulation
@@ -4397,7 +4397,7 @@
dux_dxl_np1,duz_dzl_np1,duz_dxl_np1,dux_dzl_np1,hprime_xx,hprimewgll_xx, &
hprime_zz,hprimewgll_zz,wxgll,wzgll,inv_tau_sigma_nu1, &
phi_nu1,inv_tau_sigma_nu2,phi_nu2,Mu_nu1,Mu_nu2,N_SLS, &
- deltat,coord,add_Bielak_conditions, x0_source, z0_source, &
+ deltat,coord,add_Bielak_conditions, x_source(1), z_source(1), &
A_plane, B_plane, C_plane, angleforce_refl, c_inc, c_refl, time_offset, f0(1),&
v0x_left(1,it),v0z_left(1,it),v0x_right(1,it),v0z_right(1,it),v0x_bot(1,it),v0z_bot(1,it), &
t0x_left(1,it),t0z_left(1,it),t0x_right(1,it),t0z_right(1,it),t0x_bot(1,it),t0z_bot(1,it), &
More information about the CIG-COMMITS
mailing list