Property writer

CMakeLists.txt

@@ -479,6 +479,7 @@ add_library(
         src/writer/seismogram.cpp
         src/writer/wavefield.cpp
         src/writer/kernel.cpp
+        src/writer/property.cpp
 )
 
 target_link_libraries(
@@ -530,6 +531,7 @@ add_library(
         src/parameter_parser/writer/wavefield.cpp
         src/parameter_parser/writer/plot_wavefield.cpp
         src/parameter_parser/writer/kernel.cpp
+        src/parameter_parser/writer/property.cpp
 )
 
 target_link_libraries(

include/parameter_parser/setup.hpp

@@ -11,6 +11,7 @@
 #include "specfem_setup.hpp"
 #include "time_scheme/interface.hpp"
 #include "writer/kernel.hpp"
+#include "writer/property.hpp"
 #include "writer/plot_wavefield.hpp"
 #include "writer/seismogram.hpp"
 #include "writer/wavefield.hpp"
@@ -169,6 +170,15 @@ class setup {
     }
   }
 
+  std::shared_ptr<specfem::writer::writer>
+  instantiate_property_writer(const specfem::compute::assembly &assembly) const {
+    if (this->property) {
+      return this->property->instantiate_property_writer(assembly);
+    } else {
+      return nullptr;
+    }
+  }
+
   std::shared_ptr<specfem::writer::writer>
   instantiate_kernel_writer(const specfem::compute::assembly &assembly) const {
     if (this->kernel) {
@@ -220,6 +230,7 @@ class setup {
       plot_wavefield; ///< Pointer to
                       ///< plot_wavefield object
   std::unique_ptr<specfem::runtime_configuration::kernel> kernel;
+  std::unique_ptr<specfem::runtime_configuration::property> property;
   std::unique_ptr<specfem::runtime_configuration::database_configuration>
       databases; ///< Get database filenames
   std::unique_ptr<specfem::runtime_configuration::solver::solver>

include/parameter_parser/writer/kernel.hpp

@@ -1,5 +1,4 @@
-#ifndef _SPECFEM_RUNTIME_CONFIGURATION_KERNEL_HPP
-#define _SPECFEM_RUNTIME_CONFIGURATION_KERNEL_HPP
+#pragma once
 
 #include "compute/assembly/assembly.hpp"
 #include "reader/reader.hpp"
@@ -31,5 +30,3 @@ class kernel {
 };
 } // namespace runtime_configuration
 } // namespace specfem
-
-#endif /* _SPECFEM_RUNTIME_CONFIGURATION_KERNEL_HPP */

include/parameter_parser/writer/property.hpp

@@ -0,0 +1,25 @@
+#pragma once
+
+#include "compute/assembly/assembly.hpp"
+#include "reader/reader.hpp"
+#include "writer/writer.hpp"
+#include "yaml-cpp/yaml.h"
+
+namespace specfem {
+namespace runtime_configuration {
+class property {
+public:
+  property(const std::string output_format, const std::string output_folder)
+      : output_format(output_format), output_folder(output_folder) {}
+
+  property(const YAML::Node &Node);
+
+  std::shared_ptr<specfem::writer::writer>
+  instantiate_property_writer(const specfem::compute::assembly &assembly) const;
+
+private:
+  std::string output_format;                 ///< format of output file
+  std::string output_folder;                 ///< Path to output folder
+};
+} // namespace runtime_configuration
+} // namespace specfem

include/writer/kernel.tpp

@@ -167,7 +167,7 @@ void specfem::writer::kernel<OutputLibrary>::write() {
           const specfem::point::index<specfem::dimension::type::dim2> index(
               ispec, iz, ix);
           specfem::point::kernels<specfem::dimension::type::dim2,
-                                  specfem::element::medium_tag::elastic,
+                                  specfem::element::medium_tag::acoustic,
                                   specfem::element::property_tag::isotropic,
                                   false>
               point_kernels;

include/writer/property.hpp

@@ -0,0 +1,43 @@
+#pragma once
+
+#include "compute/interface.hpp"
+#include "enumerations/interface.hpp"
+#include "writer/writer.hpp"
+
+namespace specfem {
+namespace writer {
+/**
+ * @brief Writer to model property data to disk
+ *
+ * @tparam OutputLibrary Library to use for output (HDF5, ASCII, etc.)
+ */
+template <typename OutputLibrary> class property : public writer {
+public:
+  /**
+   * @name Constructors
+   *
+   */
+  ///@{
+  /**
+   * @brief Construct a writer object
+   *
+   * @param assembly SPECFEM++ assembly
+   * @param output_folder Path to output location (will be an .h5 file if using
+   * HDF5, and a folder if using ASCII)
+   */
+  property(const specfem::compute::assembly &assembly,
+         const std::string output_folder);
+
+  /**
+   * @brief write the property data to disk
+   *
+   */
+  void write() override;
+
+private:
+  std::string output_folder;         ///< Path to output folder
+  specfem::compute::mesh mesh;       ///< Mesh object
+  specfem::compute::properties properties; ///< Properties object
+};
+} // namespace writer
+} // namespace specfem

include/writer/property.tpp

@@ -0,0 +1,181 @@
+#pragma once
+
+#include "compute/assembly/assembly.hpp"
+#include "enumerations/dimension.hpp"
+#include "enumerations/medium.hpp"
+#include "kokkos_abstractions.h"
+#include "point/properties.hpp"
+#include "writer/property.hpp"
+#include <Kokkos_Core.hpp>
+
+template <typename OutputLibrary>
+specfem::writer::property<OutputLibrary>::property(
+    const specfem::compute::assembly &assembly, const std::string output_folder)
+    : output_folder(output_folder), mesh(assembly.mesh),
+      properties(assembly.properties) {}
+
+template <typename OutputLibrary>
+void specfem::writer::property<OutputLibrary>::write() {
+
+  using DomainView =
+      Kokkos::View<type_real ***, Kokkos::LayoutLeft, Kokkos::HostSpace>;
+
+  properties.copy_to_host();
+
+  typename OutputLibrary::File file(output_folder + "/Properties");
+
+  const int nspec = mesh.points.nspec;
+  const int ngllz = mesh.points.ngllz;
+  const int ngllx = mesh.points.ngllx;
+
+  int n_elastic_isotropic;
+  int n_elastic_anisotropic;
+  int n_acoustic;
+
+  {
+    typename OutputLibrary::Group elastic = file.createGroup("/ElasticIsotropic");
+
+    const auto element_indices = properties.get_elements_on_host(
+        specfem::element::medium_tag::elastic,
+        specfem::element::property_tag::isotropic);
+    n_elastic_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 rho("rho", n_elastic_isotropic, ngllz, ngllx);
+    DomainView mu("mu", n_elastic_isotropic, ngllz, ngllx);
+    DomainView lambda("lambda", n_elastic_isotropic, ngllz, ngllx);
+
+    for (int i = 0; i < n_elastic_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);
+          const specfem::point::index<specfem::dimension::type::dim2> index(
+              ispec, iz, ix);
+          specfem::point::properties<specfem::dimension::type::dim2,
+                                  specfem::element::medium_tag::elastic,
+                                  specfem::element::property_tag::isotropic,
+                                  false>
+              point_properties;
+
+          specfem::compute::load_on_host(index, properties, point_properties);
+
+          rho(i, iz, ix) = point_properties.rho;
+          mu(i, iz, ix) = point_properties.mu;
+          lambda(i, iz, ix) = point_properties.lambda;
+        }
+      }
+    }
+
+    elastic.createDataset("X", x).write();
+    elastic.createDataset("Z", z).write();
+    elastic.createDataset("rho", rho).write();
+    elastic.createDataset("mu", mu).write();
+    elastic.createDataset("lambda", lambda).write();
+  }
+
+  {
+    typename OutputLibrary::Group elastic = file.createGroup("/ElasticAnisotropic");
+
+    const auto element_indices = properties.get_elements_on_host(
+        specfem::element::medium_tag::elastic,
+        specfem::element::property_tag::anisotropic);
+    n_elastic_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 rho("rho", n_elastic_anisotropic, ngllz, ngllx);
+    DomainView c11("c11", n_elastic_anisotropic, ngllz, ngllx);
+    DomainView c13("c13", n_elastic_anisotropic, ngllz, ngllx);
+    DomainView c15("c15", n_elastic_anisotropic, ngllz, ngllx);
+    DomainView c33("c33", n_elastic_anisotropic, ngllz, ngllx);
+    DomainView c35("c35", n_elastic_anisotropic, ngllz, ngllx);
+    DomainView c55("c55", n_elastic_anisotropic, ngllz, ngllx);
+
+    for (int i = 0; i < n_elastic_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);
+          const specfem::point::index<specfem::dimension::type::dim2> index(
+              ispec, iz, ix);
+          specfem::point::properties<specfem::dimension::type::dim2,
+                                  specfem::element::medium_tag::elastic,
+                                  specfem::element::property_tag::anisotropic,
+                                  false>
+              point_properties;
+
+          specfem::compute::load_on_host(index, properties, point_properties);
+
+          rho(i, iz, ix) = point_properties.rho;
+          c11(i, iz, ix) = point_properties.c11;
+          c13(i, iz, ix) = point_properties.c13;
+          c15(i, iz, ix) = point_properties.c15;
+          c33(i, iz, ix) = point_properties.c33;
+          c35(i, iz, ix) = point_properties.c35;
+          c55(i, iz, ix) = point_properties.c55;
+        }
+      }
+    }
+
+    elastic.createDataset("X", x).write();
+    elastic.createDataset("Z", z).write();
+    elastic.createDataset("rho", rho).write();
+    elastic.createDataset("c11", c11).write();
+    elastic.createDataset("c13", c13).write();
+    elastic.createDataset("c15", c15).write();
+    elastic.createDataset("c33", c33).write();
+    elastic.createDataset("c35", c35).write();
+    elastic.createDataset("c55", c55).write();
+  }
+
+  {
+    typename OutputLibrary::Group acoustic = file.createGroup("/Acoustic");
+
+    const auto element_indices = properties.get_elements_on_host(specfem::element::medium_tag::acoustic);
+    n_acoustic = element_indices.size();
+
+    DomainView x("xcoordinates_acoustic", n_acoustic, ngllz, ngllx);
+    DomainView z("zcoordinates_acoustic", n_acoustic, ngllz, ngllx);
+
+    DomainView rho("rho", n_acoustic, ngllz, ngllx);
+    DomainView kappa("kappa", n_acoustic, ngllz, ngllx);
+
+    for (int i = 0; i < n_acoustic; 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);
+          const specfem::point::index<specfem::dimension::type::dim2> index(
+              ispec, iz, ix);
+          specfem::point::properties<specfem::dimension::type::dim2,
+                                  specfem::element::medium_tag::acoustic,
+                                  specfem::element::property_tag::isotropic,
+                                  false>
+              point_properties;
+
+          specfem::compute::load_on_host(index, properties, point_properties);
+
+          rho(i, iz, ix) = 1.0 / point_properties.rho_inverse;
+          kappa(i, iz, ix) = point_properties.kappa;
+        }
+      }
+    }
+
+    acoustic.createDataset("X", x).write();
+    acoustic.createDataset("Z", z).write();
+    acoustic.createDataset("rho", rho).write();
+    acoustic.createDataset("kappa", kappa).write();
+  }
+
+  assert(n_elastic_isotropic + n_elastic_anisotropic + n_acoustic == nspec);
+
+  std::cout << "Properties written to " << output_folder << "/Properties"
+            << std::endl;
+}

