Wokring reading of the mass matrix and the materials

fortran/meshfem3d/generate_databases/save_arrays_solver.F90

@@ -165,7 +165,7 @@ subroutine save_arrays_solver_mesh()
   integer :: ier,i,itest
   character(len=MAX_STRING_LEN) :: filename
   ! Check actual size
-  integer :: reclen
+  ! integer :: reclen
 
   ! selects routine for file i/o format
   if (ADIOS_FOR_MESH) then
@@ -220,7 +220,7 @@ subroutine save_arrays_solver_mesh()
   write(IOUT) ystore_unique
   write(IOUT) zstore_unique
 
-  write(*,*) "xstore reclength", reclen
+  ! write(*,*) "xstore reclength", reclen
 
   write(IOUT) irregular_element_number
   write(IOUT) xix_regular
@@ -279,11 +279,14 @@ subroutine save_arrays_solver_mesh()
 
   ! this array is needed for acoustic simulations but also for elastic simulations with CPML,
   ! thus we allocate it and read it in all cases (whether the simulation is acoustic, elastic, or acoustic/elastic)
-
   call save_global_arrays(nspec, rhostore)
 
   ! write(IOUT) rhostore
 
+  ! write test value
+  itest = 9998
+  write(IOUT) itest
+
   ! elastic
   if (ELASTIC_SIMULATION) then
     write(IOUT) rmass
@@ -297,6 +300,10 @@ subroutine save_arrays_solver_mesh()
     ! write(IOUT) rho_vs
   endif
 
+  ! Write a test value
+  itest = 9997
+  write(IOUT) itest
+
   ! poroelastic
   if (POROELASTIC_SIMULATION) then
     write(IOUT) rmass_solid_poroelastic
@@ -323,6 +330,10 @@ subroutine save_arrays_solver_mesh()
     ! write(IOUT) rho_vsI
   endif
 
+  ! write test value
+  itest = 9996
+  write(IOUT) itest
+
   ! @Lucas & @Congyue need to uncomment this when implementing PML
 
   ! C-PML absorbing boundary conditions

include/mesh/dim3/element_types/element_types.hpp

@@ -36,7 +36,8 @@ template <> struct element_types<specfem::dimension::type::dim3> {
 
   void print(const int ispec) const;
 
-  template <specfem::element::medium_tag MediumTag> void print() const;
+  template <specfem::element::medium_tag MediumTag>
+  void print(const int i) const;
 
   View1D<bool> ispec_is_elastic;     ///< Elastic spectral elements
   View1D<bool> ispec_is_acoustic;    ///< Acoustic spectral elements

include/mesh/dim3/mapping/mapping.hpp

@@ -1,6 +1,7 @@
 #pragma once
 
 #include "enumerations/dimension.hpp"
+#include "enumerations/medium.hpp"
 #include "specfem_setup.hpp"
 #include <Kokkos_Core.hpp>
 
@@ -57,7 +58,10 @@ template <> struct mapping<specfem::dimension::type::dim3> {
 
   void print() const;
 
-  void print(int ispec) const;
+  void print(const int ispec) const;
+
+  template <specfem::element::medium_tag MediumTag>
+  void print(const int i) const;
 };
 
 } // namespace mesh

include/mesh/dim3/mass_matrix/mass_matrix.hpp

@@ -46,7 +46,6 @@ template <> struct mass_matrix<specfem::dimension::type::dim3> {
     }
   };
 
-private:
   View1D<type_real> elastic;
   View1D<type_real> acoustic;
   View1D<type_real> ocean_load;

src/IO/mesh/impl/fortran/dim3/mesh.cpp

