supports ASDF and SU format outputs together with WRITE_SEISMOGRAMS_B…

…Y_MAIN

src/cuda/write_seismograms_cuda.cu

@@ -537,8 +537,10 @@ void FC_FUNC_(compute_seismograms_cuda,
     realw * seismo_temp;
     if (mp->save_seismograms_p){
       // EB EB We need to reorganize data to match host array shape :
-      // if NB_RUNS_ACOUSTIC_GPU = 1 from fortran shape (1,nrec_local,NTSTEP_BETWEEN_OUTPUT_SEISMOS) to (NDIM,nrec_local,NTSTEP_BETWEEN_OUTPUT_SEISMOS)
-      // if NB_RUNS_ACOUSTIC_GPU > 1 from fortran shape (NB_RUNS_ACOUSTIC_GPU,nrec_local,NTSTEP_BETWEEN_OUTPUT_SEISMOS) to (NDIM,nrec_local*NB_RUNS_ACOUSTIC_GPU,NTSTEP_BETWEEN_OUTPUT_SEISMOS)
+      // if NB_RUNS_ACOUSTIC_GPU = 1:
+      //   from fortran shape (1,nrec_local,NTSTEP_BETWEEN_OUTPUT_SEISMOS) to (NDIM,nrec_local,NTSTEP_BETWEEN_OUTPUT_SEISMOS)
+      // if NB_RUNS_ACOUSTIC_GPU > 1:
+      //    from fortran shape (NB_RUNS_ACOUSTIC_GPU,nrec_local,NTSTEP_BETWEEN_OUTPUT_SEISMOS) to (NDIM,nrec_local*NB_RUNS_ACOUSTIC_GPU,NTSTEP_BETWEEN_OUTPUT_SEISMOS)
       seismo_temp = (realw*)malloc(size * NB_RUNS_ACOUSTIC_GPU * sizeof(realw));
       print_CUDA_error_if_any(cudaMemcpy(seismo_temp,mp->d_seismograms_p,
                                          size * NB_RUNS_ACOUSTIC_GPU * sizeof(realw),cudaMemcpyDeviceToHost),72004);

src/meshfem3D/save_databases_adios.F90

@@ -44,7 +44,9 @@ subroutine save_databases_adios(LOCAL_PATH,sizeprocs, &
                                 NMATERIALS,material_properties, &
                                 nspec_CPML,CPML_to_spec,CPML_regions,is_CPML)
 
-  use constants, only: MAX_STRING_LEN,IDOMAIN_ACOUSTIC,IDOMAIN_ELASTIC,ADIOS_TRANSPORT_METHOD, &
+  use constants, only: MAX_STRING_LEN, &
+    IDOMAIN_ACOUSTIC,IDOMAIN_ELASTIC,IDOMAIN_POROELASTIC, &
+    ADIOS_TRANSPORT_METHOD, &
     NGLLX,NGLLY,NGLLZ,NDIM,myrank
 
   use constants_meshfem3D, only: NGLLX_M,NGLLY_M,NGLLZ_M,NUMBER_OF_MATERIAL_PROPERTIES

src/specfem3D/asdf_data.f90

@@ -45,7 +45,7 @@ module asdf_data
 
     !size info
     integer :: nreceivers
-    integer :: nrec_local
+    integer :: nrec_store
 
     !seismic record info
     real, pointer :: receiver_lat(:), receiver_lo(:)

src/specfem3D/compute_seismograms.f90

@@ -36,12 +36,14 @@ subroutine compute_seismograms()
     hxir_store,hetar_store,hgammar_store,number_receiver_global,nrec_local, &
     nu_source,nu_rec,Mxx,Myy,Mzz,Mxy,Mxz,Myz,tshift_src,hdur_Gaussian, &
     hprime_xx,hprime_yy,hprime_zz, &
-    seismograms_d,seismograms_v,seismograms_a,seismograms_p, &
     hpxir_store,hpetar_store,hpgammar_store,seismograms_eps, &
     Mxx_der,Myy_der,Mzz_der,Mxy_der,Mxz_der,Myz_der,sloc_der, &
     USE_TRICK_FOR_BETTER_PRESSURE, &
     SAVE_SEISMOGRAMS_DISPLACEMENT,SAVE_SEISMOGRAMS_VELOCITY,SAVE_SEISMOGRAMS_ACCELERATION,SAVE_SEISMOGRAMS_PRESSURE
 
+  ! seismograms
+  use specfem_par, only: seismograms_d,seismograms_v,seismograms_a,seismograms_p
+
   use specfem_par_acoustic, only: ispec_is_acoustic,potential_acoustic,potential_dot_acoustic,potential_dot_dot_acoustic, &
     b_potential_acoustic,b_potential_dot_acoustic,b_potential_dot_dot_acoustic
   use specfem_par_elastic, only: ispec_is_elastic,displ,veloc,accel, &

src/specfem3D/initialize_simulation.F90

@@ -419,17 +419,28 @@ subroutine initialize_simulation_check()
 
   ! safety check
   if (NB_RUNS_ACOUSTIC_GPU > 1) then
-
-   if (NUMBER_OF_SIMULTANEOUS_RUNS > 1) stop 'NB_RUNS_ACOUSTIC_GPU > 1 not compatible with NUMBER_OF_SIMULTANEOUS_RUNS > 1 yet'
-   if (SIMULATION_TYPE /= 1) stop 'NB_RUNS_ACOUSTIC_GPU > 1 not compatible with SIMULATION_TYPE /= 1 yet'
-   if (STACEY_ABSORBING_CONDITIONS) stop 'NB_RUNS_ACOUSTIC_GPU > 1 not compatible with STACEY_ABSORBING_CONDITIONS yet'
-   if (.not. SAVE_SEISMOGRAMS_PRESSURE ) stop 'NB_RUNS_ACOUSTIC_GPU > 1 not compatible with elastic wavefield seismograms yet'
-   if (.not. GPU_MODE ) stop 'NB_RUNS_ACOUSTIC_GPU > 1 only applies with GPU_MODE'
-   if (INVERSE_FWI_FULL_PROBLEM) stop 'NB_RUNS_ACOUSTIC_GPU > 1 not compatible with INVERSE_FWI_FULL_PROBLEM yet'
-   if (myrank == 1) stop 'NB_RUNS_ACOUSTIC_GPU > 1 not compatible with MPI mode yet'
-
+    if (NUMBER_OF_SIMULTANEOUS_RUNS > 1) &
+      stop 'NB_RUNS_ACOUSTIC_GPU > 1 not compatible with NUMBER_OF_SIMULTANEOUS_RUNS > 1 yet'
+    if (SIMULATION_TYPE /= 1) &
+      stop 'NB_RUNS_ACOUSTIC_GPU > 1 not compatible with SIMULATION_TYPE /= 1 yet'
+    if (STACEY_ABSORBING_CONDITIONS) &
+      stop 'NB_RUNS_ACOUSTIC_GPU > 1 not compatible with STACEY_ABSORBING_CONDITIONS yet'
+    if (SAVE_SEISMOGRAMS_DISPLACEMENT .or. SAVE_SEISMOGRAMS_VELOCITY .or. SAVE_SEISMOGRAMS_ACCELERATION) &
+      stop 'Invalid seismogram output for NB_RUNS_ACOUSTIC_GPU > 1, only pressure output implemented yet'
+    if (.not. SAVE_SEISMOGRAMS_PRESSURE ) &
+      stop 'NB_RUNS_ACOUSTIC_GPU > 1 not compatible with elastic wavefield seismograms yet'
+    if (.not. GPU_MODE ) &
+      stop 'NB_RUNS_ACOUSTIC_GPU > 1 only applies with GPU_MODE'
+    if (INVERSE_FWI_FULL_PROBLEM) &
+      stop 'NB_RUNS_ACOUSTIC_GPU > 1 not compatible with INVERSE_FWI_FULL_PROBLEM yet'
+    if (myrank == 1) &
+      stop 'NB_RUNS_ACOUSTIC_GPU > 1 not compatible with MPI mode yet'
   endif
 
+  ! file output
+  if (SU_FORMAT .and. ASDF_FORMAT) &
+    stop 'Please choose either SU_FORMAT or ASDF_FORMAT, both outputs together are not implemented yet...'
+
   end subroutine initialize_simulation_check
 
 !

src/specfem3D/noise_tomography.f90

@@ -39,10 +39,6 @@
 !
 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ATTENTION !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
-! =============================================================================================================
-! =============================================================================================================
-! =============================================================================================================
-
 module user_noise_distribution
 
 !daniel: TODO -- setting USE_PIERO_DISTRIBUTION = .true. will produce errors
@@ -56,8 +52,8 @@ module user_noise_distribution
 ! this subroutine must be modified by USERS for their own noise distribution
 
   subroutine noise_distribution_direction(xcoord_in,ycoord_in,zcoord_in, &
-                  normal_x_noise_out,normal_y_noise_out,normal_z_noise_out, &
-                  mask_noise_out)
+                                          normal_x_noise_out,normal_y_noise_out,normal_z_noise_out, &
+                                          mask_noise_out)
 
   use constants
 
@@ -92,8 +88,8 @@ end subroutine noise_distribution_direction
 !
 ! USERS: need to modify this subroutine for their own noise characteristics
   subroutine noise_distribution_direction_d(xcoord_in,ycoord_in,zcoord_in, &
-                                           normal_x_noise_out,normal_y_noise_out,normal_z_noise_out, &
-                                           mask_noise_out)
+                                            normal_x_noise_out,normal_y_noise_out,normal_z_noise_out, &
+                                            mask_noise_out)
   use constants
 
   implicit none
