Refactor of the reading function using wrappers for the try catch blocks

CMakeLists.txt

@@ -245,6 +245,7 @@ add_library(
         src/IO/mesh/impl/fortran/dim3/read_parameters.cpp
         src/IO/mesh/impl/fortran/dim3/read_coordinates.cpp
         src/IO/mesh/impl/fortran/dim3/read_partial_derivatives.cpp
+        src/IO/mesh/impl/fortran/dim3/utilities.cpp
 )
 
 if (NOT HDF5_CXX_BUILD)
@@ -350,13 +351,16 @@ add_library(
         src/mesh/dim2/tags/tags.cpp
         # 3-D
         src/mesh/dim3/mesh.cpp
+        src/mesh/dim3/boundaries/absorbing_boundary.cpp
+        src/mesh/dim3/boundaries/free_surface.cpp
         src/mesh/dim3/element_types/element_types.cpp
         src/mesh/dim3/parameters/parameters.cpp
         src/mesh/dim3/parameters/parameters.cpp
         src/mesh/dim3/mapping/mapping.cpp
         src/mesh/dim3/materials/materials.cpp
         src/mesh/dim3/coordinates/coordinates.cpp
         src/mesh/dim3/partial_derivatives/partial_derivatives.cpp
+
 )
 
 target_link_libraries(

fortran/meshfem3d/generate_databases/save_arrays_solver.F90

@@ -437,7 +437,7 @@ subroutine save_arrays_solver_mesh()
   endif
 
   ! stamp for checking i/o so far
-  itest = 9998
+  itest = 9995
   write(IOUT) itest
 
   ! boundaries
@@ -448,6 +448,14 @@ subroutine save_arrays_solver_mesh()
   write(IOUT) NSPEC2D_BOTTOM
   write(IOUT) NSPEC2D_TOP
 
+  write(*,*) 'nspec2D_xmin: ', nspec2D_xmin
+  write(*,*) 'nspec2D_xmax: ', nspec2D_xmax
+  write(*,*) 'nspec2D_ymin: ', nspec2D_ymin
+  write(*,*) 'nspec2D_ymax: ', nspec2D_ymax
+  write(*,*) 'NSPEC2D_BOTTOM: ', NSPEC2D_BOTTOM
+  write(*,*) 'NSPEC2D_TOP: ', NSPEC2D_TOP
+
+
   if (nspec2D_xmin > 0) write(IOUT) ibelm_xmin
   if (nspec2D_xmax > 0) write(IOUT) ibelm_xmax
   if (nspec2D_ymin > 0) write(IOUT) ibelm_ymin

fortran/meshfem3d/generate_databases/save_parameters.f90

@@ -8,7 +8,7 @@ subroutine save_parameters()
     IOUT, myrank, sizeprocs, &
     NDIM, NSPEC_AB, NGLOB_AB, &
     NGLLX, NGLLY, NGLLZ, nfaces_surface, num_neighbors_all, NGLLSQUARE, &
-    nspec2D_xmin, nspec2D_xmin, nspec2D_ymin, nspec2D_ymin, &
+    nspec2D_xmin, nspec2D_xmax, nspec2D_ymin, nspec2D_ymax, &
     NSPEC2D_BOTTOM, NSPEC2D_TOP, &
     ACOUSTIC_SIMULATION, ELASTIC_SIMULATION, POROELASTIC_SIMULATION, &
     ANISOTROPY, &

include/IO/mesh/impl/fortran/dim3/interface.hpp

@@ -92,6 +92,10 @@ void read_index_array(std::ifstream &stream, View3D<T> &array);
 template <typename T>
 void read_index_array(std::ifstream &stream, View4D<T> &array);
 
+void check_read_test_value(std::ifstream &stream, int test_value);
+
+void check_values(std::string message, int value, int expected);
+
 } // namespace dim3
 } // namespace fortran
 } // namespace impl

include/IO/mesh/impl/fortran/dim3/utilities.hpp

