Added assembly properties and kernels

PrincetonUniversity · Feb 26, 2025 · 74fc87f · 74fc87f
1 parent fe6635e
commit 74fc87f
Show file tree

Hide file tree

Showing 42 changed files with 1,487 additions and 710 deletions.
diff --git a/include/IO/property/reader.hpp b/include/IO/property/reader.hpp
@@ -1,8 +1,8 @@
 #pragma once
 
+#include "IO/reader.hpp"
 #include "compute/interface.hpp"
 #include "enumerations/interface.hpp"
-#include "IO/reader.hpp"
 
 namespace specfem {
 namespace IO {

diff --git a/include/IO/property/reader.tpp b/include/IO/property/reader.tpp
@@ -1,18 +1,21 @@
 #pragma once
 
+#include "IO/property/reader.hpp"
 #include "compute/assembly/assembly.hpp"
 #include "enumerations/dimension.hpp"
 #include "enumerations/medium.hpp"
 #include "kokkos_abstractions.h"
 #include "point/properties.hpp"
-#include "IO/property/reader.hpp"
 #include <Kokkos_Core.hpp>
 
 template <typename InputLibrary>
-specfem::IO::property_reader<InputLibrary>::property_reader(const std::string input_folder): input_folder(input_folder) {}
+specfem::IO::property_reader<InputLibrary>::property_reader(
+    const std::string input_folder)
+    : input_folder(input_folder) {}
 
 template <typename InputLibrary>
-void specfem::IO::property_reader<InputLibrary>::read(specfem::compute::assembly &assembly) {
+void specfem::IO::property_reader<InputLibrary>::read(
+    specfem::compute::assembly &assembly) {
   auto &properties = assembly.properties;
 
   using DomainView =
@@ -21,32 +24,67 @@ void specfem::IO::property_reader<InputLibrary>::read(specfem::compute::assembly
   typename InputLibrary::File file(input_folder + "/Properties");
 
   {
-    typename InputLibrary::Group elastic = file.openGroup("/ElasticIsotropic");
+    typename InputLibrary::Group elastic =
+        file.openGroup("/ElasticSVIsotropic");
+
+    elastic.openDataset("rho", properties.elastic_sv_isotropic.h_rho).read();
+    elastic.openDataset("mu", properties.elastic_sv_isotropic.h_mu).read();
+    elastic
+        .openDataset("lambdaplus2mu",
+                     properties.elastic_sv_isotropic.h_lambdaplus2mu)
+        .read();
+  }
+
+  {
+    typename InputLibrary::Group elastic =
+        file.openGroup("/ElasticSHIsotropic");
 
-    elastic.openDataset("rho", properties.elastic_isotropic.h_rho).read();
-    elastic.openDataset("mu", properties.elastic_isotropic.h_mu).read();
-    elastic.openDataset("lambdaplus2mu", properties.elastic_isotropic.h_lambdaplus2mu).read();
+    elastic.openDataset("rho", properties.elastic_sh_isotropic.h_rho).read();
+    elastic.openDataset("mu", properties.elastic_sh_isotropic.h_mu).read();
+    elastic
+        .openDataset("lambdaplus2mu",
+                     properties.elastic_sh_isotropic.h_lambdaplus2mu)
+        .read();
   }
 
   {
-    typename InputLibrary::Group elastic = file.openGroup("/ElasticAnisotropic");
-
-    elastic.openDataset("rho", properties.elastic_anisotropic.h_rho).read();
-    elastic.openDataset("c11", properties.elastic_anisotropic.h_c11).read();
-    elastic.openDataset("c13", properties.elastic_anisotropic.h_c13).read();
-    elastic.openDataset("c15", properties.elastic_anisotropic.h_c15).read();
-    elastic.openDataset("c33", properties.elastic_anisotropic.h_c33).read();
-    elastic.openDataset("c35", properties.elastic_anisotropic.h_c35).read();
-    elastic.openDataset("c55", properties.elastic_anisotropic.h_c55).read();
-    elastic.openDataset("c12", properties.elastic_anisotropic.h_c12).read();
-    elastic.openDataset("c23", properties.elastic_anisotropic.h_c23).read();
-    elastic.openDataset("c25", properties.elastic_anisotropic.h_c25).read();
+    typename InputLibrary::Group elastic =
+        file.openGroup("/ElasticSVAnisotropic");
+
+    elastic.openDataset("rho", properties.elastic_sv_anisotropic.h_rho).read();
+    elastic.openDataset("c11", properties.elastic_sv_anisotropic.h_c11).read();
+    elastic.openDataset("c13", properties.elastic_sv_anisotropic.h_c13).read();
+    elastic.openDataset("c15", properties.elastic_sv_anisotropic.h_c15).read();
+    elastic.openDataset("c33", properties.elastic_sv_anisotropic.h_c33).read();
+    elastic.openDataset("c35", properties.elastic_sv_anisotropic.h_c35).read();
+    elastic.openDataset("c55", properties.elastic_sv_anisotropic.h_c55).read();
+    elastic.openDataset("c12", properties.elastic_sv_anisotropic.h_c12).read();
+    elastic.openDataset("c23", properties.elastic_sv_anisotropic.h_c23).read();
+    elastic.openDataset("c25", properties.elastic_sv_anisotropic.h_c25).read();
+  }
+
+  {
+    typename InputLibrary::Group elastic =
+        file.openGroup("/ElasticSHAnisotropic");
+
+    elastic.openDataset("rho", properties.elastic_sh_anisotropic.h_rho).read();
+    elastic.openDataset("c11", properties.elastic_sh_anisotropic.h_c11).read();
+    elastic.openDataset("c13", properties.elastic_sh_anisotropic.h_c13).read();
+    elastic.openDataset("c15", properties.elastic_sh_anisotropic.h_c15).read();
+    elastic.openDataset("c33", properties.elastic_sh_anisotropic.h_c33).read();
+    elastic.openDataset("c35", properties.elastic_sh_anisotropic.h_c35).read();
+    elastic.openDataset("c55", properties.elastic_sh_anisotropic.h_c55).read();
+    elastic.openDataset("c12", properties.elastic_sh_anisotropic.h_c12).read();
+    elastic.openDataset("c23", properties.elastic_sh_anisotropic.h_c23).read();
+    elastic.openDataset("c25", properties.elastic_sh_anisotropic.h_c25).read();
   }
 
   {
     typename InputLibrary::Group acoustic = file.openGroup("/Acoustic");
 
-    acoustic.openDataset("rho_inverse", properties.acoustic_isotropic.h_rho_inverse).read();
+    acoustic
+        .openDataset("rho_inverse", properties.acoustic_isotropic.h_rho_inverse)
+        .read();
     acoustic.openDataset("kappa", properties.acoustic_isotropic.h_kappa).read();
   }
 

diff --git a/include/IO/property/writer.hpp b/include/IO/property/writer.hpp
@@ -1,8 +1,8 @@
 #pragma once
 
+#include "IO/writer.hpp"
 #include "compute/interface.hpp"
 #include "enumerations/interface.hpp"
-#include "IO/writer.hpp"
 
 namespace specfem {
 namespace IO {
@@ -37,5 +37,5 @@ template <typename OutputLibrary> class property_writer : public writer {
 private:
   std::string output_folder; ///< Path to output folder
 };
-} // namespace writer
+} // namespace IO
 } // namespace specfem
diff --git a/include/IO/property/writer.tpp b/include/IO/property/writer.tpp
@@ -29,22 +29,24 @@ void specfem::IO::property_writer<OutputLibrary>::write(specfem::compute::assemb
   const int ngllz = mesh.points.ngllz;
   const int ngllx = mesh.points.ngllx;
 
-  int n_elastic_isotropic;
-  int n_elastic_anisotropic;
+  int n_elastic_sv_isotropic;
+  int n_elastic_sh_isotropic;
+  int n_elastic_sv_anisotropic;
+  int n_elastic_sh_anisotropic;
   int n_acoustic;
 
   {
-    typename OutputLibrary::Group elastic = file.createGroup("/ElasticIsotropic");
+    typename OutputLibrary::Group elastic = file.createGroup("/ElasticSVIsotropic");
 
     const auto element_indices = element_types.get_elements_on_host(
         specfem::element::medium_tag::elastic_sv,
         specfem::element::property_tag::isotropic);
-    n_elastic_isotropic = element_indices.size();
+        n_elastic_sv_isotropic = element_indices.size();
 
-    DomainView x("xcoordinates_elastic_isotropic", n_elastic_isotropic, ngllz, ngllx);
-    DomainView z("zcoordinates_elastic_isotropic", n_elastic_isotropic, ngllz, ngllx);
+    DomainView x("xcoordinates_elastic_isotropic", n_elastic_sv_isotropic, ngllz, ngllx);
+    DomainView z("zcoordinates_elastic_isotropic", n_elastic_sv_isotropic, ngllz, ngllx);
 
-    for (int i = 0; i < n_elastic_isotropic; i++) {
+    for (int i = 0; i < n_elastic_sv_isotropic; i++) {
       const int ispec = element_indices(i);
       for (int iz = 0; iz < ngllz; iz++) {
         for (int ix = 0; ix < ngllx; ix++) {
@@ -57,23 +59,88 @@ void specfem::IO::property_writer<OutputLibrary>::write(specfem::compute::assemb
     elastic.createDataset("X", x).write();
     elastic.createDataset("Z", z).write();
 
-    elastic.createDataset("rho", properties.elastic_isotropic.h_rho).write();
-    elastic.createDataset("mu", properties.elastic_isotropic.h_mu).write();
-    elastic.createDataset("lambdaplus2mu", properties.elastic_isotropic.h_lambdaplus2mu).write();
+    elastic.createDataset("rho", properties.elastic_sv_isotropic.h_rho).write();
+    elastic.createDataset("mu", properties.elastic_sv_isotropic.h_mu).write();
+    elastic.createDataset("lambdaplus2mu", properties.elastic_sv_isotropic.h_lambdaplus2mu).write();
   }
 
   {
-    typename OutputLibrary::Group elastic = file.createGroup("/ElasticAnisotropic");
+    typename OutputLibrary::Group elastic = file.createGroup("/ElasticSHIsotropic");
+
+    const auto element_indices = element_types.get_elements_on_host(
+        specfem::element::medium_tag::elastic_sh,
+        specfem::element::property_tag::isotropic);
+        n_elastic_sh_isotropic = element_indices.size();
+
+    DomainView x("xcoordinates_elastic_isotropic", n_elastic_sh_isotropic, ngllz, ngllx);
+    DomainView z("zcoordinates_elastic_isotropic", n_elastic_sh_isotropic, ngllz, ngllx);
+
+    for (int i = 0; i < n_elastic_sh_isotropic; i++) {
+      const int ispec = element_indices(i);
+      for (int iz = 0; iz < ngllz; iz++) {
+        for (int ix = 0; ix < ngllx; ix++) {
+          x(i, iz, ix) = mesh.points.h_coord(0, ispec, iz, ix);
+          z(i, iz, ix) = mesh.points.h_coord(1, ispec, iz, ix);
+        }
+      }
+    }
+
+    elastic.createDataset("X", x).write();
+    elastic.createDataset("Z", z).write();
+
+    elastic.createDataset("rho", properties.elastic_sh_isotropic.h_rho).write();
+    elastic.createDataset("mu", properties.elastic_sh_isotropic.h_mu).write();
+    elastic.createDataset("lambdaplus2mu", properties.elastic_sh_isotropic.h_lambdaplus2mu).write();
+  }
+
+  {
+    typename OutputLibrary::Group elastic = file.createGroup("/ElasticSVAnisotropic");
 
     const auto element_indices = element_types.get_elements_on_host(
         specfem::element::medium_tag::elastic_sv,
         specfem::element::property_tag::anisotropic);
-    n_elastic_anisotropic = element_indices.size();
+    n_elastic_sv_anisotropic = element_indices.size();
+
+    DomainView x("xcoordinates_elastic_anisotropic", n_elastic_sv_anisotropic, ngllz, ngllx);
+    DomainView z("zcoordinates_elastic_anisotropic", n_elastic_sv_anisotropic, ngllz, ngllx);
+
+    for (int i = 0; i < n_elastic_sv_anisotropic; i++) {
+      const int ispec = element_indices(i);
+      for (int iz = 0; iz < ngllz; iz++) {
+        for (int ix = 0; ix < ngllx; ix++) {
+          x(i, iz, ix) = mesh.points.h_coord(0, ispec, iz, ix);
+          z(i, iz, ix) = mesh.points.h_coord(1, ispec, iz, ix);
+        }
+      }
+    }
+
+    elastic.createDataset("X", x).write();
+    elastic.createDataset("Z", z).write();
+
+    elastic.createDataset("rho", properties.elastic_sv_anisotropic.h_rho).write();
+    elastic.createDataset("c11", properties.elastic_sv_anisotropic.h_c11).write();
+    elastic.createDataset("c13", properties.elastic_sv_anisotropic.h_c13).write();
+    elastic.createDataset("c15", properties.elastic_sv_anisotropic.h_c15).write();
+    elastic.createDataset("c33", properties.elastic_sv_anisotropic.h_c33).write();
+    elastic.createDataset("c35", properties.elastic_sv_anisotropic.h_c35).write();
+    elastic.createDataset("c55", properties.elastic_sv_anisotropic.h_c55).write();
+    elastic.createDataset("c12", properties.elastic_sv_anisotropic.h_c12).write();
+    elastic.createDataset("c23", properties.elastic_sv_anisotropic.h_c23).write();
+    elastic.createDataset("c25", properties.elastic_sv_anisotropic.h_c25).write();
+  }
+
+  {
+    typename OutputLibrary::Group elastic = file.createGroup("/ElasticSHAnisotropic");
+
+    const auto element_indices = element_types.get_elements_on_host(
+        specfem::element::medium_tag::elastic_sh,
+        specfem::element::property_tag::anisotropic);
+    n_elastic_sh_anisotropic = element_indices.size();
 
-    DomainView x("xcoordinates_elastic_anisotropic", n_elastic_anisotropic, ngllz, ngllx);
-    DomainView z("zcoordinates_elastic_anisotropic", n_elastic_anisotropic, ngllz, ngllx);
+    DomainView x("xcoordinates_elastic_anisotropic", n_elastic_sh_anisotropic, ngllz, ngllx);
+    DomainView z("zcoordinates_elastic_anisotropic", n_elastic_sh_anisotropic, ngllz, ngllx);
 
-    for (int i = 0; i < n_elastic_anisotropic; i++) {
+    for (int i = 0; i < n_elastic_sh_anisotropic; i++) {
       const int ispec = element_indices(i);
       for (int iz = 0; iz < ngllz; iz++) {
         for (int ix = 0; ix < ngllx; ix++) {
@@ -86,16 +153,16 @@ void specfem::IO::property_writer<OutputLibrary>::write(specfem::compute::assemb
     elastic.createDataset("X", x).write();
     elastic.createDataset("Z", z).write();
 
-    elastic.createDataset("rho", properties.elastic_anisotropic.h_rho).write();
-    elastic.createDataset("c11", properties.elastic_anisotropic.h_c11).write();
-    elastic.createDataset("c13", properties.elastic_anisotropic.h_c13).write();
-    elastic.createDataset("c15", properties.elastic_anisotropic.h_c15).write();
-    elastic.createDataset("c33", properties.elastic_anisotropic.h_c33).write();
-    elastic.createDataset("c35", properties.elastic_anisotropic.h_c35).write();
-    elastic.createDataset("c55", properties.elastic_anisotropic.h_c55).write();
-    elastic.createDataset("c12", properties.elastic_anisotropic.h_c12).write();
-    elastic.createDataset("c23", properties.elastic_anisotropic.h_c23).write();
-    elastic.createDataset("c25", properties.elastic_anisotropic.h_c25).write();
+    elastic.createDataset("rho", properties.elastic_sh_anisotropic.h_rho).write();
+    elastic.createDataset("c11", properties.elastic_sh_anisotropic.h_c11).write();
+    elastic.createDataset("c13", properties.elastic_sh_anisotropic.h_c13).write();
+    elastic.createDataset("c15", properties.elastic_sh_anisotropic.h_c15).write();
+    elastic.createDataset("c33", properties.elastic_sh_anisotropic.h_c33).write();
+    elastic.createDataset("c35", properties.elastic_sh_anisotropic.h_c35).write();
+    elastic.createDataset("c55", properties.elastic_sh_anisotropic.h_c55).write();
+    elastic.createDataset("c12", properties.elastic_sh_anisotropic.h_c12).write();
+    elastic.createDataset("c23", properties.elastic_sh_anisotropic.h_c23).write();
+    elastic.createDataset("c25", properties.elastic_sh_anisotropic.h_c25).write();
   }
 
   {
@@ -124,7 +191,7 @@ void specfem::IO::property_writer<OutputLibrary>::write(specfem::compute::assemb
     acoustic.createDataset("kappa", properties.acoustic_isotropic.h_kappa).write();
   }
 
-  assert(n_elastic_isotropic + n_elastic_anisotropic + n_acoustic == nspec);
+  assert(n_elastic_sv_isotropic + n_elastic_sv_anisotropic + n_elastic_sh_isotropic + n_elastic_sh_anisotropic + n_acoustic == nspec);
 
   std::cout << "Properties written to " << output_folder << "/Properties"
             << std::endl;

diff --git a/include/IO/seismogram/reader.hpp b/include/IO/seismogram/reader.hpp
@@ -1,8 +1,8 @@
 #ifndef SPECFEM_READER_SEISMOGRAM_HPP
 #define SPECFEM_READER_SEISMOGRAM_HPP
 
-#include "enumerations/specfem_enums.hpp"
 #include "IO/reader.hpp"
+#include "enumerations/specfem_enums.hpp"
 
 namespace specfem {
 namespace forcing_function {
@@ -17,13 +17,13 @@ class seismogram_reader {
 public:
   seismogram_reader(){};
   seismogram_reader(const char *filename,
-             const specfem::enums::seismogram::format type,
-             specfem::kokkos::HostView2d<type_real> source_time_function)
+                    const specfem::enums::seismogram::format type,
+                    specfem::kokkos::HostView2d<type_real> source_time_function)
       : filename(filename), type(type),
         source_time_function(source_time_function) {}
   seismogram_reader(const std::string &filename,
-             const specfem::enums::seismogram::format type,
-             specfem::kokkos::HostView2d<type_real> source_time_function)
+                    const specfem::enums::seismogram::format type,
+                    specfem::kokkos::HostView2d<type_real> source_time_function)
       : filename(filename), type(type),
         source_time_function(source_time_function) {}
   void read();

diff --git a/include/IO/wavefield/writer.hpp b/include/IO/wavefield/writer.hpp
@@ -1,8 +1,8 @@
 #pragma once
 
+#include "IO/writer.hpp"
 #include "compute/interface.hpp"
 #include "enumerations/interface.hpp"
-#include "IO/writer.hpp"
 
 namespace specfem {
 namespace IO {

diff --git a/include/IO/writer.hpp b/include/IO/writer.hpp
@@ -23,5 +23,5 @@ class writer {
   virtual void write(specfem::compute::assembly &assembly){};
 };
 
-} // namespace writer
+} // namespace IO
 } // namespace specfem
diff --git a/include/compute/element_types/element_types.hpp b/include/compute/element_types/element_types.hpp
@@ -112,6 +112,7 @@ struct element_types {
   CALL_MACRO_FOR_ALL_MEDIUM_TAGS(MEDIUM_TAG_VARIABLES,
                                  WHERE(DIMENSION_TAG_DIM2)
                                      WHERE(MEDIUM_TAG_ELASTIC_SV,
+                                           MEDIUM_TAG_ELASTIC_SH,
                                            MEDIUM_TAG_ACOUSTIC))
 
 #undef MEDIUM_TAG_VARIABLES
@@ -128,6 +129,7 @@ struct element_types {
   CALL_MACRO_FOR_ALL_MATERIAL_SYSTEMS(MATERIAL_SYSTEMS_VARIABLE_NAMES,
                                       WHERE(DIMENSION_TAG_DIM2)
                                           WHERE(MEDIUM_TAG_ELASTIC_SV,
+                                                MEDIUM_TAG_ELASTIC_SH,
                                                 MEDIUM_TAG_ACOUSTIC)
                                               WHERE(PROPERTY_TAG_ISOTROPIC,
                                                     PROPERTY_TAG_ANISOTROPIC))
@@ -146,7 +148,8 @@ struct element_types {
   CALL_MACRO_FOR_ALL_ELEMENT_TYPES(
       ELEMENT_TYPES_VARIABLE_NAMES,
       WHERE(DIMENSION_TAG_DIM2)
-          WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ACOUSTIC)
+          WHERE(MEDIUM_TAG_ELASTIC_SV, MEDIUM_TAG_ELASTIC_SH,
+                MEDIUM_TAG_ACOUSTIC)
               WHERE(PROPERTY_TAG_ISOTROPIC, PROPERTY_TAG_ANISOTROPIC)
                   WHERE(BOUNDARY_TAG_NONE, BOUNDARY_TAG_ACOUSTIC_FREE_SURFACE,
                         BOUNDARY_TAG_STACEY,