[cig-commits] [commit] devel, master: added kernel example for Tromp2005 (23b1dfd)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Wed Jun 18 15:18:46 PDT 2014
Repository : https://github.com/geodynamics/specfem2d
On branches: devel,master
Link : https://github.com/geodynamics/specfem2d/compare/fc67e6fd7ad890705b2b72b4b3c509accb22249e...e9ca46c40131588d89d7b0883250bc6584ce6b4c
>---------------------------------------------------------------
commit 23b1dfdf7de9e0b75beb69241d1acf8a82b2b381
Author: Carl Tape <carltape at gi.alaska.edu>
Date: Tue Apr 12 00:05:16 2011 +0000
added kernel example for Tromp2005
>---------------------------------------------------------------
23b1dfdf7de9e0b75beb69241d1acf8a82b2b381
Tape2007_kernel/adj_seismogram_Tape2007.f90 | 197 ++++++++++-----------
{Tromp2005 => Tromp2005_kernel}/Par_file_Tromp2005 | 4 +-
{Tape2007_kernel => Tromp2005_kernel}/README | 6 +-
{Tromp2005 => Tromp2005_kernel}/SOURCE_Tromp2005 | 0
Tromp2005_kernel/adj_seismogram_Tromp2005.f90 | 176 ++++++++++++++++++
.../interfaces_Tromp2005.dat | 0
Tromp2005_kernel/kernel_Tromp2005_kernel_PSV.pdf | Bin 0 -> 15185 bytes
{Tromp2005 => Tromp2005_kernel}/process.sh | 0
.../process_kernel.sh | 0
9 files changed, 278 insertions(+), 105 deletions(-)
diff --git a/Tape2007_kernel/adj_seismogram_Tape2007.f90 b/Tape2007_kernel/adj_seismogram_Tape2007.f90
index 0a03b9c..64e0522 100644
--- a/Tape2007_kernel/adj_seismogram_Tape2007.f90
+++ b/Tape2007_kernel/adj_seismogram_Tape2007.f90
@@ -42,88 +42,85 @@
!
!========================================================================
- 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
-!
-! rm -rf xadj_seismogram ; gfortran adj_seismogram.f90 -o xadj_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
+ 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, 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"/)
-! 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'
+ ! 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
+ ! KEY: 'absolute' time interval for window
+ tstart(1) = 95.d0 + t0
+ tend(1) = 130.d0 + t0
-! chose the component for the adjoint source (adj_comp = 1:X, 2:Y, 3:Z)
- adj_comp = 2
-!!!!
+ ! chose the component for the adjoint source (adj_comp = 1:X, 2:Y, 3:Z)
+ adj_comp = 2
- do irec =1,nrec
+ do irec =1,nrec
- do icomp = 1, NDIMr
+ 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
+ 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
+
+ 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)
+ 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
+ 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
+ do icomp = 1, NDIM
- print*,comp(icomp)
+ print*,comp(icomp)
- filename = 'OUTPUT_FILES/'//trim(station_name(irec))//'.AA.'// comp(icomp) // '.adj'
- open(unit = 11, file = trim(filename))
+ 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)
@@ -131,49 +128,49 @@
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_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)
+ 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
+
+ ! 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)
+
+ 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
+end program adj_seismogram
diff --git a/Tromp2005/Par_file_Tromp2005 b/Tromp2005_kernel/Par_file_Tromp2005
similarity index 97%
copy from Tromp2005/Par_file_Tromp2005
copy to Tromp2005_kernel/Par_file_Tromp2005
index 86b134c..83461c2 100644
--- a/Tromp2005/Par_file_Tromp2005
+++ b/Tromp2005_kernel/Par_file_Tromp2005
@@ -3,7 +3,7 @@ title = Tromp-Tape-Liu (GJI 2005)
# 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
@@ -63,7 +63,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/Tromp2005_kernel/README
similarity index 82%
copy from Tape2007_kernel/README
copy to Tromp2005_kernel/README
index 73b8e56..a01b012 100644
--- a/Tape2007_kernel/README
+++ b/Tromp2005_kernel/README
@@ -2,7 +2,7 @@
README
----------------------------------------------------------------------
-Kernel example for Tape-Liu-Tromp (GJI 2007).
+Kernel example for Tromp-Tape-Liu (GJI 2005).
TO RUN:
@@ -17,11 +17,11 @@ TO RUN:
> make all
4. run mesher and solver for forward wavefield:
- > cd EXAMPLES/Tape2007_kernel/
+ > cd EXAMPLES/Tromp2005_kernel/
> ./process.sh
5. compute adjoint source:
- > rm -rf xadj_seismogram ; gfortran adj_seismogram_Tape2007.f90 -o xadj_seismogram
+ > rm -rf xadj_seismogram ; gfortran adj_seismogram_Tromp2005.f90 -o xadj_seismogram
> xadj_seismogram
6. change Par_file with save_forward = .false. and SIMULATION_TYPE = 2
diff --git a/Tromp2005/SOURCE_Tromp2005 b/Tromp2005_kernel/SOURCE_Tromp2005
similarity index 100%
copy from Tromp2005/SOURCE_Tromp2005
copy to Tromp2005_kernel/SOURCE_Tromp2005
diff --git a/Tromp2005_kernel/adj_seismogram_Tromp2005.f90 b/Tromp2005_kernel/adj_seismogram_Tromp2005.f90
new file mode 100644
index 0000000..bec1fa5
--- /dev/null
+++ b/Tromp2005_kernel/adj_seismogram_Tromp2005.f90
@@ -0,0 +1,176 @@
+
+!========================================================================
+!
+! 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 = 8.0 ! labeled as 'time delay'
+ double precision, parameter :: deltat = 0.02
+ 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) = 27.d0 + t0
+ tend(1) = 32.d0 + t0
+
+ ! chose the component for the adjoint source (adj_comp = 1:X, 2:Y, 3:Z)
+ adj_comp = 1
+
+ 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
+
+ ! 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
+ 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/Tromp2005/interfaces_Tromp2005.dat b/Tromp2005_kernel/interfaces_Tromp2005.dat
similarity index 100%
copy from Tromp2005/interfaces_Tromp2005.dat
copy to Tromp2005_kernel/interfaces_Tromp2005.dat
diff --git a/Tromp2005_kernel/kernel_Tromp2005_kernel_PSV.pdf b/Tromp2005_kernel/kernel_Tromp2005_kernel_PSV.pdf
new file mode 100644
index 0000000..40defa2
Binary files /dev/null and b/Tromp2005_kernel/kernel_Tromp2005_kernel_PSV.pdf differ
diff --git a/Tromp2005/process.sh b/Tromp2005_kernel/process.sh
similarity index 100%
copy from Tromp2005/process.sh
copy to Tromp2005_kernel/process.sh
diff --git a/Tape2007_kernel/process_kernel.sh b/Tromp2005_kernel/process_kernel.sh
similarity index 100%
copy from Tape2007_kernel/process_kernel.sh
copy to Tromp2005_kernel/process_kernel.sh
More information about the CIG-COMMITS
mailing list