@@ -1,5 +1,7 @@
+#pragma once
 #include "IO/fortranio/interface.hpp"
 #include "IO/mesh/impl/fortran/dim3/interface.hpp"
+#include "IO/mesh/impl/fortran/dim3/utilities.hpp"
 #include "enumerations/dimension.hpp"
 #include "mesh/mesh.hpp"
 #include "specfem_mpi/interface.hpp"
@@ -34,7 +36,7 @@ void specfem::IO::mesh::impl::fortran::dim3::read_array(std::ifstream &stream,
                   << e.what() << "(" << __FILE__ << ":" << __LINE__ << ")";
     throw std::runtime_error(error_message.str());
   }
-}
+};
 
 template <typename T>
 void specfem::IO::mesh::impl::fortran::dim3::read_array(std::ifstream &stream,
@@ -283,3 +285,77 @@ void specfem::IO::mesh::impl::fortran::dim3::read_index_array(
     throw std::runtime_error(error_message.str());
   }
 }
+
+namespace specfem {
+namespace IO {
+namespace mesh {
+namespace impl {
+namespace fortran {
+namespace dim3 {
+
+// Primary template for try-catch wrapper
+template <typename... Args>
+auto try_read_line(const std::string &message, Args &&...args)
+    -> decltype(specfem::IO::fortran_read_line(std::forward<Args>(args)...)) {
+  try {
+    return specfem::IO::fortran_read_line(std::forward<Args>(args)...);
+  } catch (const std::exception &e) {
+    std::ostringstream error_message;
+    error_message << "Exception caught in try_read_line(" << message
+                  << ", ...): " << e.what();
+    throw std::runtime_error(error_message.str());
+  } catch (...) {
+    std::ostringstream error_message;
+    error_message << "Unknown exception caught in " << message;
+    throw std::runtime_error(error_message.str());
+  }
+}
+
+// Primary template for try-catch wrapper
+template <typename... Args>
+auto try_read_array(const std::string &message, Args &&...args)
+    -> decltype(specfem::IO::mesh::impl::fortran::dim3::read_array(
+        std::forward<Args>(args)...)) {
+  try {
+    return specfem::IO::mesh::impl::fortran::dim3::read_array(
+        std::forward<Args>(args)...);
+  } catch (const std::exception &e) {
+    std::ostringstream error_message;
+    error_message << "Exception caught in try_read_array('" << message
+                  << "', ...): " << e.what();
+    throw std::runtime_error(error_message.str());
+  } catch (...) {
+    std::ostringstream error_message;
+    error_message << "Unknown exception caught in try_read_array('" << message
+                  << "', ...).";
+    throw std::runtime_error(error_message.str());
+  }
+}
+
+// Primary template for try-catch wrapper
+template <typename... Args>
+auto try_read_index_array(const std::string &message, Args &&...args)
+    -> decltype(specfem::IO::mesh::impl::fortran::dim3::read_index_array(
+        std::forward<Args>(args)...)) {
+  try {
+    return specfem::IO::mesh::impl::fortran::dim3::read_index_array(
+        std::forward<Args>(args)...);
+  } catch (const std::exception &e) {
+    std::ostringstream error_message;
+    error_message << "Exception caught in try_read_index('" << message
+                  << "', ...): " << e.what();
+    throw std::runtime_error(error_message.str());
+  } catch (...) {
+    std::ostringstream error_message;
+    error_message << "Unknown exception caught in try_read_index('" << message
+                  << "', ...).";
+    throw std::runtime_error(error_message.str());
+  }
+}
+
+} // namespace dim3
+} // namespace fortran
+} // namespace impl
+} // namespace mesh
+} // namespace IO
+} // namespace specfem

include/mesh/dim3/boundaries/absorbing_boundary.hpp

