[cig-commits] r18976 - in seismo/2D/SPECFEM2D/trunk: EXAMPLES/noise_uniform src/specfem2D
rmodrak at geodynamics.org
rmodrak at geodynamics.org
Tue Sep 27 11:38:04 PDT 2011
Author: rmodrak
Date: 2011-09-27 11:38:03 -0700 (Tue, 27 Sep 2011)
New Revision: 18976
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/src/specfem2D/noise_tomography.f90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
Log:
for noise tomography, adds flags "save_everywhere" and "output_wavefields_noise"
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-27 17:54:55 UTC (rev 18975)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_1 2011-09-27 18:38:03 UTC (rev 18976)
@@ -22,7 +22,7 @@
p_sv = .false. # set the type of calculation (P-SV or SH/membrane waves)
# time step parameters
-nt = 2000 # total number of time steps
+nt = 2500 # total number of time steps
deltat = 4.d-2 # duration of a time step
# source parameters
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-27 17:54:55 UTC (rev 18975)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_2 2011-09-27 18:38:03 UTC (rev 18976)
@@ -22,7 +22,7 @@
p_sv = .false. # set the type of calculation (P-SV or SH/membrane waves)
# time step parameters
-nt = 2000 # total number of time steps
+nt = 2500 # total number of time steps
deltat = 4.d-2 # duration of a time step
# source parameters
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-27 17:54:55 UTC (rev 18975)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_3 2011-09-27 18:38:03 UTC (rev 18976)
@@ -22,7 +22,7 @@
p_sv = .false. # set the type of calculation (P-SV or SH/membrane waves)
# time step parameters
-nt = 2000 # total number of time steps
+nt = 2500 # total number of time steps
deltat = 4.d-2 # duration of a time step
# source parameters
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/noise_tomography.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/noise_tomography.f90 2011-09-27 17:54:55 UTC (rev 18975)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/noise_tomography.f90 2011-09-27 18:38:03 UTC (rev 18976)
@@ -81,7 +81,7 @@
! =============================================================================================================
-! read noise parameters and check for consitency
+! read noise parameters and check for consistency
subroutine read_parameters_noise(NOISE_TOMOGRAPHY,SIMULATION_TYPE,SAVE_FORWARD, &
any_acoustic,any_poroelastic,p_sv,irec_master, &
Mxx,Mxz,Mzz,factor,NSOURCES, &
@@ -413,22 +413,32 @@
! =============================================================================================================
! save a snapshot of the "generating wavefield" eta that will be used to drive
! the "ensemble forward wavefield"
- subroutine save_surface_movie_noise(p_sv,it,nglob,displ_elastic)
+ subroutine save_surface_movie_noise(NOISE_TOMOGRAPHY,p_sv,it,nglob,displ_elastic)
implicit none
include "constants.h"
!input paramters
+ integer :: NOISE_TOMOGRAPHY
logical :: p_sv
integer :: it, nglob
real(kind=CUSTOM_REAL), dimension(3,nglob) :: displ_elastic
!local parameters
character(len=60) file_out_noise
+
+ if (NOISE_TOMOGRAPHY == 1) then
+ write(file_out_noise,"('OUTPUT_FILES/NOISE_TOMOGRAPHY/eta_',i6.6)") it
- write(file_out_noise,"('eta_',i6.6)") it
- open(unit=500,file='OUTPUT_FILES/NOISE_TOMOGRAPHY/'//trim(file_out_noise),status='unknown',form='unformatted',action='write')
+ elseif (NOISE_TOMOGRAPHY == 2) then
+ write(file_out_noise,"('OUTPUT_FILES/NOISE_TOMOGRAPHY/phi_',i6.6)") it
+ else
+ call exit_mpi('Bad value of NOISE_TOMOGRAPHY in save_surface_movie_noise.')
+
+ endif
+
+ open(unit=500,file=trim(file_out_noise),status='unknown',form='unformatted',action='write')
if(p_sv) then
write(500) displ_elastic(1,:)
write(500) displ_elastic(3,:)
@@ -439,15 +449,14 @@
end subroutine save_surface_movie_noise
! =============================================================================================================
-! auxillary routine
- subroutine snapshot_all(ncol,nglob,filename,array_all)
+ subroutine snapshots_noise(ncol,nglob,filename,array_all)
implicit none
include "constants.h"
!input paramters
integer :: ncol,nglob
- character(len=100) filename
+ character(len=512) filename
real(kind=CUSTOM_REAL), dimension(ncol,nglob) :: array_all
@@ -468,7 +477,7 @@
close(504)
- end subroutine snapshot_all
+ end subroutine snapshots_noise
! =============================================================================================================
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2011-09-27 17:54:55 UTC (rev 18975)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2011-09-27 18:38:03 UTC (rev 18976)
@@ -799,18 +799,23 @@
!>SU_FORMAT
!<NOISE_TOMOGRAPHY
- ! NOISE_TOMOGRAPHY = 0 - turn noise tomography subroutines off, resulting in
- ! an earthquake simulation rather than a noise simulation
+ ! NOISE_TOMOGRAPHY = 0 - turn noise tomography subroutines off; setting
+ ! NOISE_TOMOGRAPHY equal to 0, in other words, results in an earthquake
+ ! simulation rather than a noise simulation
- ! NOISE_TOMOGRAPHY = 1 - compute "generating wavefield" and store the result
+ ! NOISE_TOMOGRAPHY = 1 - compute "generating" wavefield and store the result;
+ ! this stored wavefield is then used to compute the "ensemble forward"
+ ! wavefield in the next noise simulation
- ! NOISE_TOMOGRAPHY = 2 - compute "ensemble forward wavefield"; if an adjoint
- ! simulation is planned, users need to manually extract cross-correlograms
+ ! NOISE_TOMOGRAPHY = 2 - compute "ensemble forward" wavefield and store the
+ ! result; if an adjoint simulation is planned, users need to store
+ ! seismograms (actually, cross-correlograms) for later processing
- ! NOISE_TOMOGRAPHY = 3 - carry out adjoint simulation; for noise tomography
- ! applications, users need to supply adjoint source(s) based on cross-
- ! -correlograms from previous simulation
+ ! NOISE_TOMOGRAPHY = 3 - carry out adjoint simulation; users need to supply
+ ! adjoint sources constructed from cross-correlograms computed during the
+ ! "ensemble forward" step
+
! For an explanation of terms and concepts in noise tomography, see "Tromp et
! al., 2011, Noise Cross-Correlation Sensitivity Kernels, Geophysical Journal
! International"
@@ -821,15 +826,28 @@
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 dimensions depending on SH/P_SV, use separate arrays for x,y,z
+
+ ! The following three arrays are used to hold snapshots of the generating
+ ! wavefield or of the ensemble forward wavefield, depending on the type of
+ ! noise simulation specified. In some cases, the entire generating wavefield
+ ! or ensemble forward wavefield needs to be saved for all times steps. Since
+ ! the disk space required to do this is usually quite large, separate arrays
+ ! are used for x,y,z to avoid having empty dimensions (one empty dimension in
+ ! the case of P-SV or two empty dimensions in the case of SH).
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) :: noise_file
+ ! For writing noise wavefields
+ logical :: output_wavefields_noise = .true.
+ integer :: noise_output_ncol
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: noise_output_array
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: noise_output_dim_5
+ character(len=512) :: noise_output_file
+ ! For noise tomography only - specify whether to reconstruct ensemble forward
+ ! wavefield by saving everywhere or by saving only at the boundaries (the
+ ! latter usually much faster but prone to artefacts)
+ logical :: save_everywhere = .true.
!>NOISE_TOMOGRAPHY
@@ -3605,8 +3623,9 @@
call compute_source_array_noise(p_sv,NSTEP,deltat,nglob,ibool,ispec_noise, &
xi_noise,gamma_noise,xigll,zigll, &
time_function_noise,source_array_noise)
- !write out local coordinates of mesh
- open(unit=504,file='OUTPUT_FILES/elem',status='unknown',action='write')
+
+ !write out coordinates of mesh
+ open(unit=504,file='OUTPUT_FILES/mesh_spec',status='unknown',action='write')
do ispec = 1, nspec
do j = 1, NGLLZ
do i = 1, NGLLX
@@ -3618,6 +3637,13 @@
enddo
close(504)
+ open(unit=504,file='OUTPUT_FILES/mesh_glob',status='unknown',action='write')
+ do iglob = 1, nglob
+ write(504,'(1pe11.3,1pe11.3,2i3,i7)') &
+ coord(1,iglob), coord(2,iglob), iglob
+ enddo
+ close(504)
+
!write out spatial distribution of noise sources
call create_mask_noise(nglob,coord,mask_noise)
open(unit=504,file='OUTPUT_FILES/mask_noise',status='unknown',action='write')
@@ -3662,13 +3688,16 @@
call create_mask_noise(nglob,coord,mask_noise)
elseif (NOISE_TOMOGRAPHY == 3) then
- call create_mask_noise(nglob,coord,mask_noise)
- !prepare array that will hold wavefield snapshots
- ncol = 5
- allocate(noise_all(ncol,nglob))
- allocate(rho_klglob(nglob))
+ if (output_wavefields_noise) then
+ call create_mask_noise(nglob,coord,mask_noise)
+ !prepare array that will hold wavefield snapshots
+ noise_output_ncol = 5
+ allocate(noise_output_array(noise_output_ncol,nglob))
+ allocate(noise_output_dim_5(nglob))
+ endif
+
endif
!>NOISE_TOMOGRAPHY
@@ -5037,10 +5066,11 @@
surface_movie_z_noise,mask_noise,jacobian,wxgll,wzgll)
elseif (NOISE_TOMOGRAPHY == 3) then
- call add_surface_movie_noise(p_sv,it,NSTEP,nspec,nglob,ibool,b_accel_elastic, &
- surface_movie_x_noise,surface_movie_y_noise, &
- surface_movie_z_noise,mask_noise,jacobian,wxgll,wzgll)
-
+ if (.not. save_everywhere) then
+ call add_surface_movie_noise(p_sv,it,NSTEP,nspec,nglob,ibool,b_accel_elastic, &
+ surface_movie_x_noise,surface_movie_y_noise, &
+ surface_movie_z_noise,mask_noise,jacobian,wxgll,wzgll)
+ endif
endif
!>NOISE_TOMOGRAPHY
@@ -5905,6 +5935,33 @@
endif ! if(it == 1 .and. SIMULATION_TYPE == 2)
+!<NOISE_TOMOGRAPHY
+
+ if ( NOISE_TOMOGRAPHY == 1 ) then
+ call save_surface_movie_noise(NOISE_TOMOGRAPHY,p_sv,it,nglob,displ_elastic)
+
+ elseif ( NOISE_TOMOGRAPHY == 2 .and. save_everywhere ) then
+ call save_surface_movie_noise(NOISE_TOMOGRAPHY,p_sv,it,nglob,displ_elastic)
+
+ elseif ( NOISE_TOMOGRAPHY == 3 .and. save_everywhere ) then
+ write(noise_output_file,"('phi_',i6.6)") NSTEP-it+1
+ open(unit=500,file='OUTPUT_FILES/NOISE_TOMOGRAPHY/'//trim(noise_output_file), &
+ status='old',form='unformatted',action='read',iostat=ios)
+ if( ios /= 0) write(*,*) 'Error retrieving ensemble forward wavefield.'
+ if(p_sv) then
+ read(500) b_displ_elastic(1,:)
+ read(500) b_displ_elastic(3,:)
+ else
+ read(500) b_displ_elastic(2,:)
+ endif
+ close(500)
+
+ endif
+
+
+
+!>NOISE_TOMOGRAPHY
+
! ********************************************************************************************
! kernels calculation
! ********************************************************************************************
@@ -6426,13 +6483,7 @@
endif ! if(SIMULATION_TYPE == 2)
-!<NOISE_TOMOGRAPHY
- if ( NOISE_TOMOGRAPHY == 1 ) then
- call save_surface_movie_noise(p_sv,it,nglob,displ_elastic)
- endif
-!>NOISE_TOMOGRAPHY
-
!
!---- display results at given time steps
!
@@ -6518,22 +6569,28 @@
!<NOISE_TOMOGRAPHY
- if (NOISE_TOMOGRAPHY == 1) then
+ if ( NOISE_TOMOGRAPHY == 3 .and. output_wavefields_noise ) then
- elseif (NOISE_TOMOGRAPHY == 2) then
+ !load ensemble foward source
+ write(noise_output_file,"('phi_',i6.6)") NSTEP-it+1
+ open(unit=500,file='OUTPUT_FILES/NOISE_TOMOGRAPHY/'//trim(noise_output_file), &
+ status='old',form='unformatted',action='read',iostat=ios)
+ if( ios /= 0) write(*,*) 'Error preparing noise output.'
+ read(500) surface_movie_y_noise
+ close(500)
- elseif (NOISE_TOMOGRAPHY == 3) then
- call elem_to_glob(nspec,nglob,ibool,rho_kl,rho_klglob)
+ !load product of fwd, adj wavefields
+ call elem_to_glob(nspec,nglob,ibool,rho_kl,noise_output_dim_5)
- noise_all(1,:) = surface_movie_y_noise
- noise_all(2,:) = b_displ_elastic(2,:)
- noise_all(3,:) = accel_elastic(2,:)
- noise_all(4,:) = rho_k
- noise_all(5,:) = rho_klglob
+ !write text file
+ noise_output_array(1,:) = surface_movie_y_noise(:) * mask_noise(:)
+ noise_output_array(2,:) = b_displ_elastic(2,:)
+ noise_output_array(3,:) = accel_elastic(2,:)
+ noise_output_array(4,:) = rho_k(:)
+ noise_output_array(5,:) = noise_output_dim_5(:)
+ write(noise_output_file,"('OUTPUT_FILES/snapshot_all_',i6.6)") it
+ call snapshots_noise(noise_output_ncol,nglob,noise_output_file,noise_output_array)
- write(noise_file,"('OUTPUT_FILES/snapshot_all_',i6.6)") it
- call snapshot_all(ncol,nglob,noise_file,noise_all)
-
endif
!>NOISE_TOMOGRAPHY
More information about the CIG-COMMITS
mailing list