Issue 1593 part 3 - Reads attenuation database

core/specfem/io/mesh/impl/fortran/dim2/mesh.cpp

@@ -193,7 +193,7 @@ specfem::mesh::mesh<specfem::dimension::type::dim2> specfem::io::read_2d_mesh(
        MEDIUM_TAG(ELASTIC_PSV, ELASTIC_SH, ACOUSTIC, POROELASTIC, ELASTIC_PSV_T,
                   ELECTROMAGNETIC_TE),
        PROPERTY_TAG(ISOTROPIC, ANISOTROPIC, ISOTROPIC_COSSERAT),
-       ATTENUATION_TAG(NONE)),
+       ATTENUATION_TAG(NONE, CONSTANT_ISOTROPIC)),
       {
         total_materials_read += mesh.materials
                                     .get_container<_medium_tag_, _property_tag_,

core/specfem/io/mesh/impl/fortran/dim3/read_materials.cpp

@@ -52,7 +52,16 @@ specfem::io::mesh::impl::fortran::dim3::read_materials(std::ifstream &stream,
                 "materials.");
           }
 
-          if ((std::abs(Qmu - 9999.0) < 1e-6) || (std::abs(Qmu) < 1e-6)) {
+          if (!((std::abs(Qmu - 9999.0) < 1e-6) || (std::abs(Qmu) < 1e-6))) {
+            std::ostringstream error_message;
+            error_message
+                << "Qmu should be set to 9999 or 0 for acoustic materials. "
+                << "Found Qmu = " << Qmu << " for material index " << imat
+                << "." << "[" << __FILE__ << ":" << __LINE__ << "]\n";
+            throw std::runtime_error(error_message.str());
+          }
+
+          if ((std::abs(Qkappa - 9999.0) < 1e-6) || (std::abs(Qkappa) < 1e-6)) {
 
             specfem::medium_container::material<
                 specfem::dimension::type::dim3,
@@ -66,8 +75,18 @@ specfem::io::mesh::impl::fortran::dim3::read_materials(std::ifstream &stream,
                                 specfem::element::attenuation_tag::none, index,
                                 imat });
           } else {
-            throw std::runtime_error(
-                "Attenuation not yet supported for acoustic materials in 3D");
+            specfem::medium_container::material<
+                specfem::dimension::type::dim3,
+                specfem::element::medium_tag::acoustic,
+                specfem::element::property_tag::isotropic,
+                specfem::element::attenuation_tag::constant_isotropic>
+                material(rho, vp, Qkappa, static_cast<type_real>(0.0));
+            const int index = materials.add_material(material);
+            mapping.push_back(
+                { specfem::element::medium_tag::acoustic,
+                  specfem::element::property_tag::isotropic,
+                  specfem::element::attenuation_tag::constant_isotropic, index,
+                  imat });
           }
         } else if (vs > 0.0) {
           // Isotropic elastic material
@@ -92,8 +111,18 @@ specfem::io::mesh::impl::fortran::dim3::read_materials(std::ifstream &stream,
                                 specfem::element::attenuation_tag::none, index,
                                 imat });
           } else {
-            throw std::runtime_error(
-                "Attenuation not yet supported for elastic materials in 3D");
+            specfem::medium_container::material<
+                specfem::dimension::type::dim3,
+                specfem::element::medium_tag::elastic,
+                specfem::element::property_tag::isotropic,
+                specfem::element::attenuation_tag::constant_isotropic>
+                material(rho, vs, vp, Qmu, Qkappa, static_cast<type_real>(0.0));
+            const int index = materials.add_material(material);
+            mapping.push_back(
+                { specfem::element::medium_tag::elastic,
+                  specfem::element::property_tag::isotropic,
+                  specfem::element::attenuation_tag::constant_isotropic, index,
+                  imat });
           }
 
         } else {

core/specfem/macros/material_iterators.hpp

@@ -198,14 +198,10 @@
        ATTENUATION_TAG_CONSTANT_ISOTROPIC))(                                   \
       (DIMENSION_TAG_DIM2, MEDIUM_TAG_ELASTIC_SH, PROPERTY_TAG_ANISOTROPIC,    \
        ATTENUATION_TAG_CONSTANT_ISOTROPIC))(                                   \
-      (DIMENSION_TAG_DIM2, MEDIUM_TAG_ELASTIC_PSV_T,                           \
-       PROPERTY_TAG_ISOTROPIC_COSSERAT, ATTENUATION_TAG_CONSTANT_ISOTROPIC))(  \
       (DIMENSION_TAG_DIM2, MEDIUM_TAG_ACOUSTIC, PROPERTY_TAG_ISOTROPIC,        \
        ATTENUATION_TAG_CONSTANT_ISOTROPIC))(                                   \
       (DIMENSION_TAG_DIM2, MEDIUM_TAG_POROELASTIC, PROPERTY_TAG_ISOTROPIC,     \
-       ATTENUATION_TAG_CONSTANT_ISOTROPIC))(                                   \
-      (DIMENSION_TAG_DIM2, MEDIUM_TAG_ELECTROMAGNETIC_TE,                      \
-       PROPERTY_TAG_ISOTROPIC, ATTENUATION_TAG_CONSTANT_ISOTROPIC))
+       ATTENUATION_TAG_CONSTANT_ISOTROPIC))
 
 #define MATERIAL_SYSTEMS_DIM3                                                  \
   ((DIMENSION_TAG_DIM3, MEDIUM_TAG_ELASTIC, PROPERTY_TAG_ISOTROPIC,            \

core/specfem/medium/dim2/acoustic/isotropic/material.hpp

@@ -19,6 +19,8 @@ struct AttenuationValues<
 public:
   type_real Qkappa; ///< Attenuation factor for bulk modulus
 
+  AttenuationValues() = default;
+
   AttenuationValues(const type_real &Qkappa) : Qkappa(Qkappa) {
     if (this->Qkappa <= 0.0) {
       throw std::runtime_error(

core/specfem/medium/dim2/elastic/isotropic/material.hpp

@@ -22,6 +22,8 @@ struct AttenuationValues<
   type_real Qkappa; ///< Attenuation factor for bulk modulus
   type_real Qmu;    ///< Attenuation factor for shear modulus
 
+  AttenuationValues() = default;
+
   AttenuationValues(const type_real &Qkappa, const type_real &Qmu)
       : Qkappa(Qkappa), Qmu(Qmu) {
     if (this->Qkappa <= 0.0 || this->Qmu <= 0.0) {

core/specfem/medium/dim2/poroelastic/isotropic/material.hpp

@@ -27,6 +27,8 @@ struct AttenuationValues<
     }
   };
 
+  AttenuationValues() = default;
+
   bool operator==(const AttenuationValues &other) const {
     return (std::abs(this->Qmu - other.Qmu) < 1e-6);
   }

core/specfem/mesh/dim2/materials/materials.cpp

@@ -0,0 +1,40 @@
+#include "materials.hpp"
+#include "enumerations/interface.hpp"
+#include "specfem/logger.hpp"
+#include "specfem/macros.hpp"
+
+void specfem::mesh::materials<specfem::dimension::type::dim2>::print() const {
+  std::ostringstream message;
+  message << "Total number of materials: " << this->n_materials << "\n";
+  FOR_EACH_IN_PRODUCT(
+      (DIMENSION_TAG(DIM2),
+       MEDIUM_TAG(ELASTIC_PSV, ELASTIC_SH, ACOUSTIC, POROELASTIC, ELASTIC_PSV_T,
+                  ELECTROMAGNETIC_TE),
+       PROPERTY_TAG(ISOTROPIC, ANISOTROPIC, ISOTROPIC_COSSERAT),
+       ATTENUATION_TAG(NONE, CONSTANT_ISOTROPIC)),
+      CAPTURE() {
+        const auto &material_container =
+            this->get_container<_medium_tag_, _property_tag_,
+                                _attenuation_tag_>();
+
+        if (material_container.n_materials > 0) {
+
+          message << "Material Type: \n"
+                  << "\t Medium Tag: "
+                  << specfem::element::to_string(_medium_tag_) << "\n"
+                  << "\tProperty Tag: "
+                  << specfem::element::to_string(_property_tag_) << "\n"
+                  << "\tAttenuation Tag: "
+                  << specfem::element::to_string(_attenuation_tag_) << "\n";
+
+          for (int i = 0; i < material_container.n_materials; ++i) {
+            message << "Material Index: " << i << "\n";
+            message << material_container.element_materials[i].print();
+          }
+        }
+      })
+
+  specfem::Logger::info(message.str());
+
+  return;
+}

core/specfem/mesh/dim2/materials/materials.hpp

@@ -52,7 +52,7 @@ template <> struct materials<specfem::dimension::type::dim2> {
             specfem::element::property_tag property,
             specfem::element::attenuation_tag attenuation>
   struct material {
-    int n_materials; ///< Number of elements
+    int n_materials = 0; ///< Number of different materials
     std::vector<
         specfem::medium_container::material<dimension_tag, type, property,
                                             attenuation> >
@@ -65,7 +65,7 @@ template <> struct materials<specfem::dimension::type::dim2> {
                  dimension_tag, type, property, attenuation> > &l_material);
   };
 
-  int n_materials; ///< Total number of different materials
+  int n_materials = 0; ///< Total number of different materials
   Kokkos::View<material_specification *, Kokkos::HostSpace>
       material_index_mapping; ///< Mapping of spectral element to material
                               ///< properties
@@ -75,7 +75,7 @@ template <> struct materials<specfem::dimension::type::dim2> {
        MEDIUM_TAG(ELASTIC_PSV, ELASTIC_SH, ACOUSTIC, POROELASTIC, ELASTIC_PSV_T,
                   ELECTROMAGNETIC_TE),
        PROPERTY_TAG(ISOTROPIC, ANISOTROPIC, ISOTROPIC_COSSERAT),
-       ATTENUATION_TAG(NONE)),
+       ATTENUATION_TAG(NONE, CONSTANT_ISOTROPIC)),
       DECLARE(((specfem::mesh::materials, (_DIMENSION_TAG_), ::material,
                 (_MEDIUM_TAG_, _PROPERTY_TAG_, _ATTENUATION_TAG_)),
                material)))
@@ -93,16 +93,27 @@ template <> struct materials<specfem::dimension::type::dim2> {
    * @brief Constructor used to allocate views
    *
    * @param nspec Number of spectral elements
-   * @param ngnod Number of control nodes per spectral element
    */
-  materials(const int nspec, const int numat)
-      : n_materials(numat),
-        material_index_mapping("specfem::mesh::material_index_mapping", nspec) {
-        };
+  materials(const int nspec)
+      : n_materials(0), material_index_mapping(
+                            "specfem::mesh::material_index_mapping", nspec) {};
 
   ///@}
 
 public:
+  template <specfem::element::medium_tag MediumTag,
+            specfem::element::property_tag PropertyTag,
+            specfem::element::attenuation_tag AttenuationTag>
+  int add_material(
+      const specfem::medium_container::material<
+          dimension_tag, MediumTag, PropertyTag, AttenuationTag> &material) {
+    this->n_materials += 1;
+    auto &material_container =
+        this->get_container<MediumTag, PropertyTag, AttenuationTag>();
+    material_container.element_materials.push_back(material);
+    material_container.n_materials += 1;
+    return material_container.n_materials - 1;
+  }
   /**
    * @brief Material material at spectral element index
    *
@@ -122,7 +133,7 @@ template <> struct materials<specfem::dimension::type::dim2> {
          MEDIUM_TAG(ELASTIC_PSV, ELASTIC_SH, ACOUSTIC, POROELASTIC,
                     ELASTIC_PSV_T, ELECTROMAGNETIC_TE),
          PROPERTY_TAG(ISOTROPIC, ANISOTROPIC, ISOTROPIC_COSSERAT),
-         ATTENUATION_TAG(NONE)),
+         ATTENUATION_TAG(NONE, CONSTANT_ISOTROPIC)),
         CAPTURE(material) {
           if constexpr (MediumTag == _medium_tag_ &&
                         PropertyTag == _property_tag_ &&
@@ -132,8 +143,6 @@ template <> struct materials<specfem::dimension::type::dim2> {
         })
 
     Kokkos::abort("Invalid material type detected in material specification");
-
-    return {};
   }
 
   /**
@@ -155,7 +164,31 @@ template <> struct materials<specfem::dimension::type::dim2> {
          MEDIUM_TAG(ELASTIC_PSV, ELASTIC_SH, ACOUSTIC, POROELASTIC,
                     ELASTIC_PSV_T, ELECTROMAGNETIC_TE),
          PROPERTY_TAG(ISOTROPIC, ANISOTROPIC, ISOTROPIC_COSSERAT),
-         ATTENUATION_TAG(NONE)),
+         ATTENUATION_TAG(NONE, CONSTANT_ISOTROPIC)),
+        CAPTURE(material) {
+          if constexpr (_medium_tag_ == MediumTag &&
+                        _property_tag_ == PropertyTag &&
+                        _attenuation_tag_ == AttenuationTag) {
+            return _material_;
+          }
+        })
+
+    Kokkos::abort("Invalid material type detected in material specification");
+  }
+
+  template <specfem::element::medium_tag MediumTag,
+            specfem::element::property_tag PropertyTag,
+            specfem::element::attenuation_tag AttenuationTag>
+  const specfem::mesh::materials<dimension_tag>::material<
+      MediumTag, PropertyTag, AttenuationTag> &
+  get_container() const {
+
+    FOR_EACH_IN_PRODUCT(
+        (DIMENSION_TAG(DIM2),
+         MEDIUM_TAG(ELASTIC_PSV, ELASTIC_SH, ACOUSTIC, POROELASTIC,
+                    ELASTIC_PSV_T, ELECTROMAGNETIC_TE),
+         PROPERTY_TAG(ISOTROPIC, ANISOTROPIC, ISOTROPIC_COSSERAT),
+         ATTENUATION_TAG(NONE, CONSTANT_ISOTROPIC)),
         CAPTURE(material) {
           if constexpr (_medium_tag_ == MediumTag &&
                         _property_tag_ == PropertyTag &&
@@ -199,6 +232,8 @@ template <> struct materials<specfem::dimension::type::dim2> {
     return std::make_tuple(material_specification.type,
                            material_specification.property);
   }
+
+  void print() const;
 };
 
 } // namespace mesh

core/specfem/mesh/dim3/materials/materials.hpp

@@ -120,7 +120,7 @@ template <> struct materials<specfem::dimension::type::dim3> {
             specfem::element::attenuation_tag attenuation>
   struct material {
     /** @brief Number of materials stored in this container */
-    int n_materials;
+    int n_materials = 0;
 
     /** @brief Storage for material objects of this type/property combination */
     std::vector<
@@ -147,7 +147,7 @@ template <> struct materials<specfem::dimension::type::dim3> {
   /** @{ */
 
   /** @brief Total number of materials across all type containers */
-  int n_materials;
+  int n_materials = 0;
 
   int nspec; ///< Total number of spectral elements in the mesh
 
@@ -191,7 +191,8 @@ template <> struct materials<specfem::dimension::type::dim3> {
    */
   FOR_EACH_IN_PRODUCT(
       (DIMENSION_TAG(DIM3), MEDIUM_TAG(ACOUSTIC, ELASTIC),
-       PROPERTY_TAG(ISOTROPIC, ANISOTROPIC), ATTENUATION_TAG(NONE)),
+       PROPERTY_TAG(ISOTROPIC, ANISOTROPIC),
+       ATTENUATION_TAG(NONE, CONSTANT_ISOTROPIC)),
       DECLARE(((specfem::mesh::materials, (_DIMENSION_TAG_), ::material,
                 (_MEDIUM_TAG_, _PROPERTY_TAG_, _ATTENUATION_TAG_)),
                material)))
@@ -265,7 +266,8 @@ template <> struct materials<specfem::dimension::type::dim3> {
 
     FOR_EACH_IN_PRODUCT(
         (DIMENSION_TAG(DIM3), MEDIUM_TAG(ACOUSTIC, ELASTIC),
-         PROPERTY_TAG(ISOTROPIC, ANISOTROPIC), ATTENUATION_TAG(NONE)),
+         PROPERTY_TAG(ISOTROPIC, ANISOTROPIC),
+         ATTENUATION_TAG(NONE, CONSTANT_ISOTROPIC)),
         CAPTURE(material) {
           if constexpr (MediumTag == _medium_tag_ &&
                         PropertyTag == _property_tag_ &&
@@ -324,7 +326,7 @@ template <> struct materials<specfem::dimension::type::dim3> {
 
     FOR_EACH_IN_PRODUCT((DIMENSION_TAG(DIM3), MEDIUM_TAG(ACOUSTIC, ELASTIC),
                          PROPERTY_TAG(ISOTROPIC, ANISOTROPIC),
-                         ATTENUATION_TAG(NONE)),
+                         ATTENUATION_TAG(NONE, CONSTANT_ISOTROPIC)),
                         CAPTURE(material) {
                           if constexpr (_medium_tag_ == MediumTag &&
                                         _property_tag_ == PropertyTag &&

tests/unit-tests/data/dim2/fluid_solid_mesh_flat_ocean_bottom_attenuation/STATIONS

@@ -0,0 +1,11 @@
+S0001    AA         1450.0000000         2400.0000000       0.0         0.0
+S0002    AA         1462.5000000         2500.0000000       0.0         0.0
+S0003    AA         1475.0000000         2600.0000000       0.0         0.0
+S0004    AA         1487.5000000         2700.0000000       0.0         0.0
+S0005    AA         1500.0000000         2800.0000000       0.0         0.0
+S0006    AA         1512.5000000         2900.0000000       0.0         0.0
+S0007    AA         1525.0000000         3000.0000000       0.0         0.0
+S0008    AA         1537.5000000         3100.0000000       0.0         0.0
+S0009    AA         1550.0000000         3200.0000000       0.0         0.0
+S0010    AA         1562.5000000         3300.0000000       0.0         0.0
+S0011    AA         1575.0000000         3400.0000000       0.0         0.0

tests/unit-tests/data/dim2/fluid_solid_mesh_flat_ocean_bottom_attenuation/Snakefile

@@ -0,0 +1,6 @@
+include: "../../meshfem2d.smk"
+
+rule all:
+    input:
+        "<cwd>/database.bin",
+    localrule: True