[cig-commits] r18228 - in seismo/2D/SPECFEM2D/trunk: EXAMPLES EXAMPLES/Tape2007_kernel UTILS/visualization src/shared
carltape at geodynamics.org
carltape at geodynamics.org
Mon Apr 11 16:33:17 PDT 2011
Author: carltape
Date: 2011-04-11 16:33:17 -0700 (Mon, 11 Apr 2011)
New Revision: 18228
Added:
seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/
seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/Par_file_Tape2007_onerec
seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/README
seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/SOURCE_001
seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/adj_seismogram_Tape2007.f90
seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/interfaces_Tape2007.dat
seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/kernel_Tape2007_kernel_onerec.pdf
seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/process.sh
seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/process_kernel.sh
Modified:
seismo/2D/SPECFEM2D/trunk/UTILS/visualization/plot_wavefield.pl
seismo/2D/SPECFEM2D/trunk/src/shared/adj_seismogram.f90
Log:
added kernel example for Tape2007
Added: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/Par_file_Tape2007_onerec
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/Par_file_Tape2007_onerec (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/Par_file_Tape2007_onerec 2011-04-11 23:33:17 UTC (rev 18228)
@@ -0,0 +1,116 @@
+# title of job
+title = Tape-Liu-Tromp (GJI 2007)
+
+# forward or adjoint simulation
+SIMULATION_TYPE = 1 # 1 = forward, 2 = adjoint + kernels
+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
+
+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
+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
+Q0 = 1 # quality factor for viscous attenuation
+freq0 = 10 # frequency for viscous attenuation
+p_sv = .false. # set the type of calculation (P-SV or SH/membrane waves)
+
+# time step parameters
+nt = 3000 # total number of time steps
+deltat = 6.0d-2 # duration of a time step
+
+# source parameters
+NSOURCES = 1 # number of sources [source info read in CMTSOLUTION file]
+force_normal_to_surface = .false. # angleforce normal to surface (external mesh and curve file needed)
+
+# constants for attenuation
+N_SLS = 2 # number of standard linear solids 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
+rec_normal_to_surface = .false. # base anglerec normal to surface (external mesh and curve file needed)
+
+# first receiver line
+nrec = 1 # number of receivers
+xdeb = 180081.41 # first receiver x in meters
+zdeb = 388768.71 # first receiver z in meters
+xfin = 450000. # last receiver x in meters (ignored if only one receiver)
+zfin = 10000. # last receiver z in meters (ignored if only one receiver)
+enreg_surf_same_vertical = .false. # receivers inside the medium or at the surface
+
+
+# display parameters
+NTSTEP_BETWEEN_OUTPUT_INFO = 400 # 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
+cutsnaps = 1. # minimum amplitude in % for snapshots
+meshvect = .false. # display mesh on vector plots or not
+modelvect = .false. # display velocity model on vector plots
+boundvect = .true. # display boundary conditions on plots
+interpol = .true. # interpolation of the display or not
+pointsdisp = 6 # points for interpolation of display (set to 1 for lower-left corner only)
+subsamp = 1 # subsampling of color snapshots
+sizemax_arrows = 1.d0 # maximum size of arrows on vector plots in cm
+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)
+
+# velocity and density models
+nbmodels = 1 # nb of different models
+# define models as
+# 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 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 2600.d0 5800.d0 3500.0d0 0 0 10.d0 10.d0 0 0 0 0 0 0
+
+# external mesh or not
+read_external_mesh = .false.
+
+# absorbing boundary active or not
+absorbing_conditions = .true.
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR EXTERNAL MESHING
+
+# data concerning mesh, when generated using third-party app (more info in README)
+# (see also absorbing_conditions above)
+mesh_file = ./DATA/ice_water_rock_1D/ice_water_rock_1D_elements # file containing elements
+nodes_coords_file = ./DATA/ice_water_rock_1D/ice_water_rock_1D_nodes # file containing node coordinates
+materials_file = ./DATA/ice_water_rock_1D/ice_water_rock_1D_material # file containing material index for each element
+free_surface_file = ./DATA/ice_water_rock_1D/ice_water_rock_1D_surface_free # file containing free surface
+absorbing_surface_file = ./DATA/ice_water_rock_1D/ice_water_rock_1D_surface_absorb # file containing absorbing surface
+tangential_detection_curve_file = ./DATA/courbe_eros_nodes # file containing the curve delimiting the velocity model
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR INTERNAL MESHING
+
+# file containing interfaces for internal mesh
+interfacesfile = ../interfaces_Tape2007.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 = 480000.d0 # abscissa of right side of the model
+nx = 40 # number of elements along X
+
+# absorbing boundary parameters (see absorbing_conditions above)
+absorbbottom = .true.
+absorbright = .true.
+absorbtop = .true.
+absorbleft = .true.
+
+# define the different regions of the model in the (nx,nz) spectral element mes
+nbregions = 1 # nb of regions and model number for each
+1 40 1 40 1
Added: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/README
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/README (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/README 2011-04-11 23:33:17 UTC (rev 18228)
@@ -0,0 +1,35 @@
+----------------------------------------------------------------------
+README
+----------------------------------------------------------------------
+
+Kernel example for Tape-Liu-Tromp (GJI 2007).
+
+TO RUN:
+
+0. Read the user manual in SPECFEM2D/doc/manual_SPECFEM2D.pdf
+
+1. in SPECFEM2D root directory, configure, e.g.,
+ > ./configure FC=gfortran
+
+2. modify setup/constants.h file, setting USER_TO = 48.0 s
+
+3. compile:
+ > make all
+
+4. run mesher and solver for forward wavefield:
+ > cd EXAMPLES/Tape2007_kernel/
+ > ./process.sh
+
+5. compute adjoint source:
+ > rm -rf xadj_seismogram ; gfortran adj_seismogram_Tape2007.f90 -o xadj_seismogram
+ > xadj_seismogram
+
+6. change Par_file with save_forward = .false. and SIMULATION_TYPE = 2
+
+7. run adjoint simulation:
+ > ./process_kernel.sh
+
+ optional: try plotting the kernels using the script
+ SPECFEM2D/UTILS/visualization/plot_wavefield.pl
+
+---------------------------
Added: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/SOURCE_001
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/SOURCE_001 (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/SOURCE_001 2011-04-11 23:33:17 UTC (rev 18228)
@@ -0,0 +1,13 @@
+# 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_surf = .false. # source inside the medium or at the surface
+xs = 362079.06 # source location x in meters
+zs = 61219.23 # source location z in meters
+source_type = 1 # elastic force or acoustic pressure = 1 or moment tensor = 2
+time_function_type = 2 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4, Heaviside = 5
+f0 = 8.400e-02 # dominant source frequency (Hz) if not Dirac or Heaviside
+tshift = 0.000e+00 # time shift when multi sources (if one source, must be zero)
+angleforce = 0.00 # angle of the source (for a force only)
+Mxx = 1.000000 # Mxx component (for a moment tensor source only)
+Mzz = -1.000000 # Mzz component (for a moment tensor source only)
+Mxz = 0.000000 # Mxz component (for a moment tensor source only)
+factor = 1.000e+10 # amplification factor
Added: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/adj_seismogram_Tape2007.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/adj_seismogram_Tape2007.f90 (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/adj_seismogram_Tape2007.f90 2011-04-11 23:33:17 UTC (rev 18228)
@@ -0,0 +1,179 @@
+
+!========================================================================
+!
+! S P E C F E M 2 D Version 6.1
+! ------------------------------
+!
+! Copyright Universite de Pau, CNRS and INRIA, France,
+! and Princeton University / California Institute of Technology, USA.
+! Contributors: Dimitri Komatitsch, dimitri DOT komatitsch aT univ-pau DOT fr
+! Nicolas Le Goff, nicolas DOT legoff aT univ-pau DOT fr
+! Roland Martin, roland DOT martin aT univ-pau DOT fr
+! Christina Morency, cmorency aT princeton DOT edu
+!
+! This software is a computer program whose purpose is to solve
+! the two-dimensional viscoelastic anisotropic or poroelastic wave equation
+! using a spectral-element method (SEM).
+!
+! This software is governed by the CeCILL license under French law and
+! abiding by the rules of distribution of free software. You can use,
+! modify and/or redistribute the software under the terms of the CeCILL
+! license as circulated by CEA, CNRS and INRIA at the following URL
+! "http://www.cecill.info".
+!
+! As a counterpart to the access to the source code and rights to copy,
+! modify and redistribute granted by the license, users are provided only
+! with a limited warranty and the software's author, the holder of the
+! economic rights, and the successive licensors have only limited
+! liability.
+!
+! In this respect, the user's attention is drawn to the risks associated
+! with loading, using, modifying and/or developing or reproducing the
+! software by the user in light of its specific status of free software,
+! that may mean that it is complicated to manipulate, and that also
+! therefore means that it is reserved for developers and experienced
+! professionals having in-depth computer knowledge. Users are therefore
+! encouraged to load and test the software's suitability as regards their
+! requirements in conditions enabling the security of their systems and/or
+! data to be ensured and, more generally, to use and operate it in the
+! same conditions as regards security.
+!
+! The full text of the license is available in file "LICENSE".
+!
+!========================================================================
+
+ program adj_seismogram
+
+! This program cuts a certain portion of the seismograms and convert it
+! into the adjoint source for generating banana-dougnut kernels
+!
+! rm -rf xadj_seismogram ; gfortran adj_seismogram.f90 -o xadj_seismogram
+!
+
+ implicit none
+
+ integer, parameter :: NSTEP = 3000
+ integer, parameter :: nrec = 1
+ double precision, parameter :: t0 = 48.0 ! labeled as 'time delay'
+ double precision, parameter :: deltat = 0.06
+ double precision, parameter :: EPS = 1.d-40
+
+ integer :: itime,icomp,istart,iend,nlen,irec,NDIM,NDIMr,adj_comp
+ double precision :: time,tstart(nrec),tend(nrec)
+ character(len=150), dimension(nrec) :: station_name
+ double precision, dimension(NSTEP) :: time_window
+ double precision :: seism(NSTEP,3),Nnorm,seism_win(NSTEP)
+ double precision :: seism_veloc(NSTEP),seism_accel(NSTEP),ft_bar(NSTEP)
+ character(len=3) :: compr(2),comp(3)
+ character(len=150) :: filename
+
+ NDIM=3
+ comp = (/"BHX","BHY","BHZ"/)
+
+! number of components
+! NDIMr=2 ! P-SV
+ NDIMr=1 ! SH (membrane)
+! compr = (/"BHX","BHZ"/) ! P-SV
+ compr = (/"BHY","tmp"/) ! SH (membrane)
+! list of stations
+ station_name(1) = 'S0001'
+
+! KEY: 'absolute' time interval for window
+ tstart(1) = 135.d0
+ tend(1) = 180.d0
+
+! chose the component for the adjoint source (adj_comp = 1:X, 2:Y, 3:Z)
+ adj_comp = 2
+!!!!
+
+ do irec =1,nrec
+
+ do icomp = 1, NDIMr
+
+ filename = 'OUTPUT_FILES/'//trim(station_name(irec))//'.AA.'// compr(icomp) // '.semd'
+ open(unit = 10, file = trim(filename))
+
+ do itime = 1,NSTEP
+ read(10,*) time , seism(itime,icomp)
+ enddo
+
+ enddo
+
+ if(NDIMr==2)then
+ seism(:,3) = seism(:,2)
+ seism(:,2) = 0.d0
+ else
+ seism(:,2) = seism(:,1)
+ seism(:,1) = 0.d0
+ seism(:,3) = 0.d0
+ endif
+
+ close(10)
+
+
+ istart = max(floor(tstart(irec)/deltat),1)
+ iend = min(floor(tend(irec)/deltat),NSTEP)
+ print*,'istart =',istart, 'iend =', iend
+ print*,'tstart =',istart*deltat, 'tend =', iend*deltat
+ if(istart >= iend) stop 'check istart,iend'
+ nlen = iend - istart +1
+
+ do icomp = 1, NDIM
+
+ print*,comp(icomp)
+
+ filename = 'OUTPUT_FILES/'//trim(station_name(irec))//'.AA.'// comp(icomp) // '.adj'
+ open(unit = 11, file = trim(filename))
+
+ time_window(:) = 0.d0
+ seism_win(:) = seism(:,icomp)
+ seism_veloc(:) = 0.d0
+ seism_accel(:) = 0.d0
+
+ do itime =istart,iend
+! time_window(itime) = 1.d0 - cos(pi*(itime-1)/NSTEP+1)**10 ! cosine window
+ time_window(itime) = 1.d0 - (2* (dble(itime) - istart)/(iend-istart) -1.d0)**2 ! Welch window
+ enddo
+
+ do itime = 2,NSTEP-1
+ seism_veloc(itime) = (seism_win(itime+1) - seism_win(itime-1))/(2*deltat)
+ enddo
+ seism_veloc(1) = (seism_win(2) - seism_win(1))/deltat
+ seism_veloc(NSTEP) = (seism_win(NSTEP) - seism_win(NSTEP-1))/deltat
+
+ do itime = 2,NSTEP-1
+ seism_accel(itime) = (seism_veloc(itime+1) - seism_veloc(itime-1))/(2*deltat)
+ enddo
+ seism_accel(1) = (seism_veloc(2) - seism_veloc(1))/deltat
+ seism_accel(NSTEP) = (seism_veloc(NSTEP) - seism_veloc(NSTEP-1))/deltat
+
+ Nnorm = deltat * sum(time_window(:) * seism_win(:) * seism_accel(:))
+! Nnorm = deltat * sum(time_window(:) * seism_veloc(:) * seism_veloc(:))
+! cross-correlation traveltime adjoint source
+ if(abs(Nnorm) > EPS) then
+! ft_bar(:) = - seism_veloc(:) * time_window(:) / Nnorm
+ ft_bar(:) = seism_veloc(:) * time_window(:) / Nnorm
+ print*,'Norm =', Nnorm
+ else
+ print *, 'norm < EPS for file '
+ print*,'Norm =', Nnorm
+ ft_bar(:) = 0.d0
+ endif
+
+ do itime =1,NSTEP
+ if(icomp == adj_comp) then
+ write(11,*) (itime-1)*deltat - t0, ft_bar(itime)
+ else
+ write(11,*) (itime-1)*deltat - t0, 0.d0
+ endif
+ enddo
+
+ enddo
+ close(11)
+
+ enddo
+ print*,'*************************'
+ print*,'The input files (S****.AA.BHX/BHY/BHZ.adj) needed to run the adjoint simulation are in OUTPUT_FILES'
+ print*,'*************************'
+
+ end program adj_seismogram
Added: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/interfaces_Tape2007.dat
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/interfaces_Tape2007.dat (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/interfaces_Tape2007.dat 2011-04-11 23:33:17 UTC (rev 18228)
@@ -0,0 +1,18 @@
+# number of interfaces
+ 2
+#
+# for each interface below, we give the number of points and then x,z for each point
+#
+# interface number 1 (bottom of the mesh)
+ 2
+ 0 0
+ 480000 0
+# interface number 5 (topography, top of the mesh)
+ 2
+ 0 480000
+ 480000 480000
+#
+# for each layer, we give the number of spectral elements in the vertical direction
+#
+# layer number 1
+ 40
\ No newline at end of file
Added: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/kernel_Tape2007_kernel_onerec.pdf
===================================================================
(Binary files differ)
Property changes on: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/kernel_Tape2007_kernel_onerec.pdf
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/process.sh
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/process.sh (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/process.sh 2011-04-11 23:33:17 UTC (rev 18228)
@@ -0,0 +1,64 @@
+#!/bin/bash
+#
+# script runs mesher and solver (in serial)
+# using this example setup
+#
+
+echo "running example: `date`"
+currentdir=`pwd`
+
+echo
+echo "(will take about 1 minute)"
+echo
+
+# sets up directory structure in current example directoy
+echo
+echo " setting up example..."
+echo
+
+mkdir -p OUTPUT_FILES
+mkdir -p DATA
+
+# sets up local DATA/ directory
+cd DATA/
+ln -s ../Par_file_Tape2007_onerec Par_file
+ln -s ../SOURCE_001 SOURCE
+cd ../
+
+# cleans output files
+rm -rf OUTPUT_FILES/*
+
+# compiles executables in root directory
+cd ../../
+make > tmp.log
+cd $currentdir
+
+# links executables
+rm -f xmeshfem2D xspecfem2D
+ln -s ../../bin/xmeshfem2D
+ln -s ../../bin/xspecfem2D
+
+# stores setup
+cp DATA/Par_file OUTPUT_FILES/
+cp DATA/SOURCE OUTPUT_FILES/
+cp DATA/SOURCE_xz.dat OUTPUT_FILES/
+cp DATA/STATIONS OUTPUT_FILES/
+cp DATA/STATIONS_target OUTPUT_FILES/
+
+# runs database generation
+echo
+echo " running mesher..."
+echo
+./xmeshfem2D > OUTPUT_FILES/output_mesher.txt
+
+# runs simulation
+echo
+echo " running solver..."
+echo
+./xspecfem2D > OUTPUT_FILES/output_solver.txt
+
+echo
+echo "see results in directory: OUTPUT_FILES/"
+echo
+echo "done"
+echo `date`
Property changes on: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/process.sh
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/process_kernel.sh
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/process_kernel.sh (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/process_kernel.sh 2011-04-11 23:33:17 UTC (rev 18228)
@@ -0,0 +1,24 @@
+#!/bin/bash
+#
+# script runs solver in serial
+#
+
+echo "running example: `date`"
+currentdir=`pwd`
+
+echo
+echo "(will take about 1 minute)"
+echo
+
+# runs simulation
+echo
+echo " running mesher and solver..."
+echo
+./xmeshfem2D > OUTPUT_FILES/output_mesher_kernel.txt
+./xspecfem2D > OUTPUT_FILES/output_solver_kernel.txt
+
+echo
+echo "see results in directory: OUTPUT_FILES/"
+echo
+echo "done"
+echo `date`
Property changes on: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Tape2007_kernel/process_kernel.sh
___________________________________________________________________
Name: svn:executable
+ *
Modified: seismo/2D/SPECFEM2D/trunk/UTILS/visualization/plot_wavefield.pl
===================================================================
--- seismo/2D/SPECFEM2D/trunk/UTILS/visualization/plot_wavefield.pl 2011-04-11 21:39:22 UTC (rev 18227)
+++ seismo/2D/SPECFEM2D/trunk/UTILS/visualization/plot_wavefield.pl 2011-04-11 23:33:17 UTC (rev 18228)
@@ -13,13 +13,16 @@
# 1. must have output_wavefield_snapshot = .true. in Par_file to generate wavefield snapshots
# 2. boundaries between material regions are NOT plotted here
#
-# EXAMPLES:
+# WAVEFIELD EXAMPLES:
# plot_wavefield.pl 100/1600/400 200/800/200 0/4/0/4 1/0.1/1/0.1 -3/-3/-3 6/6/6 0.0/0.001 1/0/1/1/1 1.7/1/1 1/0/1/120 M2_UPPA PSV
# plot_wavefield.pl 400/2800/400 400/2000/400 0/200/0/80 50/10/40/10 -3/-3/-3 4/4/4 -8.0/0.02 1/0/1/1/1 3.0/1/0 1/0/1/200 Tromp2005 PSV
# plot_wavefield.pl 400/2800/400 400/2000/400 0/200/0/80 50/10/40/10 -3/-3/-3 10/10/10 -8.0/0.02 0/1/0/0/1 3.0/1/0 1/0/1/200 Tromp2005 SH
# plot_wavefield.pl 400/4800/400 400/2800/800 0/480/0/480 120/20/120/20 -3/-3/-3 6/6/6 -48.0/0.06 0/1/0/0/1 1.7/1/0 1/0/1/120 Tape2007 onerec
#
+# KERNEL EXAMPLES:
+# plot_wavefield.pl 400/3000/400 400/2800/800 0/480/0/480 120/20/120/20 -8/-8/-8 1/1/1 -48.0/0.06 0/1/0/0/0 2.0/1/0 1/0/1/120 Tape2007_kernel onerec
#
+#
#----------------------------------------------
# USER INPUT
@@ -42,11 +45,9 @@
$numf = ($pend - $pfirst)/$pint + 1;
# directory with data files
-#$idir1 = "/data/svn/seismo/2D/SPECFEM2D_work/OUTPUT_FILES"; # if running from the default directory
-$idir1 = "/data/svn/seismo/2D/SPECFEM2D_work/EXAMPLES/$tlab/OUTPUT_FILES"; # if running from an examples directory
-if($itype==0) {$idir1 = "/data/svn/seismo/2D/SPECFEM2D_work/OUTPUT_FILES_${tlab}";}
-#$idir1 = "${idir1}_${ttag}";
-#$idir1 = "${idir1}_${tlab}";
+$bdir = "/data/svn/seismo/2D/SPECFEM2D_work";
+#$idir1 = "$bdir/OUTPUT_FILES"; # if running from the default directory
+$idir1 = "$bdir/EXAMPLES/$tlab/OUTPUT_FILES"; # if running from an examples directory
if (not -e $idir1) {die("check if idir1 $idir1 exist or not\n");}
# plot the color frames or not
@@ -598,7 +599,7 @@
if($iheader==1) {
# plot title and GMT header
$plabel = "plot_wavefield.pl";
- $ux = -$xwid;
+ $ux = 0;
$uy = $ywid + 0.3;
$Utag = "-U/$ux/$uy/$plabel"; # GMT header
$shift = "-X0i -Y0.7i";
Modified: seismo/2D/SPECFEM2D/trunk/src/shared/adj_seismogram.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/shared/adj_seismogram.f90 2011-04-11 21:39:22 UTC (rev 18227)
+++ seismo/2D/SPECFEM2D/trunk/src/shared/adj_seismogram.f90 2011-04-11 23:33:17 UTC (rev 18228)
@@ -42,133 +42,135 @@
!
!========================================================================
- program adj_seismogram
+program adj_seismogram
-! This program cuts a certain portion of the seismograms and convert it
-! into the adjoint source for generating banana-dougnut kernels
+ ! This program cuts a certain portion of the seismograms and convert it
+ ! into the adjoint source for generating banana-dougnut kernels
+ ! rm -rf xadj_seismogram ; gfortran adj_seismogram.f90 -o xadj_seismogram
- implicit none
+ implicit none
- integer, parameter :: NSTEP = 3000
- integer, parameter :: nrec = 1
- double precision, parameter :: t0 = 12
- double precision, parameter :: deltat = 6d-2
- double precision, parameter :: EPS = 1.d-40
+ integer, parameter :: NSTEP = 3000
+ integer, parameter :: nrec = 1
+ double precision, parameter :: t0 = 48.0 ! labeled as 'time delay'
+ double precision, parameter :: deltat = 0.06
+ double precision, parameter :: EPS = 1.d-40
- integer :: itime,icomp,istart,iend,nlen,irec,NDIM,NDIMr,adj_comp
- double precision :: time,tstart(nrec),tend(nrec)
- character(len=150), dimension(nrec) :: station_name
- double precision, dimension(NSTEP) :: time_window
- double precision :: seism(NSTEP,3),Nnorm,seism_win(NSTEP)
- double precision :: seism_veloc(NSTEP),seism_accel(NSTEP),ft_bar(NSTEP)
- character(len=3) :: compr(2),comp(3)
- character(len=150) :: filename
+ integer :: itime,icomp,istart,iend,nlen,irec,NDIM,NDIMr,adj_comp
+ double precision :: time,tstart(nrec),tend(nrec)
+ character(len=150), dimension(nrec) :: station_name
+ double precision, dimension(NSTEP) :: time_window
+ double precision :: seism(NSTEP,3),Nnorm,seism_win(NSTEP)
+ double precision :: seism_veloc(NSTEP),seism_accel(NSTEP),ft_bar(NSTEP)
+ character(len=3) :: compr(2),comp(3)
+ character(len=150) :: filename
- NDIM=3
- comp = (/"BHX","BHY","BHZ"/)
+ NDIM=3
+ comp = (/"BHX","BHY","BHZ"/)
-! which calculation: P-SV (use (1)) or SH (membrane) (use (2)) waves
- NDIMr=2 !(1)
-! NDIMr=1 !(2)
-! list of stations
- station_name(1) = 'S0001'
- tstart(1) = 100d0 + t0
- tend(1) = 120d0 + t0
-! which calculation: P-SV (use (1)) or SH (membrane) (use (2)) waves
- compr = (/"BHX","BHZ"/) !(1)
-! compr = (/"BHY","dummy"/) !(2)
-! chose the component for the adjoint source (adj_comp = 1: X, 2:Y, 3:Z)
- adj_comp = 1
-!!!!
+ ! number of components
+ !NDIMr=2 ! P-SV
+ NDIMr=1 ! SH (membrane)
+ !compr = (/"BHX","BHZ"/) ! P-SV
+ compr = (/"BHY","tmp"/) ! SH (membrane)
+ ! list of stations
+ station_name(1) = 'S0001'
- do irec =1,nrec
+ ! KEY: 'absolute' time interval for window
+ tstart(1) = 135.d0 + t0
+ tend(1) = 180.d0 + t0
- do icomp = 1, NDIMr
+ ! chose the component for the adjoint source (adj_comp = 1:X, 2:Y, 3:Z)
+ adj_comp = 2
- filename = 'OUTPUT_FILES/'//trim(station_name(irec))//'.AA.'// compr(icomp) // '.semd'
- open(unit = 10, file = trim(filename))
+ do irec =1,nrec
- do itime = 1,NSTEP
- read(10,*) time , seism(itime,icomp)
- enddo
+ do icomp = 1, NDIMr
+ filename = 'OUTPUT_FILES/'//trim(station_name(irec))//'.AA.'// compr(icomp) // '.semd'
+ open(unit = 10, file = trim(filename))
+
+ do itime = 1,NSTEP
+ read(10,*) time , seism(itime,icomp)
enddo
- if(NDIMr==2)then
- seism(:,3) = seism(:,2)
- seism(:,2) = 0.d0
- else
- seism(:,2) = seism(:,1)
- seism(:,1) = 0.d0
- seism(:,3) = 0.d0
- endif
+ enddo
- close(10)
+ if(NDIMr==2)then
+ seism(:,3) = seism(:,2)
+ seism(:,2) = 0.d0
+ else
+ seism(:,2) = seism(:,1)
+ seism(:,1) = 0.d0
+ seism(:,3) = 0.d0
+ endif
+ close(10)
- istart = max(floor(tstart(irec)/deltat),1)
- iend = min(floor(tend(irec)/deltat),NSTEP)
- print*,'istart =',istart, 'iend =', iend
- print*,'tstart =',istart*deltat, 'tend =', iend*deltat
- if(istart >= iend) stop 'check istart,iend'
- nlen = iend - istart +1
- do icomp = 1, NDIM
+ istart = max(floor(tstart(irec)/deltat),1)
+ iend = min(floor(tend(irec)/deltat),NSTEP)
+ print*,'istart =',istart, 'iend =', iend
+ print*,'tstart =',istart*deltat, 'tend =', iend*deltat
+ if(istart >= iend) stop 'check istart,iend'
+ nlen = iend - istart +1
- print*,comp(icomp)
+ do icomp = 1, NDIM
- filename = 'OUTPUT_FILES/'//trim(station_name(irec))//'.AA.'// comp(icomp) // '.adj'
- open(unit = 11, file = trim(filename))
+ print*,comp(icomp)
+ filename = 'OUTPUT_FILES/'//trim(station_name(irec))//'.AA.'// comp(icomp) // '.adj'
+ open(unit = 11, file = trim(filename))
+
time_window(:) = 0.d0
seism_win(:) = seism(:,icomp)
seism_veloc(:) = 0.d0
seism_accel(:) = 0.d0
do itime =istart,iend
-! time_window(itime) = 1.d0 - cos(pi*(itime-1)/NSTEP+1)**10 ! cosine window
- time_window(itime) = 1.d0 - (2* (dble(itime) - istart)/(iend-istart) -1.d0)**2 ! Welch window
+ !time_window(itime) = 1.d0 - cos(pi*(itime-1)/NSTEP+1)**10 ! cosine window
+ time_window(itime) = 1.d0 - (2* (dble(itime) - istart)/(iend-istart) -1.d0)**2 ! Welch window
enddo
- do itime = 2,NSTEP-1
- seism_veloc(itime) = (seism_win(itime+1) - seism_win(itime-1))/(2*deltat)
- enddo
- seism_veloc(1) = (seism_win(2) - seism_win(1))/deltat
- seism_veloc(NSTEP) = (seism_win(NSTEP) - seism_win(NSTEP-1))/deltat
+ do itime = 2,NSTEP-1
+ seism_veloc(itime) = (seism_win(itime+1) - seism_win(itime-1))/(2*deltat)
+ enddo
+ seism_veloc(1) = (seism_win(2) - seism_win(1))/deltat
+ seism_veloc(NSTEP) = (seism_win(NSTEP) - seism_win(NSTEP-1))/deltat
- do itime = 2,NSTEP-1
- seism_accel(itime) = (seism_veloc(itime+1) - seism_veloc(itime-1))/(2*deltat)
- enddo
- seism_accel(1) = (seism_veloc(2) - seism_veloc(1))/deltat
- seism_accel(NSTEP) = (seism_veloc(NSTEP) - seism_veloc(NSTEP-1))/deltat
+ do itime = 2,NSTEP-1
+ seism_accel(itime) = (seism_veloc(itime+1) - seism_veloc(itime-1))/(2*deltat)
+ enddo
+ seism_accel(1) = (seism_veloc(2) - seism_veloc(1))/deltat
+ seism_accel(NSTEP) = (seism_veloc(NSTEP) - seism_veloc(NSTEP-1))/deltat
- Nnorm = deltat * sum(time_window(:) * seism_win(:) * seism_accel(:))
-! Nnorm = deltat * sum(time_window(:) * seism_veloc(:) * seism_veloc(:))
-! cross-correlation traveltime adjoint source
- if(abs(Nnorm) > EPS) then
-! ft_bar(:) = - seism_veloc(:) * time_window(:) / Nnorm
- ft_bar(:) = seism_veloc(:) * time_window(:) / Nnorm
- print*,'Norm =', Nnorm
- else
- print *, 'norm < EPS for file '
- print*,'Norm =', Nnorm
- ft_bar(:) = 0.d0
- endif
-
- do itime =1,NSTEP
- if(icomp == adj_comp) then
- write(11,*) (itime-1)*deltat - t0, ft_bar(itime)
+ ! cross-correlation traveltime adjoint source
+ Nnorm = deltat * sum(time_window(:) * seism_win(:) * seism_accel(:))
+ !Nnorm = deltat * sum(time_window(:) * seism_veloc(:) * seism_veloc(:))
+ if(abs(Nnorm) > EPS) then
+ !ft_bar(:) = -seism_veloc(:) * time_window(:) / Nnorm
+ ft_bar(:) = seism_veloc(:) * time_window(:) / Nnorm
+ print*,'Norm =', Nnorm
else
- write(11,*) (itime-1)*deltat - t0, 0.d0
+ print *, 'norm < EPS for file '
+ print*,'Norm =', Nnorm
+ ft_bar(:) = 0.d0
endif
- enddo
+ do itime =1,NSTEP
+ if(icomp == adj_comp) then
+ write(11,*) (itime-1)*deltat - t0, ft_bar(itime)
+ else
+ write(11,*) (itime-1)*deltat - t0, 0.d0
+ endif
enddo
- close(11)
- enddo
- print*,'*************************'
- print*,'The input files (S****.AA.BHX/BHY/BHZ.adj) needed to run the adjoint simulation are in OUTPUT_FILES'
- print*,'*************************'
+ enddo
+ close(11)
- end program adj_seismogram
+ enddo
+ print*,'*************************'
+ print*,'The input files (S****.AA.BHX/BHY/BHZ.adj) needed to run the adjoint simulation are in OUTPUT_FILES'
+ print*,'*************************'
+
+end program adj_seismogram
More information about the CIG-COMMITS
mailing list