@@ -148,9 +144,9 @@ subroutine noise_distribution_dir_non_uni(xcoord_in,ycoord_in,zcoord_in, &
 
   ! coordinates "x/y/zcoord_in" actually contain r theta phi, therefore convert back to x y z
   ! call rthetaphi_2_xyz(xcoord,ycoord,zcoord, xcoord_in,ycoord_in,zcoord_in)
-  xcoord=xcoord_in
-  ycoord=ycoord_in
-  zcoord=zcoord_in
+  xcoord = xcoord_in
+  ycoord = ycoord_in
+  zcoord = zcoord_in
 
   !PB NOT UNIF DISTRIBUTION OF NOISE ON THE SURFACE OF A SPHERE
   !PB lon lat colat ARE IN RADIANS SINCE ARE OBTAINED FROM Cartesian COORDINATES
@@ -227,11 +223,9 @@ end subroutine noise_distribution_dir_non_uni
 
 end module user_noise_distribution
 
-
-
-! =============================================================================================================
-! =============================================================================================================
+!
 ! =============================================================================================================
+!
 
 ! read parameters
   subroutine read_parameters_noise(nrec,NSTEP,nmovie_points, &
@@ -367,14 +361,15 @@ subroutine read_parameters_noise(nrec,NSTEP,nmovie_points, &
 
   end subroutine read_parameters_noise
 
+!
 ! =============================================================================================================
-! =============================================================================================================
-! =============================================================================================================
+!
 
-! check for consistency of the parameters
   subroutine check_parameters_noise(NOISE_TOMOGRAPHY,SIMULATION_TYPE,SAVE_FORWARD, &
                                     LOCAL_PATH,NSPEC_TOP,NSTEP)
 
+! check for consistency of the parameters
+
   use constants, only: CUSTOM_REAL,NDIM,NGLLSQUARE,MAX_STRING_LEN,IOUT_NOISE,OUTPUT_FILES,myrank
 
   implicit none
@@ -494,8 +489,9 @@ subroutine compute_arrays_source_noise(xi_noise,eta_noise,gamma_noise, &
   double precision, dimension(NDIM,NDIM) :: nu_single  ! rotation matrix at the main receiver
   ! output parameters
   real(kind=CUSTOM_REAL) :: noise_sourcearray(NDIM,NGLLX,NGLLY,NGLLZ,NSTEP)
+
   ! local parameters
-  integer itime, i, j, k, ier, nlines
+  integer :: itime, i, j, k, ier, nlines
   real(kind=CUSTOM_REAL) :: junk
   real(kind=CUSTOM_REAL) :: noise_src(NSTEP),noise_src_u(NDIM,NSTEP)
   double precision, dimension(NDIM) :: nu_main       ! component direction chosen at the main receiver
@@ -512,12 +508,13 @@ subroutine compute_arrays_source_noise(xi_noise,eta_noise,gamma_noise, &
   ! main receiver component direction, \nu_main
   filename = trim(OUTPUT_FILES)//'/..//NOISE_TOMOGRAPHY/nu_main'
   open(unit=IIN_NOISE,file=trim(filename),status='old',action='read',iostat=ier)
-  if (ier /= 0 .and. myrank == 0) &
+  if (ier /= 0 .and. myrank == 0) then
     call exit_MPI(myrank, &
       'file '//trim(filename)//' does NOT exist! nu_main is the component direction (ENZ) for main receiver')
+  endif
 
-  do itime = 1,3
-    read(IIN_NOISE,*,iostat=ier) nu_main(itime)
+  do i = 1,3
+    read(IIN_NOISE,*,iostat=ier) nu_main(i)
     if (ier /= 0 .and. myrank == 0) &
       call exit_MPI(myrank, &
         'file '//trim(filename)//' has wrong length, the vector should have three components (ENZ)')
@@ -543,8 +540,9 @@ subroutine compute_arrays_source_noise(xi_noise,eta_noise,gamma_noise, &
   ! noise file (source time function)
   filename = trim(OUTPUT_FILES)//'/..//NOISE_TOMOGRAPHY/S_squared'
   open(unit=IIN_NOISE,file=trim(filename),status='old',action='read',iostat=ier)
-  if (ier /= 0 .and. myrank == 0)  &
+  if (ier /= 0 .and. myrank == 0) then
     call exit_MPI(myrank, 'file '//trim(filename)//' does NOT exist! This file should have been generated using Matlab scripts')
+  endif
 
   ! counts line reads noise source S(t)
   nlines = 0

src/specfem3D/read_mesh_databases.F90

@@ -178,8 +178,9 @@ subroutine read_mesh_databases()
 
   ! elastic simulation
   if (ELASTIC_SIMULATION) then
-
-    if (NB_RUNS_ACOUSTIC_GPU > 1) stop 'NB_RUNS_ACOUSTIC_GPU > 1 not compatible with elastic or coupled simulations'
+    ! checks
+    if (NB_RUNS_ACOUSTIC_GPU > 1) &
+      call exit_mpi(myrank,'NB_RUNS_ACOUSTIC_GPU > 1 not compatible with elastic or coupled simulations')
 
     ! displacement,velocity,acceleration
     allocate(displ(NDIM,NGLOB_AB),stat=ier)
@@ -346,7 +347,10 @@ subroutine read_mesh_databases()
   ! poroelastic
   if (POROELASTIC_SIMULATION) then
     ! checks
-    if (GPU_MODE) call exit_mpi(myrank,'POROELASTICITY not supported by GPU mode yet...')
+    if (GPU_MODE) &
+      call exit_mpi(myrank,'POROELASTICITY not supported by GPU mode yet...')
+    if (NB_RUNS_ACOUSTIC_GPU > 1) &
+      call exit_mpi(myrank,'NB_RUNS_ACOUSTIC_GPU > 1 not compatible with poroelastic or coupled simulations')
 
     ! displacement,velocity,acceleration for the solid (s) & fluid (w) phases
     allocate(displs_poroelastic(NDIM,NGLOB_AB),stat=ier)

src/specfem3D/setup_sources_receivers.f90

@@ -289,9 +289,10 @@ end subroutine setup_sources
 !
 !-------------------------------------------------------------------------------------------------
 !
-!EB EB When NB_RUNS_ACOUSTIC_GPU > 1, the file SOURCE_FILE actually contains the sources for all the runs.
-!This routine is intended to get the array that contains the run number of each source described in SOURCE_FILE.
-!The line i of the file run_number_of_the_source contains the run number \in [ 0;NB_RUNS_ACOUSTIC_GPU-1] of the source i
+! When NB_RUNS_ACOUSTIC_GPU > 1, the file SOURCE_FILE actually contains the sources for all the runs.
+! This routine is intended to get the array that contains the run number of each source described in SOURCE_FILE.
+! The line i of the file run_number_of_the_source contains the run number \in [ 0;NB_RUNS_ACOUSTIC_GPU-1] of the source i
+
   subroutine get_run_number_of_the_source()
 
   use constants
@@ -305,11 +306,11 @@ subroutine get_run_number_of_the_source()
   if (NB_RUNS_ACOUSTIC_GPU == 1) then
     run_number_of_the_source(:) = 0
   else
-
+    ! reads file DATA/run_number_of_the_source
     filename = IN_DATA_FILES(1:len_trim(IN_DATA_FILES))//'run_number_of_the_source'
-    open(unit=IIN,file=filename,status='old',action='read',iostat=ier)
+    open(unit=IIN,file=trim(filename),status='old',action='read',iostat=ier)
     if (ier /= 0) then
-      print *,'Error opening file: ',filename
+      print *,'Error opening file: ',trim(filename)
       stop 'Error opening run_number_of_the_source file'
     endif
 
@@ -323,7 +324,7 @@ subroutine get_run_number_of_the_source()
     if (icounter /= NSOURCES) stop 'Error total number of lines in run_number_of_the_source file is not equal to NSOURCES'
 
     ! Fills the array run_number_of_the_source
-    open(unit=IIN,file=filename,status='old',action='read')
+    open(unit=IIN,file=trim(filename),status='old',action='read')
     ! reads run number for each source
     do isource = 1,NSOURCES
       read(IIN,"(a)") string
@@ -1108,7 +1109,7 @@ subroutine setup_receivers_precompute_intp()
   implicit none
 
   ! local parameters
-  integer :: irec,irec_local,isource,ier
+  integer :: irec,irec_local,isource,ier,nrec_store
   integer,dimension(0:NPROC-1) :: tmp_rec_local_all
   integer :: maxrec,maxproc(1)
   double precision :: sizeval
@@ -1470,7 +1471,18 @@ subroutine setup_receivers_precompute_intp()
   if (ASDF_FORMAT) then
     if (.not. (SIMULATION_TYPE == 3 .and. (.not. SAVE_SEISMOGRAMS_IN_ADJOINT_RUN)) ) then
       ! initializes the ASDF data structure by allocating arrays
-      call init_asdf_data(nrec_local)
+      if (WRITE_SEISMOGRAMS_BY_MAIN) then
+        if (myrank == 0) then
+          ! main process holds all seismograms
+          nrec_store = nrec * NB_RUNS_ACOUSTIC_GPU
+        else
+          nrec_store = 0
+        endif
+      else
+        ! each process writes out its local receivers
+        nrec_store = nrec_local * NB_RUNS_ACOUSTIC_GPU
+      endif
+      call init_asdf_data(nrec_store)
       call synchronize_all()
     endif
   endif