Tests are passing with the fixed reading of the specfem2d mesh

Author: lsawade

CMakeLists.txt

@@ -178,6 +178,7 @@ else(ENABLE_DOUBLE_PRECISION)
     set(TYPE_REAL "float")
 endif(ENABLE_DOUBLE_PRECISION)
 
+# Configure the setup headers
 configure_file(${CMAKE_CURRENT_SOURCE_DIR}/setup/specfem_setup.hpp.in
                ${CMAKE_BINARY_DIR}/include/specfem_setup.hpp)
 configure_file(${CMAKE_CURRENT_SOURCE_DIR}/setup/constants.hpp.in

include/IO/fortranio/fortran_io.hpp

@@ -14,7 +14,9 @@ void fortran_IO(std::ifstream &stream, int &buffer_length);
 void fortran_read_value(bool *value, std::ifstream &stream, int &buffer_length);
 void fortran_read_value(std::string *value, std::ifstream &stream,
                         int &buffer_length);
-void fortran_read_value(type_real *value, std::ifstream &stream,
+void fortran_read_value(float *value, std::ifstream &stream,
+                        int &buffer_length);
+void fortran_read_value(double *value, std::ifstream &stream,
                         int &buffer_length);
 void fortran_read_value(int *value, std::ifstream &stream, int &buffer_length);
 } // namespace IO

src/IO/fortranio/fortran_io.cpp

@@ -41,7 +41,24 @@ void specfem::IO::fortran_read_value(int *value, std::ifstream &stream,
   return;
 }
 
-void specfem::IO::fortran_read_value(type_real *value, std::ifstream &stream,
+void specfem::IO::fortran_read_value(float *value, std::ifstream &stream,
+                                     int &buffer_length) {
+
+  float *temp;
+  buffer_length -= ffloat;
+  char *ivalue = new char[ffloat];
+  if (buffer_length < 0) {
+    std::cout << "buffer_length: " << buffer_length << std::endl;
+    throw std::runtime_error("Error reading fortran file");
+  }
+  stream.read(ivalue, ffloat);
+  temp = reinterpret_cast<float *>(ivalue);
+  *value = *temp;
+  delete[] ivalue;
+  return;
+}
+
+void specfem::IO::fortran_read_value(double *value, std::ifstream &stream,
                                      int &buffer_length) {
 
   double *temp;
@@ -53,7 +70,7 @@ void specfem::IO::fortran_read_value(type_real *value, std::ifstream &stream,
   }
   stream.read(ivalue, fdouble);
   temp = reinterpret_cast<double *>(ivalue);
-  *value = static_cast<type_real>(*temp);
+  *value = *temp;
   delete[] ivalue;
   return;
 }

src/IO/mesh/impl/fortran/dim2/read_elements.cpp

@@ -29,7 +29,7 @@ specfem::IO::mesh::impl::fortran::dim2::read_axial_elements(
 specfem::mesh::elements::tangential_elements<specfem::dimension::type::dim2>
 specfem::IO::mesh::impl::fortran::dim2::read_tangential_elements(
     std::ifstream &stream, const int nnodes_tangential_curve) {
-  type_real xread, yread;
+  double xread, yread;
 
   auto tangential_elements = specfem::mesh::elements::tangential_elements<
       specfem::dimension::type::dim2>(nnodes_tangential_curve);
@@ -41,8 +41,8 @@ specfem::IO::mesh::impl::fortran::dim2::read_tangential_elements(
   if (nnodes_tangential_curve > 0) {
     for (int inum = 0; inum < nnodes_tangential_curve; inum++) {
       specfem::IO::fortran_read_line(stream, &xread, &yread);
-      tangential_elements.x(inum) = xread;
-      tangential_elements.y(inum) = yread;
+      tangential_elements.x(inum) = static_cast<type_real>(xread);
+      tangential_elements.y(inum) = static_cast<type_real>(yread);
     }
   } else {
     tangential_elements.force_normal_to_surface = false;

src/IO/mesh/impl/fortran/dim2/read_material_properties.cpp

@@ -16,7 +16,7 @@ constexpr auto acoustic = specfem::element::medium_tag::acoustic;
 
 struct input_holder {
   // Struct to hold temporary variables read from database file
-  type_real val0, val1, val2, val3, val4, val5, val6, val7, val8, val9, val10,
+  double val0, val1, val2, val3, val4, val5, val6, val7, val8, val9, val10,
       val11, val12;
   int n, indic;
 };
@@ -90,11 +90,12 @@ std::vector<specfem::mesh::materials::material_specification> read_materials(
 
       // Acoustic Material
       if (read_values.val2 == 0) {
-        const type_real density = read_values.val0;
-        const type_real cp = read_values.val1;
-        const type_real compaction_grad = read_values.val3;
-        const type_real Qkappa = read_values.val5;
-        const type_real Qmu = read_values.val6;
+        const type_real density = static_cast<type_real>(read_values.val0);
+        const type_real cp = static_cast<type_real>(read_values.val1);
+        const type_real compaction_grad =
+            static_cast<type_real>(read_values.val3);
+        const type_real Qkappa = static_cast<type_real>(read_values.val5);
+        const type_real Qmu = static_cast<type_real>(read_values.val6);
 
         specfem::medium::material<acoustic, isotropic>
             acoustic_isotropic_holder(density, cp, Qkappa, Qmu,
@@ -113,12 +114,13 @@ std::vector<specfem::mesh::materials::material_specification> read_materials(
 
       } else {
 
-        const type_real density = read_values.val0;
-        const type_real cp = read_values.val1;
-        const type_real cs = read_values.val2;
-        const type_real compaction_grad = read_values.val3;
-        const type_real Qkappa = read_values.val5;
-        const type_real Qmu = read_values.val6;
+        const type_real density = static_cast<type_real>(read_values.val0);
+        const type_real cp = static_cast<type_real>(read_values.val1);
+        const type_real cs = static_cast<type_real>(read_values.val2);
+        const type_real compaction_grad =
+            static_cast<type_real>(read_values.val3);
+        const type_real Qkappa = static_cast<type_real>(read_values.val5);
+        const type_real Qmu = static_cast<type_real>(read_values.val6);
 
         specfem::medium::material<elastic, isotropic> elastic_isotropic_holder(
             density, cs, cp, Qkappa, Qmu, compaction_grad);
@@ -137,18 +139,18 @@ std::vector<specfem::mesh::materials::material_specification> read_materials(
     }
     // Ansotropic material
     else if (read_values.indic == 2) {
-      const type_real density = read_values.val0;
-      const type_real c11 = read_values.val1;
-      const type_real c13 = read_values.val2;
-      const type_real c15 = read_values.val3;
-      const type_real c33 = read_values.val4;
-      const type_real c35 = read_values.val5;
-      const type_real c55 = read_values.val6;
-      const type_real c12 = read_values.val7;
-      const type_real c23 = read_values.val8;
-      const type_real c25 = read_values.val9;
-      const type_real Qkappa = read_values.val11;
-      const type_real Qmu = read_values.val12;
+      const type_real density = static_cast<type_real>(read_values.val0);
+      const type_real c11 = static_cast<type_real>(read_values.val1);
+      const type_real c13 = static_cast<type_real>(read_values.val2);
+      const type_real c15 = static_cast<type_real>(read_values.val3);
+      const type_real c33 = static_cast<type_real>(read_values.val4);
+      const type_real c35 = static_cast<type_real>(read_values.val5);
+      const type_real c55 = static_cast<type_real>(read_values.val6);
+      const type_real c12 = static_cast<type_real>(read_values.val7);
+      const type_real c23 = static_cast<type_real>(read_values.val8);
+      const type_real c25 = static_cast<type_real>(read_values.val9);
+      const type_real Qkappa = static_cast<type_real>(read_values.val11);
+      const type_real Qmu = static_cast<type_real>(read_values.val12);
 
       specfem::medium::material<elastic, anisotropic>
           elastic_anisotropic_holder(density, c11, c13, c15, c33, c35, c55, c12,

src/IO/mesh/impl/fortran/dim2/read_mesh_database.cpp

@@ -15,7 +15,7 @@ specfem::IO::mesh::impl::fortran::dim2::read_mesh_database_header(
   // This subroutine reads header values of the database which are skipped
   std::string dummy_s;
   int dummy_i, dummy_i1, dummy_i2;
-  type_real dummy_d, dummy_d1;
+  double dummy_d, dummy_d1; // there is no type_real in the meshfem fortran code
   bool dummy_b, dummy_b1, dummy_b2, dummy_b3;
   int nspec, npgeo, nproc;
 
@@ -150,7 +150,7 @@ specfem::IO::mesh::impl::fortran::dim2::read_coorg_elements(
 
   int ipoin = 0;
 
-  type_real coorgi, coorgj;
+  double coorgi, coorgj;
   specfem::kokkos::HostView2d<type_real> coorg("specfem::mesh::coorg", ndim,
                                                npgeo);
 
@@ -161,8 +161,8 @@ specfem::IO::mesh::impl::fortran::dim2::read_coorg_elements(
     }
     // coorg stores the x,z for every control point
     // coorg([0, 2), i) = [x, z]
-    coorg(0, ipoin - 1) = coorgi;
-    coorg(1, ipoin - 1) = coorgj;
+    coorg(0, ipoin - 1) = static_cast<type_real>(coorgi);
+    coorg(1, ipoin - 1) = static_cast<type_real>(coorgj);
   }
 
   return coorg;
@@ -173,7 +173,7 @@ specfem::IO::mesh::impl::fortran::dim2::read_mesh_database_attenuation(
     std::ifstream &stream, const specfem::MPI::MPI *mpi) {
 
   int n_sls;
-  type_real attenuation_f0_reference;
+  double attenuation_f0_reference;
   bool read_velocities_at_f0;
   specfem::IO::fortran_read_line(stream, &n_sls, &attenuation_f0_reference,
                                  &read_velocities_at_f0);
@@ -184,6 +184,7 @@ specfem::IO::mesh::impl::fortran::dim2::read_mesh_database_attenuation(
   //                              "because it is used to assign some arrays");
   //   }
 
-  return std::make_tuple(n_sls, attenuation_f0_reference,
+  return std::make_tuple(n_sls,
+                         static_cast<type_real>(attenuation_f0_reference),
                          read_velocities_at_f0);
 }

tests/unit-tests/compute/partial_derivatives/compute_partial_derivatives_tests.cpp

@@ -84,13 +84,13 @@ TEST(COMPUTE_TESTS, compute_partial_derivatives) {
   const int ngllz = compute_mesh.quadratures.gll.N;
   const int ngllx = compute_mesh.quadratures.gll.N;
 
-  specfem::testing::array3d<type_real, Kokkos::LayoutRight> xix_ref(
+  specfem::testing::array3d<double, Kokkos::LayoutRight> xix_ref(
       test_config.xix_file, nspec, ngllz, ngllx);
-  specfem::testing::array3d<type_real, Kokkos::LayoutRight> gammax_ref(
+  specfem::testing::array3d<double, Kokkos::LayoutRight> gammax_ref(
       test_config.gammax_file, nspec, ngllz, ngllx);
-  specfem::testing::array3d<type_real, Kokkos::LayoutRight> gammaz_ref(
+  specfem::testing::array3d<double, Kokkos::LayoutRight> gammaz_ref(
       test_config.gammaz_file, nspec, ngllz, ngllx);
-  specfem::testing::array3d<type_real, Kokkos::LayoutRight> jacobian_ref(
+  specfem::testing::array3d<double, Kokkos::LayoutRight> jacobian_ref(
       test_config.jacobian_file, nspec, ngllz, ngllx);
 
   for (int ix = 0; ix < ngllx; ++ix) {

tests/unit-tests/domain/rmass_inverse_tests.cpp

@@ -166,7 +166,7 @@ TEST(DOMAIN_TESTS, rmass_inverse) {
                 << std::endl;
 
       if (Test.database.elastic_mass_matrix != "NULL") {
-        specfem::testing::array2d<type_real, Kokkos::LayoutRight>
+        specfem::testing::array2d<double, Kokkos::LayoutRight>
             h_mass_matrix_global(Test.database.elastic_mass_matrix, nglob, 2);
 
         specfem::testing::array3d<int, Kokkos::LayoutRight> index_mapping(
@@ -226,7 +226,7 @@ TEST(DOMAIN_TESTS, rmass_inverse) {
       }
 
       if (Test.database.acoustic_mass_matrix != "NULL") {
-        specfem::testing::array2d<type_real, Kokkos::LayoutRight>
+        specfem::testing::array2d<double, Kokkos::LayoutRight>
             h_mass_matrix_global(Test.database.acoustic_mass_matrix, nglob, 1);
 
         specfem::testing::array3d<int, Kokkos::LayoutRight> index_mapping(

tests/unit-tests/fortran_io/fortranio_tests.cpp

@@ -16,7 +16,7 @@ TEST(iotests, fortran_io) {
   int ival;
   bool bval;
   std::string sval;
-  type_real dval;
+  double dval;
   std::vector<int> vval(100, 0);
 
   stream.open(filename);