Merge pull request #1648 from danielpeter/devel

fixes hdf5 i/o for ascii seismogram output; updates testing

.github/scripts/run_install.sh

@@ -15,6 +15,25 @@ if [[ $? -ne 0 ]]; then exit 1; fi
 # fortran/openMPI compiler
 sudo apt-get install -yq --no-install-recommends gfortran g++ openmpi-bin libopenmpi-dev
 
+# parallel hdf5
+if [[ "${TEST}" == *"with-hdf5"* ]]; then
+  echo
+  echo "additional installation: ${TEST}"
+  echo
+  sudo apt-get install -yq --no-install-recommends libhdf5-mpi-dev
+  ## checks installation paths
+  #echo
+  #dpkg -L libhdf5-mpi-dev
+  #echo
+  #dpkg -L libhdf5-openmpi-dev
+  #echo
+  #echo "hdf5 module paths:"
+  #find /usr/ -iname 'hdf5.mod'
+  #echo "hdf5 library paths:"
+  #find /usr/ -iname 'libhdf5hl_fortran*'
+  #echo
+fi
+
 # checks exit code
 if [[ $? -ne 0 ]]; then exit 1; fi
 echo

.github/scripts/run_tests.sh

@@ -76,7 +76,8 @@ if [ "$TESTDIR" == "EXAMPLES/applications/meshfem3D_examples/socal1D/" ]; then
     rm -f REF_SEIS; ln -s REF_SEIS.1d_cascadia REF_SEIS
   else
     # default
-    continue
+    # just continue
+    :
   fi
 fi
 # coupling FK
@@ -104,6 +105,17 @@ if [ "$TESTDIR" == "EXAMPLES/applications/meshfem3D_examples/sep_bathymetry/" ];
   sed -i "s:^NSTEP .*:NSTEP    = 1000:" DATA/Par_file
 fi
 
+# hdf5 i/o example
+if [[ "${TEST}" == *"with-hdf5"* ]]; then
+  echo
+  echo "test run: ${TEST}"
+  echo
+  sed -i "s:^HDF5_ENABLED .*:HDF5_ENABLED    = .true.:" DATA/Par_file
+  sed -i "s:^HDF5_FOR_MOVIES .*:HDF5_FOR_MOVIES    = .true.:" DATA/Par_file
+  sed -i "s:^HDF5_IO_NODES .*:HDF5_IO_NODES    = 1:" DATA/Par_file
+  # replaces run script
+  cp -v run_this_example_HDF5_IO_server.sh run_this_example.sh
+fi
 
 # default script
 ./run_this_example.sh

.github/workflows/CI.yml

@@ -682,4 +682,31 @@ jobs:
       run: ./.github/scripts/run_tests.sh
       shell: bash
 
+  linuxTest_18:
+    name: Test run example 18 - socal1D hdf5 i/o
+    runs-on: ubuntu-latest
+    needs: [linuxCheck]
+
+    steps:
+    - uses: actions/checkout@v3
+
+    - name: Install packages
+      env:
+        TEST: with-hdf5
+      run: ./.github/scripts/run_install.sh
+      shell: bash
+
+    - name: Run build
+      env:
+        TESTFLAGS: --with-mpi --with-hdf5 HDF5_INC=/usr/include/hdf5/openmpi/ HDF5_LIBS=-L/usr/lib/x86_64-linux-gnu/hdf5/openmpi
+      run: ./.github/scripts/run_build.sh
+      shell: bash
+
+    - name: Run test
+      env:
+        TEST: with-hdf5
+        TESTDIR: EXAMPLES/applications/meshfem3D_examples/socal1D/
+        TESTID: 0
+      run: ./.github/scripts/run_tests.sh
+      shell: bash
 

.travis.yml

@@ -141,9 +141,12 @@ env:
     # waterlayered_poroelastic example
     - TESTID=35 TESTCOV=0 TESTDIR=28 TESTFLAGS="--enable-vectorization --with-mpi" CUDA=false
 
-    # inversion example
+    # inversion example w/ coverage
     - TESTID=36 TESTCOV=1 TESTDIR=29 TESTFLAGS="--with-mpi" CUDA=false
 
+    # LTS example
+    - TESTID=37 TESTCOV=0 TESTDIR=30 TESTFLAGS="--enable-vectorization --with-mpi" CUDA=false
+
 
 # additional ARM tests
 jobs:

.travis/run_tests.sh

@@ -48,6 +48,7 @@ case "$TESTDIR" in
   27) dir=EXAMPLES/applications/Gmsh_simple_box_hex27/ ;;
   28) dir=EXAMPLES/applications/waterlayered_poroelastic/ ;;
   29) dir=EXAMPLES/applications/inversion_examples/fwi_test_acoustic/ ;;
