Merge pull request #124 from PrincetonUniversity/issue-123

Added tests for checking misfit kernels compute container

include/compute/fields/data_access.tpp

@@ -400,7 +400,7 @@ template <
     typename std::enable_if_t<
         ViewType::isPointFieldType && !ViewType::simd::using_simd, int> = 0>
 KOKKOS_FORCEINLINE_FUNCTION void
-impl_load_on_device(const specfem::point::assembly_index &index,
+impl_load_on_device(const specfem::point::assembly_index<false> &index,
                     const WavefieldType &field, ViewType &point_field) {
   constexpr static auto MediumType = ViewType::medium_tag;
   const int iglob = index.iglob;
@@ -452,7 +452,7 @@ template <
     typename WavefieldType, typename ViewType,
     typename std::enable_if_t<
         ViewType::isPointFieldType && !ViewType::simd::using_simd, int> = 0>
-inline void impl_load_on_host(const specfem::point::assembly_index &index,
+inline void impl_load_on_host(const specfem::point::assembly_index<false> &index,
                        const WavefieldType &field, ViewType &point_field) {
 
   constexpr static auto MediumType = ViewType::medium_tag;
@@ -860,7 +860,7 @@ template <
     typename std::enable_if_t<
         ViewType::isPointFieldType && !ViewType::simd::using_simd, int> = 0>
 KOKKOS_FORCEINLINE_FUNCTION void
-impl_store_on_device(const specfem::point::assembly_index &index,
+impl_store_on_device(const specfem::point::assembly_index<false> &index,
                      const ViewType &point_field, const WavefieldType &field) {
 
   constexpr static auto MediumType = ViewType::medium_tag;
@@ -918,7 +918,7 @@ template <
     typename WavefieldType, typename ViewType,
     typename std::enable_if_t<
         ViewType::isPointFieldType && !ViewType::simd::using_simd, int> = 0>
-inline void impl_store_on_host(const specfem::point::assembly_index &index,
+inline void impl_store_on_host(const specfem::point::assembly_index<false> &index,
                         const ViewType &point_field,
                         const WavefieldType &field) {
 
@@ -1366,7 +1366,7 @@ template <
     typename std::enable_if_t<
         ViewType::isPointFieldType && !ViewType::simd::using_simd, int> = 0>
 KOKKOS_FORCEINLINE_FUNCTION void
-impl_add_on_device(const specfem::point::assembly_index &index,
+impl_add_on_device(const specfem::point::assembly_index<false> &index,
                    const ViewType &point_field, const WavefieldType &field) {
 
   constexpr static auto MediumType = ViewType::medium_tag;
@@ -1423,7 +1423,7 @@ template <
     typename WavefieldType, typename ViewType,
     typename std::enable_if_t<
         ViewType::isPointFieldType && !ViewType::simd::using_simd, int> = 0>
-inline void impl_add_on_host(const specfem::point::assembly_index &index,
+inline void impl_add_on_host(const specfem::point::assembly_index<false> &index,
                       const ViewType &point_field, const WavefieldType &field) {
 
   constexpr static auto MediumType = ViewType::medium_tag;

include/compute/kernels/impl/material_kernels.hpp

@@ -20,18 +20,19 @@ class material_kernels : public kernels_container<type, property> {
 
   material_kernels() = default;
 
-  material_kernels(const int nspec, const int n_element, const int ngllz,
-                   const int ngllx,
-                   const specfem::compute::properties &properties,
-                   specfem::kokkos::HostView1d<int> property_index_mapping)
+  material_kernels(
+      const int nspec, const int n_element, const int ngllz, const int ngllx,
+      const specfem::compute::mesh_to_compute_mapping &mapping,
+      const specfem::mesh::tags &tags,
+      const specfem::kokkos::HostView1d<int> property_index_mapping)
       : specfem::compute::impl::kernels::kernels_container<value_type,
                                                            property_type>(
             n_element, ngllz, ngllx) {
     int count = 0;
     for (int ispec = 0; ispec < nspec; ++ispec) {
-      const auto medium_tag = properties.h_element_types(ispec);
-      const auto property_tag = properties.h_element_property(ispec);
-      if ((medium_tag == type) && (property_tag == property)) {
+      const int ispec_mesh = mapping.compute_to_mesh(ispec);
+      const auto &tag = tags.tags_container(ispec_mesh);
+      if ((tag.medium_tag == type) && (tag.property_tag == property)) {
         property_index_mapping(ispec) = count;
         count++;
       }

include/compute/kernels/kernels.hpp

@@ -63,15 +63,17 @@ struct kernels {
   kernels() = default;
 
   /**
-   * @brief Construct a new kernels object from assembly information
+   * @brief Construct a new kernels object
    *
    * @param nspec Total number of spectral elements
    * @param ngllz Number of quadrature points in z dimension
    * @param ngllx Number of quadrature points in x dimension
-   * @param properties Material properties. Used to access element tags.
+   * @param mapping mesh to compute mapping
+   * @param tags Tags for every element in spectral element mesh
    */
   kernels(const int nspec, const int ngllz, const int ngllx,
-          const specfem::compute::properties &properties);
+          const specfem::compute::mesh_to_compute_mapping &mapping,
+          const specfem::mesh::tags &tags);
   ///@}
 
   /**
@@ -80,10 +82,19 @@ struct kernels {
    */
   void copy_to_host() {
     Kokkos::deep_copy(h_element_types, element_types);
+    Kokkos::deep_copy(h_element_property, element_property);
     Kokkos::deep_copy(h_property_index_mapping, property_index_mapping);
     elastic_isotropic.copy_to_host();
     acoustic_isotropic.copy_to_host();
   }
+
+  void copy_to_device() {
+    Kokkos::deep_copy(element_types, h_element_types);
+    Kokkos::deep_copy(element_property, h_element_property);
+    Kokkos::deep_copy(property_index_mapping, h_property_index_mapping);
+    elastic_isotropic.copy_to_device();
+    acoustic_isotropic.copy_to_device();
+  }
 };
 
 /**

include/parallel_configuration/chunk_config.hpp

@@ -1,6 +1,7 @@
 #pragma once
 
 #include "constants.hpp"
+#include "enumerations/dimension.hpp"
 #include <Kokkos_Core.hpp>
 
 namespace specfem {

include/point/assembly_index.hpp

@@ -5,14 +5,21 @@
 namespace specfem {
 namespace point {
 
+/**
+ * @brief Struct to store the assembled index for a quadrature point
+ *
+ * @tparam using_simd Flag to indicate if this is a simd index
+ */
+template <bool using_simd = false> struct assembly_index;
+
 /**
  * @brief Struct to store the assembled index for a quadrature point
  *
  * This struct stores a 1D index that corresponds to a global numbering of the
  * quadrature point within the mesh.
  *
  */
-struct assembly_index {
+template <> struct assembly_index<false> {
   int iglob; ///< Global index number of the quadrature point
 
   /**
@@ -45,7 +52,7 @@ struct assembly_index {
  * using SIMD instructions.
  *
  */
-struct simd_assembly_index {
+template <> struct assembly_index<true> {
   int number_points; ///< Number of points in the SIMD vector
   int iglob;         ///< Global index number of the quadrature point
 
@@ -67,7 +74,7 @@ struct simd_assembly_index {
    *
    */
   KOKKOS_FUNCTION
-  simd_assembly_index() = default;
+  assembly_index() = default;
 
   /**
    * @brief Constructor with values
@@ -76,9 +83,15 @@ struct simd_assembly_index {
    * @param number_points Number of points in the SIMD vector
    */
   KOKKOS_FUNCTION
-  simd_assembly_index(const int &iglob, const int &number_points)
+  assembly_index(const int &iglob, const int &number_points)
       : number_points(number_points), iglob(iglob) {}
   ///@}
 };
+
+/**
+ * @brief Type alias for the SIMD assembly index
+ *
+ */
+using simd_assembly_index = assembly_index<true>;
 } // namespace point
 } // namespace specfem

include/point/coordinates.hpp

@@ -83,14 +83,16 @@ template <> struct global_coordinates<specfem::dimension::type::dim2> {
  *
  * @tparam DimensionType Dimension of the element where the quadrature point is
  * located
+ * @tparam using_simd Flag to indicate if this is a simd index
  */
-template <specfem::dimension::type DimensionType> struct index;
+template <specfem::dimension::type DimensionType, bool using_simd = false>
+struct index;
 
 /**
  * @brief Template specialization for 2D elements
  *
  */
-template <> struct index<specfem::dimension::type::dim2> {
+template <> struct index<specfem::dimension::type::dim2, false> {
   int ispec; ///< Index of the spectral element
   int iz;    ///< Index of the quadrature point in the z direction within the
              ///< spectral element
@@ -121,25 +123,13 @@ template <> struct index<specfem::dimension::type::dim2> {
       : ispec(ispec), iz(iz), ix(ix) {}
 };
 
-/**
- * @brief Struct to store the SIMD indices associated with a quadrature point
- *
- * SIMD indices are used to access data into SIMD vector. Generally,
- * you'd pass the SIMD index to a @c load_on_device or @c load_on_host function
- * to load data from the device or host memory into the SIMD vector.
- *
- * @tparam DimensionType Dimension of the element where the quadrature point is
- * located
- */
-template <specfem::dimension::type DimensionType> struct simd_index;
-
 /**
  * @brief Template specialization for 2D elements
  *
  * @copydoc simd_index
  *
  */
-template <> struct simd_index<specfem::dimension::type::dim2> {
+template <> struct index<specfem::dimension::type::dim2, true> {
   int ispec; ///< Index associated with the spectral element at the start
              ///< of the SIMD vector
   int number_elements; ///< Number of elements stored in the SIMD vector
@@ -156,7 +146,7 @@ template <> struct simd_index<specfem::dimension::type::dim2> {
    *
    */
   KOKKOS_FUNCTION
-  simd_index() = default;
+  index() = default;
 
   /**
    * @brief Construct a new simd index object
@@ -169,8 +159,8 @@ template <> struct simd_index<specfem::dimension::type::dim2> {
    * spectral element
    */
   KOKKOS_FUNCTION
-  simd_index(const int &ispec, const int &number_elements, const int &iz,
-             const int &ix)
+  index(const int &ispec, const int &number_elements, const int &iz,
+        const int &ix)
       : ispec(ispec), number_elements(number_elements), iz(iz), ix(ix) {}
 
   /**
@@ -186,6 +176,15 @@ template <> struct simd_index<specfem::dimension::type::dim2> {
   }
 };
 
+/**
+ * @brief Alias for the simd index
+ *
+ * @tparam DimensionType Dimension of the element where the quadrature point is
+ * located
+ */
+template <specfem::dimension::type DimensionType>
+using simd_index = index<DimensionType, true>;
+
 /**
  * @brief Distance between two global coordinates
  *

include/point/kernels.hpp

@@ -81,6 +81,15 @@ struct kernels<specfem::dimension::type::dim2,
   KOKKOS_FUNCTION
   kernels() = default;
 
+  /**
+   * @brief single value constructor
+   *
+   */
+  KOKKOS_FUNCTION
+  kernels(const value_type value)
+      : rho(value), mu(value), kappa(value), alpha(value), beta(value),
+        rhop(value) {}
+
   /**
    * @brief Constructor
    *
@@ -95,6 +104,34 @@ struct kernels<specfem::dimension::type::dim2,
   kernels(const value_type rho, const value_type mu, const value_type kappa,
           const value_type rhop, const value_type alpha, const value_type beta)
       : rho(rho), mu(mu), kappa(kappa), rhop(rhop), alpha(alpha), beta(beta) {}
+
+  ///@}
+
+  /**
+   * @brief Equality operator
+   *
+   */
+  KOKKOS_FUNCTION
+  bool operator==(const kernels &rhs) const {
+    return rho == rhs.rho && mu == rhs.mu && kappa == rhs.kappa &&
+           rhop == rhs.rhop && alpha == rhs.alpha && beta == rhs.beta;
+  }
+
+  /**
+   * @brief Inequality operator
+   *
+   */
+  KOKKOS_FUNCTION
+  bool operator!=(const kernels &rhs) const { return !(*this == rhs); }
+
+  KOKKOS_FUNCTION
+  bool operator==(const value_type value) {
+    return rho == value && mu == value && kappa == value && rhop == value &&
+           alpha == value && beta == value;
+  }
+
+  KOKKOS_FUNCTION
+  bool operator!=(const value_type value) { return !(*this == value); }
 };
 
 /**
@@ -155,6 +192,14 @@ struct kernels<specfem::dimension::type::dim2,
   KOKKOS_FUNCTION
   kernels() = default;
 
+  /**
+   * @brief single value constructor
+   *
+   */
+  KOKKOS_FUNCTION
+  kernels(const value_type value)
+      : rho(value), kappa(value), rhop(value), alpha(value) {}
+
   /**
    * @brief Constructor
    *
@@ -168,6 +213,31 @@ struct kernels<specfem::dimension::type::dim2,
     alpha = static_cast<type_real>(2.0) * kappa;
   }
   ///@}
+
+  /**
+   * @brief Equality operator
+   *
+   */
+  KOKKOS_FUNCTION
+  bool operator==(const kernels &rhs) const {
+    return rho == rhs.rho && kappa == rhs.kappa && rhop == rhs.rhop &&
+           alpha == rhs.alpha;
+  }
+
+  /**
+   * @brief Inequality operator
+   *
+   */
+  KOKKOS_FUNCTION
+  bool operator!=(const kernels &rhs) const { return !(*this == rhs); }
+
+  KOKKOS_FUNCTION
+  bool operator==(const value_type value) {
+    return rho == value && kappa == value && rhop == value && alpha == value;
+  }
+
+  KOKKOS_FUNCTION
+  bool operator!=(const value_type value) { return !(*this == value); }
 };
 
 } // namespace point

include/policies/chunk.hpp

@@ -1,5 +1,7 @@
 #pragma once
 
+#include "enumerations/dimension.hpp"
+#include "point/coordinates.hpp"
 #include <Kokkos_Core.hpp>
 #include <string>
 #include <type_traits>

include/policies/range.hpp

@@ -19,10 +19,11 @@ template <bool UseSIMD> struct range_index_type;
  *
  */
 template <> struct range_index_type<false> {
-  specfem::point::assembly_index index; ///< Assembly index
+  specfem::point::assembly_index<false> index; ///< Assembly index
 
   KOKKOS_INLINE_FUNCTION
-  range_index_type(const specfem::point::assembly_index index) : index(index) {}
+  range_index_type(const specfem::point::assembly_index<false> index)
+      : index(index) {}
 };
 
 /**
@@ -72,7 +73,7 @@ template <typename SIMD> struct range {
   KOKKOS_INLINE_FUNCTION
   impl::range_index_type<false> operator()(const int i, std::false_type) const {
     return impl::range_index_type<false>(
-        specfem::point::assembly_index{ starting_index });
+        specfem::point::assembly_index<false>{ starting_index });
   }
 
   KOKKOS_INLINE_FUNCTION