src/parameter_parser/setup.cpp

@@ -120,6 +120,15 @@ specfem::runtime_configuration::setup::setup(const std::string &parameter_file,
           this->wavefield = nullptr;
         }
 
+        if (const YAML::Node &n_property = n_writer["model"]) {
+          at_least_one_writer = true;
+          this->property =
+              std::make_unique<specfem::runtime_configuration::property>(
+                  n_property);
+        } else {
+          this->property = nullptr;
+        }
+
         if (const YAML::Node &n_plotter = n_writer["display"]) {
           if ((n_plotter["simulation-field"] &&
                n_plotter["simulation-field"].as<std::string>() != "forward")) {

src/parameter_parser/writer/property.cpp

@@ -0,0 +1,59 @@
+#include "parameter_parser/writer/property.hpp"
+#include "IO/ASCII/ASCII.hpp"
+#include "IO/HDF5/HDF5.hpp"
+#include "writer/interface.hpp"
+#include "writer/property.hpp"
+#include <boost/filesystem.hpp>
+
+specfem::runtime_configuration::property::property(
+    const YAML::Node &Node) {
+
+  const std::string output_format = [&]() -> std::string {
+    if (Node["format"]) {
+      return Node["format"].as<std::string>();
+    } else {
+      return "ASCII";
+    }
+  }();
+
+  const std::string output_folder = [&]() -> std::string {
+    if (Node["directory"]) {
+      return Node["directory"].as<std::string>();
+    } else {
+      return boost::filesystem::current_path().string();
+    }
+  }();
+
+  if (!boost::filesystem::is_directory(
+          boost::filesystem::path(output_folder))) {
+    std::ostringstream message;
+    message << "Output folder : " << output_folder << " does not exist.";
+    throw std::runtime_error(message.str());
+  }
+
+  *this = specfem::runtime_configuration::property(output_format, output_folder);
+
+  return;
+}
+
+std::shared_ptr<specfem::writer::writer>
+specfem::runtime_configuration::property::instantiate_property_writer(
+    const specfem::compute::assembly &assembly) const {
+
+  const std::shared_ptr<specfem::writer::writer> writer =
+      [&]() -> std::shared_ptr<specfem::writer::writer> {
+      if (this->output_format == "HDF5") {
+        return std::make_shared<
+            specfem::writer::property<specfem::IO::HDF5<specfem::IO::write> > >(
+            assembly, this->output_folder);
+      } else if (this->output_format == "ASCII") {
+        return std::make_shared<
+            specfem::writer::property<specfem::IO::ASCII<specfem::IO::write> > >(
+            assembly, this->output_folder);
+      } else {
+        throw std::runtime_error("Unknown model format");
+      }
+  }();
+
+  return writer;
+}