+  30) dir=EXAMPLES/applications/LTS_homogeneous_halfspace_HEX8/ ;;
   *) dir=EXAMPLES/applications/homogeneous_halfspace/ ;;
 esac
 
@@ -223,6 +224,10 @@ else
     sed -i "s/Niter .*/Niter       : 1/" DATA/inverse_problem/inversion_fwi.dat
     sed -i "s/NSTEP .*/NSTEP         : 100/" DATA/inverse_problem/acquisition.dat
   fi
+  # LTS example
+  if [ "$TESTID" == "37" ]; then
+    sed -i "s:^NSTEP .*:NSTEP    = 1000:" DATA/Par_file
+  fi
 
   # coverage run
   if [ "$TESTCOV" == "1" ]; then
@@ -346,6 +351,23 @@ if [ "$TESTCOV" == "1" ] && [ "$TESTID" == "0" ]; then
 fi
 echo -en 'travis_fold:end:coverage.kernel\\r'
 
+## LTS example
+echo 'Coverage...' && echo -en 'travis_fold:start:coverage.LTS\\r'
+if [ "$TESTCOV" == "1" ] && [ "$TESTID" == "0" ]; then
+  ##
+  ## testing LTS simulation
+  ##
+  echo "##################################################################"
+  echo "EXAMPLES/applications/LTS_homogeneous_halfspace_HEX8/"
+  echo
+  cd EXAMPLES/applications/LTS_homogeneous_halfspace_HEX8/
+  sed -i "s:^NSTEP .*:NSTEP    = 50:" DATA/Par_file
+  ./run_this_example.sh
+  if [[ $? -ne 0 ]]; then exit 1; fi
+  cd $WORKDIR
+fi
+echo -en 'travis_fold:end:coverage.LTS\\r'
+
 #echo 'Coverage...' && echo -en 'travis_fold:start:coverage.acoustic\\r'
 #if [ "$TESTCOV" == "1" ] && [ "$TESTID" == "1" ]; then
 #  ##

doc/USER_MANUAL/copyright_and_version.tex

@@ -22,7 +22,7 @@ \chapter*{Copyright}
 
 version 4.1, November 2023:
 Rahul Garg, Hom Nath Gharti, Masaru Nagaso, Elif Oral, Daniel Peter.
-HDF5 file IO support; local time stepping; strain seismograms; movie output for stresses; 
+HDF5 file IO support; local time stepping; strain seismograms; movie output for stresses;
 reuse feature for external STF; various updates and improvements.\newline
 
 

flags.guess

@@ -108,7 +108,7 @@ case $my_FC in
 # with some versions of gfortran/gcc you make get errors about unused functions when compiling with the options below;
 # if so, either change -Wunused to -Wunused -Werror=no-unused-function, or remove -Wunused, or remove -Werror
 #
-        DEF_FFLAGS="-std=f2003 -fimplicit-none -fmax-errors=10 -pedantic -pedantic-errors -Waliasing -Wampersand -Wcharacter-truncation -Wline-truncation -Wsurprising -Wno-tabs -Wunderflow -ffpe-trap=invalid,zero,overflow -Wunused" # -mcmodel=medium
+        DEF_FFLAGS="-std=f2008 -fimplicit-none -fmax-errors=10 -pedantic -pedantic-errors -Waliasing -Wampersand -Wcharacter-truncation -Wline-truncation -Wsurprising -Wno-tabs -Wunderflow -ffpe-trap=invalid,zero,overflow -Wunused" # -mcmodel=medium
         OPT_FFLAGS="-O3 -finline-functions"
         DEBUG_FFLAGS="-g -O0 -ggdb -fbacktrace -fbounds-check -frange-check -Werror"
         # useful to track loss of accuracy because of automatic double to single precision conversion:  -Wconversion  (this may generate many warnings...)

src/shared/parallel.f90

@@ -87,7 +87,7 @@ subroutine init_mpi()
 
   ! initialize the MPI communicator and start the NPROCTOT MPI processes.
   call MPI_INIT(ier)
-  if (ier /= 0 ) stop 'Error initializing MPI'
+  if (ier /= 0) stop 'Error initializing MPI'
 
   ! initialize status size
   my_status_size   = MPI_STATUS_SIZE
@@ -212,7 +212,7 @@ subroutine synchronize_all()
 
   ! synchronizes MPI processes
   call MPI_BARRIER(my_local_mpi_comm_world,ier)
-  if (ier /= 0 ) stop 'Error synchronize MPI processes'
+  if (ier /= 0) stop 'Error synchronize MPI processes'
 
   end subroutine synchronize_all
 
@@ -234,7 +234,7 @@ subroutine synchronize_all_comm(comm)
 
   ! synchronizes MPI processes
   call MPI_BARRIER(comm,ier)
