[cig-commits] r18899 - in seismo/2D/SPECFEM2D/trunk: EXAMPLES/noise_uniform src/specfem2D
rmodrak at geodynamics.org
rmodrak at geodynamics.org
Mon Sep 12 19:20:10 PDT 2011
Author: rmodrak
Date: 2011-09-12 19:20:10 -0700 (Mon, 12 Sep 2011)
New Revision: 18899
Added:
seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/contribution_1
seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/contribution_2
Modified:
seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_1
seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_2
seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_3
seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/STATIONS_target_noise
seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/adj_cc.f90
seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/clean
seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/process.sh
seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/uniform.dat
seismo/2D/SPECFEM2D/trunk/src/specfem2D/noise_tomography.f90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
Log:
changed processing script for EXAMPLES/noise_uniform so that it computes both kernel contributions instead of just one
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_1
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_1 2011-09-12 18:58:40 UTC (rev 18898)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_1 2011-09-13 02:20:10 UTC (rev 18899)
@@ -13,7 +13,7 @@
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
+assign_external_model = .false. # define external earth model or not
READ_EXTERNAL_SEP_FILE = .false. # Read external SEP file from DATA/model_velocity.dat_input, or use routine
TURN_ATTENUATION_ON = .false. # turn attenuation on or off for solid medium
TURN_VISCATTENUATION_ON = .false. # turn viscous attenuation on or off
@@ -22,8 +22,8 @@
p_sv = .false. # set the type of calculation (P-SV or SH/membrane waves)
# time step parameters
-nt = 6000 # total number of time steps
-deltat = 1.d-3 # duration of a time step
+nt = 4000 # total number of time steps
+deltat = 4.d-2 # duration of a time step
# source parameters
NSOURCES = 1 # number of sources [source info read in CMTSOLUTION file]
@@ -49,9 +49,9 @@
enreg_surf_same_vertical = .false. # receivers inside the medium or at the surface
# display parameters
-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
+NTSTEP_BETWEEN_OUTPUT_INFO = 250 # display frequency in time steps
+output_postscript_snapshot = .false. # output Postscript snapshot 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
@@ -69,12 +69,12 @@
# 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 2300.d0 2450.d0 1425.d0 0 0 9999 9999 0 0 0 0 0 0
+1 1 3000.d0 5196.d0 3000.d0 0 0 1e6 1e6 0 0 0 0 0 0
# external mesh or not
read_external_mesh = .false.
@@ -101,9 +101,9 @@
interfacesfile = uniform.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
+xmin = 0.d0 # abscissa of left side of the model
+xmax = 300000.d0 # abscissa of right side of the model
+nx = 100 # number of elements along X
# absorbing boundary parameters (see absorbing_conditions above)
absorbbottom = .true.
@@ -113,4 +113,4 @@
# define the different regions of the model in the (nx,nz) spectral element mesh
nbregions = 1 # nb of regions and model number for each
-1 80 1 60 1
+1 100 1 80 1
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_2
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_2 2011-09-12 18:58:40 UTC (rev 18898)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_2 2011-09-13 02:20:10 UTC (rev 18899)
@@ -13,7 +13,7 @@
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
+assign_external_model = .false. # define external earth model or not
READ_EXTERNAL_SEP_FILE = .false. # Read external SEP file from DATA/model_velocity.dat_input, or use routine
TURN_ATTENUATION_ON = .false. # turn attenuation on or off for solid medium
TURN_VISCATTENUATION_ON = .false. # turn viscous attenuation on or off
@@ -22,8 +22,8 @@
p_sv = .false. # set the type of calculation (P-SV or SH/membrane waves)
# time step parameters
-nt = 6000 # total number of time steps
-deltat = 1.d-3 # duration of a time step
+nt = 4000 # total number of time steps
+deltat = 4.d-2 # duration of a time step
# source parameters
NSOURCES = 1 # number of sources [source info read in CMTSOLUTION file]
@@ -49,9 +49,9 @@
enreg_surf_same_vertical = .false. # receivers inside the medium or at the surface
# display parameters
-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
+NTSTEP_BETWEEN_OUTPUT_INFO = 250 # display frequency in time steps
+output_postscript_snapshot = .false. # output Postscript snapshot 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
@@ -69,12 +69,12 @@
# 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 2300.d0 2450.d0 1425.d0 0 0 9999 9999 0 0 0 0 0 0
+1 1 3000.d0 5196.d0 3000.d0 0 0 1e6 1e6 0 0 0 0 0 0
# external mesh or not
read_external_mesh = .false.
@@ -101,9 +101,9 @@
interfacesfile = uniform.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
+xmin = 0.d0 # abscissa of left side of the model
+xmax = 300000.d0 # abscissa of right side of the model
+nx = 100 # number of elements along X
# absorbing boundary parameters (see absorbing_conditions above)
absorbbottom = .true.
@@ -113,4 +113,4 @@
# define the different regions of the model in the (nx,nz) spectral element mesh
nbregions = 1 # nb of regions and model number for each
-1 80 1 60 1
+1 100 1 80 1
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_3
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_3 2011-09-12 18:58:40 UTC (rev 18898)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_3 2011-09-13 02:20:10 UTC (rev 18899)
@@ -13,7 +13,7 @@
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
+assign_external_model = .false. # define external earth model or not
READ_EXTERNAL_SEP_FILE = .false. # Read external SEP file from DATA/model_velocity.dat_input, or use routine
TURN_ATTENUATION_ON = .false. # turn attenuation on or off for solid medium
TURN_VISCATTENUATION_ON = .false. # turn viscous attenuation on or off
@@ -22,8 +22,8 @@
p_sv = .false. # set the type of calculation (P-SV or SH/membrane waves)
# time step parameters
-nt = 6000 # total number of time steps
-deltat = 1.d-3 # duration of a time step
+nt = 4000 # total number of time steps
+deltat = 4.d-2 # duration of a time step
# source parameters
NSOURCES = 1 # number of sources [source info read in CMTSOLUTION file]
@@ -49,9 +49,9 @@
enreg_surf_same_vertical = .false. # receivers inside the medium or at the surface
# display parameters
-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
+NTSTEP_BETWEEN_OUTPUT_INFO = 250 # display frequency in time steps
+output_postscript_snapshot = .false. # output Postscript snapshot 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
@@ -69,12 +69,12 @@
# 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 2300.d0 2450.d0 1425.d0 0 0 9999 9999 0 0 0 0 0 0
+1 1 3000.d0 5196.d0 3000.d0 0 0 1e6 1e6 0 0 0 0 0 0
# external mesh or not
read_external_mesh = .false.
@@ -101,9 +101,9 @@
interfacesfile = uniform.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
+xmin = 0.d0 # abscissa of left side of the model
+xmax = 300000.d0 # abscissa of right side of the model
+nx = 100 # number of elements along X
# absorbing boundary parameters (see absorbing_conditions above)
absorbbottom = .true.
@@ -113,4 +113,4 @@
# define the different regions of the model in the (nx,nz) spectral element mesh
nbregions = 1 # nb of regions and model number for each
-1 80 1 60 1
+1 100 1 80 1
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/STATIONS_target_noise
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/STATIONS_target_noise 2011-09-12 18:58:40 UTC (rev 18898)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/STATIONS_target_noise 2011-09-13 02:20:10 UTC (rev 18899)
@@ -1,3 +1,3 @@
-S0001 AA 1000.0000000 1500.0000000 0.0 0.0
-S0002 AA 2000.0000000 1500.0000000 0.0 0.0
-S0003 AA 3000.0000000 1500.0000000 0.0 0.0
+S0001 AA 100000.0000000 120000.0000000 0.0 0.0
+S0002 AA 150000.0000000 120000.0000000 0.0 0.0
+S0003 AA 200000.0000000 120000.0000000 0.0 0.0
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/adj_cc.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/adj_cc.f90 2011-09-12 18:58:40 UTC (rev 18898)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/adj_cc.f90 2011-09-13 02:20:10 UTC (rev 18899)
@@ -2,23 +2,32 @@
implicit none
-! flags
+!choose whether to differentiate traces 0, 1, or 2 times
+integer, parameter :: differentiate = 1
+
+!choose whether to bandpass filter
logical, parameter :: use_filtering = .true.
+
+!choose exactly one of the following window options
+logical, parameter :: use_negative_branch = .true.
logical, parameter :: use_positive_branch = .false.
logical, parameter :: use_custom_window = .false.
+
+!choose whether to time reverse (carried out subsequent to all other processing)
logical, parameter :: reverse = .true.
-integer, parameter :: differentiate = 1
+
+
! FILTERING PARAMETERS
real freq_low,freq_high
-data freq_low / 1.d-2 /
-data freq_high / 1.5d1 /
+data freq_low / 2.d-4 /
+data freq_high / 5.d-1 /
-! CUSTOM WINDOW PARAMETERS
-real :: w_delay, w_width, w_tukey
-data w_delay / 1.d0 /
-data w_width / 1.d2 /
-data w_tukey / 0.4 /
+! WINDOW PARAMETERS
+real :: t_begin, t_end, w_tukey
+data t_begin / 45.d0 /
+data t_end / 65.d0 /
+data w_tukey / 0.4 /
!see explanation below
! time variables
@@ -41,18 +50,10 @@
real :: alpha, beta
-! EXPLANATION OF CUSTOM WINDOW PARAMETERS
+! EXPLANATION OF WINDOW PARAMETERS
-!To select the desired branch of the cross-correlogram, we employ a Tukey window. A Tukey taper is just a variant of a cosine taper. We use three control parameters
+!To select the desired branch of the cross-correlogram, we employ a Tukey window. A Tukey taper is just a variant of a cosine taper.
-!W_DELAY controls the time offset of the window
-!W_WIDTH controls the width of the window (i.e., the total time range over which the window has non-zero support)
-!W_TUKEY controls the sharpness of the drop-off
-
-!In noise tomography applications, W_DELAY should be roughly equal to the surface wave travel time from the one receiver to the other.
-
-!Checks on W_WIDTH are carried out to make sure that the window makes sense and lies within a single branch of the cross-correlogram. If the the window falls outside these bounds, it will be adjusted.
-
!W_TUKEY is a number between 0 and 1, 0 being pure boxcar and 1 being pure cosine
@@ -89,8 +90,7 @@
!!!!!!!!!! READ INPUT !!!!!!!!!!!!!!!!!!!!
open(unit=1001,file=trim(file_in),status='old',action='read')
do it = 1, nt
- if (.not. reverse) read(1001,*) t(it), seismo_1(it)
- if (reverse) read(1001,*) t(it), seismo_1(nt-it+1)
+ read(1001,*) t(it), seismo_1(it)
end do
close(1001)
@@ -131,26 +131,21 @@
!!!!!!!!!! WINDOW !!!!!!!!!!!!!!!!!!!!
if (use_custom_window) then
- it_off = floor(w_delay/dt)
- it_wdt = 2*floor(w_width/(2.*dt))
-else
- it_off = int(nthalf/2)
- it_wdt = int(nthalf) + 1
-endif
-alpha = w_tukey
-
-if (use_positive_branch) then
- write(*,*) 'Choosing positive branch'
- it_begin = nthalf + it_off - it_wdt/2
- it_end = nthalf + it_off + it_wdt/2
- if (it_begin < nthalf) it_begin = nthalf
+ it_begin = floor((t_begin - t(1))/dt)
+ it_end = ceiling((t_end - t(1))/dt)
+ if (it_begin < 1) it_begin = 1
if (it_end > nt) it_end = nt
-else
+elseif (use_positive_branch) then
+ write(*,*) 'Choosing positive branch'
+ it_begin = nthalf+1
+ it_end = nt
+elseif (use_negative_branch) then
write(*,*) 'Choosing negative branch'
- it_begin = nthalf - it_off - it_wdt/2
- it_end = nthalf - it_off + it_wdt/2
- if (it_begin < 1) it_begin = 1
- if (it_end > nthalf) it_end = nthalf
+ it_begin = 1
+ it_end = nthalf
+else
+ write(*,*) 'Must select one of the following: positive_branch, &
+ negative_branch, custom_window.'
endif
write(*,'(a,2f10.3)') ' Time range: ', t(1), t(nt)
@@ -180,7 +175,6 @@
!!!!!!!!!! NORMALIZE !!!!!!!!!!!!!!!!!!!!
seismo_adj = - seismo_4/(DOT_PRODUCT(seismo_4,seismo_4)*dt)
-
!!!!!!!!!! WRITE ADJOINT SOURCE !!!!!!!!!!!!!!!!!!!!
open(unit=1002,file=trim(file_in)//'.adj',status='unknown',iostat=ios)
if (ios /= 0) write(*,*) 'Error opening output file.'
@@ -189,7 +183,8 @@
write(*,*) 'Writing to file: '//trim(file_in)//'.adj'
do it = 1,nt
- write(1002,'(f16.12,1pe18.8)'), t(it), seismo_adj(it)
+ if (.not. reverse) write(1002,'(f16.12,1pe18.8)'), t(it), seismo_adj(it)
+ if (reverse) write(1002,'(f16.12,1pe18.8)'), t(it), seismo_adj(nt+1-it)
end do
close(1002)
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/clean
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/clean 2011-09-12 18:58:40 UTC (rev 18898)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/clean 2011-09-13 02:20:10 UTC (rev 18899)
@@ -5,7 +5,7 @@
echo "Cleaning up."
-rm -rf DATA SEM NOISE_TOMOGRAPHY OUTPUT_FILES OUTPUT_ALL &> /dev/null
+rm -rf DATA SEM NOISE_TOMOGRAPHY OUTPUT_FILES OUTPUT_ALL SNAPSHOTS &> /dev/null
rm xmeshfem2D xspecfem2D &> /dev/null
rm adj_run &> /dev/null
rm log.* &> /dev/null
Added: seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/contribution_1
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/contribution_1 (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/contribution_1 2011-09-13 02:20:10 UTC (rev 18899)
@@ -0,0 +1,75 @@
+#!/bin/sh
+
+
+RUN_DIR=$PWD
+
+
+#prepare directories
+rm -rf DATA SEM NOISE_TOMOGRAPHY OUTPUT_FILES OUTPUT_ALL SNAPSHOTS
+mkdir -p DATA SEM NOISE_TOMOGRAPHY OUTPUT_FILES OUTPUT_ALL SNAPSHOTS DATA/NOISE_TOMOGRAPHY OUTPUT_FILES/NOISE_TOMOGRAPHY
+
+
+#prepare files
+cp SOURCE_noise DATA/SOURCE
+cp STATIONS_target_noise DATA/STATIONS_target
+cp uniform.dat DATA/
+echo 1 > DATA/NOISE_TOMOGRAPHY/irec_master
+if [ -f S_squared ]; then cp S_squared DATA/NOISE_TOMOGRAPHY/; fi
+
+
+#compile
+rm xmshefem2D
+rm xspecfem2D
+cd ../..
+make
+cd $RUN_DIR
+ln -s ../../bin/xmeshfem2D .
+ln -s ../../bin/xspecfem2D .
+
+
+#simulation 1
+cp Par_file_noise_1 DATA/Par_file
+./xmeshfem2D; ./xspecfem2D
+mkdir OUTPUT_ALL/step_1
+mv OUTPUT_FILES/image* OUTPUT_ALL/step_1
+mv OUTPUT_FILES/*.semd OUTPUT_ALL/step_1
+mv DATA/Par_file OUTPUT_ALL/step_1
+#mv OUTPUT_FILES/snapshot* SNAPSHOTS/
+
+
+#simulation 2
+cp Par_file_noise_2 DATA/Par_file
+./xmeshfem2D; ./xspecfem2D
+mkdir OUTPUT_ALL/step_2
+mv OUTPUT_FILES/image* OUTPUT_ALL/step_2
+mv OUTPUT_FILES/*.semd OUTPUT_ALL/step_2
+mv DATA/Par_file OUTPUT_ALL/step_2
+#mv OUTPUT_FILES/snapshot* SNAPSHOTS/
+
+
+#prepare adjoint source
+ADJ_CODE=adj_cc.f90
+sed -i'.bak' 's/reverse = .[a-z]*./reverse = .false./' $ADJ_CODE
+gfortran $ADJ_CODE -o adj_run
+cp OUTPUT_ALL/step_2/*.semd SEM/
+./adj_run SEM/S0003.AA.BXY.semd
+
+cd SEM
+rename '.semd' '' *.adj
+awk '{printf(" %20.10f %20.10f\n",$1,0.)}' < S0003.AA.BXY.adj > ZEROS
+cp ZEROS S0001.AA.BXX.adj; cp ZEROS S0001.AA.BXY.adj; cp ZEROS S0001.AA.BXZ.adj
+cp ZEROS S0002.AA.BXX.adj; cp ZEROS S0002.AA.BXY.adj; cp ZEROS S0002.AA.BXZ.adj
+cp ZEROS S0003.AA.BXX.adj; cp ZEROS S0003.AA.BXZ.adj
+cd ..
+
+
+#simulation 3
+cp Par_file_noise_3 DATA/Par_file
+./xmeshfem2D; ./xspecfem2D
+mkdir OUTPUT_ALL/step_3
+mv OUTPUT_FILES/image* OUTPUT_ALL/step_3
+mv SEM/*Y.adj OUTPUT_ALL/step_3
+mv OUTPUT_FILES/proc* OUTPUT_ALL/step_3
+mv DATA/Par_file OUTPUT_ALL/step_3
+mv OUTPUT_FILES/snapshot* SNAPSHOTS/
+
Property changes on: seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/contribution_1
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/contribution_2
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/contribution_2 (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/contribution_2 2011-09-13 02:20:10 UTC (rev 18899)
@@ -0,0 +1,75 @@
+#!/bin/sh
+
+
+RUN_DIR=$PWD
+
+
+#prepare directories
+rm -rf DATA SEM NOISE_TOMOGRAPHY OUTPUT_FILES OUTPUT_ALL SNAPSHOTS
+mkdir -p DATA SEM NOISE_TOMOGRAPHY OUTPUT_FILES OUTPUT_ALL SNAPSHOTS DATA/NOISE_TOMOGRAPHY OUTPUT_FILES/NOISE_TOMOGRAPHY
+
+
+#prepare files
+cp SOURCE_noise DATA/SOURCE
+cp STATIONS_target_noise DATA/STATIONS_target
+cp uniform.dat DATA/
+echo 3 > DATA/NOISE_TOMOGRAPHY/irec_master
+if [ -f S_squared ]; then cp S_squared DATA/NOISE_TOMOGRAPHY/; fi
+
+
+#compile
+rm xmshefem2D
+rm xspecfem2D
+cd ../..
+make
+cd $RUN_DIR
+ln -s ../../bin/xmeshfem2D .
+ln -s ../../bin/xspecfem2D .
+
+
+#simulation 1
+cp Par_file_noise_1 DATA/Par_file
+./xmeshfem2D; ./xspecfem2D
+mkdir OUTPUT_ALL/step_1
+mv OUTPUT_FILES/image* OUTPUT_ALL/step_1
+mv OUTPUT_FILES/*.semd OUTPUT_ALL/step_1
+mv DATA/Par_file OUTPUT_ALL/step_1
+#mv OUTPUT_FILES/snapshot* SNAPSHOTS/
+
+
+#simulation 2
+cp Par_file_noise_2 DATA/Par_file
+./xmeshfem2D; ./xspecfem2D
+mkdir OUTPUT_ALL/step_2
+mv OUTPUT_FILES/image* OUTPUT_ALL/step_2
+mv OUTPUT_FILES/*.semd OUTPUT_ALL/step_2
+mv DATA/Par_file OUTPUT_ALL/step_2
+#mv OUTPUT_FILES/snapshot* SNAPSHOTS/
+
+
+#prepare adjoint source
+ADJ_CODE=adj_cc.f90
+sed -i'.bak' 's/reverse = .[a-z]*./reverse = .true./' $ADJ_CODE
+gfortran $ADJ_CODE -o adj_run
+cp OUTPUT_ALL/step_2/*.semd SEM/
+./adj_run SEM/S0001.AA.BXY.semd
+
+cd SEM
+rename '.semd' '' *.adj
+awk '{printf(" %20.10f %20.10f\n",$1,0.)}' < S0001.AA.BXY.adj > ZEROS
+cp ZEROS S0001.AA.BXX.adj; cp ZEROS S0001.AA.BXZ.adj
+cp ZEROS S0002.AA.BXX.adj; cp ZEROS S0002.AA.BXY.adj; cp ZEROS S0002.AA.BXZ.adj
+cp ZEROS S0003.AA.BXX.adj; cp ZEROS S0003.AA.BXY.adj; cp ZEROS S0003.AA.BXZ.adj
+cd ..
+
+
+#simulation 3
+cp Par_file_noise_3 DATA/Par_file
+./xmeshfem2D; ./xspecfem2D
+mkdir OUTPUT_ALL/step_3
+mv OUTPUT_FILES/image* OUTPUT_ALL/step_3
+mv SEM/*Y.adj OUTPUT_ALL/step_3
+mv OUTPUT_FILES/proc* OUTPUT_ALL/step_3
+mv DATA/Par_file OUTPUT_ALL/step_3
+mv OUTPUT_FILES/snapshot* SNAPSHOTS/
+
Property changes on: seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/contribution_2
___________________________________________________________________
Name: svn:executable
+ *
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/process.sh
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/process.sh 2011-09-12 18:58:40 UTC (rev 18898)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/process.sh 2011-09-13 02:20:10 UTC (rev 18899)
@@ -1,65 +1,10 @@
-#!/bin/sh
+#!/bin/bash
-RUN_DIR=$PWD
+./contribution_1
+mv OUTPUT_ALL OUTPUT_ALL_1
+mv SNAPSHOT SNAPSHOT_1
-# prepare directories
-rm -rf DATA SEM NOISE_TOMOGRAPHY OUTPUT_FILES OUTPUT_ALL
-mkdir -p DATA SEM NOISE_TOMOGRAPHY OUTPUT_FILES OUTPUT_ALL DATA/NOISE_TOMOGRAPHY OUTPUT_FILES/NOISE_TOMOGRAPHY
+./contribution_2
+mv OUTPUT_ALL OUTPUT_ALL_2
+mv SNAPSHOT SNAPSHOT_2
-
-# prepare files
-cp SOURCE_noise DATA/SOURCE
-cp STATIONS_target_noise DATA/STATIONS_target
-cp uniform.dat DATA/
-echo 1 > DATA/NOISE_TOMOGRAPHY/irec_master
-if [ -f S_squared ]; then cp S_squared DATA/NOISE_TOMOGRAPHY/; fi
-
-# compile
-cd ../..
-make
-cd $RUN_DIR
-ln -s ../../bin/xmeshfem2D .
-ln -s ../../bin/xspecfem2D .
-
-#simulation 1
-cp Par_file_noise_1 DATA/Par_file
-./xmeshfem2D; ./xspecfem2D
-mkdir OUTPUT_ALL/step_1
-mv OUTPUT_FILES/image* OUTPUT_ALL/step_1
-mv OUTPUT_FILES/*.semd OUTPUT_ALL/step_1
-mv DATA/Par_file OUTPUT_ALL/step_1
-
-
-#simulation 2
-cp Par_file_noise_2 DATA/Par_file
-./xmeshfem2D; ./xspecfem2D
-mkdir OUTPUT_ALL/step_2
-
-ADJ_CODE=adj_cc.f90
-gfortran $ADJ_CODE -o adj_run
-cp OUTPUT_FILES/*.semd SEM
-./adj_run SEM/S0003.AA.BXY.semd
-
-cd SEM
-rename '.semd' '' *.adj
-awk '{printf(" %20.10f %20.10f\n",$1,0.)}' < S0003.AA.BXY.adj > ZEROS
-cp ZEROS S0001.AA.BXX.adj; cp ZEROS S0001.AA.BXY.adj; cp ZEROS S0001.AA.BXZ.adj
-cp ZEROS S0002.AA.BXX.adj; cp ZEROS S0002.AA.BXY.adj; cp ZEROS S0002.AA.BXZ.adj
-cp ZEROS S0003.AA.BXX.adj; cp ZEROS S0003.AA.BXZ.adj
-
-cd ..
-
-mv OUTPUT_FILES/image* OUTPUT_ALL/step_2
-mv OUTPUT_FILES/*.semd OUTPUT_ALL/step_2
-mv DATA/Par_file OUTPUT_ALL/step_2
-
-
-#simulation 3
-cp Par_file_noise_3 DATA/Par_file
-./xmeshfem2D; ./xspecfem2D
-mkdir OUTPUT_ALL/step_3
-mv OUTPUT_FILES/image* OUTPUT_ALL/step_3
-mv SEM/*Y.adj OUTPUT_ALL/step_3
-mv OUTPUT_FILES/proc* OUTPUT_ALL/step_3
-mv DATA/Par_file OUTPUT_ALL/step_3
-
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/uniform.dat
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/uniform.dat 2011-09-12 18:58:40 UTC (rev 18898)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/uniform.dat 2011-09-13 02:20:10 UTC (rev 18899)
@@ -10,17 +10,17 @@
#
2
0 0
- 4000 0
+ 300000 0
#
# interface number 2
#
2
- 0 3000
- 4000 3000
+ 0 240000
+ 300000 240000
#
# for each layer, we give the number of spectral elements in the vertical direction
#
#
# layer number 1
#
- 60
+ 80
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/noise_tomography.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/noise_tomography.f90 2011-09-12 18:58:40 UTC (rev 18898)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/noise_tomography.f90 2011-09-13 02:20:10 UTC (rev 18899)
@@ -199,11 +199,8 @@
real(kind=CUSTOM_REAL), dimension(NGLLZ) :: zigll
real(kind=CUSTOM_REAL), dimension(NGLLX) :: hxi, hpxi
real(kind=CUSTOM_REAL), dimension(NGLLZ) :: hgamma, hpgamma
- real(kind=CUSTOM_REAL) :: factor_noise, f0, aval, t0
+ real(kind=CUSTOM_REAL) :: factor_noise, aval, t0
- integer, parameter :: time_function_type = 1
-
-
! ---------------------------------------------------------------------------------
! A NOTE ABOUT TIME FUNCTIONS FOR NOISE SIMULATIONS
!
@@ -229,10 +226,14 @@
!
! ----------------------------------------------------------------------------------
+ !the following values are chosen to reproduce the time function from Fig 2a of
+ !"Tromp et al., 2010, Noise Cross-Correlation Sensitivity Kernels"
+
+ integer, parameter :: time_function_type = 4
+
time_function_noise(:) = 0._CUSTOM_REAL
t0 = ((NSTEP-1)/2.)*deltat
- f0 = 15.0d0
- aval = PI*PI*f0*f0
+ aval = 0.6d0
factor_noise = 1.d3
@@ -443,26 +444,102 @@
end subroutine save_surface_movie_noise
! =============================================================================================================
-! auxillary routine for writing out the array "mask_noise"; uses a format similar to
-! the one currently used for writing kernels
- subroutine write_mask_noise(nglob,coord,mask_noise)
+! auxillary routine
+ subroutine snapshot_all(ncol,nglob,filename,array_all)
implicit none
include "constants.h"
!input paramters
+ integer :: ncol,nglob
+ character(len=100) filename
+
+ real(kind=CUSTOM_REAL), dimension(ncol,nglob) :: array_all
+
+ !local parameters
+ integer :: i,iglob
+ real(kind=CUSTOM_REAL), dimension(ncol) :: row
+
+ open(unit=504,file=filename,status='unknown',action='write')
+
+ do iglob = 1,nglob
+
+ row(:) = array_all(:,iglob)
+
+ do i = 1,ncol-1
+ if (abs(row(i)) < 1.e-32) row(i) = 0.
+ write(unit=504,fmt='(1pe15.5e3)',advance='no') row(i)
+ enddo
+ if (abs(row(ncol)) < 1.e-32) row(ncol) = 0.
+ write(unit=504,fmt='(1pe15.5e3)') row(ncol)
+
+ enddo
+
+ close(504)
+
+
+ end subroutine snapshot_all
+
+
+! =============================================================================================================
+! auxillary routine
+ subroutine elem_to_glob(nspec,nglob,ibool,array_elem,array_glob)
+
+ implicit none
+ include "constants.h"
+
+ !input paramters
+ integer :: nspec, nglob
+ integer :: ibool(NGLLX,NGLLZ,nspec)
+
+ real(kind=CUSTOM_REAL), dimension(nglob) :: array_glob
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: array_elem
+
+ !local parameters
+ integer :: i,j,iglob,ispec
+ real(kind=CUSTOM_REAL) :: xx,zz
+
+ do ispec = 1, nspec
+ do j = 1, NGLLZ
+ do i = 1, NGLLX
+ iglob = ibool(i,j,ispec)
+ array_glob(iglob) = array_elem(i,j,ispec)
+ enddo
+ enddo
+ enddo
+
+ end subroutine elem_to_glob
+
+
+! =============================================================================================================
+! auxillary routine
+ subroutine snapshot_glob(nglob,coord,filename,array1)
+
+ implicit none
+ include "constants.h"
+
+ !input paramters
integer :: nglob
+ character(len=100) filename
+
real(kind=CUSTOM_REAL), dimension(2,nglob) :: coord
- real(kind=CUSTOM_REAL), dimension(nglob) :: mask_noise
+ real(kind=CUSTOM_REAL), dimension(nglob) :: array1
!local parameters
integer :: iglob
+ real(kind=CUSTOM_REAL) :: xx,zz
- open(unit=504,file='OUTPUT_FILES/mask_noise',status='unknown',action='write')
- do iglob = 1, nglob
- write(504,*) coord(1,iglob), coord(2,iglob), mask_noise(iglob)
- enddo
+
+ open(unit=504,file=filename,status='unknown',action='write')
+
+ do iglob = 1, nglob
+ xx = coord(1,iglob)
+ zz = coord(2,iglob)
+ write(504,*) xx, zz, array1(iglob)
+ enddo
+
close(504)
+ end subroutine snapshot_glob
- end subroutine write_mask_noise
+
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2011-09-12 18:58:40 UTC (rev 18898)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2011-09-13 02:20:10 UTC (rev 18899)
@@ -824,11 +824,17 @@
real(kind=CUSTOM_REAL), dimension(:), allocatable :: time_function_noise
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: source_array_noise
real(kind=CUSTOM_REAL), dimension(:), allocatable :: mask_noise
- !to avoid empty arrays depending on SH/P_SV, use separate arrays for x,y,z
+ !to avoid empty dimensions depending on SH/P_SV, use separate arrays for x,y,z
real(kind=CUSTOM_REAL), dimension(:), allocatable :: &
surface_movie_x_noise, surface_movie_y_noise, surface_movie_z_noise
+ integer :: ncol
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: rho_klglob
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: noise_all
+ character(len=100) :: fnoise
+
+
!>NOISE_TOMOGRAPHY
@@ -3603,6 +3609,11 @@
allocate(surface_movie_y_noise(nglob))
allocate(surface_movie_z_noise(nglob))
+ ncol = 6
+ allocate(rho_klglob(nglob))
+ allocate(noise_all(ncol,nglob))
+
+
!read in parameters for noise tomography
call read_parameters_noise(NOISE_TOMOGRAPHY,SIMULATION_TYPE,SAVE_FORWARD, &
any_acoustic,any_poroelastic,p_sv,irec_master, &
@@ -6471,6 +6482,34 @@
endif
+!<NOISE_TOMOGRAPHY
+
+ if (NOISE_TOMOGRAPHY == 1) then
+ !write(fnoise,"('OUTPUT_FILES/snapshot_eta_',i6.6)") it
+ !call snapshot_glob(nglob,coord,fnoise,accel_elastic(2,:))
+
+ elseif (NOISE_TOMOGRAPHY == 2) then
+ !write(fnoise,"('OUTPUT_FILES/snapshot_phi_',i6.6)") it
+ !call snapshot_glob(nglob,coord,fnoise,displ_elastic(2,:))
+
+ elseif (NOISE_TOMOGRAPHY == 3) then
+ call elem_to_glob(nspec,nglob,ibool,rho_kl,rho_klglob)
+
+ noise_all(1,:) = coord(1,:)
+ noise_all(2,:) = coord(2,:)
+ noise_all(3,:) = b_displ_elastic(2,:)
+ noise_all(4,:) = accel_elastic(2,:)
+ noise_all(5,:) = rho_k
+ noise_all(6,:) = rho_klglob
+
+ write(fnoise,"('OUTPUT_FILES/snapshot_all_',i6.6)") it
+ call snapshot_all(ncol,nglob,fnoise,noise_all)
+
+ endif
+
+!>NOISE_TOMOGRAPHY
+
+
!
!---- PostScript display
!
More information about the CIG-COMMITS
mailing list