[cig-commits] [commit] devel, master: commit one example fluid_solid kernel example provided by Luo Yang (d703590)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Wed Jun 18 15:23:03 PDT 2014
Repository : https://github.com/geodynamics/specfem2d
On branches: devel,master
Link : https://github.com/geodynamics/specfem2d/compare/fc67e6fd7ad890705b2b72b4b3c509accb22249e...e9ca46c40131588d89d7b0883250bc6584ce6b4c
>---------------------------------------------------------------
commit d703590eb2c940709045e450942a98e0f851bd0a
Author: Xie Zhinan <xiezhinan1984 at gmail.com>
Date: Fri Apr 5 08:49:39 2013 +0000
commit one example fluid_solid kernel example provided by Luo Yang
>---------------------------------------------------------------
d703590eb2c940709045e450942a98e0f851bd0a
.../Par_file | 75 ++++++++++-----------
LuoYang_fluid_soild_kernel/README | 31 +++++++++
.../SOURCE | 10 +--
LuoYang_fluid_soild_kernel/adj_source.f90 | 25 +++++++
.../interfaces.dat | 0
LuoYang_fluid_soild_kernel/plot_kernel | 41 +++++++++++
.../process.sh | 5 +-
.../process_kernel.sh | 0
LuoYang_fluid_soild_kernel/rho_kappa kernel.png | Bin 0 -> 25806 bytes
9 files changed, 142 insertions(+), 45 deletions(-)
diff --git a/salt_dome_CUBIT_mesh/Stacey_homogeneous/Par_file b/LuoYang_fluid_soild_kernel/Par_file
similarity index 77%
copy from salt_dome_CUBIT_mesh/Stacey_homogeneous/Par_file
copy to LuoYang_fluid_soild_kernel/Par_file
index f802d78..442e239 100644
--- a/salt_dome_CUBIT_mesh/Stacey_homogeneous/Par_file
+++ b/LuoYang_fluid_soild_kernel/Par_file
@@ -1,17 +1,17 @@
# title of job
-title = Salt dome model
+title = Test
# forward or adjoint simulation
-SIMULATION_TYPE = 1 # 1 = forward, 3 = adjoint + kernels
+SIMULATION_TYPE = 1 # 1 = forward, 2 = adjoint + kernels
NOISE_TOMOGRAPHY = 0 # 0 = earthquake simulation, 1/2/3 = noise simulation
-SAVE_FORWARD = .false. # save the last frame, needed for adjoint simulation
+SAVE_FORWARD = .true. # 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
PERFORM_CUTHILL_MCKEE = .false. # perform inverse Cuthill-McKee (1969) optimization/permutation for mesh numbering (can be very costly and not very useful)
-ngnod = 4 # number of control nodes per element (4 or 9)
+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
@@ -23,8 +23,8 @@ freq0 = 10 # frequency for viscous attenua
p_sv = .true. # set the type of calculation (P-SV or SH/membrane waves)
# time step parameters
-nt = 10000 # total number of time steps
-deltat = 4.25e-4 # duration of a time step (see section "How to choose the time step" of the manual for how to do this) (for the choice of deltat please refer to section 4.5 of the user manual)
+nt = 1600 # total number of time steps
+deltat = 1.1d-3 # duration of a time step (see section "How to choose the time step" of the manual for how to do this)
USER_T0 = 0.0d0 # use this t0 as earliest starting time rather than the automatically calculated one
time_stepping_scheme = 1 # 1 = Newmark (2nd order), 2 = LDDRK4-6 (4th-order 6-stage low storage Runge-Kutta), 3 = classical RK4 4th-order 4-stage Runge-Kutta
@@ -37,7 +37,7 @@ N_SLS = 2 # number of standard li
f0_attenuation = 5.196152422706633 # (Hz) relevant only if source is a Dirac or a Heaviside, otherwise it is f0 the dominant frequency of the source in the CMTSOLUTION file
# receiver set parameters for seismograms
-seismotype = 1 # record 1=displ 2=veloc 3=accel 4=pressure 5=curl of displ 6=the fluid potential
+seismotype = 6 # record 1=displ 2=veloc 3=accel 4=pressure 5=curl of displ 6=the fluid potential
NSTEP_BETWEEN_OUTPUT_SEISMOS = 5000000 # every how many time steps we save the seismograms (costly, do not use a very small value; if you use a very large value that is larger than the total number of time steps of the run, the seismograms will automatically be saved once at the end of the run anyway)
save_ASCII_seismograms = .true. # save seismograms in ASCII format or not
save_binary_seismograms_single = .true. # save seismograms in single precision binary format or not (can be used jointly with ASCII above to save both)
@@ -50,28 +50,28 @@ anglerec = 0.d0 # angle to rotate components at
rec_normal_to_surface = .false. # base anglerec normal to surface (external mesh and curve file needed)
# first receiver set (repeat these 6 lines and adjust nreceiversets accordingly)
-nrec = 1 # number of receivers
+nrec = 1 # number of receivers
xdeb = 1000. # first receiver x in meters
-zdeb = 1100. # first receiver z in meters
+zdeb = 1500. # first receiver z in meters
xfin = 1000. # last receiver x in meters (ignored if onlyone receiver)
-zfin = 1100. # last receiver z in meters (ignored if onlyone receiver)
+zfin = 1500. # last receiver z in meters (ignored if onlyone receiver)
enreg_surf_same_vertical = .false. # receivers inside the medium or at the surface
# display parameters
-NSTEP_BETWEEN_OUTPUT_INFO = 500 # every how many time steps we display information about the simulation (costly, do not use a very small value)
+NSTEP_BETWEEN_OUTPUT_INFO = 100 # every how many time steps we display information about the simulation (costly, do not use a very small value)
####
-NSTEP_BETWEEN_OUTPUT_IMAGES = 500 # every how many time steps we draw JPEG or PostScript pictures of the simulation (costly, do not use a very small value)
+NSTEP_BETWEEN_OUTPUT_IMAGES = 100 # every how many time steps we draw JPEG or PostScript pictures of the simulation (costly, do not use a very small value)
cutsnaps = 1. # minimum amplitude kept in % for the JPEG and PostScript snapshots; amplitudes below that are muted
#### for JPEG color images ####
output_color_image = .true. # output JPEG color image of the results every NSTEP_BETWEEN_OUTPUT_IMAGES time steps or not
-imagetype_JPEG = 2 # display 1=displ_Ux 2=displ_Uz 3=displ_norm 4=veloc_Vx 5=veloc_Vz 6=veloc_norm 7=accel_Ax 8=accel_Az 9=accel_norm 10=pressure
-factor_subsample_image = 1 # factor to spatially subsample color JPEG images output by the code (useful for very large models)
+imagetype_JPEG = 10 # display 1=displ_Ux 2=displ_Uz 3=displ_norm 4=veloc_Vx 5=veloc_Vz 6=veloc_norm 7=accel_Ax 8=accel_Az 9=accel_norm 10=pressure
+factor_subsample_image = 1.0 # (double precision) factor to subsample color images output by the code (useful for very large models)
POWER_DISPLAY_COLOR = 0.30d0 # non linear display to enhance small amplitudes in JPEG color images
DRAW_SOURCES_AND_RECEIVERS = .true. # display sources as orange crosses and receivers as green squares in JPEG images or not
DRAW_WATER_IN_BLUE = .true. # display acoustic layers as constant blue in JPEG images, because they likely correspond to water in the case of ocean acoustics or in the case of offshore oil industry experiments (if off, display them as greyscale, as for elastic or poroelastic elements, for instance for acoustic-only oil industry models of solid media)
USE_SNAPSHOT_NUMBER_IN_FILENAME = .false. # use snapshot number in the file name of JPEG color snapshots instead of the time step (for instance to create movies in an easier way later)
#### for PostScript snapshots ####
-output_postscript_snapshot = .true. # output Postscript snapshot of the results every NSTEP_BETWEEN_OUTPUT_IMAGES time steps or not
+output_postscript_snapshot = .false. # output Postscript snapshot of the results every NSTEP_BETWEEN_OUTPUT_IMAGES time steps or not
imagetype_postscript = 1 # display 1=displ vector 2=veloc vector 3=accel vector; small arrows are displayed for the vectors
meshvect = .true. # display mesh on PostScript plots or not
modelvect = .false. # display velocity model on PostScript plots or not
@@ -92,7 +92,7 @@ output_grid_ASCII = .false. # dump the grid in an ASCII tex
output_energy = .false. # compute and output total acoustic and elastic energy curves (slows down the code significantly)
# velocity and density models
-nbmodels = 1 # nb of different models
+nbmodels = 4 # 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
# II: (model_number 2 rho c11 c13 c15 c33 c35 c55 0 0 0 0 0 0) or
@@ -108,17 +108,20 @@ nbmodels = 1 # nb of different models
# sigma_xz = c15*dux_dx + c55*(duz_dx + dux_dz) + c35*duz_dz
# where the notations are for instance duz_dx = d(Uz) / dx.
#
-1 1 2000.d0 3500.d0 2000.d0 0 0 9999 9999 0 0 0 0 0 0
+1 1 2700.d0 3000.d0 1732.051d0 0 0 9999 9999 0 0 0 0 0 0
+2 1 2500.d0 2700.d0 0 0 0 9999 9999 0 0 0 0 0 0
+3 1 2200.d0 2500.d0 1443.375d0 0 0 9999 9999 0 0 0 0 0 0
+4 1 2200.d0 2200.d0 1343.375d0 0 0 9999 9999 0 0 0 0 0 0
# external mesh or not
-read_external_mesh = .true.
+read_external_mesh = .false.
# absorbing boundary active or not
PML_BOUNDARY_CONDITIONS = .false.
NELEM_PML_THICKNESS = 3
ROTATE_PML_ACTIVATE = .false.
ROTATE_PML_ANGLE = 30.
-STACEY_ABSORBING_CONDITIONS = .true.
+STACEY_ABSORBING_CONDITIONS = .false.
ADD_SPRING_TO_STACEY = .false.
# for horizontal periodic conditions: detect common points between left and right edges
@@ -135,11 +138,11 @@ PERIODIC_DETECT_TOL = 3.3334d-6
# data concerning mesh, when generated using third-party app (more info in README)
# (see also absorbing_conditions above)
-mesh_file = ./DATA/modelY1_mesh_file # file containing the mesh
-nodes_coords_file = ./DATA/modelY1_nodes_coords_file # file containing the nodes coordinates
-materials_file = ./DATA/modelY1_materials_file # file containing the material number for each element
-free_surface_file = ./DATA/modelY1_free_surface_file # file containing the free surface
-absorbing_surface_file = ./DATA/modelY1_absorbing_surface_file # file containing the absorbing surface
+mesh_file = ./DATA/Mesh_canyon/canyon_mesh_file # file containing the mesh
+nodes_coords_file = ./DATA/Mesh_canyon/canyon_nodes_coords_file # file containing the nodes coordinates
+materials_file = ./DATA/Mesh_canyon/canyon_materials_file # file containing the material number for each element
+free_surface_file = ./DATA/Mesh_canyon/canyon_free_surface_file # file containing the free surface
+absorbing_surface_file = ./DATA/Mesh_canyon/canyon_absorbing_surface_file # file containing the absorbing surface
CPML_element_file = Elements_CPML_list # file containing the CPML element numbers
tangential_detection_curve_file = ./DATA/courbe_eros_nodes # file containing the curve delimiting the velocity model
@@ -147,25 +150,21 @@ tangential_detection_curve_file = ./DATA/courbe_eros_nodes # file containing the
# PARAMETERS FOR INTERNAL MESHING
# file containing interfaces for internal mesh
-interfacesfile = interfaces_model5.dat
+interfacesfile = ./interfaces.dat
# 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
-nx = 80 # number of elements along X
+xmax = 5000.d0 # abscissa of right side of the model
+nx = 100 # number of elements along X
# absorbing boundary parameters (see absorbing_conditions above)
-absorbbottom = .true.
-absorbright = .true.
+absorbbottom = .false.
+absorbright = .false.
absorbtop = .false.
-absorbleft = .true.
+absorbleft = .false.
# define the different regions of the model in the (nx,nz) spectral element mesh
-nbregions = 7 # nb of regions and model number for each
-1 50 1 40 1
-1 50 41 50 2
-1 50 51 70 3
-1 50 71 80 4
-1 50 81 90 5
-1 50 91 100 6
-1 50 101 130 7
+nbregions = 3 # nb of regions and model number for each
+ 1 40 1 60 2
+41 60 1 60 1
+61 100 1 60 2
diff --git a/LuoYang_fluid_soild_kernel/README b/LuoYang_fluid_soild_kernel/README
new file mode 100644
index 0000000..61ab526
--- /dev/null
+++ b/LuoYang_fluid_soild_kernel/README
@@ -0,0 +1,31 @@
+----------------------------------------------------------------------
+README
+----------------------------------------------------------------------
+
+Kernel example provided by LuoYang
+
+TO RUN:
+
+0. Read the user manual in SPECFEM2D/doc/manual_SPECFEM2D.pdf
+
+1. in SPECFEM2D root directory, configure, e.g.,
+ > ./configure FC=gfortran
+ make, e.g. make all
+
+2. run mesher and solver for forward wavefield:
+ > cd EXAMPLES/LuoYang_fluid_soild_kernel/
+ > ./process.sh
+
+3. compute adjoint source:
+ > rm -rf xadj_source ; gfortran adj_source.f90 -o xadj_source
+ > xadj_source
+
+4. change Par_file with save_forward = .false. and SIMULATION_TYPE = 3
+
+5. run adjoint simulation:
+ > ./process_kernel.sh
+
+ optional: try plotting the kernels using the script
+ plot_kernel with command gnuplot plot_kernel
+
+---------------------------
diff --git a/canyon/SOURCE_canyon b/LuoYang_fluid_soild_kernel/SOURCE
similarity index 81%
copy from canyon/SOURCE_canyon
copy to LuoYang_fluid_soild_kernel/SOURCE
index d82f4c5..528ca38 100644
--- a/canyon/SOURCE_canyon
+++ b/LuoYang_fluid_soild_kernel/SOURCE
@@ -1,16 +1,16 @@
# The components of a moment tensor source must be given in N.m, not in dyne.cm as in the DATA/CMTSOLUTION source file of the 3D version of the code.
# source 1
source_surf = .false. # source inside the medium or at the surface
-xs = 8. # source location x in meters
-zs = 9. # source location z in meters
+xs = 4000. # source location x in meters
+zs = 1500. # source location z in meters
# source type: elastic force or acoustic pressure = 1 or moment tensor = 2;
# for a plane wave including converted and reflected waves at the free surface, P wave = 1, S wave = 2, Rayleigh wave = 3
# for a plane wave without converted nor reflected waves at the free surface, i.e. with the incident wave only, P wave = 4, S wave = 5
-source_type = 2
+source_type = 1
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
+f0 = 10.0 # dominant source frequency (Hz) if not Dirac or Heaviside
tshift = 0.0 # time shift when multi sources (if one source, must be zero)
-anglesource = 30. # angle of the source (for a force only); for a plane wave, this is the incidence angle; for moment tensor sources this is unused
+anglesource = 0.0 # angle of the source (for a force only); for a plane wave, this is the incidence angle; for moment tensor sources this is unused
Mxx = 1. # Mxx component (for a moment tensor source only)
Mzz = 1. # Mzz component (for a moment tensor source only)
Mxz = 0. # Mxz component (for a moment tensor source only)
diff --git a/LuoYang_fluid_soild_kernel/adj_source.f90 b/LuoYang_fluid_soild_kernel/adj_source.f90
new file mode 100644
index 0000000..0f756a8
--- /dev/null
+++ b/LuoYang_fluid_soild_kernel/adj_source.f90
@@ -0,0 +1,25 @@
+implicit none
+
+double precision :: time,pressure
+integer,parameter :: NSTEP = 1600
+integer :: i
+
+open(1,file="OUTPUT_FILES/S0001.AA.PRE.semp",status="old")
+
+open(10,file="SEM/S0001.AA.BXX.adj",status="unknown")
+
+open(20,file="SEM/S0001.AA.BXY.adj",status="unknown")
+
+open(30,file="SEM/S0001.AA.BXZ.adj",status="unknown")
+
+do i=1,NSTEP
+
+ read(1,*)time,pressure
+ if(time>1.274)pressure=0.d0
+ write(10,*)time,pressure
+ write(20,*)time,pressure
+ write(30,*)time,pressure
+
+enddo
+
+end
diff --git a/M2_UPPA/interfaces_M2_UPPA_flat.dat b/LuoYang_fluid_soild_kernel/interfaces.dat
similarity index 100%
copy from M2_UPPA/interfaces_M2_UPPA_flat.dat
copy to LuoYang_fluid_soild_kernel/interfaces.dat
diff --git a/LuoYang_fluid_soild_kernel/plot_kernel b/LuoYang_fluid_soild_kernel/plot_kernel
new file mode 100644
index 0000000..1d513b0
--- /dev/null
+++ b/LuoYang_fluid_soild_kernel/plot_kernel
@@ -0,0 +1,41 @@
+set term png
+set output "rho_kappa kernel.png"
+
+set palette defined ( -8e+16 "red", 0 "yellow", 8e+16 "blue")
+set pm3d map
+unset xtics
+unset ytics
+unset key
+unset grid
+set samples 2
+set cbrange [-8e+16:8e+16]
+set isosamples 2
+
+set multiplot
+set size 0.5,0.5
+set origin 0,0
+set title "rho kernel_acoustic only"
+splot "OUTPUT_FILES/proc000000_rho_kappa_kernel.dat" using 1:2:(($1>=2000.0 && $1<=3000.0)? 1/0 : $3) w points palette ps 0.02 pt 5
+
+set cbrange [-8e+16:8e+16]
+unset cbtics
+set size 0.5,0.5
+set origin 0.5,0
+set title "kappa kernel_acoustic only"
+splot "OUTPUT_FILES/proc000000_rho_kappa_kernel.dat" using 1:2:(($1>=2000.0 &&$1<=3000.0)? 1/0 : $4) w points palette ps 0.02 pt 5
+
+set size 0.5,0.5
+set origin 0,0.5
+set title "rho kernel_acoustic and elastic"
+splot "OUTPUT_FILES/proc000000_rho_kappa_kernel.dat" using 1:2:(($1>=2000.0 && $1<=3000.0)? 1/0 : $3) w points palette ps 0.02 pt 5,"OUTPUT_FILES/proc000000_rho_kappa_mu_kernel.dat" using 1:2:(($1>=2000.0 && $1<=3000.0)? $3 : 1/0) w points palette ps 0.02 pt 5
+
+set cbrange [-8e+16:8e+16]
+
+set size 0.5,0.5
+set origin 0.5,0.5
+set title "kappa kernel_acoustic and elastic"
+splot "OUTPUT_FILES/proc000000_rho_kappa_kernel.dat" using 1:2:(($1>=2000.0 &&$1<=3000.0)? 1/0 : $4) w points palette ps 0.02 pt 5,"OUTPUT_FILES/proc000000_rho_kappa_mu_kernel.dat" using 1:2:(($1>=2000.0 && $1<=3000.0)? $4 : 1/0) w points palette ps 0.02 pt 5
+
+unset multiplot
+
+
diff --git a/Tape2007_kernel/process.sh b/LuoYang_fluid_soild_kernel/process.sh
similarity index 92%
copy from Tape2007_kernel/process.sh
copy to LuoYang_fluid_soild_kernel/process.sh
index e0373de..b692396 100755
--- a/Tape2007_kernel/process.sh
+++ b/LuoYang_fluid_soild_kernel/process.sh
@@ -22,8 +22,9 @@ mkdir -p SEM
# sets up local DATA/ directory
cd DATA/
-ln -s ../Par_file_Tape2007_onerec Par_file
-ln -s ../SOURCE_001 SOURCE
+ln -s ../Par_file Par_file
+ln -s ../SOURCE SOURCE
+cp ../interfaces.dat interfaces.dat
cd ../
# cleans output files
diff --git a/Tape2007_kernel/process_kernel.sh b/LuoYang_fluid_soild_kernel/process_kernel.sh
similarity index 100%
copy from Tape2007_kernel/process_kernel.sh
copy to LuoYang_fluid_soild_kernel/process_kernel.sh
diff --git a/LuoYang_fluid_soild_kernel/rho_kappa kernel.png b/LuoYang_fluid_soild_kernel/rho_kappa kernel.png
new file mode 100644
index 0000000..831a1d0
Binary files /dev/null and b/LuoYang_fluid_soild_kernel/rho_kappa kernel.png differ
More information about the CIG-COMMITS
mailing list