C++20 Concepts

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" = [

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 } {}
 

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 {

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,16 +200,15 @@ 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,
                             const std::pair<ED1, ED2>&         energy_dists,
                             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);

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 };

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 } {}
 

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,

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)

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);

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");

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();

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
 

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>;
@@ -696,17 +697,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 +975,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) {

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>;

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);