[cig-commits] [commit] devel, master: added kernel example for Tape2007 (b6a1334)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Wed Jun 18 15:18:44 PDT 2014
Repository : https://github.com/geodynamics/specfem2d
On branches: devel,master
Link : https://github.com/geodynamics/specfem2d/compare/fc67e6fd7ad890705b2b72b4b3c509accb22249e...e9ca46c40131588d89d7b0883250bc6584ce6b4c
>---------------------------------------------------------------
commit b6a13342c09600ab9596193a8c01beabf3f16760
Author: Carl Tape <carltape at gi.alaska.edu>
Date: Mon Apr 11 23:33:17 2011 +0000
added kernel example for Tape2007
>---------------------------------------------------------------
b6a13342c09600ab9596193a8c01beabf3f16760
.../Par_file_Tape2007_onerec | 4 +-
Tape2007_kernel/README | 35 ++++
{Tape2007 => Tape2007_kernel}/SOURCE_001 | 0
Tape2007_kernel/adj_seismogram_Tape2007.f90 | 179 +++++++++++++++++++++
.../interfaces_Tape2007.dat | 0
Tape2007_kernel/kernel_Tape2007_kernel_onerec.pdf | Bin 0 -> 9849 bytes
{Tape2007 => Tape2007_kernel}/process.sh | 0
Tape2007_kernel/process_kernel.sh | 24 +++
8 files changed, 240 insertions(+), 2 deletions(-)
diff --git a/Tape2007/Par_file_Tape2007_onerec b/Tape2007_kernel/Par_file_Tape2007_onerec
similarity index 97%
copy from Tape2007/Par_file_Tape2007_onerec
copy to Tape2007_kernel/Par_file_Tape2007_onerec
index 1ac4c9c..14a518d 100644
--- a/Tape2007/Par_file_Tape2007_onerec
+++ b/Tape2007_kernel/Par_file_Tape2007_onerec
@@ -3,7 +3,7 @@ title = Tape-Liu-Tromp (GJI 2007)
# forward or adjoint simulation
SIMULATION_TYPE = 1 # 1 = forward, 2 = adjoint + kernels
-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
@@ -64,7 +64,7 @@ sizemax_arrows = 1.d0 # maximum size of arrows on vec
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 = .true. # 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
diff --git a/Tape2007_kernel/README b/Tape2007_kernel/README
new file mode 100644
index 0000000..73b8e56
--- /dev/null
+++ b/Tape2007_kernel/README
@@ -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
+
+---------------------------
diff --git a/Tape2007/SOURCE_001 b/Tape2007_kernel/SOURCE_001
similarity index 100%
copy from Tape2007/SOURCE_001
copy to Tape2007_kernel/SOURCE_001
diff --git a/Tape2007_kernel/adj_seismogram_Tape2007.f90 b/Tape2007_kernel/adj_seismogram_Tape2007.f90
new file mode 100644
index 0000000..0a03b9c
--- /dev/null
+++ b/Tape2007_kernel/adj_seismogram_Tape2007.f90
@@ -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
diff --git a/Tape2007/interfaces_Tape2007.dat b/Tape2007_kernel/interfaces_Tape2007.dat
similarity index 100%
copy from Tape2007/interfaces_Tape2007.dat
copy to Tape2007_kernel/interfaces_Tape2007.dat
diff --git a/Tape2007_kernel/kernel_Tape2007_kernel_onerec.pdf b/Tape2007_kernel/kernel_Tape2007_kernel_onerec.pdf
new file mode 100644
index 0000000..23d6fb7
Binary files /dev/null and b/Tape2007_kernel/kernel_Tape2007_kernel_onerec.pdf differ
diff --git a/Tape2007/process.sh b/Tape2007_kernel/process.sh
similarity index 100%
copy from Tape2007/process.sh
copy to Tape2007_kernel/process.sh
diff --git a/Tape2007_kernel/process_kernel.sh b/Tape2007_kernel/process_kernel.sh
new file mode 100755
index 0000000..0b89c41
--- /dev/null
+++ b/Tape2007_kernel/process_kernel.sh
@@ -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`
More information about the CIG-COMMITS
mailing list