@@ -57,6 +57,19 @@ template <> struct absorbing_boundary<specfem::dimension::type::dim3> {
   bool acoustic = false;      ///< Flag for acoustic simulation
   bool elastic = false;       ///< Flag for elastic simulation
 
+  int nspec2D_xmin, nspec2D_xmax; /// Number of elements on the boundaries
+  int nspec2D_ymin, nspec2D_ymax;
+  int NSPEC2D_BOTTOM, NSPEC2D_TOP;
+
+  Kokkos::View<int *, Kokkos::HostSpace> ibelm_xmin,
+      ibelm_xmax; ///< Spectral
+                  ///< element
+                  ///< index for
+                  ///< elements on
+                  ///< the boundary
+  Kokkos::View<int *, Kokkos::HostSpace> ibelm_ymin, ibelm_ymax;
+  Kokkos::View<int *, Kokkos::HostSpace> ibelm_bottom, ibelm_top;
+
   Kokkos::View<int *, Kokkos::HostSpace> ispec; ///< Spectral element index for
                                                 ///< elements on the boundary
   Kokkos::View<int ***, Kokkos::HostSpace> ijk; ///< Which edge of the element
@@ -81,16 +94,24 @@ template <> struct absorbing_boundary<specfem::dimension::type::dim3> {
 
   absorbing_boundary(const int nglob, const int num_abs_boundary_faces,
                      const int ngllsquare, const bool acoustic,
-                     const bool elastic)
+                     const bool elastic, const int nspec2D_xmin,
+                     const int nspec2D_xmax, const int nspec2D_ymin,
+                     const int nspec2D_ymax, const int NSPEC2D_BOTTOM,
+                     const int NSPEC2D_TOP)
       : nelements(num_abs_boundary_faces), ngllsquare(ngllsquare),
         num_abs_boundary_faces(num_abs_boundary_faces), acoustic(acoustic),
-        elastic(elastic) {
+        elastic(elastic), nspec2D_xmin(nspec2D_xmin),
+        nspec2D_xmax(nspec2D_xmax), nspec2D_ymin(nspec2D_ymin),
+        nspec2D_ymax(nspec2D_ymax), NSPEC2D_BOTTOM(NSPEC2D_BOTTOM),
+        NSPEC2D_TOP(NSPEC2D_TOP) {
+
     ispec = Kokkos::View<int *, Kokkos::HostSpace>("ispec", nelements);
     ijk = Kokkos::View<int ***, Kokkos::HostSpace>("ijk", nelements, 3,
                                                    ngllsquare);
     jacobian2Dw = Kokkos::View<type_real **, Kokkos::HostSpace>(
         "jacobian2Dw", nelements, ngllsquare);
     // ndim=3
+
     normal = Kokkos::View<type_real ***, Kokkos::HostSpace>("normal", nelements,
                                                             3, ngllsquare);
 
@@ -101,8 +122,37 @@ template <> struct absorbing_boundary<specfem::dimension::type::dim3> {
       mass_acoustic =
           stacey_mass<specfem::element::medium_tag::acoustic>(nglob);
     }
+
+    // Create the boundary view elements only if the number of elements is
+    // greater than 0 on said boundary
+    if (nspec2D_xmin > 0) {
+      ibelm_xmin =
+          Kokkos::View<int *, Kokkos::HostSpace>("ibelm_xmin", nspec2D_xmin);
+    }
+    if (nspec2D_xmax > 0) {
+      ibelm_xmax =
+          Kokkos::View<int *, Kokkos::HostSpace>("ibelm_xmax", nspec2D_xmax);
+    }
+    if (nspec2D_ymin > 0) {
+      ibelm_ymin =
+          Kokkos::View<int *, Kokkos::HostSpace>("ibelm_ymin", nspec2D_ymin);
+    }
+    if (nspec2D_ymax > 0) {
+      ibelm_ymax =
+          Kokkos::View<int *, Kokkos::HostSpace>("ibelm_ymax", nspec2D_ymax);
+    }
+    if (NSPEC2D_BOTTOM > 0) {
+      ibelm_bottom = Kokkos::View<int *, Kokkos::HostSpace>("ibelm_bottom",
+                                                            NSPEC2D_BOTTOM);
+    }
+    if (NSPEC2D_TOP > 0) {
+      ibelm_top =
+          Kokkos::View<int *, Kokkos::HostSpace>("ibelm_top", NSPEC2D_TOP);
+    }
   }
   ///@}
+
+  void print() const;
 };
 
 } // namespace mesh

