Implemented compute fields

- Added fields within the compute namespace
	- This follows the idea that all simulation data should be owned by compute namespace
	- The operating kernels can have a reference to this data

- Advantages of new implementation
	- Lower memory footprint through compressed field storage
	- Lower RangePolicy kernel time since kernel will only iterate though non-zero nglobs

CMakeLists.txt

@@ -249,6 +249,7 @@ add_library(
         src/compute/compute_receivers.cpp
         src/compute/coupled_interfaces.cpp
         src/compute/compute_boundaries.cpp
+        src/compute/compute_fields.cpp
 )
 
 target_link_libraries(

include/compute/fields/fields.hpp

@@ -0,0 +1,33 @@
+#ifndef _COMPUTE_FIELDS_FIELDS_HPP_
+#define _COMPUTE_FIELDS_FIELDS_HPP_
+
+#include "compute/compute_mesh.hpp"
+#include "compute/properties/interface.hpp"
+#include "enumerations/simulation.hpp"
+#include "enumerations/specfem_enums.hpp"
+#include "simulation_field.hpp"
+#include "simulation_field.tpp"
+
+namespace specfem {
+namespace compute {
+struct fields {
+
+  using forward_type = specfem::enums::simulation::forward;
+
+  fields(const specfem::compute::mesh &mesh,
+         const specfem::compute::properties &properties);
+
+  template <typename simulation>
+  KOKKOS_INLINE_FUNCTION specfem::compute::simulation_field<simulation>
+  get_simulation_field() {
+    if constexpr (std::is_same_v<simulation, forward_type>) {
+      return forward;
+    }
+  }
+
+  specfem::compute::simulation_field<forward_type> forward;
+};
+} // namespace compute
+} // namespace specfem
+
+#endif /* _COMPUTE_FIELDS_FIELDS_HPP_ */

include/compute/fields/impl/field_impl.hpp

@@ -0,0 +1,49 @@
+#ifndef _COMPUTE_FIELDS_IMPL_FIELD_IMPL_HPP_
+#define _COMPUTE_FIELDS_IMPL_FIELD_IMPL_HPP_
+
+#include "compute/compute_mesh.hpp"
+#include "compute/properties/interface.hpp"
+#include "kokkos_abstractions.h"
+#include <Kokkos_Core.hpp>
+
+namespace specfem {
+namespace compute {
+namespace impl {
+template <typename medium> class field_impl {
+public:
+  using medium_type = medium;
+
+  constexpr static int components = medium::components;
+
+  field_impl() = default;
+
+  field_impl(const specfem::compute::mesh &mesh,
+             const specfem::compute::properties &properties,
+             Kokkos::View<int * [specfem::enums::element::ntypes],
+                          Kokkos::LayoutLeft, specfem::kokkos::HostMemSpace>
+                 assembly_index_mapping);
+
+  field_impl(const int nglob, const int nspec, const int ngllz,
+             const int ngllx);
+
+  int nglob;
+  int nspec;
+  int ngllz;
+  int ngllx;
+  specfem::kokkos::DeviceView3d<type_real> index_mapping;
+  specfem::kokkos::HostMirror3d<type_real> h_index_mapping;
+  specfem::kokkos::DeviceView2d<type_real, Kokkos::LayoutLeft> field;
+  specfem::kokkos::HostMirror2d<type_real, Kokkos::LayoutLeft> h_field;
+  specfem::kokkos::DeviceView2d<type_real, Kokkos::LayoutLeft> field_dot;
+  specfem::kokkos::HostMirror2d<type_real, Kokkos::LayoutLeft> h_field_dot;
+  specfem::kokkos::DeviceView2d<type_real, Kokkos::LayoutLeft> field_dot_dot;
+  specfem::kokkos::HostMirror2d<type_real, Kokkos::LayoutLeft> h_field_dot_dot;
+  specfem::kokkos::DeviceView2d<type_real, Kokkos::LayoutLeft> mass_inverse;
+  specfem::kokkos::HostMirror2d<type_real, Kokkos::LayoutLeft> h_mass_inverse;
+};
+} // namespace impl
+
+} // namespace compute
+} // namespace specfem
+
+#endif /* _COMPUTE_FIELDS_IMPL_FIELD_IMPL_HPP_ */

include/compute/fields/impl/field_impl.tpp

@@ -0,0 +1,225 @@
+#ifndef _COMPUTE_FIELDS_IMPL_FIELD_IMPL_TPP_
+#define _COMPUTE_FIELDS_IMPL_FIELD_IMPL_TPP_
+
+#include "compute/fields/impl/field_impl.hpp"
+#include "kokkos_abstractions.h"
+#include <Kokkos_Core.hpp>
+
+namespace {
+int compute_nglob(const specfem::kokkos::HostView3d<int> index_mapping) {
+  const int nspec = index_mapping.extent(0);
+  const int ngllz = index_mapping.extent(1);
+  const int ngllx = index_mapping.extent(2);
+
+  int nglob;
+  Kokkos::parallel_reduce(
+      "specfem::utils::compute_nglob",
+      specfem::kokkos::HostMDrange<3>({ 0, 0, 0 }, { nspec, ngllz, ngllx }),
+      KOKKOS_LAMBDA(const int ispec, const int iz, const int ix, int &l_nglob) {
+        l_nglob = l_nglob > index_mapping(ispec, iz, ix)
+                      ? l_nglob
+                      : index_mapping(ispec, iz, ix);
+      },
+      Kokkos::Max<int>(nglob));
+
+  return nglob + 1;
+}
+} // namespace
+
+template <typename medium>
+specfem::compute::impl::field_impl<medium>::field_impl(const int nglob,
+                                                 const int nspec,
+                                                 const int ngllz,
+                                                 const int ngllx)
+    : nglob(nglob), nspec(nspec),
+      index_mapping("specfem::compute::fields::index_mapping", nspec, ngllz,
+                    ngllx),
+      h_index_mapping(Kokkos::create_mirror_view(index_mapping)),
+      field("specfem::compute::fields::field", nglob, medium::components),
+      h_field(Kokkos::create_mirror_view(field)),
+      field_dot("specfem::compute::fields::field_dot", nglob,
+                medium::components),
+      h_field_dot(Kokkos::create_mirror_view(field_dot)),
+      field_dot_dot("specfem::compute::fields::field_dot_dot", nglob,
+                    medium::components),
+      h_field_dot_dot(Kokkos::create_mirror_view(field_dot_dot)),
+      mass_inverse("specfem::compute::fields::mass_inverse", nglob, nglob),
+      h_mass_inverse(Kokkos::create_mirror_view(mass_inverse)) {}
+
+template <typename medium>
+specfem::compute::impl::field_impl<medium>::field_impl(
+    const specfem::compute::mesh &mesh,
+    const specfem::compute::properties &properties,
+    Kokkos::View<int * [specfem::enums::element::ntypes], Kokkos::LayoutLeft,
+                 specfem::kokkos::HostMemSpace>
+        assembly_index_mapping) {
+
+  const auto index_mapping = mesh.points.index_mapping;
+  const auto element_type = properties.h_element_types;
+  const int nspec = mesh.points.nspec;
+  const int ngllz = mesh.points.ngllz;
+  const int ngllx = mesh.points.ngllx;
+
+  // Count the total number of distinct global indices for the medium
+  int count = 0;
+
+  for (int ispec = 0; ispec < nspec; ++ispec) {
+    // increase the count only if current element is of the medium type
+    if (element_type(ispec) == medium::value) {
+      for (int iz = 0; iz < ngllz; ++iz) {
+        for (int ix = 0; ix < ngllx; ++ix) {
+          const int index = index_mapping(ispec, iz, ix);
+          // increase the count only if the global index is not already counted
+          /// static_cast<int>(medium::value) is the index of the medium in the
+          /// enum class
+          if (assembly_index_mapping(index, static_cast<int>(medium::value)) ==
+              -1) {
+            assembly_index_mapping(index, static_cast<int>(medium::value)) =
+                count;
+            count++;
+          }
+        }
+      }
+    }
+  }
+
+  nglob = count;
+  field = specfem::kokkos::DeviceView2d<type_real, Kokkos::LayoutLeft>(
+      "specfem::compute::fields::field", nglob, medium::components);
+  h_field = specfem::kokkos::HostMirror2d<type_real, Kokkos::LayoutLeft>(
+      Kokkos::create_mirror_view(field));
+  field_dot = specfem::kokkos::DeviceView2d<type_real, Kokkos::LayoutLeft>(
+      "specfem::compute::fields::field_dot", nglob, medium::components);
+  h_field_dot = specfem::kokkos::HostMirror2d<type_real, Kokkos::LayoutLeft>(
+      Kokkos::create_mirror_view(field_dot));
+  field_dot_dot = specfem::kokkos::DeviceView2d<type_real, Kokkos::LayoutLeft>(
+      "specfem::compute::fields::field_dot_dot", nglob, medium::components);
+
+  Kokkos::parallel_for(
+      "specfem::compute::fields::field_impl::initialize_field",
+      specfem::kokkos::HostRange(0, nglob), KOKKOS_LAMBDA(const int &iglob) {
+        for (int icomp = 0; icomp < medium::components; ++icomp) {
+          h_field(iglob, icomp) = 0.0;
+          h_field_dot(iglob, icomp) = 0.0;
+          h_field_dot_dot(iglob, icomp) = 0.0;
+          h_mass_inverse(iglob, icomp) = 0.0;
+        }
+      });
+
+  Kokkos::deep_copy(field, h_field);
+  Kokkos::deep_copy(field_dot, h_field_dot);
+  Kokkos::deep_copy(field_dot_dot, h_field_dot_dot);
+  Kokkos::deep_copy(mass_inverse, h_mass_inverse);
+
+  return;
+}
+
+#endif /* _COMPUTE_FIELDS_IMPL_FIELD_IMPL_TPP_ */
+
+
+// template <typename medium>
+//   KOKKOS_INLINE_FUNCTION type_real &specfem::compute::(const int &iglob, const int &icomp) {
+//     if constexpr (std::is_same_v<medium, elastic_type>) {
+//       int index =
+//           assembly_index_mapping(iglob, static_cast<int>(medium::value));
+//       return elastic.field(index, icomp);
+//     } else if constexpr (std::is_same_v<medium, acoustic_type>) {
+//       int index =
+//           assembly_index_mapping(iglob, static_cast<int>(medium::value));
+//       return acoustic.field(index, icomp);
+//     }
+//   }
+
+//   template <typename medium>
+//   inline type_real &h_field(const int &iglob, const int &icomp) {
+//     if constexpr (std::is_same_v<medium, elastic_type>) {
+//       int index =
+//           h_assembly_index_mapping(iglob, static_cast<int>(medium::value));
+//       return elastic.h_field(index, icomp);
+//     } else if constexpr (std::is_same_v<medium, acoustic_type>) {
+//       int index =
+//           h_assembly_index_mapping(iglob, static_cast<int>(medium::value));
+//       return acoustic.h_field(index, icomp);
+//     }
+//   }
+
+//   template <typename medium>
+//   KOKKOS_INLINE_FUNCTION type_real &field_dot(const int &iglob,
+//                                               const int &icomp) {
+//     if constexpr (std::is_same_v<medium, elastic_type>) {
+//       int index =
+//           assembly_index_mapping(iglob, static_cast<int>(medium::value));
+//       return elastic.field_dot(index, icomp);
+//     } else if constexpr (std::is_same_v<medium, acoustic_type>) {
+//       int index =
+//           assembly_index_mapping(iglob, static_cast<int>(medium::value));
+//       return acoustic.field_dot(index, icomp);
+//     }
+//   }
+
+//   template <typename medium>
+//   inline type_real &h_field_dot(const int &iglob, const int &icomp) {
+//     if constexpr (std::is_same_v<medium, elastic_type>) {
+//       int index =
+//           h_assembly_index_mapping(iglob, static_cast<int>(medium::value));
+//       return elastic.h_field_dot(index, icomp);
+//     } else if constexpr (std::is_same_v<medium, acoustic_type>) {
+//       int index =
+//           h_assembly_index_mapping(iglob, static_cast<int>(medium::value));
+//       return acoustic.h_field_dot(index, icomp);
+//     }
+//   }
+
+//   template <typename medium>
+//   KOKKOS_INLINE_FUNCTION type_real &field_dot_dot(const int &iglob,
+//                                                   const int &icomp) {
+//     if constexpr (std::is_same_v<medium, elastic_type>) {
+//       int index =
+//           assembly_index_mapping(iglob, static_cast<int>(medium::value));
+//       return elastic.field_dot_dot(index, icomp);
+//     } else if constexpr (std::is_same_v<medium, acoustic_type>) {
+//       int index =
+//           assembly_index_mapping(iglob, static_cast<int>(medium::value));
+//       return acoustic.field_dot_dot(index, icomp);
+//     }
+//   }
+
+//   template <typename medium>
+//   inline type_real &h_field_dot_dot(const int &iglob, const int &icomp) {
+//     if constexpr (std::is_same_v<medium, elastic_type>) {
+//       int index =
+//           h_assembly_index_mapping(iglob, static_cast<int>(medium::value));
+//       return elastic.h_field_dot_dot(index, icomp);
+//     } else if constexpr (std::is_same_v<medium, acoustic_type>) {
+//       int index =
+//           h_assembly_index_mapping(iglob, static_cast<int>(medium::value));
+//       return acoustic.h_field_dot_dot(index, icomp);
+//     }
+//   }
+
+//   template <typename medium>
+//   KOKKOS_INLINE_FUNCTION type_real &mass_inverse(const int &iglob,
+//                                                  const int &icomp) {
+//     if constexpr (std::is_same_v<medium, elastic_type>) {
+//       int index =
+//           assembly_index_mapping(iglob, static_cast<int>(medium::value));
+//       return elastic.mass_inverse(index, icomp);
+//     } else if constexpr (std::is_same_v<medium, acoustic_type>) {
+//       int index =
+//           assembly_index_mapping(iglob, static_cast<int>(medium::value));
+//       return acoustic.mass_inverse(index, icomp);
+//     }
+//   }
+
+//   template <typename medium>
+//   inline type_real &h_mass_inverse(const int &iglob, const int &icomp) {
+//     if constexpr (std::is_same_v<medium, elastic_type>) {
+//       int index =
+//           h_assembly_index_mapping(iglob, static_cast<int>(medium::value));
+//       return elastic.h_mass_inverse(index, icomp);
+//     } else if constexpr (std::is_same_v<medium, acoustic_type>) {
+//       int index =
+//           h_assembly_index_mapping(iglob, static_cast<int>(medium::value));
+//       return acoustic.h_mass_inverse(index, icomp);
+//     }
+//   }

include/compute/fields/simulation_field.hpp

@@ -0,0 +1,65 @@
+#ifndef _COMPUTE_FIELDS_SIMULATION_FIELD_HPP_
+#define _COMPUTE_FIELDS_SIMULATION_FIELD_HPP_
+
+#include "compute/fields/impl/field_impl.hpp"
+#include "enumerations/medium.hpp"
+#include "enumerations/simulation.hpp"
+#include "enumerations/specfem_enums.hpp"
+#include "kokkos_abstractions.h"
+#include "specfem_setup.hpp"
+#include <Kokkos_Core.hpp>
+
+namespace specfem {
+namespace compute {
+
+template <typename simulation> struct simulation_field {
+
+  using simulation_type = simulation;
+  using elastic_type = specfem::enums::element::medium::elastic;
+  using acoustic_type = specfem::enums::element::medium::acoustic;
+
+  simulation_field(const specfem::compute::mesh &mesh,
+                   const specfem::compute::properties &properties);
+
+  template <typename medium>
+  KOKKOS_INLINE_FUNCTION type_real &field(const int &iglob, const int &icomp);
+
+  template <typename medium>
+  inline type_real &h_field(const int &iglob, const int &icomp);
+
+  template <typename medium>
+  KOKKOS_INLINE_FUNCTION type_real &field_dot(const int &iglob,
+                                              const int &icomp);
+
+  template <typename medium>
+  inline type_real &h_field_dot(const int &iglob, const int &icomp);
+
+  template <typename medium>
+  KOKKOS_INLINE_FUNCTION type_real &field_dot_dot(const int &iglob,
+                                                  const int &icomp);
+
+  template <typename medium>
+  inline type_real &h_field_dot_dot(const int &iglob, const int &icomp);
+
+  template <typename medium>
+  KOKKOS_INLINE_FUNCTION type_real &mass_inverse(const int &iglob,
+                                                 const int &icomp);
+
+  template <typename medium>
+  inline type_real &h_mass_inverse(const int &iglob, const int &icomp);
+
+  int nglob = 0;
+  Kokkos::View<int * [specfem::enums::element::ntypes], Kokkos::LayoutLeft,
+               specfem::kokkos::DevMemSpace>
+      assembly_index_mapping;
+  Kokkos::View<int * [specfem::enums::element::ntypes], Kokkos::LayoutLeft,
+               specfem::kokkos::HostMemSpace>
+      h_assembly_index_mapping;
+  specfem::compute::impl::field_impl<elastic_type> elastic;
+  specfem::compute::impl::field_impl<acoustic_type> acoustic;
+};
+
+} // namespace compute
+} // namespace specfem
+
+#endif /* _COMPUTE_FIELDS_SIMULATION_FIELD_HPP_ */