[cig-commits] r17162 - seismo/3D/SPECFEM3D/trunk
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Wed Sep 1 19:41:23 PDT 2010
Author: danielpeter
Date: 2010-09-01 19:41:23 -0700 (Wed, 01 Sep 2010)
New Revision: 17162
Modified:
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_stacey_acoustic.f90
seismo/3D/SPECFEM3D/trunk/compute_stacey_elastic.f90
seismo/3D/SPECFEM3D/trunk/constants.h.in
seismo/3D/SPECFEM3D/trunk/iterate_time.f90
seismo/3D/SPECFEM3D/trunk/locate_source.f90
seismo/3D/SPECFEM3D/trunk/prepare_timerun.f90
seismo/3D/SPECFEM3D/trunk/setup_sources_receivers.f90
seismo/3D/SPECFEM3D/trunk/todo_list_please_dont_remove.txt
seismo/3D/SPECFEM3D/trunk/write_seismograms.f90
Log:
updates reconstruction of forward wavefield in adjoint simulations; minor bug fix in sr.vtk creation for point force sources
Modified: seismo/3D/SPECFEM3D/trunk/compute_add_sources_acoustic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/compute_add_sources_acoustic.f90 2010-09-02 00:38:38 UTC (rev 17161)
+++ seismo/3D/SPECFEM3D/trunk/compute_add_sources_acoustic.f90 2010-09-02 02:41:23 UTC (rev 17162)
@@ -182,26 +182,28 @@
endif
! NOTE: adjoint sources and backward wavefield timing:
-! idea is to start with the backward field b_potential.. at time (T)
+! idea is to start with the backward field b_potential,.. at time (T)
! and convolve with the adjoint field at time (T-t)
!
! backward/reconstructed wavefields:
-! time for b_potential..( it ) corresponds to (NSTEP - it - 1 )*DT - t0 ...
-! since we start with saved wavefields b_potential..( 0 ) = potential..( NSTEP ) which correspond
-! to a time (NSTEP - 1)*DT - t0
+! time for b_potential( it ) would correspond to (NSTEP - it - 1 )*DT - t0
+! if we read in saved wavefields b_potential() before Newark time scheme
! (see sources for simulation_type 1 and seismograms)
-! now, at the beginning of the time loop, the numerical Newark time scheme updates
-! the wavefields, that is b_potential..( it=1) corresponds now to time (NSTEP -1 - 1)*DT - t0
+! since at the beginning of the time loop, the numerical Newark time scheme updates
+! the wavefields, that is b_potential( it=1) would correspond to time (NSTEP -1 - 1)*DT - t0
!
-! let's define the start time t to (1-1)*DT - t0 = -t0, and the end time T to (NSTEP-1)*DT - t0
-! these are the start and end times of all seismograms
-!
+! b_potential is now read in after Newark time scheme:
+! we read the backward/reconstructed wavefield at the end of the first time loop,
+! such that b_potential(it=1) corresponds to -t0 + (NSTEP-1)*DT.
+! assuming that until that end the backward/reconstructed wavefield and adjoint fields
+! have a zero contribution to adjoint kernels.
+! thus the correct indexing is NSTEP - it + 1, instead of NSTEP - it
+!
! adjoint wavefields:
! since the adjoint source traces were derived from the seismograms,
! it follows that for the adjoint wavefield, the time equivalent to ( T - t ) uses the time-reversed
! adjoint source traces which start at -t0 and end at time (NSTEP-1)*DT - t0
-! for it=1: (NSTEP -1 - 1)*DT - t0 for backward wavefields corresponds to time T-1
-! and time (T-1) corresponds now to index (NSTEP -1) in the adjoint source array
+! for step it=1: (NSTEP -it + 1)*DT - t0 for backward wavefields corresponds to time T
! adjoint simulations
if (SIMULATION_TYPE == 2 .or. SIMULATION_TYPE == 3) then
@@ -220,9 +222,11 @@
do i=1,NGLLX
iglob = ibool(i,j,k,ispec)
- ! beware, for acoustic medium, source is: pressure divided by Kappa of the fluid
+ ! 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
potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) &
- - adj_sourcearrays(irec_local,NSTEP-it,1,i,j,k) / kappastore(i,j,k,ispec)
+ - adj_sourcearrays(irec_local,NSTEP-it+1,1,i,j,k) / kappastore(i,j,k,ispec)
enddo
enddo
enddo
@@ -231,6 +235,10 @@
endif ! it
endif
+! note: b_potential() is read in after Newark time scheme, thus
+! b_potential(it=1) corresponds to -t0 + (NSTEP-1)*DT.
+! thus indexing is NSTEP - it , instead of NSTEP - it - 1
+
! adjoint simulations
if (SIMULATION_TYPE == 3) then
! adds acoustic sources
@@ -267,7 +275,7 @@
! we use nu_source(:,3) here because we want a source normal to the surface.
! This is the expression of a Ricker; should be changed according maybe to the Par_file.
- stf_used = FACTOR_FORCE_SOURCE * comp_source_time_function_rickr(dble(NSTEP-it-1)*DT-t0-t_cmt(isource),f0)
+ stf_used = FACTOR_FORCE_SOURCE * comp_source_time_function_rickr(dble(NSTEP-it)*DT-t0-t_cmt(isource),f0)
! beware, for acoustic medium, source is: pressure divided by Kappa of the fluid
! the sign is negative because pressure p = - Chi_dot_dot therefore we need
@@ -283,7 +291,7 @@
else
! gaussian source time
- stf = comp_source_time_function_gauss(dble(NSTEP-it-1)*DT-t0-t_cmt(isource),hdur_gaussian(isource))
+ stf = comp_source_time_function_gauss(dble(NSTEP-it)*DT-t0-t_cmt(isource),hdur_gaussian(isource))
! distinguishes between single and double precision for reals
if(CUSTOM_REAL == SIZE_REAL) then
Modified: seismo/3D/SPECFEM3D/trunk/compute_add_sources_elastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/compute_add_sources_elastic.f90 2010-09-02 00:38:38 UTC (rev 17161)
+++ seismo/3D/SPECFEM3D/trunk/compute_add_sources_elastic.f90 2010-09-02 02:41:23 UTC (rev 17162)
@@ -108,6 +108,7 @@
ispec_selected_source(isource))
f0 = hdur(isource) !! using hdur as a FREQUENCY just to avoid changing CMTSOLUTION file format
+ ! note: for a Ricker source time function, a start time ~1.2 * main_period is a good choice
t0 = 1.2d0/f0
!if (it == 1 .and. myrank == 0) then
@@ -117,11 +118,13 @@
!endif
! This is the expression of a Ricker; should be changed according maybe to the Par_file.
- stf_used = FACTOR_FORCE_SOURCE * comp_source_time_function_rickr(dble(it-1)*DT-t0-t_cmt(isource),f0)
+ stf_used = FACTOR_FORCE_SOURCE * comp_source_time_function_rickr(dble(it-1)*DT-t0-t_cmt(isource),f0)
- ! we use nu_source(:,3) here because we want a source normal to the surface (z-direction).
+ ! we use a force in a single direction along one of the components:
+ ! x/y/z or E/N/Z-direction would correspond to 1/2/3 = COMPONENT_FORCE_SOURCE
+ ! e.g. nu_source(:,3) here would be a source normal to the surface (z-direction).
accel(:,iglob) = accel(:,iglob) &
- + sngl( nu_source(:,3,isource) ) * stf_used
+ + sngl( nu_source(:,COMPONENT_FORCE_SOURCE,isource) ) * stf_used
else
@@ -160,23 +163,24 @@
! and convolve with the adjoint field at time (T-t)
!
! backward/reconstructed wavefields:
-! time for b_displ( it ) corresponds to (NSTEP - it - 1 )*DT - t0 ...
-! since we start with saved wavefields b_displ( 0 ) = displ( NSTEP ) which correspond
-! to a time (NSTEP - 1)*DT - t0
+! time for b_displ( it ) would correspond to (NSTEP - it - 1 )*DT - t0
+! if we read in saved wavefields b_displ() before Newark time scheme
! (see sources for simulation_type 1 and seismograms)
-! now, at the beginning of the time loop, the numerical Newark time scheme updates
-! the wavefields, that is b_displ( it=1) corresponds now to time (NSTEP -1 - 1)*DT - t0
+! since at the beginning of the time loop, the numerical Newark time scheme updates
+! the wavefields, that is b_displ( it=1) would correspond to time (NSTEP -1 - 1)*DT - t0
!
-! let's define the start time t to (1-1)*DT - t0 = -t0, and the end time T to (NSTEP-1)*DT - t0
-! these are the start and end times of all seismograms
-!
+! b_displ is now read in after Newark time scheme:
+! we read the backward/reconstructed wavefield at the end of the first time loop,
+! such that b_displ(it=1) corresponds to -t0 + (NSTEP-1)*DT.
+! assuming that until that end the backward/reconstructed wavefield and adjoint fields
+! have a zero contribution to adjoint kernels.
+! thus the correct indexing is NSTEP - it + 1, instead of NSTEP - it
+!
! adjoint wavefields:
! since the adjoint source traces were derived from the seismograms,
! it follows that for the adjoint wavefield, the time equivalent to ( T - t ) uses the time-reversed
! adjoint source traces which start at -t0 and end at time (NSTEP-1)*DT - t0
-! for it=1: (NSTEP -1 - 1)*DT - t0 for backward wavefields corresponds to time T-1
-! and time (T-1) corresponds now to index (NSTEP -1) in the adjoint source array
-
+! for step it=1: (NSTEP -it + 1)*DT - t0 for backward wavefields corresponds to time T
! adjoint simulations
if (SIMULATION_TYPE == 2 .or. SIMULATION_TYPE == 3) then
@@ -194,7 +198,7 @@
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,:,i,j,k)
+ accel(:,iglob) = accel(:,iglob) + adj_sourcearrays(irec_local,NSTEP-it+1,:,i,j,k)
enddo
enddo
enddo
@@ -205,6 +209,10 @@
endif !adjoint
+! note: b_displ() is read in after Newark time scheme, thus
+! b_displ(it=1) corresponds to -t0 + (NSTEP-1)*DT.
+! thus indexing is NSTEP - it , instead of NSTEP - it - 1
+
! adjoint simulations
if (SIMULATION_TYPE == 3) then
@@ -238,19 +246,18 @@
!endif
! This is the expression of a Ricker; should be changed according maybe to the Par_file.
- stf_used = FACTOR_FORCE_SOURCE * comp_source_time_function_rickr(dble(NSTEP-it-1)*DT-t0-t_cmt(isource),f0)
+ stf_used = FACTOR_FORCE_SOURCE * comp_source_time_function_rickr(dble(NSTEP-it)*DT-t0-t_cmt(isource),f0)
- ! we use nu_source(:,3) here because we want a source normal to the surface.
+ ! e.g. we use nu_source(:,3) here if we want a source normal to the surface.
! note: time step is now at NSTEP-it
b_accel(:,iglob) = b_accel(:,iglob) &
- + sngl( nu_source(:,3,isource) ) * stf_used
-
+ + sngl( nu_source(:,COMPONENT_FORCE_SOURCE,isource) ) * stf_used
else
- ! see note above: time step corresponds now to NSTEP-it-1
+ ! see note above: time step corresponds now to NSTEP-it
! (also compare to it-1 for forward simulation)
- stf = comp_source_time_function(dble(NSTEP-it-1)*DT-t0-t_cmt(isource),hdur_gaussian(isource))
+ stf = comp_source_time_function(dble(NSTEP-it)*DT-t0-t_cmt(isource),hdur_gaussian(isource))
! distinguish between single and double precision for reals
if(CUSTOM_REAL == SIZE_REAL) then
Modified: seismo/3D/SPECFEM3D/trunk/compute_arrays_source.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/compute_arrays_source.f90 2010-09-02 00:38:38 UTC (rev 17161)
+++ seismo/3D/SPECFEM3D/trunk/compute_arrays_source.f90 2010-09-02 02:41:23 UTC (rev 17162)
@@ -196,7 +196,7 @@
integer icomp, itime, i, j, k, ios
double precision :: junk
- ! note: should have some order as orientation in write_seismograms_to_file()
+ ! 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
@@ -213,7 +213,11 @@
! 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, &
Modified: seismo/3D/SPECFEM3D/trunk/compute_stacey_acoustic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/compute_stacey_acoustic.f90 2010-09-02 00:38:38 UTC (rev 17161)
+++ seismo/3D/SPECFEM3D/trunk/compute_stacey_acoustic.f90 2010-09-02 02:41:23 UTC (rev 17162)
@@ -74,6 +74,7 @@
! adjoint simulations:
if (SIMULATION_TYPE == 3 .and. num_abs_boundary_faces > 0) then
+ ! the index NSTEP-it+1 is valid if b_displ is read in after the Newark scheme
read(IOABS_AC,rec=NSTEP-it+1) reclen1,b_absorb_potential,reclen2
if (reclen1 /= b_reclen_potential .or. reclen1 /= reclen2) &
call exit_mpi(myrank,'Error reading absorbing contribution b_absorb_potential')
Modified: seismo/3D/SPECFEM3D/trunk/compute_stacey_elastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/compute_stacey_elastic.f90 2010-09-02 00:38:38 UTC (rev 17161)
+++ seismo/3D/SPECFEM3D/trunk/compute_stacey_elastic.f90 2010-09-02 02:41:23 UTC (rev 17162)
@@ -83,6 +83,7 @@
! adjoint simulations:
if (SIMULATION_TYPE == 3 .and. num_abs_boundary_faces > 0) then
+ ! the index NSTEP-it+1 is valid if b_displ is read in after the Newark scheme
read(IOABS,rec=NSTEP-it+1) reclen1,b_absorb_field,reclen2
if (reclen1 /= b_reclen_field .or. reclen1 /= reclen2) &
call exit_mpi(myrank,'Error reading absorbing contribution b_absorb_field')
Modified: seismo/3D/SPECFEM3D/trunk/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D/trunk/constants.h.in 2010-09-02 00:38:38 UTC (rev 17161)
+++ seismo/3D/SPECFEM3D/trunk/constants.h.in 2010-09-02 02:41:23 UTC (rev 17162)
@@ -141,7 +141,8 @@
! in which the source is a vertical force, normal force, impact etc.
logical, parameter :: USE_FORCE_POINT_SOURCE = .false.
double precision, parameter :: FACTOR_FORCE_SOURCE = 1.d15
-
+ integer, parameter :: COMPONENT_FORCE_SOURCE = 3 ! takes direction in comp E/N/Z = 1/2/3
+
! the receivers can be located inside the model
logical, parameter :: RECVS_CAN_BE_BURIED_EXT_MESH = .true.
logical, parameter :: SOURCES_CAN_BE_BURIED_EXT_MESH = .true.
Modified: seismo/3D/SPECFEM3D/trunk/iterate_time.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/iterate_time.f90 2010-09-02 00:38:38 UTC (rev 17161)
+++ seismo/3D/SPECFEM3D/trunk/iterate_time.f90 2010-09-02 02:41:23 UTC (rev 17162)
@@ -81,6 +81,12 @@
! poroelastic solver
if( POROELASTIC_SIMULATION ) stop 'poroelastic simulation not implemented yet'
+
+ ! restores last time snapshot saved for backward/reconstruction of wavefields
+ ! note: this must be read in after the Newark time scheme
+ if( SIMULATION_TYPE == 3 .and. it == 1 ) then
+ call it_read_foward_arrays()
+ endif
! write the seismograms with time shift
if (nrec_local > 0) then
@@ -338,10 +344,67 @@
end subroutine it_update_displacement_scheme
-
+
!=====================================================================
+ subroutine it_read_foward_arrays()
+
+ use specfem_par
+ use specfem_par_acoustic
+ use specfem_par_elastic
+ use specfem_par_poroelastic
+ implicit none
+
+ integer :: ier
+
+! restores last time snapshot saved for backward/reconstruction of wavefields
+! note: this is done here after the Newmark time scheme, otherwise the indexing for sources
+! and adjoint sources will become more complicated
+! that is, index it for adjoint sources will match index NSTEP - 1 for backward/reconstructed wavefields
+
+ ! reads in wavefields
+ open(unit=27,file=trim(prname)//'save_forward_arrays.bin',status='old',&
+ action='read',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error: opening save_forward_arrays'
+ print*,'path: ',trim(prname)//'save_forward_arrays.bin'
+ call exit_mpi(myrank,'error open file save_forward_arrays.bin')
+ endif
+
+ if( ACOUSTIC_SIMULATION ) then
+ read(27) b_potential_acoustic
+ read(27) b_potential_dot_acoustic
+ read(27) b_potential_dot_dot_acoustic
+ endif
+
+ ! elastic wavefields
+ if( ELASTIC_SIMULATION ) then
+ read(27) b_displ
+ read(27) b_veloc
+ read(27) b_accel
+
+ ! memory variables if attenuation
+ if( ATTENUATION ) then
+ read(27) b_R_xx
+ read(27) b_R_yy
+ read(27) b_R_xy
+ read(27) b_R_xz
+ read(27) b_R_yz
+ read(27) b_epsilondev_xx
+ read(27) b_epsilondev_yy
+ read(27) b_epsilondev_xy
+ read(27) b_epsilondev_xz
+ read(27) b_epsilondev_yz
+ endif
+
+ endif
+
+ close(27)
+
+ end subroutine it_read_foward_arrays
+!=====================================================================
+
subroutine it_store_attenuation_arrays()
! resetting d/v/a/R/eps for the backward reconstruction with attenuation
Modified: seismo/3D/SPECFEM3D/trunk/locate_source.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/locate_source.f90 2010-09-02 00:38:38 UTC (rev 17161)
+++ seismo/3D/SPECFEM3D/trunk/locate_source.f90 2010-09-02 02:41:23 UTC (rev 17162)
@@ -353,9 +353,12 @@
iz_initial_guess_source = k
! store xi,eta,gamma and x,y,z of point found
+ ! note: they have range [1.0d0,NGLLX/Y/Z], used for point sources
+ ! see e.g. in compute_add_source_elastic.f90
xi_source(isource) = dble(ix_initial_guess_source)
eta_source(isource) = dble(iy_initial_guess_source)
- gamma_source(isource) = dble(iz_initial_guess_source)
+ gamma_source(isource) = dble(iz_initial_guess_source)
+
x_found_source(isource) = xstore(iglob)
y_found_source(isource) = ystore(iglob)
z_found_source(isource) = zstore(iglob)
@@ -386,9 +389,15 @@
idomain(isource) = 0
endif
- ! get normal to the face of the hexaedra if receiver is on the surface
+ ! get normal to the face of the hexahedra if receiver is on the surface
if ((.not. SOURCES_CAN_BE_BURIED_EXT_MESH) .and. &
.not. (ispec_selected_source(isource) == 0)) then
+
+ ! note: at this point, xi_source,.. are in range [1.0d0,NGLLX/Y/Z] for point sources only,
+ ! for non-point sources the range is limited to [2.0d0,NGLLX/Y/Z - 1]
+ if( .not. USE_FORCE_POINT_SOURCE ) call exit_MPI(myrank,'error locate source: no point source at surface')
+
+ ! initialize indices
pt0_ix = -1
pt0_iy = -1
pt0_iz = -1
@@ -398,8 +407,9 @@
pt2_ix = -1
pt2_iy = -1
pt2_iz = -1
- ! we get two vectors of the face (three points) to compute the normal
- if (xi_source(isource) == 1 .and. &
+
+ ! we get two vectors of the face (three points) to compute the normal
+ if (nint(xi_source(isource)) == 1 .and. &
iglob_is_surface_external_mesh(ibool(1,2,2,ispec_selected_source(isource)))) then
pt0_ix = 1
pt0_iy = NGLLY
@@ -411,7 +421,7 @@
pt2_iy = NGLLY
pt2_iz = NGLLZ
endif
- if (xi_source(isource) == NGLLX .and. &
+ if (nint(xi_source(isource)) == NGLLX .and. &
iglob_is_surface_external_mesh(ibool(NGLLX,2,2,ispec_selected_source(isource)))) then
pt0_ix = NGLLX
pt0_iy = 1
@@ -423,7 +433,7 @@
pt2_iy = 1
pt2_iz = NGLLZ
endif
- if (eta_source(isource) == 1 .and. &
+ if (nint(eta_source(isource)) == 1 .and. &
iglob_is_surface_external_mesh(ibool(2,1,2,ispec_selected_source(isource)))) then
pt0_ix = 1
pt0_iy = 1
@@ -435,7 +445,7 @@
pt2_iy = 1
pt2_iz = NGLLZ
endif
- if (eta_source(isource) == NGLLY .and. &
+ if (nint(eta_source(isource)) == NGLLY .and. &
iglob_is_surface_external_mesh(ibool(2,NGLLY,2,ispec_selected_source(isource)))) then
pt0_ix = NGLLX
pt0_iy = NGLLY
@@ -447,7 +457,7 @@
pt2_iy = NGLLY
pt2_iz = NGLLZ
endif
- if (gamma_source(isource) == 1 .and. &
+ if (nint(gamma_source(isource)) == 1 .and. &
iglob_is_surface_external_mesh(ibool(2,2,1,ispec_selected_source(isource)))) then
pt0_ix = NGLLX
pt0_iy = 1
@@ -459,7 +469,7 @@
pt2_iy = NGLLY
pt2_iz = 1
endif
- if (gamma_source(isource) == NGLLZ .and. &
+ if (nint(gamma_source(isource)) == NGLLZ .and. &
iglob_is_surface_external_mesh(ibool(2,2,NGLLZ,ispec_selected_source(isource)))) then
pt0_ix = 1
pt0_iy = 1
@@ -475,7 +485,7 @@
if (pt0_ix<0 .or.pt0_iy<0 .or. pt0_iz<0 .or. &
pt1_ix<0 .or. pt1_iy<0 .or. pt1_iz<0 .or. &
pt2_ix<0 .or. pt2_iy<0 .or. pt2_iz<0) then
- stop 'error in computing normal for sources.'
+ call exit_mpi(myrank,'error in computing normal for sources.')
endif
u_vector(1) = xstore(ibool(pt1_ix,pt1_iy,pt1_iz,ispec_selected_source(isource))) &
@@ -519,15 +529,17 @@
nu_source(3,2,isource) = v_vector(3)
nu_source(3,3,isource) = w_vector(3)
- endif ! of if (.not. RECEIVERS_CAN_BE_BURIED_EXT_MESH)
+ endif ! of if (.not. SOURCES_CAN_BE_BURIED_EXT_MESH)
! *******************************************
! find the best (xi,eta,gamma) for the source
! *******************************************
+ ! for point sources, the location will be exactly at a GLL point
+ ! otherwise this tries to find best location
if(.not. USE_FORCE_POINT_SOURCE) then
- ! use initial guess in xi, eta and gamma
+ ! uses actual location interpolators, in range [-1,1]
xi = xigll(ix_initial_guess_source)
eta = yigll(iy_initial_guess_source)
gamma = zigll(iz_initial_guess_source)
@@ -746,12 +758,16 @@
write(IMAIN,*)
if(USE_FORCE_POINT_SOURCE) then
- write(IMAIN,*) ' xi coordinate of source in that element: ',nint(xi_source(isource))
- write(IMAIN,*) ' eta coordinate of source in that element: ',nint(eta_source(isource))
- write(IMAIN,*) ' gamma coordinate of source in that element: ',nint(gamma_source(isource))
+ write(IMAIN,*) ' i index of source in that element: ',nint(xi_source(isource))
+ write(IMAIN,*) ' j index of source in that element: ',nint(eta_source(isource))
+ write(IMAIN,*) ' k index of source in that element: ',nint(gamma_source(isource))
+ write(IMAIN,*)
+ write(IMAIN,*) ' component direction: ',COMPONENT_FORCE_SOURCE
+ write(IMAIN,*)
write(IMAIN,*) ' nu1 = ',nu_source(1,:,isource)
write(IMAIN,*) ' nu2 = ',nu_source(2,:,isource)
write(IMAIN,*) ' nu3 = ',nu_source(3,:,isource)
+ write(IMAIN,*)
write(IMAIN,*) ' at (x,y,z) coordinates = ',x_found_source(isource),y_found_source(isource),z_found_source(isource)
! prints frequency content for point forces
@@ -761,11 +777,10 @@
write(IMAIN,*) ' lambda_S at dominant frequency = ',3000./sqrt(3.)/f0
write(IMAIN,*) ' lambda_S at highest significant frequency = ',3000./sqrt(3.)/(2.5*f0)
write(IMAIN,*) ' t0 = ',t0_ricker,'t_cmt = ',t_cmt(isource)
-
else
write(IMAIN,*) ' xi coordinate of source in that element: ',xi_source(isource)
write(IMAIN,*) ' eta coordinate of source in that element: ',eta_source(isource)
- write(IMAIN,*) ' gamma coordinate of source in that element: ',gamma_source(isource)
+ write(IMAIN,*) ' gamma coordinate of source in that element: ',gamma_source(isource)
endif
! add message if source is a Heaviside
@@ -776,7 +791,11 @@
endif
write(IMAIN,*)
- write(IMAIN,*) ' half duration: ',hdur(isource),' seconds'
+ if(USE_FORCE_POINT_SOURCE) then
+ write(IMAIN,*) ' half duration -> frequency: ',hdur(isource),' seconds**(-1)'
+ else
+ write(IMAIN,*) ' half duration: ',hdur(isource),' seconds'
+ endif
write(IMAIN,*) ' time shift: ',t_cmt(isource),' seconds'
write(IMAIN,*)
Modified: seismo/3D/SPECFEM3D/trunk/prepare_timerun.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/prepare_timerun.f90 2010-09-02 00:38:38 UTC (rev 17161)
+++ seismo/3D/SPECFEM3D/trunk/prepare_timerun.f90 2010-09-02 02:41:23 UTC (rev 17162)
@@ -119,7 +119,7 @@
potential_dot_acoustic(:) = 0._CUSTOM_REAL
potential_dot_dot_acoustic(:) = 0._CUSTOM_REAL
! put negligible initial value to avoid very slow underflow trapping
- if(FIX_UNDERFLOW_PROBLEM) potential_dot_dot_acoustic(:) = VERYSMALLVAL
+ if(FIX_UNDERFLOW_PROBLEM) potential_acoustic(:) = VERYSMALLVAL
endif
! initialize elastic arrays to zero/verysmallvall
@@ -437,62 +437,48 @@
b_gammaval(:,:) = b_deltat / 2._CUSTOM_REAL + b_deltat**2*tauinv(:,:) / 6._CUSTOM_REAL &
+ b_deltat**3*tauinv(:,:)**2 / 24._CUSTOM_REAL
endif
-
-! kernel calculation, reads in last frame
- if (SIMULATION_TYPE == 3) then
- ! reads in wavefields
- open(unit=27,file=trim(prname)//'save_forward_arrays.bin',status='old',&
- action='read',form='unformatted',iostat=ier)
- if( ier /= 0 ) then
- print*,'error: opening save_forward_arrays'
- print*,'path: ',trim(prname)//'save_forward_arrays.bin'
- call exit_mpi(myrank,'error open file save_forward_arrays.bin')
- endif
- if( ACOUSTIC_SIMULATION ) then
- read(27) b_potential_acoustic
- read(27) b_potential_dot_acoustic
- read(27) b_potential_dot_dot_acoustic
- endif
-
- ! elastic wavefields
- if( ELASTIC_SIMULATION ) then
- read(27) b_displ
- read(27) b_veloc
- read(27) b_accel
-
- ! memory variables if attenuation
- if( ATTENUATION ) then
- read(27) b_R_xx
- read(27) b_R_yy
- read(27) b_R_xy
- read(27) b_R_xz
- read(27) b_R_yz
- read(27) b_epsilondev_xx
- read(27) b_epsilondev_yy
- read(27) b_epsilondev_xy
- read(27) b_epsilondev_xz
- read(27) b_epsilondev_yz
- endif
-
- endif
-
- close(27)
- endif
-
-! initializes adjoint kernels
+! initializes adjoint kernels and reconstructed/backward wavefields
if (SIMULATION_TYPE == 3) then
! elastic domain
if( ELASTIC_SIMULATION ) then
rho_kl(:,:,:,:) = 0._CUSTOM_REAL
mu_kl(:,:,:,:) = 0._CUSTOM_REAL
kappa_kl(:,:,:,:) = 0._CUSTOM_REAL
+
+ ! reconstructed/backward elastic wavefields
+ b_displ = 0._CUSTOM_REAL
+ b_veloc = 0._CUSTOM_REAL
+ b_accel = 0._CUSTOM_REAL
+ if(FIX_UNDERFLOW_PROBLEM) b_displ = VERYSMALLVAL
+
+ ! memory variables if attenuation
+ if( ATTENUATION ) then
+ b_R_xx = 0._CUSTOM_REAL
+ b_R_yy = 0._CUSTOM_REAL
+ b_R_xy = 0._CUSTOM_REAL
+ b_R_xz = 0._CUSTOM_REAL
+ b_R_yz = 0._CUSTOM_REAL
+ b_epsilondev_xx = 0._CUSTOM_REAL
+ b_epsilondev_yy = 0._CUSTOM_REAL
+ b_epsilondev_xy = 0._CUSTOM_REAL
+ b_epsilondev_xz = 0._CUSTOM_REAL
+ b_epsilondev_yz = 0._CUSTOM_REAL
+ endif
+
endif
! acoustic domain
if( ACOUSTIC_SIMULATION ) then
rho_ac_kl(:,:,:,:) = 0._CUSTOM_REAL
kappa_ac_kl(:,:,:,:) = 0._CUSTOM_REAL
+
+ ! reconstructed/backward acoustic potentials
+ b_potential_acoustic = 0._CUSTOM_REAL
+ b_potential_dot_acoustic = 0._CUSTOM_REAL
+ b_potential_dot_dot_acoustic = 0._CUSTOM_REAL
+ if(FIX_UNDERFLOW_PROBLEM) b_potential_acoustic = VERYSMALLVAL
+
endif
endif
Modified: seismo/3D/SPECFEM3D/trunk/setup_sources_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/setup_sources_receivers.f90 2010-09-02 00:38:38 UTC (rev 17161)
+++ seismo/3D/SPECFEM3D/trunk/setup_sources_receivers.f90 2010-09-02 02:41:23 UTC (rev 17162)
@@ -449,7 +449,7 @@
integer :: isource,ispec
integer :: irec,irec_local
integer :: icomp,itime,nadj_files_found,nadj_files_found_tot,ier
- character(len=3),dimension(NDIM) :: comp = (/ "BHN", "BHE", "BHZ" /)
+ character(len=3),dimension(NDIM) :: comp = (/ "BHE", "BHN", "BHZ" /)
character(len=150) :: filename
@@ -491,7 +491,7 @@
! (and getting rid of 1/sqrt(2) factor from scalar moment tensor definition above)
factor_source = factor_source * sqrt(2.0) / sqrt(3.0)
- ! source array interpolated on all element gll points
+ ! source array interpolated on all element gll points (only used for non point sources)
call compute_arrays_source_acoustic(xi_source(isource),eta_source(isource),gamma_source(isource),&
sourcearray,xigll,yigll,zigll,factor_source)
endif
@@ -711,9 +711,17 @@
if( myrank == 0 ) then
! get the 3-D shape functions
- xil = xi_source(isource)
- etal = eta_source(isource)
- gammal = gamma_source(isource)
+ if( USE_FORCE_POINT_SOURCE ) then
+ ! note: we switch xi,eta,gamma range to be [-1,1]
+ ! uses initial guess in xi, eta and gamma
+ xil = xigll(nint(xi_source(isource)))
+ etal = yigll(nint(eta_source(isource)))
+ gammal = zigll(nint(gamma_source(isource)))
+ else
+ xil = xi_source(isource)
+ etal = eta_source(isource)
+ gammal = gamma_source(isource)
+ endif
call get_shape3D_single(myrank,shape3D,xil,etal,gammal)
! interpolates source locations
Modified: seismo/3D/SPECFEM3D/trunk/todo_list_please_dont_remove.txt
===================================================================
--- seismo/3D/SPECFEM3D/trunk/todo_list_please_dont_remove.txt 2010-09-02 00:38:38 UTC (rev 17161)
+++ seismo/3D/SPECFEM3D/trunk/todo_list_please_dont_remove.txt 2010-09-02 02:41:23 UTC (rev 17162)
@@ -2,6 +2,12 @@
To-do list for SPECFEM3D, by Dimitri Komatitsch
-----------------------------------------------
+NOTE (by Daniel): the 3D code here uses a potential of (density * displacement),
+ identical to the 2D code (see also the note in compute_forces_acoustic.f90).
+ thus, the remark below is partly obsolete, having
+ acoustic-acoustic discontinuities in your 3D model is fine.
+
+
Date: Fri, 23 Jul 2010 01:30:16 +0200
From: Dimitri Komatitsch <dimitri.komatitsch at univ-pau.fr>
To: Anne Sieminski <anne.sieminski at obs.ujf-grenoble.fr>
@@ -52,13 +58,14 @@
Dimitri.
+
On 06/23/2010 05:42 PM, Anne Sieminski wrote:
> Bonjour Dimitri,
>
> Je vois que SPECFEM2D-6.0.0 (disponible sur le site de CIG) permet de
-> faire des calculs dans des modèles acoustiques. La même option est-elle
+> faire des calculs dans des modËles acoustiques. La mÃme option est-elle
> en projet pour SPECFEM3D (ancien SPECFEM3D_BASIN) ? J'ai un peu perdu le
-> fil de toutes les nouveautés... Désolée !
+> fil de toutes les nouveautÈs... DÈsolÈe !
>
> Merci,
> Anne
@@ -115,7 +122,7 @@
- should add support for parMETIS or PT-Scotch (i.e. decompose the mesh in parallel rather than serial; this will be required for very large meshes, too large to fit on a serial machine)
-- etudier comment lire un maillage CUBIT stocke au format natif de CUBIT (Exodus, base sur netCDF) depuis SPECFEM3D. On pourrait utiliser la commande "ncdump" dans CUBIT si necessaire d'après Emanuele Casarotti. Une autre option serait d'utiliser le format de stockage ABAQUS dans CUBIT.
+- etudier comment lire un maillage CUBIT stocke au format natif de CUBIT (Exodus, base sur netCDF) depuis SPECFEM3D. On pourrait utiliser la commande "ncdump" dans CUBIT si necessaire d'aprËs Emanuele Casarotti. Une autre option serait d'utiliser le format de stockage ABAQUS dans CUBIT.
by Nicolas Le Goff :
Modified: seismo/3D/SPECFEM3D/trunk/write_seismograms.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/write_seismograms.f90 2010-09-02 00:38:38 UTC (rev 17161)
+++ seismo/3D/SPECFEM3D/trunk/write_seismograms.f90 2010-09-02 02:41:23 UTC (rev 17162)
@@ -365,9 +365,9 @@
! distinguish between single and double precision for reals
! note: compare time_t with time used for source term
if(CUSTOM_REAL == SIZE_REAL) then
- time_t = sngl( dble(NSTEP-isample-1)*DT - t0 )
+ time_t = sngl( dble(NSTEP-isample)*DT - t0 )
else
- time_t = dble(NSTEP-isample-1)*DT - t0
+ time_t = dble(NSTEP-isample)*DT - t0
endif
endif
@@ -498,9 +498,9 @@
! distinguish between single and double precision for reals
! note: compare time_t with time used for source term
if(CUSTOM_REAL == SIZE_REAL) then
- time_t = sngl( dble(NSTEP-isample-1)*DT - t0 )
+ time_t = sngl( dble(NSTEP-isample)*DT - t0 )
else
- time_t = dble(NSTEP-isample-1)*DT - t0
+ time_t = dble(NSTEP-isample)*DT - t0
endif
endif
More information about the CIG-COMMITS
mailing list