-  if (ier /= 0 ) stop 'Error synchronize MPI processes for specified communicator'
+  if (ier /= 0) stop 'Error synchronize MPI processes for specified communicator'
 
   end subroutine synchronize_all_comm
 
@@ -1946,11 +1946,11 @@ subroutine get_count_i(source,itag,recv_count)
 
   ! probe for incoming message from process source
   call MPI_Probe(source,itag,my_local_mpi_comm_world,msg_status,ier)
-  if ( ier /= 0 ) stop 'error mpi_probe for message status'
+  if ( ier /= 0) stop 'error mpi_probe for message status'
 
   ! gets the size of the message
   call MPI_Get_count(msg_status,MPI_INTEGER,recv_count,ier)
-  if ( ier /= 0 ) stop 'error mpi_get_count did not receive count'
+  if ( ier /= 0) stop 'error mpi_get_count did not receive count'
 
   end subroutine get_count_i
 
@@ -1973,7 +1973,7 @@ subroutine world_size(sizeval)
   integer :: ier
 
   call MPI_COMM_SIZE(my_local_mpi_comm_world,sizeval,ier)
-  if (ier /= 0 ) stop 'Error getting MPI world size'
+  if (ier /= 0) stop 'Error getting MPI world size'
 
   end subroutine world_size
 
@@ -1994,7 +1994,7 @@ subroutine world_size_comm(sizeval,comm)
   integer :: ier
 
   call MPI_COMM_SIZE(comm,sizeval,ier)
-  if (ier /= 0 ) stop 'Error getting MPI world size'
+  if (ier /= 0) stop 'Error getting MPI world size'
 
   end subroutine world_size_comm
 
@@ -2014,7 +2014,7 @@ subroutine world_rank(rank)
   integer :: ier
 
   call MPI_COMM_RANK(my_local_mpi_comm_world,rank,ier)
-  if (ier /= 0 ) stop 'Error getting MPI rank'
+  if (ier /= 0) stop 'Error getting MPI rank'
 
   end subroutine world_rank
 
@@ -2035,7 +2035,7 @@ subroutine world_rank_comm(rank,comm)
   integer :: ier
 
   call MPI_COMM_RANK(comm,rank,ier)
-  if (ier /= 0 ) stop 'Error getting MPI rank'
+  if (ier /= 0) stop 'Error getting MPI rank'
 
   end subroutine world_rank_comm
 
@@ -2123,7 +2123,7 @@ subroutine world_comm_free(comm)
   integer :: ier
 
   call MPI_Comm_free(comm,ier)
-  if (ier /= 0 ) stop 'Error freeing MPI communicator'
+  if (ier /= 0) stop 'Error freeing MPI communicator'
 
   end subroutine world_comm_free
 
@@ -2372,7 +2372,7 @@ subroutine synchronize_inter()
 
   ! synchronizes MPI processes
   call MPI_BARRIER(my_local_mpi_comm_inter,ier)
-  if (ier /= 0 ) stop 'Error synchronize MPI inter processes'
+  if (ier /= 0) stop 'Error synchronize MPI inter processes'
 
   end subroutine synchronize_inter
 
@@ -2390,7 +2390,7 @@ subroutine world_comm_free_inter()
   integer :: ier
 
   call MPI_Comm_free(my_local_mpi_comm_inter,ier)
-  if (ier /= 0 ) stop 'Error freeing MPI inter communicator'
+  if (ier /= 0) stop 'Error freeing MPI inter communicator'
 
   end subroutine world_comm_free_inter
 
@@ -2518,7 +2518,7 @@ subroutine irecvv_cr_inter(recvbuf, recvcount, dest, recvtag, req)
   integer :: ier
 
   call MPI_IRECV(recvbuf,recvcount,CUSTOM_MPI_TYPE,dest,recvtag, &
-                my_local_mpi_comm_inter,req,ier)
+                 my_local_mpi_comm_inter,req,ier)
 
   end subroutine irecvv_cr_inter
 

src/specfem3D/hdf5_io_server.F90

@@ -640,7 +640,8 @@ subroutine do_io_start_idle()
   use shared_parameters, only: NPROC,DT,NSTEP,NTSTEP_BETWEEN_FRAMES, &
     MOVIE_SURFACE,MOVIE_VOLUME,USE_HIGHRES_FOR_MOVIES,CREATE_SHAKEMAP, &
     MOVIE_VOLUME_STRESS, &
-    NSTEP,NTSTEP_BETWEEN_OUTPUT_SEISMOS
+    NSTEP,NTSTEP_BETWEEN_OUTPUT_SEISMOS, &
+    HDF5_FORMAT
 
   use specfem_par, only: nlength_seismogram
 
@@ -723,35 +724,39 @@ subroutine do_io_start_idle()
   !
 
   if (myrank == 0) then
