[cig-commits] [commit] devel: added option EXACT_UNDOING_TO_DISK, which saves the forward run to disk in order to compute exact sensitivity kernels without having to undo the forward run numerically; it will read it back from disk instead. This is for validation purposes only, since the disk files will typically require tens of terabytes (yes, tera). However this will become affordable uin a few years. (4e3f2c6)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Jul 10 08:52:46 PDT 2014
Repository : https://github.com/geodynamics/specfem3d_globe
On branch : devel
Link : https://github.com/geodynamics/specfem3d_globe/compare/b17a1fe8465ca1e8307c6f39ee901f4704f266ea...4e3f2c683cec88a4cbf79a4b91b5fec24b215868
>---------------------------------------------------------------
commit 4e3f2c683cec88a4cbf79a4b91b5fec24b215868
Author: Dimitri Komatitsch <komatitsch at lma.cnrs-mrs.fr>
Date: Thu Jul 10 17:42:47 2014 +0200
added option EXACT_UNDOING_TO_DISK, which saves the forward run to disk in order to compute exact sensitivity kernels without having to undo the forward run numerically; it will read it back from disk instead. This is for validation purposes only, since the disk files will typically require tens of terabytes (yes, tera). However this will become affordable uin a few years.
>---------------------------------------------------------------
4e3f2c683cec88a4cbf79a4b91b5fec24b215868
setup/constants.h.in | 12 +++
src/shared/read_compute_parameters.f90 | 1 +
src/shared/save_header_file.F90 | 2 +-
.../compute_forces_crust_mantle_noDev.f90 | 2 -
src/specfem3D/compute_kernels.F90 | 9 ++
src/specfem3D/initialize_simulation.f90 | 9 +-
src/specfem3D/iterate_time.F90 | 111 ++++++++++++++++++++-
7 files changed, 139 insertions(+), 7 deletions(-)
diff --git a/setup/constants.h.in b/setup/constants.h.in
index 394de9d..f99c8c2 100644
--- a/setup/constants.h.in
+++ b/setup/constants.h.in
@@ -468,6 +468,18 @@
integer, parameter :: NSTEP_FOR_BENCHMARK = 300
logical, parameter :: SET_INITIAL_FIELD_TO_1_IN_BENCH = .true.
+! compute an alpha sensitivity kernel directly by dumping the whole forward
+! run to disk instead of rebuilding it backwards from the final time step in a second stage;
+! used for debugging and validation purposes only, in order to have an exact reference for the sensitivity kernels.
+! use wisely, this requires several terabytes (yes, tera) of disk space.
+! In the future that option should become standard at some point though.
+! For now, in order to save disk space, it is implemented for the alpha kernel only
+! (because it only requires saving a scalar: the trace of epsilon) and for surface wave
+! kernels only, this way we can save the upper mantle only and ignore the rest.
+! In the future it will be easy to generalize this though.
+ logical, parameter :: EXACT_UNDOING_TO_DISK = .false.
+ ! ID of the huge file in which we dump all the time steps of the simulation
+ integer, parameter :: IFILE_FOR_EXACT_UNDOING = 244
!------------------------------------------------------
!----------- do not modify anything below -------------
diff --git a/src/shared/read_compute_parameters.f90 b/src/shared/read_compute_parameters.f90
index 1ad5e16..1a963c9 100644
--- a/src/shared/read_compute_parameters.f90
+++ b/src/shared/read_compute_parameters.f90
@@ -298,6 +298,7 @@
if( UNDO_ATTENUATION .and. NUMBER_OF_RUNS /= 1) &
stop 'NUMBER_OF_RUNS should be == 1 for now when using UNDO_ATTENUATION'
+ !! DK DK this should not be difficult to fix and test, but not done yet by lack of time
if(UNDO_ATTENUATION .and. NUMBER_OF_THIS_RUN > 1) &
stop 'we currently do not support NUMBER_OF_THIS_RUN > 1 in the case of UNDO_ATTENUATION'
diff --git a/src/shared/save_header_file.F90 b/src/shared/save_header_file.F90
index e71cb24..65247cc 100644
--- a/src/shared/save_header_file.F90
+++ b/src/shared/save_header_file.F90
@@ -382,7 +382,7 @@
write(IOUT,*) 'integer, parameter :: NGLOB_CRUST_MANTLE_OCEANS = ',NGLOB_CRUST_MANTLE_OCEANS
write(IOUT,*)
-! this to allow for code elimination by compiler in solver for performance
+! this to allow for code elimination by the compiler in the solver for performance
if(TRANSVERSE_ISOTROPY) then
write(IOUT,*) 'logical, parameter :: TRANSVERSE_ISOTROPY_VAL = .true.'
diff --git a/src/specfem3D/compute_forces_crust_mantle_noDev.f90 b/src/specfem3D/compute_forces_crust_mantle_noDev.f90
index 6280a7b..3588b98 100644
--- a/src/specfem3D/compute_forces_crust_mantle_noDev.f90
+++ b/src/specfem3D/compute_forces_crust_mantle_noDev.f90
@@ -347,8 +347,6 @@
else
! do not use transverse isotropy except if element is between d220 and Moho
-! if(.not. (TRANSVERSE_ISOTROPY_VAL .and. (idoubling(ispec)==IFLAG_220_80 .or. idoubling(ispec)==IFLAG_80_MOHO))) then
-
if( .not. ispec_is_tiso(ispec) ) then
! isotropic element
diff --git a/src/specfem3D/compute_kernels.F90 b/src/specfem3D/compute_kernels.F90
index 9fb28fc..fb8114f 100644
--- a/src/specfem3D/compute_kernels.F90
+++ b/src/specfem3D/compute_kernels.F90
@@ -38,6 +38,9 @@
! crust mantle
call compute_kernels_crust_mantle()
+ ! only information to compute the crust_mantle kernels was saved to disk
+ if(EXACT_UNDOING_TO_DISK) return
+
! outer core
call compute_kernels_outer_core(vector_displ_outer_core,vector_accel_outer_core,b_vector_displ_outer_core, &
displ_outer_core,accel_outer_core,b_displ_outer_core,b_accel_outer_core, &
@@ -87,6 +90,7 @@
integer :: i,j,k,ispec,iglob
if( .not. GPU_MODE ) then
+
! on CPU
! crust_mantle
do ispec = 1, NSPEC_CRUST_MANTLE
@@ -175,6 +179,11 @@
! isotropic kernels
+ ! if EXACT_UNDOING_TO_DISK is set, currently the rho and beta kernels below will be computed but will be equal
+ ! to zero everywhere because the backward field is then not computed, and only the alpha kernel
+ ! will be computed correctly based on the values of the forward run saved to disk in the first run
+ ! and read back at the beginning of this routine
+
do k = 1, NGLLZ
do j = 1, NGLLY
do i = 1, NGLLX
diff --git a/src/specfem3D/initialize_simulation.f90 b/src/specfem3D/initialize_simulation.f90
index 6b0af64..ccfbf9f 100644
--- a/src/specfem3D/initialize_simulation.f90
+++ b/src/specfem3D/initialize_simulation.f90
@@ -362,13 +362,16 @@
call exit_MPI(myrank,'cannot have both UNDO_ATTENUATION and PARTIAL_PHYS_DISPERSION_ONLY, please check Par_file...')
if((SIMULATION_TYPE == 1 .and. SAVE_FORWARD) .or. SIMULATION_TYPE == 3) then
+
if ( ATTENUATION_VAL) then
! checks mimic flag:
! attenuation for adjoint simulations must have PARTIAL_PHYS_DISPERSION_ONLY set by xcreate_header_file
- if(.not. UNDO_ATTENUATION )then
- if(.not. PARTIAL_PHYS_DISPERSION_ONLY) then
- call exit_MPI(myrank, &
+ if(.not. EXACT_UNDOING_TO_DISK)then
+ if(.not. UNDO_ATTENUATION )then
+ if(.not. PARTIAL_PHYS_DISPERSION_ONLY) then
+ call exit_MPI(myrank, &
'ATTENUATION for adjoint runs or SAVE_FORWARD requires UNDO_ATTENUATION or PARTIAL_PHYS_DISPERSION_ONLY')
+ endif
endif
endif
diff --git a/src/specfem3D/iterate_time.F90 b/src/specfem3D/iterate_time.F90
index 0356c81..a60d91f 100644
--- a/src/specfem3D/iterate_time.F90
+++ b/src/specfem3D/iterate_time.F90
@@ -38,6 +38,13 @@
! timing
double precision, external :: wtime
+ ! for EXACT_UNDOING_TO_DISK
+ integer :: ispec,iglob,i,j,k,counter,record_length
+ real(kind=CUSTOM_REAL) :: radius
+ integer, dimension(:), allocatable :: integer_mask_ibool_exact_undo
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: buffer_for_disk
+ character(len=150) outputname
+
!
! s t a r t t i m e i t e r a t i o n s
!
@@ -70,6 +77,69 @@
! ************* MAIN LOOP OVER THE TIME STEPS *************
! *********************************************************
+ if(EXACT_UNDOING_TO_DISK) then
+
+ if(GPU_MODE) call exit_MPI(myrank,'EXACT_UNDOING_TO_DISK not supported for GPUs')
+
+ if(UNDO_ATTENUATION) &
+ call exit_MPI(myrank,'EXACT_UNDOING_TO_DISK needs UNDO_ATTENUATION to be off because it computes the kernel directly instead')
+
+ if(SIMULATION_TYPE == 1 .and. .not. SAVE_FORWARD) &
+ call exit_MPI(myrank,'EXACT_UNDOING_TO_DISK requires SAVE_FORWARD if SIMULATION_TYPE == 1')
+
+ if(ANISOTROPIC_KL) call exit_MPI(myrank,'EXACT_UNDOING_TO_DISK requires ANISOTROPIC_KL to be turned off')
+
+!! DK DK determine the largest value of iglob that we need to save to disk,
+!! DK DK since we save the upper mantle only in the case of surface-wave kernels
+ ! crust_mantle
+ allocate(integer_mask_ibool_exact_undo(NGLOB_CRUST_MANTLE))
+ integer_mask_ibool_exact_undo(:) = -1
+
+ counter = 0
+ do ispec = 1, NSPEC_CRUST_MANTLE
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ iglob = ibool_crust_mantle(i,j,k,ispec)
+ ! xstore ystore zstore have previously been converted to r theta phi, thus xstore now stores the radius
+ radius = xstore_crust_mantle(iglob) ! <- radius r (normalized)
+ ! save that element only if it is in the upper mantle
+ if(radius >= R670 / R_EARTH) then
+ ! if this point has not yet been found before
+ if(integer_mask_ibool_exact_undo(iglob) == -1) then
+ ! create a new unique point
+ counter = counter + 1
+ integer_mask_ibool_exact_undo(iglob) = counter
+ endif
+ endif
+ enddo
+ enddo
+ enddo
+ enddo
+
+! print *,'myrank, counter, NGLOB, ratio = ',myrank, counter, NGLOB_CRUST_MANTLE, &
+! real(counter) / NGLOB_CRUST_MANTLE
+! print *,'myrank, R670, ratioR670surR_EARTH = ',myrank, R670, R670 / R_EARTH
+
+ ! allocate the buffer used to dump a single time step
+ allocate(buffer_for_disk(counter))
+
+ ! open the file in which we will dump all the time steps (in a single file)
+ write(outputname,"('huge_dumps/proc',i6.6,'_huge_dump_of_all_time_steps.bin')") myrank
+ inquire(iolength=record_length) buffer_for_disk
+ ! we write to or read from the file depending on the simulation type
+ if(SIMULATION_TYPE == 1) then
+ open(file=outputname, unit=IFILE_FOR_EXACT_UNDOING, action='write', status='unknown', &
+ form='unformatted', access='direct', recl=record_length)
+ else if(SIMULATION_TYPE == 3) then
+ open(file=outputname, unit=IFILE_FOR_EXACT_UNDOING, action='read', status='old', &
+ form='unformatted', access='direct', recl=record_length)
+ else
+ call exit_MPI(myrank,'EXACT_UNDOING_TO_DISK can only be used with SIMULATION_TYPE == 1 or SIMULATION_TYPE == 3')
+ endif
+
+ endif ! of if(EXACT_UNDOING_TO_DISK)
+
do it = it_begin,it_end
! simulation status output and stability check
@@ -96,9 +166,26 @@
enddo ! end of very big external loop on istage for all the stages of the LDDRK time scheme (only one stage if Newmark)
+ ! save the forward run to disk for the alpha kernel only
+ if(EXACT_UNDOING_TO_DISK .and. SIMULATION_TYPE == 1) then
+ do ispec = 1, NSPEC_CRUST_MANTLE
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ iglob = ibool_crust_mantle(i,j,k,ispec)
+ if(integer_mask_ibool_exact_undo(iglob) /= -1) &
+ buffer_for_disk(integer_mask_ibool_exact_undo(iglob)) = eps_trace_over_3_crust_mantle(i,j,k,ispec)
+ enddo
+ enddo
+ enddo
+ enddo
+ write(IFILE_FOR_EXACT_UNDOING,rec=it) buffer_for_disk
+ endif
+
! kernel simulations (forward and adjoint wavefields)
- if( SIMULATION_TYPE == 3 ) then
+ if(SIMULATION_TYPE == 3) then
+ if(.not. EXACT_UNDOING_TO_DISK) then
! note: we step back in time (using time steps - DT ), i.e. wavefields b_displ_..() are time-reversed here
! reconstructs forward wavefields based on last stored wavefield data
@@ -129,6 +216,25 @@
call read_forward_arrays()
endif
+ else ! of if(.not. EXACT_UNDOING_TO_DISK)
+
+ ! read the forward run from disk for the alpha kernel only
+ ! here we time revert the forward run by reading time step NSTEP - it + 1
+ read(IFILE_FOR_EXACT_UNDOING,rec=NSTEP-it+1) buffer_for_disk
+ do ispec = 1, NSPEC_CRUST_MANTLE
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ iglob = ibool_crust_mantle(i,j,k,ispec)
+ if(integer_mask_ibool_exact_undo(iglob) /= -1) &
+ b_eps_trace_over_3_crust_mantle(i,j,k,ispec) = buffer_for_disk(integer_mask_ibool_exact_undo(iglob))
+ enddo
+ enddo
+ enddo
+ enddo
+
+ endif ! of if(.not. EXACT_UNDOING_TO_DISK)
+
! adjoint simulations: kernels
call compute_kernels()
@@ -159,6 +265,9 @@
enddo ! end of main time loop
+ ! close the huge file that contains a dump of all the time steps to disk
+ if(EXACT_UNDOING_TO_DISK) close(IFILE_FOR_EXACT_UNDOING)
+
!
!---- end of time iteration loop
!
More information about the CIG-COMMITS
mailing list