[cig-commits] r17236 - seismo/3D/SPECFEM3D/trunk
yangl at geodynamics.org
yangl at geodynamics.org
Mon Oct 4 13:20:38 PDT 2010
Author: yangl
Date: 2010-10-04 13:20:38 -0700 (Mon, 04 Oct 2010)
New Revision: 17236
Modified:
seismo/3D/SPECFEM3D/trunk/check_buffers_2D.f90
seismo/3D/SPECFEM3D/trunk/combine_AVS_DX.f90
seismo/3D/SPECFEM3D/trunk/combine_vol_data.f90
seismo/3D/SPECFEM3D/trunk/compute_add_sources_acoustic.f90
seismo/3D/SPECFEM3D/trunk/compute_add_sources_elastic.f90
seismo/3D/SPECFEM3D/trunk/compute_arrays_source.f90
seismo/3D/SPECFEM3D/trunk/compute_forces_acoustic.f90
seismo/3D/SPECFEM3D/trunk/compute_forces_elastic.f90
seismo/3D/SPECFEM3D/trunk/create_header_file.f90
seismo/3D/SPECFEM3D/trunk/create_movie_shakemap_AVS_DX_GMT.f90
seismo/3D/SPECFEM3D/trunk/generate_databases.f90
seismo/3D/SPECFEM3D/trunk/initialize_simulation.f90
seismo/3D/SPECFEM3D/trunk/read_parameter_file.f90
seismo/3D/SPECFEM3D/trunk/setup_sources_receivers.f90
seismo/3D/SPECFEM3D/trunk/specfem3D_par.f90
Log:
implementing NTSTEP_BETWEEN_READ_ADJSRC
Modified: seismo/3D/SPECFEM3D/trunk/check_buffers_2D.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/check_buffers_2D.f90 2010-10-04 09:49:07 UTC (rev 17235)
+++ seismo/3D/SPECFEM3D/trunk/check_buffers_2D.f90 2010-10-04 20:20:38 UTC (rev 17236)
@@ -59,7 +59,7 @@
integer NER_SEDIM,NER_BASEMENT_SEDIM,NER_16_BASEMENT, &
NER_MOHO_16,NER_BOTTOM_MOHO,NEX_XI,NEX_ETA, &
NPROC_XI,NPROC_ETA,NTSTEP_BETWEEN_OUTPUT_SEISMOS,NSTEP,UTM_PROJECTION_ZONE,SIMULATION_TYPE
- integer NSOURCES
+ integer NSOURCES, NTSTEP_BETWEEN_READ_ADJSRC
double precision UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK
double precision DT,LATITUDE_MIN,LATITUDE_MAX,LONGITUDE_MIN,LONGITUDE_MAX,HDUR_MOVIE
@@ -107,7 +107,8 @@
MOVIE_SURFACE,MOVIE_VOLUME,CREATE_SHAKEMAP,SAVE_DISPLACEMENT, &
NTSTEP_BETWEEN_FRAMES,USE_HIGHRES_FOR_MOVIES,HDUR_MOVIE, &
SAVE_MESH_FILES,PRINT_SOURCE_TIME_FUNCTION, &
- NTSTEP_BETWEEN_OUTPUT_INFO,SUPPRESS_UTM_PROJECTION,MODEL,USE_REGULAR_MESH,SIMULATION_TYPE,SAVE_FORWARD)
+ NTSTEP_BETWEEN_OUTPUT_INFO,SUPPRESS_UTM_PROJECTION,MODEL,USE_REGULAR_MESH,SIMULATION_TYPE,SAVE_FORWARD, &
+ NTSTEP_BETWEEN_READ_ADJSRC)
! compute other parameters based upon values read
call compute_parameters(NER,NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA, &
Modified: seismo/3D/SPECFEM3D/trunk/combine_AVS_DX.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/combine_AVS_DX.f90 2010-10-04 09:49:07 UTC (rev 17235)
+++ seismo/3D/SPECFEM3D/trunk/combine_AVS_DX.f90 2010-10-04 20:20:38 UTC (rev 17236)
@@ -85,7 +85,7 @@
integer NER_SEDIM,NER_BASEMENT_SEDIM,NER_16_BASEMENT, &
NER_MOHO_16,NER_BOTTOM_MOHO,NEX_XI,NEX_ETA, &
NPROC_XI,NPROC_ETA,NTSTEP_BETWEEN_OUTPUT_SEISMOS,NSTEP,UTM_PROJECTION_ZONE,SIMULATION_TYPE
- integer NSOURCES
+ integer NSOURCES,NTSTEP_BETWEEN_READ_ADJSRC
double precision UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK
double precision DT,LATITUDE_MIN,LATITUDE_MAX,LONGITUDE_MIN,LONGITUDE_MAX,HDUR_MOVIE
@@ -139,7 +139,8 @@
MOVIE_SURFACE,MOVIE_VOLUME,CREATE_SHAKEMAP,SAVE_DISPLACEMENT, &
NTSTEP_BETWEEN_FRAMES,USE_HIGHRES_FOR_MOVIES,HDUR_MOVIE, &
SAVE_MESH_FILES,PRINT_SOURCE_TIME_FUNCTION, &
- NTSTEP_BETWEEN_OUTPUT_INFO,SUPPRESS_UTM_PROJECTION,MODEL,USE_REGULAR_MESH,SIMULATION_TYPE,SAVE_FORWARD)
+ NTSTEP_BETWEEN_OUTPUT_INFO,SUPPRESS_UTM_PROJECTION,MODEL,USE_REGULAR_MESH,SIMULATION_TYPE,SAVE_FORWARD, &
+ NTSTEP_BETWEEN_READ_ADJSRC)
if(.not. SAVE_MESH_FILES) stop 'AVS or DX files were not saved by the mesher'
Modified: seismo/3D/SPECFEM3D/trunk/combine_vol_data.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/combine_vol_data.f90 2010-10-04 09:49:07 UTC (rev 17235)
+++ seismo/3D/SPECFEM3D/trunk/combine_vol_data.f90 2010-10-04 20:20:38 UTC (rev 17236)
@@ -65,7 +65,7 @@
double precision :: HDUR_MOVIE
integer :: NPROC,NTSTEP_BETWEEN_OUTPUT_SEISMOS,NSTEP, &
UTM_PROJECTION_ZONE,SIMULATION_TYPE
- integer :: NSOURCES
+ integer :: NSOURCES,NTSTEP_BETWEEN_READ_ADJSRC
integer :: NTSTEP_BETWEEN_FRAMES,NTSTEP_BETWEEN_OUTPUT_INFO
logical :: MOVIE_SURFACE,MOVIE_VOLUME,CREATE_SHAKEMAP,SAVE_DISPLACEMENT, &
USE_HIGHRES_FOR_MOVIES,SUPPRESS_UTM_PROJECTION
@@ -149,7 +149,8 @@
MOVIE_SURFACE,MOVIE_VOLUME,CREATE_SHAKEMAP,SAVE_DISPLACEMENT, &
NTSTEP_BETWEEN_FRAMES,USE_HIGHRES_FOR_MOVIES,HDUR_MOVIE, &
SAVE_MESH_FILES,PRINT_SOURCE_TIME_FUNCTION, &
- NTSTEP_BETWEEN_OUTPUT_INFO,SIMULATION_TYPE,SAVE_FORWARD)
+ NTSTEP_BETWEEN_OUTPUT_INFO,SIMULATION_TYPE,SAVE_FORWARD, &
+ NTSTEP_BETWEEN_READ_ADJSRC)
print *, 'Slice list: '
print *, node_list(1:num_node)
Modified: seismo/3D/SPECFEM3D/trunk/compute_add_sources_acoustic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/compute_add_sources_acoustic.f90 2010-10-04 09:49:07 UTC (rev 17235)
+++ seismo/3D/SPECFEM3D/trunk/compute_add_sources_acoustic.f90 2010-10-04 20:20:38 UTC (rev 17236)
@@ -33,9 +33,12 @@
sourcearrays,kappastore,ispec_is_acoustic,&
SIMULATION_TYPE,NSTEP,NGLOB_ADJOINT, &
nrec,islice_selected_rec,ispec_selected_rec, &
- nadj_rec_local,adj_sourcearrays,b_potential_dot_dot_acoustic )
+ nadj_rec_local,adj_sourcearrays,b_potential_dot_dot_acoustic, &
+ NTSTEP_BETWEEN_READ_ADJSRC )
- use specfem_par,only: PRINT_SOURCE_TIME_FUNCTION,stf_used_total
+ use specfem_par,only: PRINT_SOURCE_TIME_FUNCTION,stf_used_total, &
+ xigll,yigll,zigll,xi_receiver,eta_receiver,gamma_receiver,&
+ station_name,network_name,adj_source_file
implicit none
include "constants.h"
@@ -72,8 +75,14 @@
integer:: nrec
integer,dimension(nrec) :: islice_selected_rec,ispec_selected_rec
integer:: nadj_rec_local
- real(kind=CUSTOM_REAL),dimension(nadj_rec_local,NSTEP,NDIM,NGLLX,NGLLY,NGLLZ):: adj_sourcearrays
real(kind=CUSTOM_REAL),dimension(NGLOB_ADJOINT):: b_potential_dot_dot_acoustic
+!<YANGL
+! real(kind=CUSTOM_REAL),dimension(nadj_rec_local,NSTEP,NDIM,NGLLX,NGLLY,NGLLZ):: adj_sourcearrays
+ real(kind=CUSTOM_REAL),dimension(nadj_rec_local,NTSTEP_BETWEEN_READ_ADJSRC,NDIM,NGLLX,NGLLY,NGLLZ):: adj_sourcearrays
+ integer :: it_sub_adj,itime,NTSTEP_BETWEEN_READ_ADJSRC
+ real(kind=CUSTOM_REAL),dimension(NTSTEP_BETWEEN_READ_ADJSRC,NDIM,NGLLX,NGLLY,NGLLZ):: adj_sourcearray
+ logical :: ibool_read_adj_arrays
+!>YANGL
! local parameters
double precision :: f0
@@ -208,6 +217,40 @@
! adjoint simulations
if (SIMULATION_TYPE == 2 .or. SIMULATION_TYPE == 3) then
+!<YANGL
+! read in adjoint sources block by block (for memory consideration)
+! e.g., in exploration experiments, both the number of receivers (nrec) and the number of time steps (NSTEP) are huge,
+! which may cause problems since we have a large array: adj_sourcearrays(nadj_rec_local,NSTEP,NDIM,NGLLX,NGLLY,NGLLZ)
+
+ ! figure out if we need to read in a chunk of the adjoint source at this timestep
+ it_sub_adj = ceiling( dble(it)/dble(NTSTEP_BETWEEN_READ_ADJSRC) ) !chunk_number
+ ibool_read_adj_arrays = (((mod(it-1,NTSTEP_BETWEEN_READ_ADJSRC) == 0)) .and. (nadj_rec_local > 0))
+
+ ! needs to read in a new chunk/block of the adjoint source
+ if(ibool_read_adj_arrays) then
+
+ irec_local = 0
+ do irec = 1, nrec
+ ! compute source arrays
+ if(myrank == islice_selected_rec(irec)) then
+ irec_local = irec_local + 1
+
+ ! reads in **sta**.**net**.**LH**.adj files
+ adj_source_file = trim(station_name(irec))//'.'//trim(network_name(irec))
+ call compute_arrays_adjoint_source(myrank,adj_source_file, &
+ xi_receiver(irec),eta_receiver(irec),gamma_receiver(irec), &
+ adj_sourcearray, xigll,yigll,zigll, &
+ it_sub_adj,NSTEP,NTSTEP_BETWEEN_READ_ADJSRC)
+ do itime = 1,NTSTEP_BETWEEN_READ_ADJSRC
+ adj_sourcearrays(irec_local,itime,:,:,:,:) = adj_sourcearray(itime,:,:,:,:)
+ enddo
+
+ endif
+ enddo
+
+ endif ! if(ibool_read_adj_arrays)
+!>YANGL
+
if( it < NSTEP ) then
! receivers act as sources
irec_local = 0
@@ -225,8 +268,12 @@
! beware, for acoustic medium, source is: pressure divided by Kappa of the fluid
! note: it takes the first component of the adj_sourcearrays
! the idea is to have e.g. a pressure source, where all 3 components would be the same
+!<YANGL
+! potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) &
+! - adj_sourcearrays(irec_local,NSTEP-it+1,1,i,j,k) / kappastore(i,j,k,ispec)
potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) &
- - adj_sourcearrays(irec_local,NSTEP-it+1,1,i,j,k) / kappastore(i,j,k,ispec)
+ - adj_sourcearrays(irec_local,NTSTEP_BETWEEN_READ_ADJSRC-mod(it-1,NTSTEP_BETWEEN_READ_ADJSRC),1,i,j,k) / kappastore(i,j,k,ispec)
+!>YANGL
enddo
enddo
enddo
Modified: seismo/3D/SPECFEM3D/trunk/compute_add_sources_elastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/compute_add_sources_elastic.f90 2010-10-04 09:49:07 UTC (rev 17235)
+++ seismo/3D/SPECFEM3D/trunk/compute_add_sources_elastic.f90 2010-10-04 20:20:38 UTC (rev 17236)
@@ -32,9 +32,12 @@
hdur,hdur_gaussian,t_cmt,dt,t0,sourcearrays, &
ispec_is_elastic,SIMULATION_TYPE,NSTEP,NGLOB_ADJOINT, &
nrec,islice_selected_rec,ispec_selected_rec, &
- nadj_rec_local,adj_sourcearrays,b_accel )
+ nadj_rec_local,adj_sourcearrays,b_accel, &
+ NTSTEP_BETWEEN_READ_ADJSRC )
- use specfem_par,only: PRINT_SOURCE_TIME_FUNCTION,stf_used_total
+ use specfem_par,only: PRINT_SOURCE_TIME_FUNCTION,stf_used_total, &
+ xigll,yigll,zigll,xi_receiver,eta_receiver,gamma_receiver,&
+ station_name,network_name,adj_source_file
implicit none
include "constants.h"
@@ -69,8 +72,14 @@
integer:: nrec
integer,dimension(nrec) :: islice_selected_rec,ispec_selected_rec
integer:: nadj_rec_local
- real(kind=CUSTOM_REAL),dimension(nadj_rec_local,NSTEP,NDIM,NGLLX,NGLLY,NGLLZ):: adj_sourcearrays
real(kind=CUSTOM_REAL),dimension(NDIM,NGLOB_ADJOINT):: b_accel
+!<YANGL
+! real(kind=CUSTOM_REAL),dimension(nadj_rec_local,NSTEP,NDIM,NGLLX,NGLLY,NGLLZ):: adj_sourcearrays
+ real(kind=CUSTOM_REAL),dimension(nadj_rec_local,NTSTEP_BETWEEN_READ_ADJSRC,NDIM,NGLLX,NGLLY,NGLLZ):: adj_sourcearrays
+ integer :: it_sub_adj,itime,NTSTEP_BETWEEN_READ_ADJSRC
+ real(kind=CUSTOM_REAL),dimension(NTSTEP_BETWEEN_READ_ADJSRC,NDIM,NGLLX,NGLLY,NGLLZ):: adj_sourcearray
+ logical :: ibool_read_adj_arrays
+!>YANGL
! local parameters
double precision :: f0
@@ -185,6 +194,40 @@
! adjoint simulations
if (SIMULATION_TYPE == 2 .or. SIMULATION_TYPE == 3) then
+!<YANGL
+! read in adjoint sources block by block (for memory consideration)
+! e.g., in exploration experiments, both the number of receivers (nrec) and the number of time steps (NSTEP) are huge,
+! which may cause problems since we have a large array: adj_sourcearrays(nadj_rec_local,NSTEP,NDIM,NGLLX,NGLLY,NGLLZ)
+
+ ! figure out if we need to read in a chunk of the adjoint source at this timestep
+ it_sub_adj = ceiling( dble(it)/dble(NTSTEP_BETWEEN_READ_ADJSRC) ) !chunk_number
+ ibool_read_adj_arrays = (((mod(it-1,NTSTEP_BETWEEN_READ_ADJSRC) == 0)) .and. (nadj_rec_local > 0))
+
+ ! needs to read in a new chunk/block of the adjoint source
+ if(ibool_read_adj_arrays) then
+
+ irec_local = 0
+ do irec = 1, nrec
+ ! compute source arrays
+ if(myrank == islice_selected_rec(irec)) then
+ irec_local = irec_local + 1
+
+ ! reads in **sta**.**net**.**LH**.adj files
+ adj_source_file = trim(station_name(irec))//'.'//trim(network_name(irec))
+ call compute_arrays_adjoint_source(myrank,adj_source_file, &
+ xi_receiver(irec),eta_receiver(irec),gamma_receiver(irec), &
+ adj_sourcearray, xigll,yigll,zigll, &
+ it_sub_adj,NSTEP,NTSTEP_BETWEEN_READ_ADJSRC)
+ do itime = 1,NTSTEP_BETWEEN_READ_ADJSRC
+ adj_sourcearrays(irec_local,itime,:,:,:,:) = adj_sourcearray(itime,:,:,:,:)
+ enddo
+
+ endif
+ enddo
+
+ endif ! if(ibool_read_adj_arrays)
+!>YANGL
+
if( it < NSTEP ) then
! receivers act as sources
@@ -198,7 +241,10 @@
do j=1,NGLLY
do i=1,NGLLX
iglob = ibool(i,j,k,ispec_selected_rec(irec))
- accel(:,iglob) = accel(:,iglob) + adj_sourcearrays(irec_local,NSTEP-it+1,:,i,j,k)
+!<YANGL
+! accel(:,iglob) = accel(:,iglob) + adj_sourcearrays(irec_local,NSTEP-it+1,:,i,j,k)
+ accel(:,iglob) = accel(:,iglob) + adj_sourcearrays(irec_local,NTSTEP_BETWEEN_READ_ADJSRC-mod(it-1,NTSTEP_BETWEEN_READ_ADJSRC),:,i,j,k)
+!>YANGL
enddo
enddo
enddo
Modified: seismo/3D/SPECFEM3D/trunk/compute_arrays_source.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/compute_arrays_source.f90 2010-10-04 09:49:07 UTC (rev 17235)
+++ seismo/3D/SPECFEM3D/trunk/compute_arrays_source.f90 2010-10-04 20:20:38 UTC (rev 17236)
@@ -164,25 +164,109 @@
end subroutine multiply_arrays_source
!=============================================================================
+!<YANGL
+! testing read in adjoint sources block by block
+!!!the original version
+!!!
+!!!subroutine compute_arrays_adjoint_source(myrank, adj_source_file, &
+!!! xi_receiver,eta_receiver,gamma_receiver, adj_sourcearray, &
+!!! xigll,yigll,zigll,NSTEP)
+!!!
+!!!
+!!! implicit none
+!!!
+!!! include 'constants.h'
+!!!
+!!!! input
+!!! integer myrank, NSTEP
+!!!
+!!! double precision xi_receiver, eta_receiver, gamma_receiver
+!!!
+!!! character(len=*) adj_source_file
+!!!
+!!!! output
+!!! real(kind=CUSTOM_REAL),dimension(NSTEP,NDIM,NGLLX,NGLLY,NGLLZ) :: adj_sourcearray
+!!!
+!!!! Gauss-Lobatto-Legendre points of integration and weights
+!!! double precision, dimension(NGLLX) :: xigll
+!!! double precision, dimension(NGLLY) :: yigll
+!!! double precision, dimension(NGLLZ) :: zigll
+!!!
+!!! double precision :: hxir(NGLLX), hpxir(NGLLX), hetar(NGLLY), hpetar(NGLLY), &
+!!! hgammar(NGLLZ), hpgammar(NGLLZ)
+!!!
+!!! real(kind=CUSTOM_REAL) :: adj_src(NSTEP,NDIM)
+!!!
+!!! integer icomp, itime, i, j, k, ios
+!!! double precision :: junk
+!!! ! note: should have same order as orientation in write_seismograms_to_file()
+!!! character(len=3),dimension(NDIM) :: comp = (/ "BHE", "BHN", "BHZ" /)
+!!! character(len=256) :: filename
+!!!
+!!! !adj_sourcearray(:,:,:,:,:) = 0.
+!!! adj_src = 0._CUSTOM_REAL
+!!!
+!!! ! loops over components
+!!! do icomp = 1, NDIM
+!!!
+!!! filename = 'SEM/'//trim(adj_source_file) // '.'// comp(icomp) // '.adj'
+!!! open(unit=IIN,file=trim(filename),status='old',action='read',iostat = ios)
+!!! if (ios /= 0) cycle ! cycles to next file
+!!! !if (ios /= 0) call exit_MPI(myrank, ' file '//trim(filename)//'does not exist')
+!!!
+!!! ! reads in adjoint source trace
+!!! do itime = 1, NSTEP
+!!!
+!!! ! things become a bit tricky because of the Newark time scheme at
+!!! ! the very beginning of the time loop. however, when we read in the backward/reconstructed
+!!! ! wavefields at the end of the first time loop, we can use the adjoint source index from 1 to NSTEP
+!!! ! (and then access it in reverse NSTEP-it+1 down to 1, for it=1,..NSTEP; see compute_add_sources*.f90).
+!!! read(IIN,*,iostat=ios) junk, adj_src(itime,icomp)
+!!! if( ios /= 0 ) &
+!!! call exit_MPI(myrank, &
+!!! 'file '//trim(filename)//' has wrong length, please check with your simulation duration')
+!!! enddo
+!!! close(IIN)
+!!!
+!!! enddo
+!!!
+!!! ! lagrange interpolators for receiver location
+!!! call lagrange_any(xi_receiver,NGLLX,xigll,hxir,hpxir)
+!!! call lagrange_any(eta_receiver,NGLLY,yigll,hetar,hpetar)
+!!! call lagrange_any(gamma_receiver,NGLLZ,zigll,hgammar,hpgammar)
+!!!
+!!! ! interpolates adjoint source onto GLL points within this element
+!!! do k = 1, NGLLZ
+!!! do j = 1, NGLLY
+!!! do i = 1, NGLLX
+!!! adj_sourcearray(:,:,i,j,k) = hxir(i) * hetar(j) * hgammar(k) * adj_src(:,:)
+!!! enddo
+!!! enddo
+!!! enddo
+!!!
+!!!end subroutine compute_arrays_adjoint_source
+
+
+!!! test version
subroutine compute_arrays_adjoint_source(myrank, adj_source_file, &
xi_receiver,eta_receiver,gamma_receiver, adj_sourcearray, &
- xigll,yigll,zigll,NSTEP)
+ xigll,yigll,zigll, &
+ it_sub_adj,NSTEP,NTSTEP_BETWEEN_READ_ADJSRC)
-
implicit none
include 'constants.h'
! input
- integer myrank, NSTEP
+ integer myrank, NSTEP, it_sub_adj, NTSTEP_BETWEEN_READ_ADJSRC
double precision xi_receiver, eta_receiver, gamma_receiver
character(len=*) adj_source_file
! output
- real(kind=CUSTOM_REAL),dimension(NSTEP,NDIM,NGLLX,NGLLY,NGLLZ) :: adj_sourcearray
+ real(kind=CUSTOM_REAL),dimension(NTSTEP_BETWEEN_READ_ADJSRC,NDIM,NGLLX,NGLLY,NGLLZ) :: adj_sourcearray
! Gauss-Lobatto-Legendre points of integration and weights
double precision, dimension(NGLLX) :: xigll
@@ -192,14 +276,18 @@
double precision :: hxir(NGLLX), hpxir(NGLLX), hetar(NGLLY), hpetar(NGLLY), &
hgammar(NGLLZ), hpgammar(NGLLZ)
- real(kind=CUSTOM_REAL) :: adj_src(NSTEP,NDIM)
+ real(kind=CUSTOM_REAL) :: adj_src(NTSTEP_BETWEEN_READ_ADJSRC,NDIM)
- integer icomp, itime, i, j, k, ios
+ integer icomp, itime, i, j, k, ios, it_start, it_end
double precision :: junk
! note: should have same order as orientation in write_seismograms_to_file()
character(len=3),dimension(NDIM) :: comp = (/ "BHE", "BHN", "BHZ" /)
character(len=256) :: filename
+ ! range of the block we need to read
+ it_start = NSTEP - it_sub_adj*NTSTEP_BETWEEN_READ_ADJSRC + 1
+ it_end = it_start + NTSTEP_BETWEEN_READ_ADJSRC - 1
+
!adj_sourcearray(:,:,:,:,:) = 0.
adj_src = 0._CUSTOM_REAL
@@ -212,17 +300,20 @@
!if (ios /= 0) call exit_MPI(myrank, ' file '//trim(filename)//'does not exist')
! reads in adjoint source trace
- do itime = 1, NSTEP
-
- ! things become a bit tricky because of the Newark time scheme at
- ! the very beginning of the time loop. however, when we read in the backward/reconstructed
- ! wavefields at the end of the first time loop, we can use the adjoint source index from 1 to NSTEP
- ! (and then access it in reverse NSTEP-it+1 down to 1, for it=1,..NSTEP; see compute_add_sources*.f90).
- read(IIN,*,iostat=ios) junk, adj_src(itime,icomp)
+ !! skip unused blocks
+ do itime = 1, it_start-1
+ read(IIN,*,iostat=ios) junk, junk
if( ios /= 0 ) &
call exit_MPI(myrank, &
- 'file '//trim(filename)//' has wrong length, please check with your simulation duration')
+ 'file '//trim(filename)//' has wrong length, please check with your simulation duration (1111)')
enddo
+ !! read the block we need
+ do itime = it_start, it_end
+ read(IIN,*,iostat=ios) junk, adj_src(itime-it_start+1,icomp)
+ if( ios /= 0 ) &
+ call exit_MPI(myrank, &
+ 'file '//trim(filename)//' has wrong length, please check with your simulation duration (2222)')
+ enddo
close(IIN)
enddo
@@ -243,7 +334,9 @@
end subroutine compute_arrays_adjoint_source
+!>YANGL
+
! =======================================================================
! compute the integrated derivatives of source parameters (M_jk and X_s)
Modified: seismo/3D/SPECFEM3D/trunk/compute_forces_acoustic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/compute_forces_acoustic.f90 2010-10-04 09:49:07 UTC (rev 17235)
+++ seismo/3D/SPECFEM3D/trunk/compute_forces_acoustic.f90 2010-10-04 20:20:38 UTC (rev 17236)
@@ -204,7 +204,8 @@
sourcearrays,kappastore,ispec_is_acoustic,&
SIMULATION_TYPE,NSTEP,NGLOB_ADJOINT, &
nrec,islice_selected_rec,ispec_selected_rec, &
- nadj_rec_local,adj_sourcearrays,b_potential_dot_dot_acoustic )
+ nadj_rec_local,adj_sourcearrays,b_potential_dot_dot_acoustic, &
+ NTSTEP_BETWEEN_READ_ADJSRC )
! assemble all the contributions between slices using MPI
if( phase_is_inner .eqv. .false. ) then
Modified: seismo/3D/SPECFEM3D/trunk/compute_forces_elastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/compute_forces_elastic.f90 2010-10-04 09:49:07 UTC (rev 17235)
+++ seismo/3D/SPECFEM3D/trunk/compute_forces_elastic.f90 2010-10-04 20:20:38 UTC (rev 17236)
@@ -213,7 +213,8 @@
hdur,hdur_gaussian,t_cmt,dt,t0,sourcearrays, &
ispec_is_elastic,SIMULATION_TYPE,NSTEP,NGLOB_ADJOINT, &
nrec,islice_selected_rec,ispec_selected_rec, &
- nadj_rec_local,adj_sourcearrays,b_accel )
+ nadj_rec_local,adj_sourcearrays,b_accel, &
+ NTSTEP_BETWEEN_READ_ADJSRC )
! assemble all the contributions between slices using MPI
if( phase_is_inner .eqv. .false. ) then
Modified: seismo/3D/SPECFEM3D/trunk/create_header_file.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/create_header_file.f90 2010-10-04 09:49:07 UTC (rev 17235)
+++ seismo/3D/SPECFEM3D/trunk/create_header_file.f90 2010-10-04 20:20:38 UTC (rev 17236)
@@ -33,7 +33,7 @@
include "constants.h"
integer NTSTEP_BETWEEN_OUTPUT_SEISMOS,NSTEP,UTM_PROJECTION_ZONE,SIMULATION_TYPE
- integer NSOURCES
+ integer NSOURCES,NTSTEP_BETWEEN_READ_ADJSRC
! parameters to be computed based upon parameters above read from file
integer NPROC
@@ -67,7 +67,8 @@
MOVIE_SURFACE,MOVIE_VOLUME,CREATE_SHAKEMAP,SAVE_DISPLACEMENT, &
NTSTEP_BETWEEN_FRAMES,USE_HIGHRES_FOR_MOVIES,HDUR_MOVIE, &
SAVE_AVS_DX_MESH_FILES,PRINT_SOURCE_TIME_FUNCTION, &
- NTSTEP_BETWEEN_OUTPUT_INFO,SIMULATION_TYPE,SAVE_FORWARD)
+ NTSTEP_BETWEEN_OUTPUT_INFO,SIMULATION_TYPE,SAVE_FORWARD, &
+ NTSTEP_BETWEEN_READ_ADJSRC)
! create include file for the solver
call save_header_file(NSPEC_AB,NGLOB_AB,NPROC, &
Modified: seismo/3D/SPECFEM3D/trunk/create_movie_shakemap_AVS_DX_GMT.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/create_movie_shakemap_AVS_DX_GMT.f90 2010-10-04 09:49:07 UTC (rev 17235)
+++ seismo/3D/SPECFEM3D/trunk/create_movie_shakemap_AVS_DX_GMT.f90 2010-10-04 20:20:38 UTC (rev 17236)
@@ -83,7 +83,7 @@
! parameters read from parameter file
integer NTSTEP_BETWEEN_OUTPUT_SEISMOS,NSTEP,UTM_PROJECTION_ZONE,SIMULATION_TYPE
- integer NSOURCES
+ integer NSOURCES,NTSTEP_BETWEEN_READ_ADJSRC
logical MOVIE_SURFACE,MOVIE_VOLUME,CREATE_SHAKEMAP,SAVE_DISPLACEMENT, &
USE_HIGHRES_FOR_MOVIES,SUPPRESS_UTM_PROJECTION
integer NTSTEP_BETWEEN_FRAMES,NTSTEP_BETWEEN_OUTPUT_INFO
@@ -132,7 +132,8 @@
MOVIE_SURFACE,MOVIE_VOLUME,CREATE_SHAKEMAP,SAVE_DISPLACEMENT, &
NTSTEP_BETWEEN_FRAMES,USE_HIGHRES_FOR_MOVIES,HDUR_MOVIE, &
SAVE_MESH_FILES,PRINT_SOURCE_TIME_FUNCTION, &
- NTSTEP_BETWEEN_OUTPUT_INFO,SIMULATION_TYPE,SAVE_FORWARD)
+ NTSTEP_BETWEEN_OUTPUT_INFO,SIMULATION_TYPE,SAVE_FORWARD, &
+ NTSTEP_BETWEEN_READ_ADJSRC)
! get the base pathname for output files
call get_value_string(OUTPUT_FILES, 'OUTPUT_FILES', 'OUTPUT_FILES')
Modified: seismo/3D/SPECFEM3D/trunk/generate_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/generate_databases.f90 2010-10-04 09:49:07 UTC (rev 17235)
+++ seismo/3D/SPECFEM3D/trunk/generate_databases.f90 2010-10-04 20:20:38 UTC (rev 17236)
@@ -226,7 +226,7 @@
logical :: MOVIE_SURFACE,MOVIE_VOLUME,CREATE_SHAKEMAP,SAVE_DISPLACEMENT, &
USE_HIGHRES_FOR_MOVIES
- integer :: NTSTEP_BETWEEN_FRAMES,NTSTEP_BETWEEN_OUTPUT_INFO
+ integer :: NTSTEP_BETWEEN_FRAMES,NTSTEP_BETWEEN_OUTPUT_INFO,NTSTEP_BETWEEN_READ_ADJSRC
character(len=256) OUTPUT_FILES,LOCAL_PATH
@@ -369,7 +369,8 @@
MOVIE_SURFACE,MOVIE_VOLUME,CREATE_SHAKEMAP,SAVE_DISPLACEMENT, &
NTSTEP_BETWEEN_FRAMES,USE_HIGHRES_FOR_MOVIES,HDUR_MOVIE, &
SAVE_MESH_FILES,PRINT_SOURCE_TIME_FUNCTION, &
- NTSTEP_BETWEEN_OUTPUT_INFO,SIMULATION_TYPE,SAVE_FORWARD)
+ NTSTEP_BETWEEN_OUTPUT_INFO,SIMULATION_TYPE,SAVE_FORWARD, &
+ NTSTEP_BETWEEN_READ_ADJSRC)
! check that the code is running with the requested nb of processes
if(sizeprocs /= NPROC) then
Modified: seismo/3D/SPECFEM3D/trunk/initialize_simulation.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/initialize_simulation.f90 2010-10-04 09:49:07 UTC (rev 17235)
+++ seismo/3D/SPECFEM3D/trunk/initialize_simulation.f90 2010-10-04 20:20:38 UTC (rev 17236)
@@ -44,7 +44,8 @@
MOVIE_SURFACE,MOVIE_VOLUME,CREATE_SHAKEMAP,SAVE_DISPLACEMENT, &
NTSTEP_BETWEEN_FRAMES,USE_HIGHRES_FOR_MOVIES,HDUR_MOVIE, &
SAVE_MESH_FILES,PRINT_SOURCE_TIME_FUNCTION, &
- NTSTEP_BETWEEN_OUTPUT_INFO,SIMULATION_TYPE,SAVE_FORWARD)
+ NTSTEP_BETWEEN_OUTPUT_INFO,SIMULATION_TYPE,SAVE_FORWARD, &
+ NTSTEP_BETWEEN_READ_ADJSRC)
! get the base pathname for output files
call get_value_string(OUTPUT_FILES, 'OUTPUT_FILES', 'OUTPUT_FILES')
Modified: seismo/3D/SPECFEM3D/trunk/read_parameter_file.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/read_parameter_file.f90 2010-10-04 09:49:07 UTC (rev 17235)
+++ seismo/3D/SPECFEM3D/trunk/read_parameter_file.f90 2010-10-04 20:20:38 UTC (rev 17236)
@@ -30,13 +30,14 @@
MOVIE_SURFACE,MOVIE_VOLUME,CREATE_SHAKEMAP,SAVE_DISPLACEMENT, &
NTSTEP_BETWEEN_FRAMES,USE_HIGHRES_FOR_MOVIES,HDUR_MOVIE, &
SAVE_MESH_FILES,PRINT_SOURCE_TIME_FUNCTION,NTSTEP_BETWEEN_OUTPUT_INFO, &
- SIMULATION_TYPE,SAVE_FORWARD )
+ SIMULATION_TYPE,SAVE_FORWARD, &
+ NTSTEP_BETWEEN_READ_ADJSRC )
implicit none
include "constants.h"
- integer NPROC,NTSTEP_BETWEEN_OUTPUT_SEISMOS,NSTEP,SIMULATION_TYPE
+ integer NPROC,NTSTEP_BETWEEN_OUTPUT_SEISMOS,NSTEP,SIMULATION_TYPE, NTSTEP_BETWEEN_READ_ADJSRC
integer NSOURCES,NTSTEP_BETWEEN_FRAMES,NTSTEP_BETWEEN_OUTPUT_INFO,UTM_PROJECTION_ZONE
double precision DT,HDUR_MOVIE
@@ -86,6 +87,13 @@
endif
call read_value_integer(NSTEP, 'solver.NSTEP')
if(err_occurred() /= 0) return
+!<YANGL
+! read in adjoint sources block by block
+! we shall later put this parameter in DATA/Par_file
+! the default value is to read in the whole trace
+! should you change the value, make sure "mod(NTSTEP_BETWEEN_READ_ADJSRC,NSTEP) == 0"
+ NTSTEP_BETWEEN_READ_ADJSRC = NSTEP
+!>YANGL
call read_value_double_precision(DT, 'solver.DT')
if(err_occurred() /= 0) return
call read_value_logical(OCEANS, 'model.OCEANS')
Modified: seismo/3D/SPECFEM3D/trunk/setup_sources_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/setup_sources_receivers.f90 2010-10-04 09:49:07 UTC (rev 17235)
+++ seismo/3D/SPECFEM3D/trunk/setup_sources_receivers.f90 2010-10-04 20:20:38 UTC (rev 17236)
@@ -550,8 +550,12 @@
endif
! reads in adjoint source traces
- allocate(adj_sourcearray(NSTEP,NDIM,NGLLX,NGLLY,NGLLZ))
- allocate(adj_sourcearrays(nadj_rec_local,NSTEP,NDIM,NGLLX,NGLLY,NGLLZ))
+!<YANGL
+! allocate(adj_sourcearray(NSTEP,NDIM,NGLLX,NGLLY,NGLLZ))
+! allocate(adj_sourcearrays(nadj_rec_local,NSTEP,NDIM,NGLLX,NGLLY,NGLLZ))
+ allocate(adj_sourcearray(NTSTEP_BETWEEN_READ_ADJSRC,NDIM,NGLLX,NGLLY,NGLLZ))
+ allocate(adj_sourcearrays(nadj_rec_local,NTSTEP_BETWEEN_READ_ADJSRC,NDIM,NGLLX,NGLLY,NGLLZ))
+!>YANGL
adj_sourcearrays = 0._CUSTOM_REAL
adj_sourcearray = 0._CUSTOM_REAL
@@ -567,7 +571,7 @@
call compute_arrays_adjoint_source(myrank, adj_source_file, &
xi_receiver(irec), eta_receiver(irec), gamma_receiver(irec), &
- adj_sourcearray, xigll,yigll,zigll,NSTEP)
+ adj_sourcearray, xigll,yigll,zigll,1,NSTEP,NTSTEP_BETWEEN_READ_ADJSRC)
adj_sourcearrays(irec_local,:,:,:,:,:) = adj_sourcearray(:,:,:,:,:)
Modified: seismo/3D/SPECFEM3D/trunk/specfem3D_par.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/specfem3D_par.f90 2010-10-04 09:49:07 UTC (rev 17235)
+++ seismo/3D/SPECFEM3D/trunk/specfem3D_par.f90 2010-10-04 20:20:38 UTC (rev 17236)
@@ -234,6 +234,8 @@
! norm of the backward displacement
real(kind=CUSTOM_REAL) b_Usolidnorm, b_Usolidnorm_all
+ ! length of reading blocks
+ integer :: NTSTEP_BETWEEN_READ_ADJSRC
end module specfem_par
More information about the CIG-COMMITS
mailing list