From 9e32980d3e76a61d841ba6bc00d781529087d700 Mon Sep 17 00:00:00 2001
From: haykh <haykh.astro@gmail.com>
Date: Fri, 19 Dec 2025 17:48:31 -0500
Subject: [PATCH 1/5] WIP concepts + requires (CPUTEST)

---
 src/archetypes/energy_dist.h       |   4 +-
 src/archetypes/field_setter.h      |  62 +++----------
 src/archetypes/particle_injector.h |  14 ++-
 src/archetypes/spatial_dist.h      |   4 +-
 src/engines/grpic.hpp              |  36 ++++----
 src/global/arch/traits.h           | 137 +++++++++++++++++++++++++++++
 src/kernels/ampere_gr.hpp          |   6 +-
 src/kernels/ampere_sr.hpp          |   6 +-
 src/kernels/aux_fields_gr.hpp      |  35 +++-----
 src/kernels/currents_deposit.hpp   |   6 +-
 src/kernels/divergences.hpp        |   4 +
 src/kernels/faraday_gr.hpp         |   4 +-
 src/kernels/faraday_sr.hpp         |   4 +-
 src/kernels/fields_bcs.hpp         |  35 ++++----
 src/kernels/fields_to_phys.hpp     |   6 +-
 src/kernels/injectors.hpp          |  22 +++--
 src/kernels/particle_moments.hpp   |   5 +-
 src/kernels/particle_pusher_gr.hpp |  21 ++++-
 src/kernels/particle_pusher_sr.hpp |   5 +-
 src/kernels/prtls_to_phys.hpp      |   5 +-
 src/kernels/reduced_stats.hpp      |  19 ++--
 21 files changed, 290 insertions(+), 150 deletions(-)

diff --git a/src/archetypes/energy_dist.h b/src/archetypes/energy_dist.h
index 56a797751..f387922d3 100644
--- a/src/archetypes/energy_dist.h
+++ b/src/archetypes/energy_dist.h
@@ -22,6 +22,7 @@
 #include "global.h"
 
 #include "arch/kokkos_aliases.h"
+#include "arch/traits.h"
 #include "utils/comparators.h"
 #include "utils/error.h"
 #include "utils/numeric.h"