@@ -271,7 +271,7 @@ specfem::IO::read_3d_mesh(const std::string mesh_parameters_file,
       specfem::mesh::element_types<specfem::dimension::type::dim3>(
           mesh.parameters.nspec);
 
-  // Read is_acoustic using read_array
+  // Read boolean is x array!
   try {
     specfem::IO::mesh::impl::fortran::dim3::read_array(
         stream, mesh.elements_types.ispec_is_acoustic);
@@ -289,6 +289,181 @@ specfem::IO::read_3d_mesh(const std::string mesh_parameters_file,
     throw std::runtime_error(error_message.str());
   }
 
+  // Print the element types
+  mesh.elements_types.print();
+  mesh.elements_types.print(0);
+
+  // Print elements first element of each category
+  try {
+    mesh.elements_types.print<specfem::element::medium_tag::acoustic>(0);
+  } catch (std::runtime_error &e) {
+    std::cout << e.what();
+  };
+  try {
+    mesh.elements_types.print<specfem::element::medium_tag::elastic>(0);
+  } catch (std::runtime_error &e) {
+    std::cout << e.what();
+  };
+  try {
+    mesh.elements_types.print<specfem::element::medium_tag::poroelastic>(0);
+  } catch (std::runtime_error &e) {
+    std::cout << e.what();
+  };
+
+  // read test value
+  try {
+    check_read_test_value(stream, 9999);
+  } catch (std::runtime_error &e) {
+    std::ostringstream error_message;
+    error_message << "Error reading test value from database file: " << e.what()
+                  << "(" << __FILE__ << ":" << __LINE__ << ")";
+    throw std::runtime_error(error_message.str());
+  }
+
+  // Intialize the mass matrix object
+  mesh.mass_matrix = specfem::mesh::mass_matrix<specfem::dimension::type::dim3>(
+      mesh.parameters.nglob, mesh.parameters.acoustic_simulation,
+      mesh.parameters.elastic_simulation,
+      mesh.parameters.poroelastic_simulation,
+      mesh.parameters.approximate_ocean_load);
+
+  // Read the acoustic mass matrix if acoustic simulation
+  if (mesh.parameters.acoustic_simulation) {
+    try {
+      specfem::IO::mesh::impl::fortran::dim3::read_array(
+          stream, mesh.mass_matrix.acoustic);
+    } catch (std::runtime_error &e) {
+      std::ostringstream error_message;
+      error_message << "Error reading acoustic mass matrix from database file: "
+                    << e.what() << "(" << __FILE__ << ":" << __LINE__ << ")";
+      throw std::runtime_error(error_message.str());
+    }
+  }
+
+  // Read the density rho
+  try {
+    specfem::IO::mesh::impl::fortran::dim3::read_array(stream,
+                                                       mesh.materials.rho);
+  } catch (std::runtime_error &e) {
+    std::ostringstream error_message;
+    error_message << "Error reading rho from database file: " << e.what() << "("
+                  << __FILE__ << ":" << __LINE__ << ")";
+    throw std::runtime_error(error_message.str());
+  }
+
+  // Read test value 9998
+  try {
+    check_read_test_value(stream, 9998);
+  } catch (std::runtime_error &e) {
+    std::ostringstream error_message;
+    error_message << "Error reading test value from database file: " << e.what()
+                  << "(" << __FILE__ << ":" << __LINE__ << ")";
+    throw std::runtime_error(error_message.str());
+  }
+
+  // Read the elastic mass matrix if elastic simulation
+  if (mesh.parameters.elastic_simulation) {
+    try {
+      specfem::IO::mesh::impl::fortran::dim3::read_array(
+          stream, mesh.mass_matrix.elastic);
+      if (mesh.parameters.approximate_ocean_load) {
+        specfem::IO::mesh::impl::fortran::dim3::read_array(
+            stream, mesh.mass_matrix.ocean_load);
+      }
+    } catch (std::runtime_error &e) {
+      std::ostringstream error_message;
+      error_message << "Error reading elastic mass matrix from database file: "
+                    << e.what() << "(" << __FILE__ << ":" << __LINE__ << ")";
+      throw std::runtime_error(error_message.str());
+    }
+
+    // Read the stacey boundary values
+    try {
+      specfem::IO::mesh::impl::fortran::dim3::read_array(stream,
+                                                         mesh.materials.rho_vp);
+    } catch (std::runtime_error &e) {
+      std::ostringstream error_message;
+      error_message << "Error reading rho_vp from database file: " << e.what()
+                    << "(" << __FILE__ << ":" << __LINE__ << ")";
+      throw std::runtime_error(error_message.str());
+    }
+
+    // Read the stacey boundary values
+    try {
+      specfem::IO::mesh::impl::fortran::dim3::read_array(stream,
+                                                         mesh.materials.rho_vs);
+    } catch (std::runtime_error &e) {
+      std::ostringstream error_message;
+      error_message << "Error reading rho_vs from database file: " << e.what()
+                    << "(" << __FILE__ << ":" << __LINE__ << ")";
+      throw std::runtime_error(error_message.str());
+    }
+  }
+
+  // Read test value 9997
+  try {
+    check_read_test_value(stream, 9997);
+  } catch (std::runtime_error &e) {
+    std::ostringstream error_message;
+    error_message << "Error reading test value from database file: " << e.what()
+                  << "(" << __FILE__ << ":" << __LINE__ << ")";
+    throw std::runtime_error(error_message.str());
+  }
+
+  // If simulation poroelastic
+  if (mesh.parameters.poroelastic_simulation) {
+    try {
+      specfem::IO::mesh::impl::fortran::dim3::read_array(
+          stream, mesh.mass_matrix.solid_poroelastic);
+      specfem::IO::mesh::impl::fortran::dim3::read_array(
+          stream, mesh.mass_matrix.fluid_poroelastic);
+    } catch (std::runtime_error &e) {
+      std::ostringstream error_message;
+      error_message
+          << "Error reading poroelastic mass matrix from database file: "
+          << e.what() << "(" << __FILE__ << ":" << __LINE__ << ")";
+      throw std::runtime_error(error_message.str());
+    }
+
+    // read the poroelastic material properties
+    try {
+      specfem::IO::mesh::impl::fortran::dim3::read_array(
+          stream, mesh.materials.poro_rho);
+      specfem::IO::mesh::impl::fortran::dim3::read_array(
+          stream, mesh.materials.poro_kappa);
+      specfem::IO::mesh::impl::fortran::dim3::read_array(
+          stream, mesh.materials.poro_eta);
+      specfem::IO::mesh::impl::fortran::dim3::read_array(
+          stream, mesh.materials.poro_tort);
+      specfem::IO::mesh::impl::fortran::dim3::read_array(
+          stream, mesh.materials.poro_perm);
+      specfem::IO::mesh::impl::fortran::dim3::read_array(
+          stream, mesh.materials.poro_phi);
+      specfem::IO::mesh::impl::fortran::dim3::read_array(
+          stream, mesh.materials.poro_rho_vpI);
+      specfem::IO::mesh::impl::fortran::dim3::read_array(
+          stream, mesh.materials.poro_rho_vpII);
+      specfem::IO::mesh::impl::fortran::dim3::read_array(
+          stream, mesh.materials.poro_rho_vsI);
+    } catch (std::runtime_error &e) {
+      std::ostringstream error_message;
+      error_message << "Error reading poroelastic material properties from "
+                       "database file: "
+                    << e.what() << "(" << __FILE__ << ":" << __LINE__ << ")";
+      throw std::runtime_error(error_message.str());
+    }
+  }
+
+  // Read test value 9996
+  try {
+    check_read_test_value(stream, 9996);
+  } catch (std::runtime_error &e) {
+    std::ostringstream error_message;
+    error_message << "Error reading test value from database file: " << e.what()
+                  << "(" << __FILE__ << ":" << __LINE__ << ")";
+    throw std::runtime_error(error_message.str());
+  }
+
   stream.close();
 
   return mesh;

src/mesh/dim3/element_types/element_types.cpp

@@ -10,14 +10,22 @@ void specfem::mesh::element_types<
   std::vector<int> ispec_acoustic;
   std::vector<int> ispec_poroelastic;
 
+  // Initialize the number of elements
+  this->nelastic = 0;
+  this->nacoustic = 0;
+  this->nporoelastic = 0;
+
   // Initialize the vectors
   for (int ispec = 0; ispec < nspec; ispec++) {
     if (ispec_is_elastic(ispec) == true) {
       ispec_elastic.push_back(ispec);
+      this->nelastic++;
     } else if (ispec_is_acoustic(ispec) == true) {
       ispec_acoustic.push_back(ispec);
+      this->nacoustic++;
     } else if (ispec_is_poroelastic(ispec) == true) {
       ispec_poroelastic.push_back(ispec);
+      this->nporoelastic++;
     } else {
       std::ostringstream message;
       message << "Error: Spectral element " << ispec
@@ -30,11 +38,6 @@ void specfem::mesh::element_types<
     this->ispec_acoustic = View1D<int>("ispec_acoustic", ispec_acoustic.size());
     this->ispec_poroelastic =
         View1D<int>("ispec_poroelastic", ispec_poroelastic.size());
-
-    // get numbers of elements
-    this->nelastic = ispec_elastic.size();
-    this->nacoustic = ispec_acoustic.size();
-    this->nporoelastic = ispec_poroelastic.size();
   };
 }
 
@@ -56,19 +59,56 @@ void specfem::mesh::element_types<specfem::dimension::type::dim3>::print(
 }
 
 template <specfem::element::medium_tag MediumTag>
-void specfem::mesh::element_types<specfem::dimension::type::dim3>::print()
-    const {
+void specfem::mesh::element_types<specfem::dimension::type::dim3>::print(
+    const int i) const {
   if (MediumTag == specfem::element::medium_tag::elastic) {
-    for (int i = 0; i < nelastic; i++) {
-      std::cout << "Elastic element " << ispec_elastic(i) << std::endl;
+    if (i >= nelastic) {
+      std::ostringstream message;
+      message << "Error: Elastic element " << i
+              << " does not exist.\nNumber "
+                 "of elastic elements is :"
+              << nelastic << ". "
+              << "(" << __FILE__ << ":" << __LINE__ << ")\n";
+      throw std::runtime_error(message.str());
     }
+    std::cout << "Elastic element " << ispec_elastic(i) << std::endl;
   } else if (MediumTag == specfem::element::medium_tag::acoustic) {
-    for (int i = 0; i < nacoustic; i++) {
-      std::cout << "Acoustic element " << ispec_acoustic(i) << std::endl;
+
+    if (i >= nacoustic) {
+      std::ostringstream message;
+      message << "Error: Acoustic element " << i
+              << " does not exist.\nNumber "
+                 " of acoustic elements: "
+              << nacoustic << ". "
+              << "(" << __FILE__ << ":" << __LINE__ << ")\n";
+      throw std::runtime_error(message.str());
     }
+
+    std::cout << "Acoustic element " << ispec_acoustic(i) << std::endl;
+
   } else if (MediumTag == specfem::element::medium_tag::poroelastic) {
-    for (int i = 0; i < nporoelastic; i++) {
-      std::cout << "Poroelastic element " << ispec_poroelastic(i) << std::endl;
+
+    if (i >= nporoelastic) {
+      std::ostringstream message;
+      message << "Error: Poroelastic element " << i
+              << " does not exist. Number "
+                 " of poroelastic elements: "
+              << nporoelastic << ". "
+              << "(" << __FILE__ << ":" << __LINE__ << ")\n";
+      throw std::runtime_error(message.str());
     }
+
+    std::cout << "Poroelastic element " << ispec_poroelastic(i) << std::endl;
   }
 }
+
+// Explicit instantiations
+template void
+specfem::mesh::element_types<specfem::dimension::type::dim3>::print<
+    specfem::element::medium_tag::elastic>(const int) const;
+template void
+specfem::mesh::element_types<specfem::dimension::type::dim3>::print<
+    specfem::element::medium_tag::acoustic>(const int) const;
+template void
+specfem::mesh::element_types<specfem::dimension::type::dim3>::print<
+    specfem::element::medium_tag::poroelastic>(const int) const;