Merge pull request #501 from PrincetonUniversity/issue-484

Issue 484/498/499/500 -- Implements reading of ispec_is_<medium>, materials, and mass matrices.

CMakeLists.txt

@@ -350,9 +350,11 @@ add_library(
         src/mesh/dim2/tags/tags.cpp
         # 3-D
         src/mesh/dim3/mesh.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
 )

docs/conf.py

@@ -78,7 +78,7 @@
 # }
 # html_logo = ''
 github_url = "https://github.com/PrincetonUniversity/SPECFEMPP"
-# html_baseurl = ''
+html_baseurl = 'https://specfem2d-kokkos.readthedocs.io/'
 
 # Add any paths that contain custom static files (such as style sheets) here,
 # relative to this directory. They are copied after the builtin static files,

docs/requirements.txt

@@ -1,3 +1,4 @@
+sphinx==8.0.0
 breathe
 furo
 sphinx-copybutton

fortran/meshfem3d/generate_databases/save_arrays_solver.F90

@@ -75,13 +75,41 @@ subroutine save_global_arrays_with_components(nspec, array)
 
   end subroutine save_global_arrays_with_components
 
+  ! subroutine that converts a logical array to an integer array and writes it
+  ! a file IOUT (see save_global_arrays)
+  subroutine save_logical_nspec_array(array)
+
+    use constants, only: IOUT
+    use generate_databases_par, only: nspec => NSPEC_AB
+
+    logical, dimension(nspec), intent(in) :: array
+
+    integer, dimension(nspec) :: array_int
+
+    integer :: i
+
+    do i = 1, nspec
+      if (array(i)) then
+        array_int(i) = 1
+      else
+        array_int(i) = 0
+      endif
+    enddo
+
+    write(IOUT) array_int
+
+  end subroutine save_logical_nspec_array
+
   end module save_arrays_module
 ! for external mesh
 
   subroutine save_arrays_solver_mesh()
 
   use constants, only: IMAIN,IOUT,myrank
-  use save_arrays_module, only: save_global_arrays, save_global_arrays_with_components
+  use save_arrays_module, only: &
+    save_global_arrays, &
+    save_global_arrays_with_components, &
+    save_logical_nspec_array
 
   use shared_parameters, only: ACOUSTIC_SIMULATION, ELASTIC_SIMULATION, POROELASTIC_SIMULATION, &
     APPROXIMATE_OCEAN_LOAD, SAVE_MESH_FILES, ANISOTROPY
@@ -136,6 +164,8 @@ subroutine save_arrays_solver_mesh()
   integer :: nglob
   integer :: ier,i,itest
   character(len=MAX_STRING_LEN) :: filename
+  ! Check actual size
+  ! integer :: reclen
 
   ! selects routine for file i/o format
   if (ADIOS_FOR_MESH) then
@@ -190,6 +220,8 @@ subroutine save_arrays_solver_mesh()
   write(IOUT) ystore_unique
   write(IOUT) zstore_unique
 
+  ! write(*,*) "xstore reclength", reclen
+
   write(IOUT) irregular_element_number
   write(IOUT) xix_regular
   write(IOUT) jacobian_regular
@@ -205,6 +237,10 @@ subroutine save_arrays_solver_mesh()
   call save_global_arrays(nspec, gammazstore)
   call save_global_arrays(nspec, jacobianstore)
 
+  ! write test value
+  itest = 10000
+  write(IOUT) itest
+
   ! write(IOUT) xixstore
   ! write(IOUT) xiystore
   ! write(IOUT) xizstore
@@ -221,10 +257,17 @@ subroutine save_arrays_solver_mesh()
   ! write(IOUT) kappastore
   ! write(IOUT) mustore
 
-  write(IOUT)
+
+  ! Converting the logical arrays to integer arrays and writing them as ints
+  ! call save_logical_nspec_array(ispec_is_acoustic)
+  ! call save_logical_nspec_array(ispec_is_elastic)
+  ! call save_logical_nspec_array(ispec_is_poroelastic)
+  ! inquire(iolength=reclen) xstore_unique
+  write(IOUT) ispec_is_acoustic
   write(IOUT) ispec_is_elastic
   write(IOUT) ispec_is_poroelastic
 
+
   ! stamp for checking i/o
   itest = 9999
   write(IOUT) itest
@@ -236,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
@@ -254,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
@@ -280,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/IO/fortranio/fortran_io.tpp

@@ -23,6 +23,13 @@ void fortran_read_value(std::vector<T> *value,
   return;
 }
 