include/mesh/dim3/boundaries/free_surface.hpp

@@ -0,0 +1,61 @@
+#pragma once
+
+#include "enumerations/dimension.hpp"
+#include "enumerations/medium.hpp"
+#include "mesh/mesh_base.hpp"
+#include "specfem_setup.hpp"
+#include <Kokkos_Core.hpp>
+
+namespace specfem {
+namespace mesh {
+
+template <> struct free_surface<specfem::dimension::type::dim3> {
+
+  constexpr static auto dimension =
+      specfem::dimension::type::dim3; ///< Dimension type
+  // const int ndim = 3; ///< Dimension
+
+  int num_free_surface_faces; ///< Number of boundaries faces
+  int nelements;              ///< Number of elements on the boundary
+  int ngllsquare;             ///< Number of GLL points squared
+
+  Kokkos::View<int *, Kokkos::HostSpace> ispec; ///< Spectral element index for
+                                                ///< elements on the boundary
+  Kokkos::View<int ***, Kokkos::HostSpace> ijk; ///< Which edge of the element
+                                                ///< is on the boundary
+  Kokkos::View<type_real **, Kokkos::HostSpace> jacobian2Dw; ///< Jacobian of
+                                                             ///< the 2D
+  Kokkos::View<type_real ***, Kokkos::HostSpace> normal; ///< Jacobian of the 2D
+
+  /**
+   * @name Constructors
+   *
+   */
+  ///@{
+  /**
+   * @brief Default constructor
+   *
+   */
+  free_surface(){};
+
+  free_surface(const int num_free_surface_faces, const int ngllsquare)
+      : nelements(num_free_surface_faces), ngllsquare(ngllsquare),
+        num_free_surface_faces(num_free_surface_faces) {
+
+    ispec = Kokkos::View<int *, Kokkos::HostSpace>("ispec", nelements);
+    ijk = Kokkos::View<int ***, Kokkos::HostSpace>("ijk", nelements, 3,
+                                                   ngllsquare);
+    jacobian2Dw = Kokkos::View<type_real **, Kokkos::HostSpace>(
+        "jacobian2Dw", nelements, ngllsquare);
+    // ndim=3
+
+    normal = Kokkos::View<type_real ***, Kokkos::HostSpace>("normal", nelements,
+                                                            3, ngllsquare);
+  }
+  ///@}
+
+  void print() const;
+};
+
+} // namespace mesh
+} // namespace specfem

include/mesh/dim3/mesh.hpp

@@ -2,6 +2,7 @@
 
 #include "boundaries/absorbing_boundary.hpp"
 #include "boundaries/boundaries.hpp"
+#include "boundaries/free_surface.hpp"
 #include "coordinates/coordinates.hpp"
 #include "element_types/element_types.hpp"
 #include "mass_matrix/mass_matrix.hpp"
@@ -56,6 +57,9 @@ template <> struct mesh<specfem::dimension::type::dim3> {
   // Struct to store the absorbing boundaries
   specfem::mesh::absorbing_boundary<dimension> absorbing_boundary;
 
+  // Struct to store the free surface
+  specfem::mesh::free_surface<dimension> free_surface;
+
   /**
    * @name Constructors
    *

include/mesh/mesh_base.hpp

@@ -52,6 +52,12 @@ template <specfem::dimension::type DimensionType> struct absorbing_boundary;
  */
 template <specfem::dimension::type DimensionType> struct acoustic_free_surface;
 
+/**
+ * @brief Struct to store acoustic free surface boundaries
+ *
+ */
+template <specfem::dimension::type DimensionType> struct free_surface;
+
 /**
  * @brief Struct to store forcing boundaries
  *