@@ -33,10 +34,9 @@ namespace arch {
   using namespace ntt;
 
   template <SimEngine::type S, class M>
+    requires traits::metric::HasD<M>
   struct EnergyDistribution {
     static constexpr auto D = M::Dim;
-    static constexpr bool is_energy_dist { true };
-    static_assert(M::is_metric, "M must be a metric class");
 
     EnergyDistribution(const M& metric) : metric { metric } {}
 
diff --git a/src/archetypes/field_setter.h b/src/archetypes/field_setter.h
index 5c5c4dbe4..56a24af15 100644
--- a/src/archetypes/field_setter.h
+++ b/src/archetypes/field_setter.h
@@ -33,6 +33,10 @@ namespace arch {
   using namespace ntt;
 
   template <class I, SimEngine::type S, class M>
+    requires traits::metric::HasD<M> &&
+             ((S == SimEngine::SRPIC && traits::metric::HasConvert<M> &&
+               traits::metric::HasTransform_i<M>) ||
+              (S == SimEngine::GRPIC && traits::metric::HasConvert_i<M>))
   class SetEMFields_kernel {
     static constexpr Dimension D = M::Dim;
     static constexpr bool defines_ex1 = traits::has_method<traits::ex1_t, I>::value;
@@ -53,7 +57,6 @@ namespace arch {
                     (defines_dx1 == defines_dx2 && defines_dx2 == defines_dx3 &&
                      defines_bx1 == defines_bx2 && defines_bx2 == defines_bx3),
                   "In GR mode, all components must be defined or none");
-    static_assert(M::is_metric, "M must be a metric class");
 
     ndfield_t<M::Dim, 6> EM;
     const I              finit;
@@ -271,69 +274,24 @@ namespace arch {
 
           if constexpr (defines_dx1 && defines_dx2 && defines_dx3) {
             { // dx1
-              vec_t<Dim::_3D> d_PU { finit.dx1({ x1_H, x2_0, x3_0 }),
-                                     finit.dx2({ x1_H, x2_0, x3_0 }),
-                                     finit.dx3({ x1_H, x2_0, x3_0 }) };
-              vec_t<Dim::_3D> d_U { ZERO };
-              metric.template transform<Idx::PU, Idx::U>({ i1_ + HALF, i2_, i3_ },
-                                                         d_PU,
-                                                         d_U);
-              EM(i1, i2, i3, em::dx1) = d_U[0];
+              EM(i1, i2, i3, em::dx1) = finit.dx1({ x1_H, x2_0, x3_0 });
             }
             { // dx2
-              vec_t<Dim::_3D> d_PU { finit.dx1({ x1_0, x2_H, x3_0 }),
-                                     finit.dx2({ x1_0, x2_H, x3_0 }),
-                                     finit.dx3({ x1_0, x2_H, x3_0 }) };
-              vec_t<Dim::_3D> d_U { ZERO };
-              metric.template transform<Idx::PU, Idx::U>({ i1_, i2_ + HALF, i3_ },
-                                                         d_PU,
-                                                         d_U);
-              EM(i1, i2, i3, em::dx2) = d_U[1];
+              EM(i1, i2, i3, em::dx2) = finit.dx2({ x1_0, x2_H, x3_0 });
             }
             { // dx3
-              vec_t<Dim::_3D> d_PU { finit.dx1({ x1_0, x2_0, x3_H }),
-                                     finit.dx2({ x1_0, x2_0, x3_H }),
-                                     finit.dx3({ x1_0, x2_0, x3_H }) };
-              vec_t<Dim::_3D> d_U { ZERO };
-              metric.template transform<Idx::PU, Idx::U>({ i1_, i2_, i3_ + HALF },
-                                                         d_PU,
-                                                         d_U);
-              EM(i1, i2, i3, em::dx3) = d_U[2];
+              EM(i1, i2, i3, em::dx3) = finit.dx3({ x1_0, x2_0, x3_H });
             }
           }
           if constexpr (defines_bx1 && defines_bx2 && defines_bx3) {
             { // bx1
-              vec_t<Dim::_3D> b_PU { finit.bx1({ x1_0, x2_H, x3_H }),
-                                     finit.bx2({ x1_0, x2_H, x3_H }),
-                                     finit.bx3({ x1_0, x2_H, x3_H }) };
-              vec_t<Dim::_3D> b_U { ZERO };
-              metric.template transform<Idx::PU, Idx::U>(
-                { i1_, i2_ + HALF, i3_ + HALF },
-                b_PU,
-                b_U);
-              EM(i1, i2, i3, em::bx1) = b_U[0];
+              EM(i1, i2, i3, em::bx1) = finit.bx1({ x1_0, x2_H, x3_H });
             }
             { // bx2
-              vec_t<Dim::_3D> b_PU { finit.bx1({ x1_H, x2_0, x3_H }),
-                                     finit.bx2({ x1_H, x2_0, x3_H }),
-                                     finit.bx3({ x1_H, x2_0, x3_H }) };
-              vec_t<Dim::_3D> b_U { ZERO };
-              metric.template transform<Idx::PU, Idx::U>(
-                { i1_ + HALF, i2_, i3_ + HALF },
-                b_PU,
-                b_U);
-              EM(i1, i2, i3, em::bx2) = b_U[1];
+              EM(i1, i2, i3, em::bx2) = finit.bx2({ x1_H, x2_0, x3_H });
             }
             { // bx3
-              vec_t<Dim::_3D> b_PU { finit.bx1({ x1_H, x2_H, x3_0 }),
-                                     finit.bx2({ x1_H, x2_H, x3_0 }),
-                                     finit.bx3({ x1_H, x2_H, x3_0 }) };
-              vec_t<Dim::_3D> b_U { ZERO };
-              metric.template transform<Idx::PU, Idx::U>(
-                { i1_ + HALF, i2_ + HALF, i3_ },
-                b_PU,
-                b_U);
-              EM(i1, i2, i3, em::bx3) = b_U[2];
+              EM(i1, i2, i3, em::bx3) = finit.bx3({ x1_H, x2_H, x3_0 });
             }
           }
         } else {
diff --git a/src/archetypes/particle_injector.h b/src/archetypes/particle_injector.h
index 79743ce65..1543e2654 100644
--- a/src/archetypes/particle_injector.h
+++ b/src/archetypes/particle_injector.h
@@ -53,6 +53,7 @@ namespace arch {
    *   - array_t<real_t*>: maximum coordinates of the region in computational coords
    */
   template <SimEngine::type S, class M>
+    requires traits::metric::HasD<M> && traits::metric::HasConvert<M>
   auto DeduceRegion(const Domain<S, M>& domain, const boundaries_t<real_t>& box)
     -> std::tuple<bool, array_t<real_t*>, array_t<real_t*>> {
     if (not domain.mesh.Intersects(box)) {
@@ -107,6 +108,7 @@ namespace arch {
    *   - array_t<real_t*>: maximum coordinates of the region in computational coords
    */
   template <SimEngine::type S, class M>
+    requires traits::metric::HasD<M>
   auto ComputeNumInject(const SimulationParams&     params,
                         const Domain<S, M>&         domain,
                         real_t                      number_density,
@@ -198,6 +200,8 @@ namespace arch {
    * @tparam ED2 Energy distribution type for species 2
    */
   template <SimEngine::type S, class M, class ED1, class ED2>
+    requires traits::metric::HasD<M> && traits::energydist::IsValid<ED1> &&
+             traits::energydist::IsValid<ED2>
   inline void InjectUniform(const SimulationParams&            params,
                             Domain<S, M>&                      domain,
                             const std::pair<spidx_t, spidx_t>& species,
@@ -205,9 +209,6 @@ namespace arch {
                             real_t                             number_density,
                             bool                        use_weights = false,
                             const boundaries_t<real_t>& box         = {}) {
-    static_assert(M::is_metric, "M must be a metric class");
-    static_assert(ED1::is_energy_dist, "ED1 must be an energy distribution class");
-    static_assert(ED2::is_energy_dist, "ED2 must be an energy distribution class");
     raise::ErrorIf((M::CoordType != Coord::Cart) && (not use_weights),
                    "Weights must be used for non-Cartesian coordinates",
                    HERE);
@@ -278,7 +279,6 @@ namespace arch {
                              spidx_t                 spidx,
                              const std::map<std::string, std::vector<real_t>>& data,
                              bool use_weights = false) {
-    static_assert(M::is_metric, "M must be a metric class");
     const auto n_inject        = data.at("ux1").size();
     auto       injector_kernel = kernel::GlobalInjector_kernel<S, M>(
       local_domain.species[spidx - 1],
@@ -309,6 +309,8 @@ namespace arch {
    * @tparam SD Spatial distribution type
    */
   template <SimEngine::type S, class M, class ED1, class ED2, class SD>
+    requires traits::metric::HasD<M> && traits::energydist::IsValid<ED1> &&
+             traits::energydist::IsValid<ED2> && traits::spatialdist::IsValid<SD>
   inline void InjectNonUniform(const SimulationParams&            params,
                                Domain<S, M>&                      domain,
                                const std::pair<spidx_t, spidx_t>& species,
@@ -317,10 +319,6 @@ namespace arch {
                                real_t                      number_density,
                                bool                        use_weights = false,
                                const boundaries_t<real_t>& box         = {}) {
-    static_assert(M::is_metric, "M must be a metric class");
-    static_assert(ED1::is_energy_dist, "ED1 must be an energy distribution class");
-    static_assert(ED2::is_energy_dist, "ED2 must be an energy distribution class");
-    static_assert(SD::is_spatial_dist, "SD must be a spatial distribution class");
     raise::ErrorIf((M::CoordType != Coord::Cart) && (not use_weights),
                    "Weights must be used for non-Cartesian coordinates",
                    HERE);
diff --git a/src/archetypes/spatial_dist.h b/src/archetypes/spatial_dist.h
index 68477208c..be42c3ba7 100644
--- a/src/archetypes/spatial_dist.h
+++ b/src/archetypes/spatial_dist.h
@@ -27,9 +27,9 @@ namespace arch {
   using namespace ntt;
 
   template <SimEngine::type S, class M>
+    requires traits::metric::HasD<M>
   struct SpatialDistribution {
-    static constexpr bool is_spatial_dist { true };
-    static_assert(M::is_metric, "M must be a metric class");
+    static constexpr auto D = M::Dim;
 
     SpatialDistribution(const M& metric) : metric { metric } {}
 
diff --git a/src/engines/grpic.hpp b/src/engines/grpic.hpp
index 35e544a99..38548469c 100644
--- a/src/engines/grpic.hpp
+++ b/src/engines/grpic.hpp
@@ -696,17 +696,17 @@ namespace ntt {
         Kokkos::parallel_for(
           "OpenBCFields",
           range,
-          kernel::bc::gr::HorizonBoundaries_kernel<M>(domain.fields.em,
-                                                      i1_min,
-                                                      tags,
-                                                      nfilter));
+          kernel::bc::gr::HorizonBoundaries_kernel<M::Dim>(domain.fields.em,
+                                                           i1_min,
+                                                           tags,
+                                                           nfilter));
         Kokkos::parallel_for(
           "OpenBCFields",
           range,
-          kernel::bc::gr::HorizonBoundaries_kernel<M>(domain.fields.em0,
-                                                      i1_min,
-                                                      tags,
-                                                      nfilter));
+          kernel::bc::gr::HorizonBoundaries_kernel<M::Dim>(domain.fields.em0,
+                                                           i1_min,
+                                                           tags,
+                                                           nfilter));
       }
     }
 
@@ -974,19 +974,19 @@ namespace ntt {
     }
 
     void TimeAverageDB(domain_t& domain) {
-      Kokkos::parallel_for("TimeAverageDB",
-                           domain.mesh.rangeActiveCells(),
-                           kernel::gr::TimeAverageDB_kernel<M>(domain.fields.em,
-                                                               domain.fields.em0,
-                                                               domain.mesh.metric));
+      Kokkos::parallel_for(
+        "TimeAverageDB",
+        domain.mesh.rangeActiveCells(),
+        kernel::gr::TimeAverageDB_kernel<M::Dim>(domain.fields.em,
+                                                 domain.fields.em0));
     }
 
     void TimeAverageJ(domain_t& domain) {
-      Kokkos::parallel_for("TimeAverageJ",
-                           domain.mesh.rangeActiveCells(),
-                           kernel::gr::TimeAverageJ_kernel<M>(domain.fields.cur,
-                                                              domain.fields.cur0,
-                                                              domain.mesh.metric));
+      Kokkos::parallel_for(
+        "TimeAverageJ",
+        domain.mesh.rangeActiveCells(),
+        kernel::gr::TimeAverageJ_kernel<M::Dim>(domain.fields.cur,
+                                                domain.fields.cur0));
     }
 
     void CurrentsDeposit(domain_t& domain) {
diff --git a/src/global/arch/traits.h b/src/global/arch/traits.h
index 269e76fc0..7e96beb0c 100644
--- a/src/global/arch/traits.h
+++ b/src/global/arch/traits.h
@@ -32,6 +32,8 @@
 #ifndef GLOBAL_ARCH_TRAITS_H
 #define GLOBAL_ARCH_TRAITS_H
 
+#include "global.h"
+
 #include <type_traits>
 #include <utility>
 
@@ -179,6 +181,141 @@ namespace traits {
   template <typename T, typename U>
   struct is_pair<std::pair<T, U>> : std::true_type {};
 
+  // c++20
+  namespace metric {
+
+    template <class M>
+    concept HasD = requires {
+      { M::Dim } -> std::convertible_to<Dimension>;
+    };
+
+    template <class M>
+    concept HasH_ij = requires(const M& m, const coord_t<M::Dim>& xi) {
+      { m.template h_<1, 1>(xi) } -> std::convertible_to<real_t>;
+    };
+
+    template <class M>
+    concept HasHij = requires(const M& m, const coord_t<M::Dim>& xi) {
+      { m.template h<1, 1>(xi) } -> std::convertible_to<real_t>;
+    };
+
+    template <class M>
+    concept HasSqrtDetH = requires(const M& m, const coord_t<M::Dim>& xi) {
+      { m.sqrt_det_h(xi) } -> std::convertible_to<real_t>;
+    };
+
+    template <class M>
+    concept HasSqrtDetHTilde = requires(const M& m, const coord_t<M::Dim>& xi) {
+      { m.sqrt_det_h_tilde(xi) } -> std::convertible_to<real_t>;
+    };
+
+    template <class M>
+    concept HasPolarArea = requires(const M& m, real_t xi_2) {
+      { m.polar_area(xi_2) } -> std::convertible_to<real_t>;
+    };
+
+    template <class M>
+    concept HasTransform_i = requires(const M&               m,
+                                      const coord_t<M::Dim>& xi,
+                                      real_t                 v_in) {
+      {
+        m.template transform<1, Idx::U, Idx::D>(xi, v_in)
+      } -> std::convertible_to<real_t>;
+    };
+
+    template <class M>
+    concept HasTransform = requires(const M&               m,
+                                    const coord_t<M::Dim>& xi,
+                                    const vec_t<Dim::_3D>& v_in,
+                                    vec_t<Dim::_3D>&       v_out) {
+      {
+        m.template transform<Idx::U, Idx::D>(xi, v_in, v_out)
+      } -> std::same_as<void>;
+    };
+
+    template <class M>
+    concept HasTransformXYZ = requires(const M&                   m,
+                                       const coord_t<M::PrtlDim>& xi,
+                                       const vec_t<Dim::_3D>&     v_in,
+                                       vec_t<Dim::_3D>&           v_out) {
+      {
+        m.template transform_xyz<Idx::XYZ, Idx::D>(xi, v_in, v_out)
+      } -> std::same_as<void>;
+    };
+
+    template <class M>
+    concept HasConvert_i = requires(const M& m, real_t x) {
+      {
+        m.template convert<1, Crd::Cd, Crd::Ph>(x)
+      } -> std::convertible_to<real_t>;
+    };
+
+    template <class M>
+    concept HasConvert = requires(const M&               m,
+                                  const coord_t<M::Dim>& x_in,
+                                  coord_t<M::Dim>&       x_out) {
+      {
+        m.template convert<Crd::Cd, Crd::Ph>(x_in, x_out)
+      } -> std::same_as<void>;
+    };
+
+    template <class M>
+    concept HasConvertXYZ = requires(const M&                   m,
+                                     const coord_t<M::PrtlDim>& x_in,
+                                     coord_t<M::PrtlDim>&       x_out) {
+      {
+        m.template convert_xyz<Crd::Cd, Crd::XYZ>(x_in, x_out)
+      } -> std::same_as<void>;
+    };
+
+    template <class M>
+    concept HasAlpha = requires(const M& m, const coord_t<M::Dim>& xi) {
+      { m.alpha(xi) } -> std::convertible_to<real_t>;
+    };
+
+    template <class M>
+    concept HasBeta1 = requires(const M& m, const coord_t<M::Dim>& xi) {
+      { m.beta1(xi) } -> std::convertible_to<real_t>;
+    };
+
+    template <class M>
+    concept HasMetricDerivatives = requires(const M& m, const coord_t<M::Dim>& xi) {
+      { m.dr_alpha(xi) } -> std::convertible_to<real_t>;
+      { m.dr_beta1(xi) } -> std::convertible_to<real_t>;
+      { m.dr_h11(xi) } -> std::convertible_to<real_t>;
+      { m.dr_h22(xi) } -> std::convertible_to<real_t>;
+      { m.dr_h33(xi) } -> std::convertible_to<real_t>;
+      { m.dr_h13(xi) } -> std::convertible_to<real_t>;
+      { m.dt_alpha(xi) } -> std::convertible_to<real_t>;
+      { m.dt_beta1(xi) } -> std::convertible_to<real_t>;
+      { m.dt_h11(xi) } -> std::convertible_to<real_t>;
+      { m.dt_h22(xi) } -> std::convertible_to<real_t>;
+      { m.dt_h33(xi) } -> std::convertible_to<real_t>;
+      { m.dt_h13(xi) } -> std::convertible_to<real_t>;
+    };
+
+  } // namespace metric
+
+  namespace energydist {
+
+    template <class ED>
+    concept IsValid = requires(const ED&             edist,
+                               const coord_t<ED::D>& x_Ph,
+                               vec_t<Dim::_3D>&      v) {
+      { edist(x_Ph, v) } -> std::same_as<void>;
+    };
+
+  } // namespace energydist
+
+  namespace spatialdist {
+
+    template <class SD>
+    concept IsValid = requires(const SD& sdist, const coord_t<SD::D>& x_Ph) {
+      { sdist(x_Ph) } -> std::convertible_to<real_t>;
+    };
+
+  } // namespace spatialdist
+
 } // namespace traits
 
 #endif // GLOBAL_ARCH_TRAITS_H
diff --git a/src/kernels/ampere_gr.hpp b/src/kernels/ampere_gr.hpp
index 327ce0cdf..d4bf569da 100644
--- a/src/kernels/ampere_gr.hpp
+++ b/src/kernels/ampere_gr.hpp
@@ -17,6 +17,7 @@
 #include "global.h"
 
 #include "arch/kokkos_aliases.h"
+#include "arch/traits.h"
 #include "utils/error.h"
 #include "utils/numeric.h"
 
@@ -29,8 +30,9 @@ namespace kernel::gr {
    * @tparam M Metric.
    */
   template <class M>
+    requires traits::metric::HasD<M> && traits::metric::HasH_ij<M> &&
+             traits::metric::HasSqrtDetH<M> && traits::metric::HasPolarArea<M>
   class Ampere_kernel {
-    static_assert(M::is_metric, "M must be a metric class");
     static constexpr auto D = M::Dim;
 
     const ndfield_t<D, 6> Din;
@@ -110,6 +112,8 @@ namespace kernel::gr {
    * @brief Add the currents to the D field with the appropriate conversion.
    */
   template <class M>
+    requires traits::metric::HasD<M> && traits::metric::HasH_ij<M> &&
+             traits::metric::HasSqrtDetH<M> && traits::metric::HasPolarArea<M>
   class CurrentsAmpere_kernel {
     static constexpr auto D = M::Dim;
 
diff --git a/src/kernels/ampere_sr.hpp b/src/kernels/ampere_sr.hpp
index dbe9d3dbf..4cb695f29 100644
--- a/src/kernels/ampere_sr.hpp
+++ b/src/kernels/ampere_sr.hpp
@@ -17,6 +17,7 @@
 #include "global.h"
 
 #include "arch/kokkos_aliases.h"
+#include "arch/traits.h"
 #include "utils/error.h"
 #include "utils/numeric.h"
 
@@ -28,8 +29,9 @@ namespace kernel::sr {
    * @tparam M Metric
    */
   template <class M>
+    requires traits::metric::HasD<M> && traits::metric::HasH_ij<M> &&
+             traits::metric::HasSqrtDetH<M> && traits::metric::HasPolarArea<M>
   class Ampere_kernel {
-    static_assert(M::is_metric, "M must be a metric class");
     static constexpr auto D = M::Dim;
 
     ndfield_t<D, 6> EB;
@@ -121,6 +123,8 @@ namespace kernel::sr {
    * @brief Add the currents to the E field with the appropriate conversion
    */
   template <class M>
+    requires traits::metric::HasD<M> && traits::metric::HasH_ij<M> &&
+             traits::metric::HasSqrtDetH<M> && traits::metric::HasPolarArea<M>
   class CurrentsAmpere_kernel {
     static constexpr auto     D     = M::Dim;
     static constexpr ncells_t i2min = N_GHOSTS;
diff --git a/src/kernels/aux_fields_gr.hpp b/src/kernels/aux_fields_gr.hpp
index 16ea11f46..ac3959e94 100644
--- a/src/kernels/aux_fields_gr.hpp
+++ b/src/kernels/aux_fields_gr.hpp
@@ -16,6 +16,7 @@
 #include "global.h"
 
 #include "arch/kokkos_aliases.h"
+#include "arch/traits.h"
 #include "utils/error.h"
 #include "utils/numeric.h"
 
@@ -27,8 +28,10 @@ namespace kernel::gr {
    * @tparam M Metric
    */
   template <class M>
+    requires traits::metric::HasD<M> && traits::metric::HasSqrtDetH<M> &&
+             traits::metric::HasSqrtDetHTilde<M> && traits::metric::HasH_ij<M> &&
+             traits::metric::HasAlpha<M> && traits::metric::HasBeta1<M>
   class ComputeAuxE_kernel {
-    static_assert(M::is_metric, "M must be a metric class");
     static constexpr auto D = M::Dim;
 
     const ndfield_t<D, 6> Df;
@@ -134,8 +137,10 @@ namespace kernel::gr {
    * @tparam M Metric
    */
   template <class M>
+    requires traits::metric::HasD<M> && traits::metric::HasSqrtDetH<M> &&
+             traits::metric::HasSqrtDetHTilde<M> && traits::metric::HasH_ij<M> &&
+             traits::metric::HasAlpha<M> && traits::metric::HasBeta1<M>
   class ComputeAuxH_kernel {
-    static_assert(M::is_metric, "M must be a metric class");
     static constexpr auto D = M::Dim;
 
     const ndfield_t<D, 6> Df;
@@ -241,22 +246,15 @@ namespace kernel::gr {
    * @brief Kernel for computing time average of B and D
    * @tparam M Metric
    */
-  template <class M>
+  template <Dimension D>
   class TimeAverageDB_kernel {
-    static_assert(M::is_metric, "M must be a metric class");
-    static constexpr auto D = M::Dim;
-
     const ndfield_t<D, 6> BDf;
     ndfield_t<D, 6>       BDf0;
-    const M               metric;
 
   public:
-    TimeAverageDB_kernel(const ndfield_t<D, 6>& BDf,
-                         const ndfield_t<D, 6>& BDf0,
-                         const M&               metric)
+    TimeAverageDB_kernel(const ndfield_t<D, 6>& BDf, const ndfield_t<D, 6>& BDf0)
       : BDf { BDf }
-      , BDf0 { BDf0 }
-      , metric { metric } {}
+      , BDf0 { BDf0 } {}
 
     Inline void operator()(index_t i1, index_t i2) const {
       if constexpr (D == Dim::_2D) {
@@ -284,22 +282,15 @@ namespace kernel::gr {
    * @brief Kernel for computing time average of J
    * @tparam M Metric
    */
-  template <class M>
+  template <Dimension D>
   class TimeAverageJ_kernel {
-    static_assert(M::is_metric, "M must be a metric class");
-    static constexpr auto D = M::Dim;
-
     ndfield_t<D, 3>       Jf;
     const ndfield_t<D, 3> Jf0;
-    const M               metric;
 
   public:
-    TimeAverageJ_kernel(const ndfield_t<D, 3>& Jf,
-                        const ndfield_t<D, 3>& Jf0,
-                        const M&               metric)
+    TimeAverageJ_kernel(const ndfield_t<D, 3>& Jf, const ndfield_t<D, 3>& Jf0)
       : Jf { Jf }
-      , Jf0 { Jf0 }
-      , metric { metric } {}
+      , Jf0 { Jf0 } {}
 
     Inline void operator()(index_t i1, index_t i2) const {
       if constexpr (D == Dim::_2D) {
diff --git a/src/kernels/currents_deposit.hpp b/src/kernels/currents_deposit.hpp
index 3bd205949..24273da65 100644
--- a/src/kernels/currents_deposit.hpp
+++ b/src/kernels/currents_deposit.hpp
@@ -14,6 +14,7 @@
 #include "global.h"
 
 #include "arch/kokkos_aliases.h"
+#include "arch/traits.h"
 #include "utils/error.h"
 #include "utils/numeric.h"
 
@@ -30,8 +31,11 @@ namespace kernel {
    * @brief Algorithm for the current deposition
    */
   template <SimEngine::type S, class M, unsigned short O = 1u>
+    requires traits::metric::HasD<M> &&
+             ((S == SimEngine::SRPIC && traits::metric::HasTransformXYZ<M>) ||
+              (S == SimEngine::GRPIC && traits::metric::HasTransform<M> &&
+               traits::metric::HasConvert_i<M> && traits::metric::HasAlpha<M>))
   class DepositCurrents_kernel {
-    static_assert(M::is_metric, "M must be a metric class");
     static_assert(O <= 11u, "Shape function order O must be <= 11");
     static constexpr auto D = M::Dim;
 
diff --git a/src/kernels/divergences.hpp b/src/kernels/divergences.hpp
index c60be564b..0635798ed 100644
--- a/src/kernels/divergences.hpp
+++ b/src/kernels/divergences.hpp
@@ -10,16 +10,20 @@
 #ifndef KERNELS_DIVERGENCES_HPP
 #define KERNELS_DIVERGENCES_HPP
 
+#include "enums.h"
 #include "global.h"
 
 #include "arch/kokkos_aliases.h"
+#include "arch/traits.h"
 #include "utils/error.h"
+#include "utils/numeric.h"
 
 namespace kernel {
   using namespace ntt;
 
   // @TODO: take care of boundaries
   template <class M, unsigned short N>
+    requires traits::metric::HasD<M> && traits::metric::HasSqrtDetH<M>
   class ComputeDivergence_kernel {
     const M metric;
 
diff --git a/src/kernels/faraday_gr.hpp b/src/kernels/faraday_gr.hpp
index b79afa460..1df53ec6d 100644
--- a/src/kernels/faraday_gr.hpp
+++ b/src/kernels/faraday_gr.hpp
@@ -16,6 +16,7 @@
 #include "global.h"
 
 #include "arch/kokkos_aliases.h"
+#include "arch/traits.h"
 #include "utils/error.h"
 #include "utils/numeric.h"
 
@@ -28,8 +29,9 @@ namespace kernel::gr {
    * @tparam M Metric
    */
   template <class M>
+    requires traits::metric::HasD<M> && traits::metric::HasH_ij<M> &&
+             traits::metric::HasSqrtDetH<M>
   class Faraday_kernel {
-    static_assert(M::is_metric, "M must be a metric class");
     static constexpr auto D = M::Dim;
 
     const ndfield_t<D, 6> Bin;
diff --git a/src/kernels/faraday_sr.hpp b/src/kernels/faraday_sr.hpp
index be3e60dd0..6d87a5bd9 100644
--- a/src/kernels/faraday_sr.hpp
+++ b/src/kernels/faraday_sr.hpp
@@ -16,6 +16,7 @@
 #include "global.h"
 
 #include "arch/kokkos_aliases.h"
+#include "arch/traits.h"
 #include "utils/error.h"
 #include "utils/numeric.h"
 
@@ -27,8 +28,9 @@ namespace kernel::sr {
    * space (diagonal metric)
    */
   template <class M>
+    requires traits::metric::HasD<M> && traits::metric::HasH_ij<M> &&
+             traits::metric::HasSqrtDetH<M>
   class Faraday_kernel {
-    static_assert(M::is_metric, "M must be a metric class");
     static constexpr auto D = M::Dim;
 
     ndfield_t<D, 6> EB;
diff --git a/src/kernels/fields_bcs.hpp b/src/kernels/fields_bcs.hpp
index 575bda879..aa252d901 100644
--- a/src/kernels/fields_bcs.hpp
+++ b/src/kernels/fields_bcs.hpp
@@ -38,8 +38,11 @@ namespace kernel::bc {
    * @note It is supposed to only be called on the active side of the absorbing edge (so sign is not needed).
    */
   template <SimEngine::type S, class I, class M, in o>
+    requires traits::metric::HasD<M> && traits::metric::HasConvert_i<M> &&
+             traits::metric::HasConvert<M> &&
+             ((S == SimEngine::SRPIC && traits::metric::HasTransform_i<M>) ||
+              S == SimEngine::GRPIC)
   struct MatchBoundaries_kernel {
-    static_assert(M::is_metric, "M must be a metric class");
     static_assert(static_cast<dim_t>(o) < static_cast<dim_t>(M::Dim),
                   "Invalid component index");
     static constexpr auto  D = M::Dim;
@@ -925,6 +928,9 @@ namespace kernel::bc {
   // };
 
   template <class I, class M, bool P, in O>
+    requires traits::metric::HasD<M> && traits::metric::HasTransform_i<M> &&
+             traits::metric::HasConvert<M> &&
+             (static_cast<dim_t>(O) < static_cast<dim_t>(M::Dim))
   struct EnforcedBoundaries_kernel {
     static constexpr Dimension D = M::Dim;
     static constexpr bool defines_ex1 = traits::has_method<traits::ex1_t, I>::value;
@@ -937,9 +943,6 @@ namespace kernel::bc {
     static_assert(defines_ex1 or defines_ex2 or defines_ex3 or defines_bx1 or
                     defines_bx2 or defines_bx3,
                   "none of the components of E or B are specified in PGEN");
-    static_assert(M::is_metric, "M must be a metric class");
-    static_assert(static_cast<dim_t>(O) < static_cast<dim_t>(M::Dim),
-                  "Invalid Orientation");
 
     ndfield_t<D, 6> Fld;
     const I         fset;
@@ -1245,17 +1248,17 @@ namespace kernel::bc {
 
   namespace gr {
 
-    template <class M>
+    template <Dimension D>
     struct HorizonBoundaries_kernel {
-      ndfield_t<M::Dim, 6> Fld;
-      const std::size_t    i1_min;
-      const bool           setE, setB;
-      const std::size_t    nfilter;
-
-      HorizonBoundaries_kernel(ndfield_t<M::Dim, 6> Fld,
-                               std::size_t          i1_min,
-                               BCTags               tags,
-                               std::size_t          nfilter)
+      ndfield_t<D, 6>   Fld;
+      const std::size_t i1_min;
+      const bool        setE, setB;
+      const std::size_t nfilter;
+
+      HorizonBoundaries_kernel(ndfield_t<D, 6> Fld,
+                               std::size_t     i1_min,
+                               BCTags          tags,
+                               std::size_t     nfilter)
         : Fld { Fld }
         , i1_min { i1_min }
         , setE { (tags & BC::Ex1 or tags & BC::Ex2 or tags & BC::Ex3) or
@@ -1265,7 +1268,7 @@ namespace kernel::bc {
         , nfilter { nfilter } {}
 
       Inline void operator()(index_t i2) const {
-        if constexpr (M::Dim == Dim::_2D) {
+        if constexpr (D == Dim::_2D) {
           if (setE) {
             for (unsigned short i = 0; i <= 2 + nfilter; ++i) {
               Fld(i1_min - N_GHOSTS + i,
@@ -1301,8 +1304,8 @@ namespace kernel::bc {
     };
 
     template <class M, idx_t i>
+      requires traits::metric::HasD<M> && traits::metric::HasConvert_i<M>
     struct AbsorbCurrents_kernel {
-      static_assert(M::is_metric, "M must be a metric class");
       static_assert(i <= static_cast<unsigned short>(M::Dim),
                     "Invalid component index");
 
diff --git a/src/kernels/fields_to_phys.hpp b/src/kernels/fields_to_phys.hpp
index a44a7b30d..daeb032e3 100644
--- a/src/kernels/fields_to_phys.hpp
+++ b/src/kernels/fields_to_phys.hpp
@@ -24,6 +24,7 @@
 #include "global.h"
 
 #include "arch/kokkos_aliases.h"
+#include "arch/traits.h"
 #include "utils/error.h"
 #include "utils/numeric.h"
 
@@ -31,10 +32,9 @@ namespace kernel {
   using namespace ntt;
 
   template <class M, int N1, int N2>
+    requires traits::metric::HasD<M> && traits::metric::HasTransform<M> &&
+             (N1 >= 3) && (N2 >= 3)
   class FieldsToPhys_kernel {
-    static_assert(M::is_metric, "M must be a metric class");
-    static_assert(N1 >= 3 && N2 >= 3, "Invalid N1 and/or N2");
-
     static constexpr auto D = M::Dim;
 
     const ndfield_t<D, N1>   Ffrom;
diff --git a/src/kernels/injectors.hpp b/src/kernels/injectors.hpp
index cddb18839..bc55b0ff2 100644
--- a/src/kernels/injectors.hpp
+++ b/src/kernels/injectors.hpp
@@ -16,6 +16,7 @@
 #include "global.h"
 
 #include "arch/kokkos_aliases.h"
+#include "arch/traits.h"
 #include "utils/error.h"
 #include "utils/numeric.h"
 
@@ -76,10 +77,11 @@ namespace kernel {
   }
 
   template <SimEngine::type S, class M, class ED1, class ED2>
+    requires traits::metric::HasD<M> && traits::metric::HasConvert<M> &&
+             ((S == SimEngine::SRPIC && traits::metric::HasTransformXYZ<M>) ||
+              (S == SimEngine::GRPIC && traits::metric::HasTransform<M>)) &&
+             traits::energydist::IsValid<ED1> && traits::energydist::IsValid<ED2>
   struct UniformInjector_kernel {
-    static_assert(ED1::is_energy_dist, "ED1 must be an energy distribution class");
-    static_assert(ED2::is_energy_dist, "ED2 must be an energy distribution class");
-    static_assert(M::is_metric, "M must be a metric class");
 
     array_t<int*>      i1s_1, i2s_1, i3s_1;
     array_t<prtldx_t*> dx1s_1, dx2s_1, dx3s_1;
@@ -291,8 +293,10 @@ namespace kernel {
   }; // struct UniformInjector_kernel
 
   template <SimEngine::type S, class M>
+    requires traits::metric::HasD<M> && traits::metric::HasConvert<M> &&
+             ((S == SimEngine::SRPIC && traits::metric::HasTransformXYZ<M>) ||
+              (S == SimEngine::GRPIC && traits::metric::HasTransform<M>))
   struct GlobalInjector_kernel {
-    static_assert(M::is_metric, "M must be a metric class");
     static constexpr auto D = M::Dim;
 
     const bool use_weights;
@@ -522,11 +526,13 @@ namespace kernel {
   }; // struct GlobalInjector_kernel
 
   template <SimEngine::type S, class M, class ED1, class ED2, class SD>
+    requires traits::metric::HasD<M> && traits::metric::HasConvert<M> &&
+             traits::metric::HasSqrtDetH<M> &&
+             ((S == SimEngine::SRPIC && traits::metric::HasTransformXYZ<M>) ||
+              (S == SimEngine::GRPIC && traits::metric::HasTransform<M>)) &&
+             traits::energydist::IsValid<ED1> &&
+             traits::energydist::IsValid<ED2> && traits::spatialdist::IsValid<SD>
   struct NonUniformInjector_kernel {
-    static_assert(M::is_metric, "M must be a metric class");
-    static_assert(ED1::is_energy_dist, "ED1 must be an energy distribution class");
-    static_assert(ED2::is_energy_dist, "ED2 must be an energy distribution class");
-    static_assert(SD::is_spatial_dist, "SD must be a spatial distribution class");
 
     const real_t ppc0;
 
diff --git a/src/kernels/particle_moments.hpp b/src/kernels/particle_moments.hpp
index a2145fa89..8fc8616aa 100644
--- a/src/kernels/particle_moments.hpp
+++ b/src/kernels/particle_moments.hpp
@@ -14,6 +14,7 @@
 #include "global.h"
 
 #include "arch/kokkos_aliases.h"
+#include "arch/traits.h"
 #include "utils/comparators.h"
 #include "utils/error.h"
 #include "utils/numeric.h"
@@ -35,8 +36,10 @@ namespace kernel {
   }
 
   template <SimEngine::type S, class M, FldsID::type F, unsigned short N>
+    requires traits::metric::HasD<M> && traits::metric::HasSqrtDetH<M> &&
+             ((S == SimEngine::SRPIC && traits::metric::HasTransformXYZ<M>) ||
+              (S == SimEngine::GRPIC && traits::metric::HasTransform<M>))
   class ParticleMoments_kernel {
-    static_assert(M::is_metric, "M must be a metric class");
     static constexpr auto D = M::Dim;
 
     static_assert(!((S == SimEngine::GRPIC) && (F == FldsID::V)),
diff --git a/src/kernels/particle_pusher_gr.hpp b/src/kernels/particle_pusher_gr.hpp
index 4432792ba..fb96660d0 100644
--- a/src/kernels/particle_pusher_gr.hpp
+++ b/src/kernels/particle_pusher_gr.hpp
@@ -18,6 +18,7 @@
 #include "global.h"
 
 #include "arch/kokkos_aliases.h"
+#include "arch/traits.h"
 #include "utils/error.h"
 #include "utils/numeric.h"
 
@@ -56,13 +57,21 @@ namespace kernel::gr {
 
   struct Massless_t {};
 
+  template <class M>
+  concept IsCompatibleWithPusherGR = traits::metric::HasD<M> &&
+                                     traits::metric::HasTransform<M> &&
+                                     traits::metric::HasHij<M> &&
+                                     traits::metric::HasAlpha<M> &&
+                                     traits::metric::HasBeta1<M> &&
+                                     traits::metric::HasMetricDerivatives<M>;
+
   /**
    * @brief Algorithm for the Particle pusher
    * @tparam M Metric
    */
   template <class M>
+    requires IsCompatibleWithPusherGR<M>
   class Pusher_kernel {
-    static_assert(M::is_metric, "M must be a metric class");
     static constexpr auto D = M::Dim;
 
   private:
@@ -97,7 +106,7 @@ namespace kernel::gr {
                   array_t<prtldx_t*>&         dx1,
                   array_t<prtldx_t*>&         dx2,
                   array_t<prtldx_t*>&         dx3,
-                  array_t<prtldx_t*>&         dx1_prev,
+                  array_t<prtldx_t*>&         dx1_pGev,
                   array_t<prtldx_t*>&         dx2_prev,
                   array_t<prtldx_t*>&         dx3_prev,
                   array_t<real_t*>&           ux1,
@@ -306,6 +315,7 @@ namespace kernel::gr {
   /* -------------------------------------------------------------------------- */
 
   template <class M>
+    requires IsCompatibleWithPusherGR<M>
   template <typename T>
   Inline void Pusher_kernel<M>::GeodesicMomentumPush(T,
                                                      const coord_t<D>&      xp,
@@ -379,6 +389,7 @@ namespace kernel::gr {
   }
 
   template <class M>
+    requires IsCompatibleWithPusherGR<M>
   template <typename T>
   Inline void Pusher_kernel<M>::GeodesicCoordinatePush(T,
                                                        const coord_t<D>& xp,
@@ -414,6 +425,7 @@ namespace kernel::gr {
   }
 
   template <class M>
+    requires IsCompatibleWithPusherGR<M>
   template <typename T>
   Inline void Pusher_kernel<M>::GeodesicFullPush(T,
                                                  const coord_t<D>&      xp,
@@ -488,6 +500,7 @@ namespace kernel::gr {
   /*                                 Phi pusher */
   /* -------------------------------------------------------------------------- */
   template <class M>
+    requires IsCompatibleWithPusherGR<M>
   template <typename T>
   Inline void Pusher_kernel<M>::UpdatePhi(T,
                                           const coord_t<D>&      xp,
@@ -511,6 +524,7 @@ namespace kernel::gr {
   }
 
   template <class M>
+    requires IsCompatibleWithPusherGR<M>
   Inline void Pusher_kernel<M>::interpolateFields(index_t          p,
                                                   vec_t<Dim::_3D>& e0,
                                                   vec_t<Dim::_3D>& b0) const {
@@ -591,6 +605,7 @@ namespace kernel::gr {
   /* ------------------------------ Photon pusher ----------------------------- */
 
   template <class M>
+    requires IsCompatibleWithPusherGR<M>
   Inline void Pusher_kernel<M>::operator()(Massless_t, index_t p) const {
     if constexpr (D == Dim::_1D) {
       raise::KernelError(HERE, "Photon pusher not implemented for 1D");
@@ -651,6 +666,7 @@ namespace kernel::gr {
   /* ------------------------- Massive particle pusher ------------------------ */
 
   template <class M>
+    requires IsCompatibleWithPusherGR<M>
   Inline void Pusher_kernel<M>::operator()(Massive_t, index_t p) const {
     if constexpr (D == Dim::_1D) {
       raise::KernelError(HERE, "Massive pusher not implemented for 1D");
@@ -743,6 +759,7 @@ namespace kernel::gr {
   // Boundary conditions
 
   template <class M>
+    requires IsCompatibleWithPusherGR<M>
   Inline void Pusher_kernel<M>::boundaryConditions(index_t p) const {
     if constexpr (D == Dim::_1D || D == Dim::_2D || D == Dim::_3D) {
       if (i1(p) < 0 && is_absorb_i1min) {
diff --git a/src/kernels/particle_pusher_sr.hpp b/src/kernels/particle_pusher_sr.hpp
index 232481450..9f9a48a86 100644
--- a/src/kernels/particle_pusher_sr.hpp
+++ b/src/kernels/particle_pusher_sr.hpp
@@ -19,6 +19,7 @@
 #include "global.h"
 
 #include "arch/kokkos_aliases.h"
+#include "arch/traits.h"
 #include "utils/error.h"
 #include "utils/numeric.h"
 
@@ -189,8 +190,10 @@ namespace kernel::sr {
    * @tparam F Additional force
    */
   template <class M, class F = NoForce_t>
+    requires traits::metric::HasD<M> && traits::metric::HasTransformXYZ<M> &&
+             traits::metric::HasConvertXYZ<M> &&
+             traits::metric::HasTransform_i<M> && traits::metric::HasConvert_i<M>
   struct Pusher_kernel {
-    static_assert(M::is_metric, "M must be a metric class");
     static constexpr auto D        = M::Dim;
     static constexpr auto ExtForce = not std::is_same<F, NoForce_t>::value;
 
diff --git a/src/kernels/prtls_to_phys.hpp b/src/kernels/prtls_to_phys.hpp
index 678fb15ab..9c6a3d8bc 100644
--- a/src/kernels/prtls_to_phys.hpp
+++ b/src/kernels/prtls_to_phys.hpp
@@ -19,6 +19,7 @@
 #include "global.h"
 
 #include "arch/kokkos_aliases.h"
+#include "arch/traits.h"
 #include "utils/error.h"
 #include "utils/numeric.h"
 
@@ -26,8 +27,10 @@ namespace kernel {
   using namespace ntt;
 
   template <SimEngine::type S, class M, bool T>
+    requires traits::metric::HasD<M> && traits::metric::HasConvert_i<M> &&
+             ((S == SimEngine::SRPIC && traits::metric::HasTransformXYZ<M>) ||
+              S == SimEngine::GRPIC && traits::metric::HasTransform<M>)
   class PrtlToPhys_kernel {
-    static_assert(M::is_metric, "M must be a metric class");
     static constexpr Dimension D = M::Dim;
 
   protected:
diff --git a/src/kernels/reduced_stats.hpp b/src/kernels/reduced_stats.hpp
index bdc7633d6..0f3e130c6 100644
--- a/src/kernels/reduced_stats.hpp
+++ b/src/kernels/reduced_stats.hpp
@@ -15,16 +15,18 @@
 #include "global.h"
 
 #include "arch/kokkos_aliases.h"
+#include "arch/traits.h"
 #include "utils/numeric.h"
 
 namespace kernel {
   using namespace ntt;
 
   template <SimEngine::type S, class M, StatsID::type F, unsigned short I = 0>
+    requires traits::metric::HasD<M> &&
+             (traits::metric::HasTransform_i<M> || I == 0) &&
+             traits::metric::HasTransform<M> && traits::metric::HasSqrtDetH<M> &&
+             (I <= 3)
   class ReducedFields_kernel {
-    static_assert(M::is_metric, "M must be a metric class");
-    static_assert(I <= 3,
-                  "I must be less than or equal to 3 for ReducedFields_kernel");
     static constexpr auto D = M::Dim;
 
     ndfield_t<D, 6> EM;
@@ -398,15 +400,14 @@ namespace kernel {
   }
 
   template <SimEngine::type S, class M, StatsID::type P>
+    requires traits::metric::HasD<M> &&
+             ((S == SimEngine::SRPIC && traits::metric::HasTransformXYZ<M>) ||
+              S == SimEngine::GRPIC && traits::metric::HasTransform<M>) &&
+             ((P == StatsID::Rho) || (P == StatsID::Charge) ||
+              (P == StatsID::N) || (P == StatsID::Npart) || (P == StatsID::T))
   class ReducedParticleMoments_kernel {
-    static_assert(M::is_metric, "M must be a metric class");
     static constexpr auto D = M::Dim;
 
-    static_assert((P == StatsID::Rho) || (P == StatsID::Charge) ||
-                    (P == StatsID::N) || (P == StatsID::Npart) ||
-                    (P == StatsID::T),
-                  "Invalid stats ID");
-
     const unsigned short     c1, c2;
     const array_t<int*>      i1, i2, i3;
     const array_t<prtldx_t*> dx1, dx2, dx3;

From 8de9280fb5424d18275c38de809dc1ea7e5a9af9 Mon Sep 17 00:00:00 2001
From: haykh <haykh.astro@gmail.com>
Date: Sat, 20 Dec 2025 00:21:03 -0500
Subject: [PATCH 2/5] compatibility bug with kokkos 5

---
 src/global/arch/kokkos_aliases.h | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/global/arch/kokkos_aliases.h b/src/global/arch/kokkos_aliases.h
index 06fff56ac..ea10ea4e1 100644
--- a/src/global/arch/kokkos_aliases.h
+++ b/src/global/arch/kokkos_aliases.h
@@ -41,7 +41,7 @@ using array_h_t = Kokkos::View<T, Kokkos::HostSpace>;
 
 // Array mirror alias of arbitrary type
 template <typename T>
-using array_mirror_t = typename array_t<T>::HostMirror;
+using array_mirror_t = typename array_t<T>::host_mirror_type;
 
 // Scatter view alias of arbitrary type
 template <typename T>
@@ -107,7 +107,7 @@ using ndfield_t = typename kokkos_aliases_hidden::ndfield_impl<D, N>::type;
 
 // D x N dimensional array (host memspace) for storing fields on ND hypercube
 template <Dimension D, unsigned short N>
-using ndfield_mirror_t = typename ndfield_t<D, N>::HostMirror;
+using ndfield_mirror_t = typename ndfield_t<D, N>::host_mirror_type;
 
 // D x N dimensional scatter array for storing fields on ND hypercubes
 namespace kokkos_aliases_hidden {

From b90cf9d5776174eb12d42ace3a88cf35770f1c88 Mon Sep 17 00:00:00 2001
From: haykh <haykh.astro@gmail.com>
Date: Sat, 20 Dec 2025 00:21:17 -0500
Subject: [PATCH 3/5] metadomain + domain concepts

---
 src/framework/domain/checkpoint.cpp     | 3 +++
 src/framework/domain/communications.cpp | 4 ++++
 src/framework/domain/domain.h           | 4 +++-
 src/framework/domain/mesh.h             | 3 ++-
 src/framework/domain/metadomain.cpp     | 9 +++++++++
 src/framework/domain/metadomain.h       | 9 +++++++--
 src/framework/domain/output.cpp         | 3 +++
 src/framework/domain/stats.cpp          | 2 ++
 src/global/arch/traits.h                | 5 +++++
 9 files changed, 38 insertions(+), 4 deletions(-)

diff --git a/src/framework/domain/checkpoint.cpp b/src/framework/domain/checkpoint.cpp
index cc4eac43e..6e10eda2e 100644
--- a/src/framework/domain/checkpoint.cpp
+++ b/src/framework/domain/checkpoint.cpp
@@ -14,6 +14,7 @@
 namespace ntt {
 
   template <SimEngine::type S, class M>
+    requires IsCompatibleWithMetadomain<M>
   void Metadomain<S, M>::InitCheckpointWriter(adios2::ADIOS*          ptr_adios,
                                               const SimulationParams& params) {
     raise::ErrorIf(ptr_adios == nullptr, "adios == nullptr", HERE);
@@ -64,6 +65,7 @@ namespace ntt {
   }
 
   template <SimEngine::type S, class M>
+    requires IsCompatibleWithMetadomain<M>
   auto Metadomain<S, M>::WriteCheckpoint(const SimulationParams& params,
                                          timestep_t              current_step,
                                          timestep_t              finished_step,
@@ -110,6 +112,7 @@ namespace ntt {
   }
 
   template <SimEngine::type S, class M>
+    requires IsCompatibleWithMetadomain<M>
   void Metadomain<S, M>::ContinueFromCheckpoint(adios2::ADIOS* ptr_adios,
                                                 const SimulationParams& params) {
     raise::ErrorIf(ptr_adios == nullptr, "adios == nullptr", HERE);
diff --git a/src/framework/domain/communications.cpp b/src/framework/domain/communications.cpp
index ae4ec83d1..66ee4e063 100644
--- a/src/framework/domain/communications.cpp
+++ b/src/framework/domain/communications.cpp
@@ -197,6 +197,7 @@ namespace ntt {
   }
 
   template <SimEngine::type S, class M>
+    requires IsCompatibleWithMetadomain<M>
   void Metadomain<S, M>::CommunicateFields(Domain<S, M>& domain, CommTags tags) {
     // const auto comm_fields = (tags & Comm::E) or (tags & Comm::B) or
     //                          (tags & Comm::J) or (tags & Comm::D) or
@@ -417,6 +418,7 @@ namespace ntt {
   }
 
   template <SimEngine::type S, class M>
+    requires IsCompatibleWithMetadomain<M>
   void Metadomain<S, M>::SynchronizeFields(Domain<S, M>&        domain,
                                            CommTags             tags,
                                            const range_tuple_t& components) {
@@ -573,6 +575,7 @@ namespace ntt {
   }
 
   template <SimEngine::type S, class M>
+    requires IsCompatibleWithMetadomain<M>
   void Metadomain<S, M>::CommunicateParticles(Domain<S, M>& domain) {
 #if defined(MPI_ENABLED)
     logger::Checkpoint("Communicating particles\n", HERE);
@@ -664,6 +667,7 @@ namespace ntt {
   }
 
   template <SimEngine::type S, class M>
+    requires IsCompatibleWithMetadomain<M>
   void Metadomain<S, M>::RemoveDeadParticles(Domain<S, M>& domain) {
     for (auto& species : domain.species) {
       species.RemoveDead();
diff --git a/src/framework/domain/domain.h b/src/framework/domain/domain.h
index b6cbb985a..c4546110b 100644
--- a/src/framework/domain/domain.h
+++ b/src/framework/domain/domain.h
@@ -43,6 +43,7 @@
 #include "global.h"
 
 #include "arch/directions.h"
+#include "arch/traits.h"
 #include "utils/formatting.h"
 #include "utils/numeric.h"
 
@@ -57,9 +58,10 @@
 #include <vector>
 
 namespace ntt {
+
   template <SimEngine::type S, class M>
+    requires traits::metric::HasD<M>
   struct Domain {
-    static_assert(M::is_metric, "template arg for Mesh class has to be a metric");
     static constexpr Dimension D { M::Dim };
 
     Mesh<M>                                 mesh;
diff --git a/src/framework/domain/mesh.h b/src/framework/domain/mesh.h
index 89f6a1b9f..3cfd9117d 100644
--- a/src/framework/domain/mesh.h
+++ b/src/framework/domain/mesh.h
@@ -17,6 +17,7 @@
 #include "global.h"
 
 #include "arch/directions.h"
+#include "arch/traits.h"
 #include "utils/comparators.h"
 #include "utils/error.h"
 #include "utils/numeric.h"
@@ -31,8 +32,8 @@
 namespace ntt {
 
   template <class M>
+    requires traits::metric::HasD<M> && traits::metric::HasConvert_i<M>
   struct Mesh : public Grid<M::Dim> {
-    static_assert(M::is_metric, "template arg for Mesh class has to be a metric");
     static constexpr bool      is_mesh { true };
     static constexpr Dimension D { M::Dim };
 
diff --git a/src/framework/domain/metadomain.cpp b/src/framework/domain/metadomain.cpp
index cbea40843..e31015c24 100644
--- a/src/framework/domain/metadomain.cpp
+++ b/src/framework/domain/metadomain.cpp
@@ -26,6 +26,7 @@
 namespace ntt {
 
   template <SimEngine::type S, class M>
+    requires IsCompatibleWithMetadomain<M>
   Metadomain<S, M>::Metadomain(unsigned int            global_ndomains,
                                const std::vector<int>& global_decomposition,
                                const std::vector<ncells_t>& global_ncells,
@@ -57,6 +58,7 @@ namespace ntt {
   }
 
   template <SimEngine::type S, class M>
+    requires IsCompatibleWithMetadomain<M>
   void Metadomain<S, M>::initialValidityCheck() const {
     // ensure everything has the correct shape
     raise::ErrorIf(g_decomposition.size() != (std::size_t)D,
@@ -94,6 +96,7 @@ namespace ntt {
   }
 
   template <SimEngine::type S, class M>
+    requires IsCompatibleWithMetadomain<M>
   void Metadomain<S, M>::createEmptyDomains() {
     /* decompose and compute cell & domain offsets ------------------------ */
     auto d_ncells = tools::Decompose(g_ndomains, g_mesh.n_active(), g_decomposition);
@@ -189,6 +192,7 @@ namespace ntt {
   }
 
   template <SimEngine::type S, class M>
+    requires IsCompatibleWithMetadomain<M>
   void Metadomain<S, M>::redefineNeighbors() {
     for (unsigned int idx { 0 }; idx < g_ndomains; ++idx) {
       // offset of the subdomain[idx]
@@ -229,6 +233,7 @@ namespace ntt {
   }
 
   template <SimEngine::type S, class M>
+    requires IsCompatibleWithMetadomain<M>
   void Metadomain<S, M>::redefineBoundaries() {
     for (unsigned int idx { 0 }; idx < g_ndomains; ++idx) {
       // offset of the subdomain[idx]
@@ -318,6 +323,7 @@ namespace ntt {
   }
 
   template <SimEngine::type S, class M>
+    requires IsCompatibleWithMetadomain<M>
   void Metadomain<S, M>::finalValidityCheck() const {
     for (unsigned int idx { 0 }; idx < g_ndomains; ++idx) {
       const auto& current_domain = g_subdomains[idx];
@@ -363,6 +369,7 @@ namespace ntt {
   }
 
   template <SimEngine::type S, class M>
+    requires IsCompatibleWithMetadomain<M>
   void Metadomain<S, M>::metricCompatibilityCheck() const {
     const auto epsilon = std::numeric_limits<real_t>::epsilon() *
                          static_cast<real_t>(100.0);
@@ -397,6 +404,7 @@ namespace ntt {
   }
 
   template <SimEngine::type S, class M>
+    requires IsCompatibleWithMetadomain<M>
   void Metadomain<S, M>::setFldsBC(const bc_in& dir, const FldsBC& new_bcs) {
     if (dir == bc_in::Mx1) {
       if constexpr (M::Dim == Dim::_1D) {
@@ -465,6 +473,7 @@ namespace ntt {
   }
 
   template <SimEngine::type S, class M>
+    requires IsCompatibleWithMetadomain<M>
   void Metadomain<S, M>::setPrtlBC(const bc_in& dir, const PrtlBC& new_bcs) {
     if (dir == bc_in::Mx1) {
       if constexpr (M::Dim == Dim::_1D) {
diff --git a/src/framework/domain/metadomain.h b/src/framework/domain/metadomain.h
index b1afcf522..bd6ec8641 100644
--- a/src/framework/domain/metadomain.h
+++ b/src/framework/domain/metadomain.h
@@ -19,6 +19,7 @@
 #include "global.h"
 
 #include "arch/kokkos_aliases.h"
+#include "arch/traits.h"
 
 #include "framework/containers/species.h"
 #include "framework/domain/domain.h"
@@ -46,10 +47,14 @@
 
 namespace ntt {
 
+  template <class M>
+  concept IsCompatibleWithMetadomain = traits::metric::HasD<M> &&
+                                       traits::metric::HasConvert<M> &&
+                                       traits::metric::HasTotVolume<M>;
+
   template <SimEngine::type S, class M>
+    requires IsCompatibleWithMetadomain<M>
   struct Metadomain {
-    static_assert(M::is_metric,
-                  "template arg for Metadomain class has to be a metric");
     static constexpr Dimension D { M::Dim };
 
     void initialValidityCheck() const;
diff --git a/src/framework/domain/output.cpp b/src/framework/domain/output.cpp
index 4d625c4d0..b4ff21235 100644
--- a/src/framework/domain/output.cpp
+++ b/src/framework/domain/output.cpp
@@ -31,6 +31,7 @@
 namespace ntt {
 
   template <SimEngine::type S, class M>
+    requires IsCompatibleWithMetadomain<M>
   void Metadomain<S, M>::InitWriter(adios2::ADIOS*          ptr_adios,
                                     const SimulationParams& params) {
     raise::ErrorIf(
@@ -263,6 +264,7 @@ namespace ntt {
   }
 
   template <SimEngine::type S, class M>
+    requires IsCompatibleWithMetadomain<M>
   void Metadomain<S, M>::CommunicateVectorPotential(unsigned short buff_idx) {
     if constexpr (M::Dim == Dim::_2D) {
       auto       local_domain = subdomain_ptr(l_subdomain_indices()[0]);
@@ -316,6 +318,7 @@ namespace ntt {
 #endif
 
   template <SimEngine::type S, class M>
+    requires IsCompatibleWithMetadomain<M>
   auto Metadomain<S, M>::Write(
     const SimulationParams&                  params,
     timestep_t                               current_step,
diff --git a/src/framework/domain/stats.cpp b/src/framework/domain/stats.cpp
index 64690d695..83d659562 100644
--- a/src/framework/domain/stats.cpp
+++ b/src/framework/domain/stats.cpp
@@ -23,6 +23,7 @@
 namespace ntt {
 
   template <SimEngine::type S, class M>
+    requires IsCompatibleWithMetadomain<M>
   void Metadomain<S, M>::InitStatsWriter(const SimulationParams& params,
                                          bool                    is_resuming) {
     raise::ErrorIf(
@@ -182,6 +183,7 @@ namespace ntt {
   }
 
   template <SimEngine::type S, class M>
+    requires IsCompatibleWithMetadomain<M>
   auto Metadomain<S, M>::WriteStats(
     const SimulationParams& params,
     timestep_t              current_step,
diff --git a/src/global/arch/traits.h b/src/global/arch/traits.h
index 7e96beb0c..6edc1d53c 100644
--- a/src/global/arch/traits.h
+++ b/src/global/arch/traits.h
@@ -214,6 +214,11 @@ namespace traits {
       { m.polar_area(xi_2) } -> std::convertible_to<real_t>;
     };
 
+    template <class M>
+    concept HasTotVolume = requires(const M& m) {
+      { m.totVolume() } -> std::convertible_to<real_t>;
+    };
+
     template <class M>
     concept HasTransform_i = requires(const M&               m,
                                       const coord_t<M::Dim>& xi,

From f9d0e801bbf7f2f95a686240eab96db64c61ed00 Mon Sep 17 00:00:00 2001
From: haykh <haykh.astro@gmail.com>
Date: Sat, 20 Dec 2025 00:40:37 -0500
Subject: [PATCH 4/5] engine+metric concepts (CPUTEST)

---
 src/archetypes/problem_generator.h      | 4 +++-
 src/archetypes/utils.h                  | 2 --
 src/engines/engine.hpp                  | 7 +++++--
 src/engines/engine_init.cpp             | 1 +
 src/engines/engine_printer.cpp          | 1 +
 src/engines/engine_run.cpp              | 1 +
 src/engines/engine_traits.h             | 2 --
 src/engines/grpic.hpp                   | 1 +
 src/engines/srpic.hpp                   | 3 +++
 src/framework/specialization_registry.h | 4 ----
 src/global/arch/traits.h                | 4 ++++
 11 files changed, 19 insertions(+), 11 deletions(-)

diff --git a/src/archetypes/problem_generator.h b/src/archetypes/problem_generator.h
index da116672d..7b9af5fad 100644
--- a/src/archetypes/problem_generator.h
+++ b/src/archetypes/problem_generator.h
@@ -23,14 +23,16 @@
 #include "enums.h"
 #include "global.h"
 
+#include "arch/traits.h"
+
 #include "framework/parameters.h"
 
 namespace arch {
   using namespace ntt;
 
   template <SimEngine::type S, class M>
+    requires traits::metric::HasD<M>
   struct ProblemGenerator {
-    static_assert(M::is_metric, "M must be a metric class");
     static constexpr bool      is_pgen { true };
     static constexpr Dimension D { M::Dim };
     static constexpr Coord     C { M::CoordType };
diff --git a/src/archetypes/utils.h b/src/archetypes/utils.h
index 3447558fb..b3eb24aeb 100644
--- a/src/archetypes/utils.h
+++ b/src/archetypes/utils.h
@@ -45,7 +45,6 @@ namespace arch {
     const std::pair<std::vector<real_t>, std::vector<real_t>>& drift_four_vels = {{ ZERO, ZERO, ZERO }, { ZERO, ZERO, ZERO }},
     bool                               use_weights = false,
     const boundaries_t<real_t>&        box         = {}) {
-    static_assert(M::is_metric, "M must be a metric class");
 
     const auto mass_1        = domain.species[species.first - 1].mass();
     const auto mass_2        = domain.species[species.second - 1].mass();
@@ -93,7 +92,6 @@ namespace arch {
     const std::pair<std::vector<real_t>, std::vector<real_t>>& drift_four_vels = {{ ZERO, ZERO, ZERO }, { ZERO, ZERO, ZERO }},
     bool                               use_weights = false,
     const boundaries_t<real_t>&        box         = {}) {
-    static_assert(M::is_metric, "M must be a metric class");
 
     InjectUniformMaxwellians<S, M>(params,
                                    domain,
diff --git a/src/engines/engine.hpp b/src/engines/engine.hpp
index 07ac82e68..88c4f20ea 100644
--- a/src/engines/engine.hpp
+++ b/src/engines/engine.hpp
@@ -49,9 +49,12 @@
 namespace ntt {
 
   template <SimEngine::type S, class M>
+  concept IsCompatibleWithEngine = traits::metric::HasD<M> &&
+                                   user::PGen<S, M>::is_pgen;
+
+  template <SimEngine::type S, class M>
+    requires IsCompatibleWithEngine<S, M>
   class Engine {
-    static_assert(M::is_metric, "template arg for Engine class has to be a metric");
-    static_assert(user::PGen<S, M>::is_pgen, "unrecognized problem generator");
 
   protected:
 #if defined(OUTPUT_ENABLED)
diff --git a/src/engines/engine_init.cpp b/src/engines/engine_init.cpp
index f98b117dd..84764e4ae 100644
--- a/src/engines/engine_init.cpp
+++ b/src/engines/engine_init.cpp
@@ -15,6 +15,7 @@
 namespace ntt {
 
   template <SimEngine::type S, class M>
+    requires IsCompatibleWithEngine<S, M>
   void Engine<S, M>::init() {
     if constexpr (pgen_is_ok) {
       m_metadomain.InitStatsWriter(m_params, is_resuming);
diff --git a/src/engines/engine_printer.cpp b/src/engines/engine_printer.cpp
index f3d6f3b38..34df741c3 100644
--- a/src/engines/engine_printer.cpp
+++ b/src/engines/engine_printer.cpp
@@ -115,6 +115,7 @@ namespace ntt {
   } // namespace
 
   template <SimEngine::type S, class M>
+    requires IsCompatibleWithEngine<S, M>
   void Engine<S, M>::print_report() const {
     const auto colored_stdout = m_params.template get<bool>(
       "diagnostics.colored_stdout");
diff --git a/src/engines/engine_run.cpp b/src/engines/engine_run.cpp
index 961e39206..aa64e0e57 100644
--- a/src/engines/engine_run.cpp
+++ b/src/engines/engine_run.cpp
@@ -11,6 +11,7 @@
 namespace ntt {
 
   template <SimEngine::type S, class M>
+    requires IsCompatibleWithEngine<S, M>
   void Engine<S, M>::run() {
     if constexpr (pgen_is_ok) {
       init();
diff --git a/src/engines/engine_traits.h b/src/engines/engine_traits.h
index 28a74bea3..d26b0addc 100644
--- a/src/engines/engine_traits.h
+++ b/src/engines/engine_traits.h
@@ -1,5 +1,3 @@
-// SPDX-License-Identifier: BSD-3-Clause
-
 #ifndef ENGINES_ENGINE_TRAITS_H
 #define ENGINES_ENGINE_TRAITS_H
 
diff --git a/src/engines/grpic.hpp b/src/engines/grpic.hpp
index 38548469c..275c70b57 100644
--- a/src/engines/grpic.hpp
+++ b/src/engines/grpic.hpp
@@ -66,6 +66,7 @@ namespace ntt {
   };
 
   template <class M>
+    requires IsCompatibleWithEngine<SimEngine::GRPIC, M>
   class GRPICEngine : public Engine<SimEngine::GRPIC, M> {
     using base_t   = Engine<SimEngine::GRPIC, M>;
     using pgen_t   = user::PGen<SimEngine::GRPIC, M>;
diff --git a/src/engines/srpic.hpp b/src/engines/srpic.hpp
index 9a4ec31fa..7c1360359 100644
--- a/src/engines/srpic.hpp
+++ b/src/engines/srpic.hpp
@@ -49,6 +49,9 @@
 namespace ntt {
 
   template <class M>
+    requires IsCompatibleWithEngine<SimEngine::SRPIC, M> &&
+             traits::metric::HasH_ij<M> && traits::metric::HasConvert_i<M> &&
+             traits::metric::HasSqrtH_ij<M>
   class SRPICEngine : public Engine<SimEngine::SRPIC, M> {
 
     using base_t   = Engine<SimEngine::SRPIC, M>;
diff --git a/src/framework/specialization_registry.h b/src/framework/specialization_registry.h
index e2628c420..3b82e3d03 100644
--- a/src/framework/specialization_registry.h
+++ b/src/framework/specialization_registry.h
@@ -1,12 +1,8 @@
-// SPDX-License-Identifier: BSD-3-Clause
-
 #ifndef FRAMEWORK_SPECIALIZATION_REGISTRY_H
 #define FRAMEWORK_SPECIALIZATION_REGISTRY_H
 
 #include "enums.h"
 
-#include "arch/traits.h"
-
 #include "metrics/kerr_schild.h"
 #include "metrics/kerr_schild_0.h"
 #include "metrics/minkowski.h"
diff --git a/src/global/arch/traits.h b/src/global/arch/traits.h
index 6edc1d53c..fd5f23cea 100644
--- a/src/global/arch/traits.h
+++ b/src/global/arch/traits.h
@@ -198,6 +198,10 @@ namespace traits {
     concept HasHij = requires(const M& m, const coord_t<M::Dim>& xi) {
       { m.template h<1, 1>(xi) } -> std::convertible_to<real_t>;
     };
+    template <class M>
+    concept HasSqrtH_ij = requires(const M& m, const coord_t<M::Dim>& xi) {
+      { m.template sqrt_h_<1, 1>(xi) } -> std::convertible_to<real_t>;
+    };
 
     template <class M>
     concept HasSqrtDetH = requires(const M& m, const coord_t<M::Dim>& xi) {

From c088ae0d5f92d8507871f8e4d2e319d0afa8c9a2 Mon Sep 17 00:00:00 2001
From: haykh <haykh.astro@gmail.com>
Date: Mon, 22 Dec 2025 15:40:53 -0500
Subject: [PATCH 5/5] amdgpu_targets & rocm libs in nix shell

---
 dev/nix/kokkos.nix | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/dev/nix/kokkos.nix b/dev/nix/kokkos.nix
index 455ced21b..d8ae115ce 100644
--- a/dev/nix/kokkos.nix
+++ b/dev/nix/kokkos.nix
@@ -10,7 +10,7 @@ let
   pversion = "5.0.0";
   compilerPkgs = {
     "HIP" = with pkgs.rocmPackages; [
-      llvm.rocm-merged-llvm
+      llvm.llvm
       rocm-core
       clr
       rocthrust
@@ -39,7 +39,7 @@ let
     "HIP" = [
       "-D Kokkos_ENABLE_HIP=ON"
       "-D Kokkos_ARCH_${getArch { }}=ON"
-      "-D AMDGPU_TARGETS=${builtins.replaceStrings [ "amd_" ] [ "" ] (pkgs.lib.toLower (getArch { }))}"
+      "-D GPU_TARGETS=${builtins.replaceStrings [ "amd_" ] [ "" ] (pkgs.lib.toLower (getArch { }))}"
       "-D CMAKE_CXX_COMPILER=hipcc"
     ];
     "CUDA" = [