-    ! get receiver info from compute nodes
-    call get_receiver_info(islice_num_rec_local)
+    ! for seismograms stored in HDF5 format
+    if (HDF5_FORMAT) then
+      ! get receiver info from compute nodes
+      call get_receiver_info(islice_num_rec_local)
 
-    ! initialize output file for seismo
-    call write_output_hdf5_seismogram_init()
+      ! initialize output file for seismo
+      call write_output_hdf5_seismogram_init()
+      if (VERBOSE) print *, "io_server: seismo init done"
 
-    ! count the number of procs having receivers (n_procs_with_rec)
-    ! and the number of receivers on each procs (islice...)
-    call count_nprocs_with_recs(islice_num_rec_local)
+      ! count the number of procs having receivers (n_procs_with_rec)
+      ! and the number of receivers on each procs (islice...)
+      call count_nprocs_with_recs(islice_num_rec_local)
 
-    ! check the seismo types to be saved
-    call count_seismo_type()
+      ! check the seismo types to be saved
+      call count_seismo_type()
 
-    ! allocate temporal arrays for seismo signals
-    call allocate_seismo_arrays(islice_num_rec_local)
+      ! allocate temporal arrays for seismo signals
+      call allocate_seismo_arrays(islice_num_rec_local)
 
-    ! initialize receive count
-    ! count the number of messages being sent
-    n_recv_msg_seismo = n_procs_with_rec * n_msg_seismo_each_proc * n_seismo_type
+      ! initialize receive count
+      ! count the number of messages being sent
+      n_recv_msg_seismo = n_procs_with_rec * n_msg_seismo_each_proc * n_seismo_type
 
-    if (NTSTEP_BETWEEN_OUTPUT_SEISMOS < NSTEP) then
-      max_seismo_out = int(NSTEP/NTSTEP_BETWEEN_OUTPUT_SEISMOS)
-      if (mod(NSTEP,NTSTEP_BETWEEN_OUTPUT_SEISMOS) /= 0) max_seismo_out = max_seismo_out+1
-    else
-      max_seismo_out = 1
-    endif
+      if (NTSTEP_BETWEEN_OUTPUT_SEISMOS < NSTEP) then
+        max_seismo_out = int(NSTEP/NTSTEP_BETWEEN_OUTPUT_SEISMOS)
+        if (mod(NSTEP,NTSTEP_BETWEEN_OUTPUT_SEISMOS) /= 0) max_seismo_out = max_seismo_out+1
+      else
+        max_seismo_out = 1
+      endif
 
-    ! receive the global id of received
-    call recv_id_rec(islice_num_rec_local)
+      ! receive the global id of received
+      call recv_id_rec(islice_num_rec_local)
+    endif
 
     !
     ! initialize surface movie

src/specfem3D/specfem3D.F90

@@ -271,7 +271,7 @@
 !
 ! version 4.1, November 2023:
 !   Rahul Garg, Hom Nath Gharti, Masaru Nagaso, Elif Oral, Daniel Peter.
-!     HDF5 file IO support; local time stepping; strain seismograms; movie output for stresses; 
+!     HDF5 file IO support; local time stepping; strain seismograms; movie output for stresses;
 !     reuse feature for external STF; various updates and improvements.
 !
 ! version 4.0, March 2023:

src/specfem3D/write_output_HDF5.F90

@@ -58,6 +58,13 @@ subroutine write_output_HDF5(all_seismograms,nrec_store,istore)
 
 ! HDF5 writes into a single file (by main process only)
 
+  ! checks if anything to do
+  if (.not. HDF5_FORMAT) return
+
+  ! safety check
+  if (.not. WRITE_SEISMOGRAMS_BY_MAIN) &
+    stop 'HDF5_FORMAT must have WRITE_SEISMOGRAMS_BY_MAIN set to .true.'
+
   ! io server
   if (HDF5_IO_NODES > 0) then
     ! only io nodes do the file output
@@ -88,10 +95,6 @@ subroutine write_output_HDF5(all_seismograms,nrec_store,istore)
     continue
   endif
 
-  ! safety check
-  if (.not. WRITE_SEISMOGRAMS_BY_MAIN) &
-    stop 'HDF5_FORMAT must have WRITE_SEISMOGRAMS_BY_MAIN set to .true.'
-
   ! initializes
   ! we only want to do this once
   if (.not. is_initialized) then
@@ -217,6 +220,9 @@ subroutine write_output_hdf5_seismogram_init()
   character(len=3),dimension(1) :: channels_press
   character(len=MAX_STRING_LEN) :: fname_h5_seismo
 
+  ! checks if anything to do
+  if (.not. HDF5_FORMAT) return
+
   ! only main process creates file
   if (myrank /= 0) return