+
+// Template specialization for std::vector<bool>
+template <>
+void fortran_read_value(std::vector<bool> *value,
+                        std::ifstream &stream,
+                        int &buffer_length);
+
 template <typename T, typename... Args>
 void fortran_IO(std::ifstream &stream, int &buffer_length, T *value,
                 Args... values) {
@@ -41,6 +48,7 @@ void fortran_read_line(std::ifstream &stream, Args... values) {
   }
 
   stream.read(reinterpret_cast<char *>(&buffer_length), fint);
+
   specfem::IO::fortran_IO(stream, buffer_length, values...);
 
   stream.read(reinterpret_cast<char *>(&buffer_length), fint);

include/IO/mesh/impl/fortran/dim2/read_material_properties.hpp

@@ -1,5 +1,6 @@
 #pragma once
 
+#include "enumerations/dimension.hpp"
 #include "mesh/mesh.hpp"
 #include "specfem_mpi/interface.hpp"
 #include <fstream>
@@ -25,9 +26,11 @@ namespace dim2 {
  * from the database file
  */
 
-specfem::mesh::materials read_material_properties(
-    std::ifstream &stream, const int numat, const int nspec,
-    const specfem::kokkos::HostView2d<int> knods, const specfem::MPI::MPI *mpi);
+specfem::mesh::materials<specfem::dimension::type::dim2>
+read_material_properties(std::ifstream &stream, const int numat,
+                         const int nspec,
+                         const specfem::kokkos::HostView2d<int> knods,
+                         const specfem::MPI::MPI *mpi);
 
 } // namespace dim2
 } // namespace fortran

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

@@ -57,11 +57,6 @@ void specfem::IO::mesh::impl::fortran::dim3::read_array(std::ifstream &stream,
           }
         }
       }
-      if ((ispec == 0) || (ispec == nspec - 1)) {
-        std::cout << "ispec=" << ispec
-                  << " array(ispec,0,0,0)=" << array(ispec, 0, 0, 0)
-                  << std::endl;
-      }
     }
 
   } catch (std::runtime_error &e) {
@@ -157,11 +152,6 @@ void specfem::IO::mesh::impl::fortran::dim3::read_index_array(
           }
         }
       }
-      if ((ispec == 0) || (ispec == nspec - 1)) {
-        std::cout << "ispec=" << ispec
-                  << " array(ispec,0,0,0)=" << array(ispec, 0, 0, 0)
-                  << std::endl;
-      }
     }
 
   } catch (std::runtime_error &e) {

include/compute/properties/properties.hpp

@@ -2,6 +2,7 @@
 
 #include "compute/element_types/element_types.hpp"
 #include "compute/impl/value_containers.hpp"
+#include "enumerations/dimension.hpp"
 #include "enumerations/specfem_enums.hpp"
 #include "kokkos_abstractions.h"
 #include "macros.hpp"
@@ -47,9 +48,11 @@ struct properties
    * @param has_gll_model Whether a GLL model is present (skip material property
    * assignment if true)
    */
-  properties(const int nspec, const int ngllz, const int ngllx,
-             const specfem::compute::element_types &element_types,
-             const specfem::mesh::materials &materials, bool has_gll_model);
+  properties(
+      const int nspec, const int ngllz, const int ngllx,
+      const specfem::compute::element_types &element_types,
+      const specfem::mesh::materials<specfem::dimension::type::dim2> &materials,
+      bool has_gll_model);
 
   ///@}
 

include/medium/material_properties.hpp

@@ -21,7 +21,8 @@ struct material_properties
   material_properties(
       const Kokkos::View<int *, Kokkos::DefaultHostExecutionSpace> elements,
       const int ngllz, const int ngllx,
-      const specfem::mesh::materials &materials, const bool has_gll_model,
+      const specfem::mesh::materials<specfem::dimension::type::dim2> &materials,
+      const bool has_gll_model,
       const specfem::kokkos::HostView1d<int> property_index_mapping)
       : specfem::medium::properties_container<type, property>(
             elements.extent(0), ngllz, ngllx) {

include/mesh/dim2/materials/materials.hpp

@@ -1,7 +1,9 @@
 #pragma once
 
+#include "enumerations/dimension.hpp"
 #include "kokkos_abstractions.h"
 #include "medium/material.hpp"
+#include "mesh/mesh_base.hpp"
 #include "specfem_mpi/interface.hpp"
 #include "specfem_setup.hpp"
 #include <variant>
@@ -12,7 +14,10 @@ namespace mesh {
  * @brief Material properties information
  *
  */
-struct materials {
+
+template <> struct materials<specfem::dimension::type::dim2> {
+  constexpr static auto dimension =
+      specfem::dimension::type::dim2; ///< Dimension type
 
   struct material_specification {
     specfem::element::medium_tag type;       ///< Type of element
@@ -59,17 +64,19 @@ struct materials {
       material_index_mapping; ///< Mapping of spectral element to material
                               ///< properties
 
-  specfem::mesh::materials::material<specfem::element::medium_tag::elastic,
-                                     specfem::element::property_tag::isotropic>
+  specfem::mesh::materials<specfem::dimension::type::dim2>::material<
+      specfem::element::medium_tag::elastic,
+      specfem::element::property_tag::isotropic>
       elastic_isotropic; ///< Elastic isotropic material properties
 
-  specfem::mesh::materials::material<
+  specfem::mesh::materials<specfem::dimension::type::dim2>::material<
       specfem::element::medium_tag::elastic,
       specfem::element::property_tag::anisotropic>
       elastic_anisotropic; ///< Elastic anisotropic material properties
 
-  specfem::mesh::materials::material<specfem::element::medium_tag::acoustic,
-                                     specfem::element::property_tag::isotropic>
+  specfem::mesh::materials<specfem::dimension::type::dim2>::material<
+      specfem::element::medium_tag::acoustic,
+      specfem::element::property_tag::isotropic>
       acoustic_isotropic; ///< Acoustic isotropic material properties
 
   /**

include/mesh/dim2/materials/materials.tpp

@@ -9,7 +9,7 @@
 
 template <specfem::element::medium_tag type,
           specfem::element::property_tag property>
-specfem::mesh::materials::material<type, property>::material(
+          specfem::mesh::materials<specfem::dimension::type::dim2>::material<type, property>::material(
     const int n_materials,
     const std::vector<specfem::medium::material<type, property> >
         &l_materials)

include/mesh/dim2/mesh.hpp

@@ -60,7 +60,8 @@ template <> struct mesh<specfem::dimension::type::dim2> {
                                                                   ///< nodes
                                                                   ///< (never
                                                                   ///< used)
-  specfem::mesh::materials materials; ///< Defines material properties
+  specfem::mesh::materials<dimension> materials; ///< Defines material
+                                                 ///< properties
 
   /**
    * @name Constructors
@@ -73,21 +74,22 @@ template <> struct mesh<specfem::dimension::type::dim2> {
    */
   mesh(){};
 
-  mesh(const int npgeo, const int nspec, const int nproc,
-       const specfem::mesh::control_nodes<specfem::dimension::type::dim2>
-           &control_nodes,
-       const specfem::mesh::parameters<specfem::dimension::type::dim2>
-           &parameters,
-       const specfem::mesh::coupled_interfaces<specfem::dimension::type::dim2>
-           &coupled_interfaces,
-       const specfem::mesh::boundaries<specfem::dimension::type::dim2>
-           &boundaries,
-       const specfem::mesh::tags<specfem::dimension::type::dim2> &tags,
-       const specfem::mesh::elements::tangential_elements<
-           specfem::dimension::type::dim2> &tangential_nodes,
-       const specfem::mesh::elements::axial_elements<
-           specfem::dimension::type::dim2> &axial_nodes,
-       const specfem::mesh::materials &materials)
+  mesh(
+      const int npgeo, const int nspec, const int nproc,
+      const specfem::mesh::control_nodes<specfem::dimension::type::dim2>
+          &control_nodes,
+      const specfem::mesh::parameters<specfem::dimension::type::dim2>
+          &parameters,
+      const specfem::mesh::coupled_interfaces<specfem::dimension::type::dim2>
+          &coupled_interfaces,
+      const specfem::mesh::boundaries<specfem::dimension::type::dim2>
+          &boundaries,
+      const specfem::mesh::tags<specfem::dimension::type::dim2> &tags,
+      const specfem::mesh::elements::tangential_elements<
+          specfem::dimension::type::dim2> &tangential_nodes,
+      const specfem::mesh::elements::axial_elements<
+          specfem::dimension::type::dim2> &axial_nodes,
+      const specfem::mesh::materials<specfem::dimension::type::dim2> &materials)
       : npgeo(npgeo), nspec(nspec), nproc(nproc), control_nodes(control_nodes),
         parameters(parameters), coupled_interfaces(coupled_interfaces),
         boundaries(boundaries), tags(tags), tangential_nodes(tangential_nodes),

include/mesh/dim2/tags/tags.hpp

@@ -41,9 +41,10 @@ template <> struct tags<specfem::dimension::type::dim2> {
    * @param materials Material properties
    * @param boundaries Boundary information
    */
-  tags(const specfem::mesh::materials &materials,
-       const specfem::mesh::boundaries<specfem::dimension::type::dim2>
-           &boundaries);
+  tags(
+      const specfem::mesh::materials<specfem::dimension::type::dim2> &materials,
+      const specfem::mesh::boundaries<specfem::dimension::type::dim2>
+          &boundaries);
   ///@}
 };
 } // namespace mesh