diff --git a/.github/workflows/cpuarch.yml b/.github/workflows/cpuarch.yml
index 672d2639..56184e03 100644
--- a/.github/workflows/cpuarch.yml
+++ b/.github/workflows/cpuarch.yml
@@ -2,6 +2,8 @@ name: CPU Compilation/Unit Tests
 
 on:
   push:
+  pull_request:
+    types: [ready_for_review] 
 
 jobs:
   check-commit:
@@ -14,6 +16,10 @@ jobs:
       - name: Check commit message
         id: check_message
         run: |
+          if [[ "${{ github.event_name }}" == "pull_request" && "${{ github.event.action }}" == "ready_for_review" ]]; then
+            echo "run_tests=true" >> "$GITHUB_OUTPUT"
+            exit 0
+          fi
           if git log -1 --pretty=%B | grep -q "CPUTEST"; then
             echo "run_tests=true" >> "$GITHUB_OUTPUT"
           else
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 260c3e87..109a1f2f 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -7,7 +7,7 @@ set(PROJECT_NAME entity)
 
 project(
   ${PROJECT_NAME}
-  VERSION 1.3.0
+  VERSION 1.4.0
   LANGUAGES CXX C)
 add_compile_options("-D ENTITY_VERSION=\"${PROJECT_VERSION}\"")
 set(hash_cmd "git diff --quiet src/ && echo $(git rev-parse HEAD) ")
@@ -62,7 +62,7 @@ set(gpu_aware_mpi
     CACHE BOOL "Enable GPU-aware MPI")
 
 # -------------------------- Compilation settings -------------------------- #
-set(CMAKE_CXX_STANDARD 17)
+set(CMAKE_CXX_STANDARD 20)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
 
@@ -162,9 +162,9 @@ if(${output})
   find_or_fetch_dependency(adios2 FALSE QUIET)
   add_compile_options("-D OUTPUT_ENABLED")
   if(${mpi})
-    set(DEPENDENCIES ${DEPENDENCIES} adios2::cxx11_mpi)
+    set(DEPENDENCIES ${DEPENDENCIES} adios2::cxx_mpi)
   else()
-    set(DEPENDENCIES ${DEPENDENCIES} adios2::cxx11)
+    set(DEPENDENCIES ${DEPENDENCIES} adios2::cxx)
   endif()
 endif()
 
diff --git a/cmake/dependencies.cmake b/cmake/dependencies.cmake
index 1f3ed3c6..6561c304 100644
--- a/cmake/dependencies.cmake
+++ b/cmake/dependencies.cmake
@@ -58,12 +58,12 @@ function(find_or_fetch_dependency package_name header_only mode)
         FetchContent_Declare(
           ${package_name}
           GIT_REPOSITORY ${${package_name}_REPOSITORY}
-          GIT_TAG 4.7.01)
+          GIT_TAG 5.0.0)
       elseif(${package_name} STREQUAL "adios2")
         FetchContent_Declare(
           ${package_name}
           GIT_REPOSITORY ${${package_name}_REPOSITORY}
-          GIT_TAG v2.10.2)
+          GIT_TAG v2.11.0)
       else()
         FetchContent_Declare(${package_name}
                              GIT_REPOSITORY ${${package_name}_REPOSITORY})
diff --git a/dev/nix/adios2.nix b/dev/nix/adios2.nix
index a3890b78..810c228d 100644
--- a/dev/nix/adios2.nix
+++ b/dev/nix/adios2.nix
@@ -6,9 +6,9 @@
 
 let
   name = "adios2";
-  version = "2.10.2";
+  version = "2.11.0";
   cmakeFlags = {
-    CMAKE_CXX_STANDARD = "17";
+    CMAKE_CXX_STANDARD = "20";
     CMAKE_CXX_EXTENSIONS = "OFF";
     CMAKE_POSITION_INDEPENDENT_CODE = "TRUE";
     BUILD_SHARED_LIBS = "ON";
@@ -30,7 +30,7 @@ stdenv.mkDerivation {
   src = pkgs.fetchgit {
     url = "https://github.com/ornladios/ADIOS2/";
     rev = "v${version}";
-    sha256 = "sha256-NVyw7xoPutXeUS87jjVv1YxJnwNGZAT4QfkBLzvQbwg=";
+    sha256 = "sha256-yHPI///17poiCEb7Luu5qfqxTWm9Nh+o9r57mZT26U0=";
   };
 
   nativeBuildInputs = with pkgs; [
diff --git a/dev/nix/kokkos.nix b/dev/nix/kokkos.nix
index 7d86e665..d8ae115c 100644
--- a/dev/nix/kokkos.nix
+++ b/dev/nix/kokkos.nix
@@ -7,10 +7,10 @@
 
 let
   name = "kokkos";
-  pversion = "4.7.01";
+  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" = [
@@ -56,7 +56,7 @@ pkgs.stdenv.mkDerivation rec {
   src = pkgs.fetchgit {
     url = "https://github.com/kokkos/kokkos/";
     rev = "${pversion}";
-    sha256 = "sha256-MgphOsKE8umgYxVQZzex+elgvDDC09JaMCoU5YXaLco=";
+    sha256 = "sha256-C4DarqnEcdF3+19TPqcM0A9bcQSkKTJkB8b7OkzC7T8=";
   };
 
   nativeBuildInputs = with pkgs; [
@@ -78,7 +78,7 @@ pkgs.stdenv.mkDerivation rec {
 
   configurePhase = ''
     cmake -B build -D CMAKE_BUILD_TYPE=Release \
-      -D CMAKE_CXX_STANDARD=17 \
+      -D CMAKE_CXX_STANDARD=20 \
       -D CMAKE_CXX_EXTENSIONS=OFF \
       -D CMAKE_POSITION_INDEPENDENT_CODE=TRUE \
       ${pkgs.lib.concatStringsSep " " cmakeExtraFlags.${gpu}} \
diff --git a/dev/nix/shell.nix b/dev/nix/shell.nix
index 0d4cc911..5ff7a389 100644
--- a/dev/nix/shell.nix
+++ b/dev/nix/shell.nix
@@ -45,8 +45,7 @@ pkgs.mkShell {
     adios2Pkg
     kokkosPkg
 
-    python312
-    python312Packages.jupyter
+    python314
 
     cmake-format
     cmake-lint
diff --git a/extern/Kokkos b/extern/Kokkos
index 1b1383c6..37f70304 160000
--- a/extern/Kokkos
+++ b/extern/Kokkos
@@ -1 +1 @@
-Subproject commit 1b1383c6001f3bfe9fe309ca923c2d786600cc79
+Subproject commit 37f70304dc3676691af88d3ac3ba50cddbfa337f
diff --git a/extern/adios2 b/extern/adios2
index a19dad6c..1ef0b579 160000
--- a/extern/adios2
+++ b/extern/adios2
@@ -1 +1 @@
-Subproject commit a19dad6cecb00319825f20fd9f455ebbab903d34
+Subproject commit 1ef0b5797aeb8a1cc1bb36ec4089eaec19a2eea0
diff --git a/extern/plog b/extern/plog
index e21baecd..7881eeac 160000
--- a/extern/plog
+++ b/extern/plog
@@ -1 +1 @@
-Subproject commit e21baecd4753f14da64ede979c5a19302618b752
+Subproject commit 7881eeacf0f31b7862f90b88e62be7ecf0dd0a45
diff --git a/minimal/CMakeLists.txt b/minimal/CMakeLists.txt
index b21dd0fe..72e30026 100644
--- a/minimal/CMakeLists.txt
+++ b/minimal/CMakeLists.txt
@@ -89,7 +89,7 @@ if("ADIOS2_NOMPI" IN_LIST MODES)
 
   find_kokkos()
   find_adios2()
-  list(APPEND libs Kokkos::kokkos adios2::cxx11)
+  list(APPEND libs Kokkos::kokkos adios2::cxx)
 
   add_executable(${exec} ${src})
 
@@ -104,7 +104,7 @@ if("ADIOS2_MPI" IN_LIST MODES)
   find_package(MPI REQUIRED)
   find_kokkos()
   find_adios2()
-  list(APPEND libs MPI::MPI_CXX Kokkos::kokkos adios2::cxx11_mpi)
+  list(APPEND libs MPI::MPI_CXX Kokkos::kokkos adios2::cxx_mpi)
 
   add_executable(${exec} ${src})
 
diff --git a/minimal/adios2.cpp b/minimal/adios2.cpp
index cd4ca3d6..c89bf383 100644
--- a/minimal/adios2.cpp
+++ b/minimal/adios2.cpp
@@ -1,6 +1,6 @@
 #include <Kokkos_Core.hpp>
 #include <adios2.h>
-#include <adios2/cxx11/KokkosView.h>
+#include <adios2/cxx/KokkosView.h>
 
 #if defined(MPI_ENABLED)
   #include <mpi.h>
diff --git a/pgens/magnetosphere/pgen.hpp b/pgens/magnetosphere/pgen.hpp
index 1afc8ccc..eab55ffb 100644
--- a/pgens/magnetosphere/pgen.hpp
+++ b/pgens/magnetosphere/pgen.hpp
@@ -52,13 +52,13 @@ namespace user {
 
   template <Dimension D>
   struct DriveFields : public InitFields<D> {
-    DriveFields(real_t             time,
+    DriveFields(simtime_t          time,
                 real_t             bsurf,
                 real_t             rstar,
                 real_t             omega,
                 const std::string& field_geometry)
       : InitFields<D> { bsurf, rstar, field_geometry }
-      , time { time }
+      , time { (real_t)time }
       , Omega { omega } {}
 
     using InitFields<D>::bx1;
@@ -113,11 +113,11 @@ namespace user {
 
     inline PGen() {}
 
-    auto AtmFields(real_t time) const -> DriveFields<D> {
+    auto AtmFields(simtime_t time) const -> DriveFields<D> {
       return DriveFields<D> { time, Bsurf, Rstar, Omega, field_geom };
     }
 
-    auto MatchFields(real_t) const -> InitFields<D> {
+    auto MatchFields(simtime_t) const -> InitFields<D> {
       return InitFields<D> { Bsurf, Rstar, field_geom };
     }
   };
diff --git a/pgens/shock/pgen.hpp b/pgens/shock/pgen.hpp
index fc579777..c3904f02 100644
--- a/pgens/shock/pgen.hpp
+++ b/pgens/shock/pgen.hpp
@@ -113,7 +113,7 @@ namespace user {
 
     inline PGen() {}
 
-    auto MatchFields(real_t time) const -> InitFields<D> {
+    auto MatchFields(simtime_t) const -> InitFields<D> {
       return init_flds;
     }
 
diff --git a/src/archetypes/energy_dist.h b/src/archetypes/energy_dist.h
index 56a79775..f387922d 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 5c5c4dbe..56a24af1 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 79743ce6..1543e265 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/problem_generator.h b/src/archetypes/problem_generator.h
index da116672..efcec90a 100644
--- a/src/archetypes/problem_generator.h
+++ b/src/archetypes/problem_generator.h
@@ -23,15 +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/spatial_dist.h b/src/archetypes/spatial_dist.h
index 68477208..be42c3ba 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/archetypes/utils.h b/src/archetypes/utils.h
index 3447558f..b3eb24ae 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 17103f1d..3117b082 100644
--- a/src/engines/engine.hpp
+++ b/src/engines/engine.hpp
@@ -36,7 +36,7 @@
 
 #if defined(OUTPUT_ENABLED)
   #include <adios2.h>
-  #include <adios2/cxx11/KokkosView.h>
+  #include <adios2/cxx/KokkosView.h>
 #endif // OUTPUT_ENABLED
 
 #if defined(MPI_ENABLED)
@@ -49,9 +49,11 @@
 namespace ntt {
 
   template <SimEngine::type S, class M>
+  concept IsCompatibleWithEngine = traits::metric::HasD<M>;
+
+  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)
@@ -83,7 +85,6 @@ namespace ntt {
     };
 
     static constexpr Dimension D { M::Dim };
-    static constexpr bool      is_engine { true };
 
     Engine(const SimulationParams& params)
       : m_params { params }
diff --git a/src/engines/engine_init.cpp b/src/engines/engine_init.cpp
index f98b117d..e822ba3f 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);
@@ -26,8 +27,7 @@ namespace ntt {
       if (not is_resuming) {
         // start a new simulation with initial conditions
         logger::Checkpoint("Loading initial conditions", HERE);
-        if constexpr (
-          traits::has_member<traits::pgen::init_flds_t, user::PGen<S, M>>::value) {
+        if constexpr (traits::pgen::HasInitFlds<user::PGen<S, M>>) {
           logger::Checkpoint("Initializing fields from problem generator", HERE);
           m_metadomain.runOnLocalDomains([&](auto& loc_dom) {
             Kokkos::parallel_for(
@@ -39,8 +39,7 @@ namespace ntt {
                 loc_dom.mesh.metric });
           });
         }
-        if constexpr (
-          traits::has_member<traits::pgen::init_prtls_t, user::PGen<S, M>>::value) {
+        if constexpr (traits::pgen::HasInitPrtls<user::PGen<S, M>, Domain<S, M>>) {
           logger::Checkpoint("Initializing particles from problem generator", HERE);
           m_metadomain.runOnLocalDomains([&](auto& loc_dom) {
             m_pgen.InitPrtls(loc_dom);
diff --git a/src/engines/engine_printer.cpp b/src/engines/engine_printer.cpp
index f3d6f3b3..34df741c 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 961e3920..7d81fe4b 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();
@@ -43,7 +44,7 @@ namespace ntt {
         });
         // poststep (if defined)
         if constexpr (
-          traits::has_method<traits::pgen::custom_poststep_t, decltype(m_pgen)>::value) {
+          traits::pgen::HasCustomPostStep<decltype(m_pgen), Domain<S, M>>) {
           timers.start("Custom");
           m_metadomain.runOnLocalDomains([&timers, this](auto& dom) {
             m_pgen.CustomPostStep(step, time, dom);
@@ -62,7 +63,7 @@ namespace ntt {
 #if defined(OUTPUT_ENABLED)
         timers.start("Output");
         if constexpr (
-          traits::has_method<traits::pgen::custom_field_output_t, decltype(m_pgen)>::value) {
+          traits::pgen::HasCustomFieldOutput<decltype(m_pgen), Domain<S, M>, M::Dim>) {
           auto lambda_custom_field_output = [&](const std::string&    name,
                                                 ndfield_t<M::Dim, 6>& buff,
                                                 index_t               idx,
@@ -81,7 +82,7 @@ namespace ntt {
           print_output &= m_metadomain.Write(m_params, step, step - 1, time, time - dt);
         }
         if constexpr (
-          traits::has_method<traits::pgen::custom_stat_t, decltype(m_pgen)>::value) {
+          traits::pgen::HasCustomStatOutput<decltype(m_pgen), Domain<S, M>>) {
           auto lambda_custom_stat = [&](const std::string&  name,
                                         timestep_t          step,
                                         simtime_t           time,
diff --git a/src/engines/engine_traits.h b/src/engines/engine_traits.h
index 28a74bea..d26b0add 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 35e544a9..275c70b5 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>;
@@ -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) {
diff --git a/src/engines/srpic.hpp b/src/engines/srpic.hpp
index 9a4ec31f..fc95e6cf 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>;
@@ -369,7 +372,7 @@ namespace ntt {
                                               : ZERO;
         // toggle to indicate whether pgen defines the external force
         bool       has_extforce = false;
-        if constexpr (traits::has_member<traits::pgen::ext_force_t, pgen_t>::value) {
+        if constexpr (traits::pgen::HasExtForce<pgen_t>) {
           has_extforce = true;
           // toggle to indicate whether the ext force applies to current species
           if (traits::has_member<traits::species_t, decltype(pgen_t::ext_force)>::value) {
@@ -442,7 +445,7 @@ namespace ntt {
                 gca_larmor_max, gca_eovrb_max, sync_coeff, comp_coeff
             ));
         } else if (not has_atmosphere and has_extforce) {
-          if constexpr (traits::has_member<traits::pgen::ext_force_t, pgen_t>::value) {
+          if constexpr (traits::pgen::HasExtForce<pgen_t>) {
             const auto force =
               kernel::sr::Force<M::PrtlDim, M::CoordType, decltype(m_pgen.ext_force), false> {
                 m_pgen.ext_force
@@ -474,7 +477,7 @@ namespace ntt {
             raise::Error("External force not implemented", HERE);
           }
         } else { // has_atmosphere and has_extforce
-          if constexpr (traits::has_member<traits::pgen::ext_force_t, pgen_t>::value) {
+          if constexpr (traits::pgen::HasExtForce<pgen_t>) {
             const auto force =
               kernel::sr::Force<M::PrtlDim, M::CoordType, decltype(m_pgen.ext_force), true> {
                 m_pgen.ext_force, {gx1, gx2, gx3}, x_surf, ds
@@ -570,8 +573,7 @@ namespace ntt {
         const auto V0    = m_params.template get<real_t>("scales.V0");
         const auto ppc0  = m_params.template get<real_t>("particles.ppc0");
         const auto coeff = -dt * q0 / (B0 * V0);
-        if constexpr (
-          traits::has_member<traits::pgen::ext_current_t, pgen_t>::value) {
+        if constexpr (traits::pgen::HasExtCurrent<pgen_t>) {
           const std::vector<real_t> xmin { domain.mesh.extent(in::x1).first,
                                            domain.mesh.extent(in::x2).first,
                                            domain.mesh.extent(in::x3).first };
@@ -720,8 +722,7 @@ namespace ntt {
       }
 
       if (dim == in::x1) {
-        if constexpr (
-          traits::has_member<traits::pgen::match_fields_t, pgen_t>::value) {
+        if constexpr (traits::pgen::HasMatchFields<pgen_t>) {
           auto match_fields = m_pgen.MatchFields(time);
           call_match_fields<decltype(match_fields), in::x1>(domain.fields.em,
                                                             domain.mesh.flds_bc(),
@@ -732,8 +733,7 @@ namespace ntt {
                                                             tags,
                                                             range_min,
                                                             range_max);
-        } else if constexpr (
-          traits::has_member<traits::pgen::match_fields_in_x1_t, pgen_t>::value) {
+        } else if constexpr (traits::pgen::HasMatchFieldsInX1<pgen_t>) {
           auto match_fields = m_pgen.MatchFieldsInX1(time);
           call_match_fields<decltype(match_fields), in::x1>(domain.fields.em,
                                                             domain.mesh.flds_bc(),
@@ -747,8 +747,7 @@ namespace ntt {
         }
       } else if (dim == in::x2) {
         if constexpr (M::Dim == Dim::_2D or M::Dim == Dim::_3D) {
-          if constexpr (
-            traits::has_member<traits::pgen::match_fields_t, pgen_t>::value) {
+          if constexpr (traits::pgen::HasMatchFields<pgen_t>) {
             auto match_fields = m_pgen.MatchFields(time);
             call_match_fields<decltype(match_fields), in::x2>(domain.fields.em,
                                                               domain.mesh.flds_bc(),
@@ -759,8 +758,7 @@ namespace ntt {
                                                               tags,
                                                               range_min,
                                                               range_max);
-          } else if constexpr (
-            traits::has_member<traits::pgen::match_fields_in_x2_t, pgen_t>::value) {
+          } else if constexpr (traits::pgen::HasMatchFieldsInX2<pgen_t>) {
             auto match_fields = m_pgen.MatchFieldsInX2(time);
             call_match_fields<decltype(match_fields), in::x2>(domain.fields.em,
                                                               domain.mesh.flds_bc(),
@@ -777,8 +775,7 @@ namespace ntt {
         }
       } else if (dim == in::x3) {
         if constexpr (M::Dim == Dim::_3D) {
-          if constexpr (
-            traits::has_member<traits::pgen::match_fields_t, pgen_t>::value) {
+          if constexpr (traits::pgen::HasMatchFields<pgen_t>) {
             auto match_fields = m_pgen.MatchFields(time);
             call_match_fields<decltype(match_fields), in::x3>(domain.fields.em,
                                                               domain.mesh.flds_bc(),
@@ -789,8 +786,7 @@ namespace ntt {
                                                               tags,
                                                               range_min,
                                                               range_max);
-          } else if constexpr (
-            traits::has_member<traits::pgen::match_fields_in_x3_t, pgen_t>::value) {
+          } else if constexpr (traits::pgen::HasMatchFieldsInX3<pgen_t>) {
             auto match_fields = m_pgen.MatchFieldsInX3(time);
             call_match_fields<decltype(match_fields), in::x3>(domain.fields.em,
                                                               domain.mesh.flds_bc(),
@@ -900,22 +896,16 @@ namespace ntt {
       if (tags & BC::B) {
         comps.push_back(normal_b_comp);
       }
-      if constexpr (traits::has_member<traits::pgen::fix_fields_t, pgen_t>::value) {
-        raise::Error("Non-const fixed fields not implemented", HERE);
-      } else if constexpr (
-        traits::has_member<traits::pgen::fix_fields_const_t, pgen_t>::value) {
+      if constexpr (traits::pgen::HasFixFieldsConst<pgen_t>) {
         for (const auto& comp : comps) {
-          auto value     = ZERO;
-          bool shouldset = false;
-          if constexpr (
-            traits::has_member<traits::pgen::fix_fields_const_t, pgen_t>::value) {
-            // if fix field function present, read from it
-            const auto newset = m_pgen.FixFieldsConst(
-              (bc_in)(sign * ((short)dim + 1)),
-              (em)comp);
-            value     = newset.first;
-            shouldset = newset.second;
-          }
+          auto       value     = ZERO;
+          bool       shouldset = false;
+          // if fix field function present, read from it
+          const auto newset    = m_pgen.FixFieldsConst(
+            (bc_in)(sign * ((short)dim + 1)),
+            (em)comp);
+          value     = newset.first;
+          shouldset = newset.second;
           if (shouldset) {
             if constexpr (M::Dim == Dim::_1D) {
               Kokkos::deep_copy(
@@ -1078,7 +1068,7 @@ namespace ntt {
       /**
        * atmosphere field boundaries
        */
-      if constexpr (traits::has_member<traits::pgen::atm_fields_t, pgen_t>::value) {
+      if constexpr (traits::pgen::HasAtmFields<pgen_t>) {
         const auto [sign, dim, xg_min, xg_max] = get_atm_extent(direction);
         const auto           dd                = static_cast<dim_t>(dim);
         boundaries_t<real_t> box;
diff --git a/src/framework/domain/checkpoint.cpp b/src/framework/domain/checkpoint.cpp
index cc4eac43..6e10eda2 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 ae4ec83d..66ee4e06 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 b6cbb985..c4546110 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 89f6a1b9..3cfd9117 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 cbea4084..e31015c2 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 a456879b..bd6ec864 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"
@@ -35,7 +36,7 @@
   #include "output/writer.h"
 
   #include <adios2.h>
-  #include <adios2/cxx11/KokkosView.h>
+  #include <adios2/cxx/KokkosView.h>
 #endif // OUTPUT_ENABLED
 
 #include <functional>
@@ -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 4d625c4d..b4ff2123 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 64690d69..83d65956 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/framework/simulation.h b/src/framework/simulation.h
index 33750030..d8696532 100644
--- a/src/framework/simulation.h
+++ b/src/framework/simulation.h
@@ -36,12 +36,9 @@ namespace ntt {
     ~Simulation();
 
     template <template <class> class E, template <Dimension> class M, Dimension D>
+      requires traits::engine::HasRun<E<M<D>>>
     inline void run() {
       using engine_t = E<M<D>>;
-      static_assert(engine_t::is_engine,
-                    "template arg for Simulation::run has to be an engine");
-      static_assert(traits::has_method<traits::run_t, engine_t>::value,
-                    "Engine must contain a ::run() method");
       try {
         engine_t engine { m_params };
         engine.run();
diff --git a/src/framework/specialization_registry.h b/src/framework/specialization_registry.h
index e2628c42..3b82e3d0 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/kokkos_aliases.h b/src/global/arch/kokkos_aliases.h
index 06fff56a..ea10ea4e 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 {
diff --git a/src/global/arch/traits.h b/src/global/arch/traits.h
index 269e76fc..e950af53 100644
--- a/src/global/arch/traits.h
+++ b/src/global/arch/traits.h
@@ -8,18 +8,6 @@
  *   - traits::bx1_t, ::bx2_t, ::bx3_t
  *   - traits::dx1_t, ::dx2_t, ::dx3_t
  *   - traits::run_t, traits::to_string_t
- *   - traits::pgen::init_flds_t
- *   - traits::pgen::ext_force_t
- *   - traits::pgen::ext_current_t
- *   - traits::pgen::atm_fields_t
- *   - traits::pgen::match_fields_const_t
- *   - traits::pgen::match_fields_t
- *   - traits::pgen::fix_fields_const_t
- *   - traits::pgen::fix_fields_t
- *   - traits::pgen::init_prtls_t
- *   - traits::pgen::custom_fields_t
- *   - traits::pgen::custom_field_output_t
- *   - traits::pgen::custom_poststep_t
  *   - traits::check_compatibility<>
  *   - traits::compatibility<>
  *   - traits::is_pair<>
@@ -32,6 +20,10 @@
 #ifndef GLOBAL_ARCH_TRAITS_H
 #define GLOBAL_ARCH_TRAITS_H
 
+#include "global.h"
+
+#include "arch/kokkos_aliases.h"
+
 #include <type_traits>
 #include <utility>
 
@@ -79,77 +71,6 @@ namespace traits {
   template <typename T>
   using dx3_t = decltype(&T::dx3);
 
-  // for simengine
-  template <typename T>
-  using run_t = decltype(&T::run);
-
-  // for params
-  template <typename T>
-  using to_string_t = decltype(&T::to_string);
-
-  // for pgen
-  namespace pgen {
-    template <typename T>
-    using init_flds_t = decltype(&T::init_flds);
-
-    template <typename T>
-    using init_prtls_t = decltype(&T::InitPrtls);
-
-    template <typename T>
-    using ext_force_t = decltype(&T::ext_force);
-
-    template <typename T>
-    using ext_current_t = decltype(&T::ext_current);
-
-    template <typename T>
-    using atm_fields_t = decltype(&T::AtmFields);
-
-    template <typename T>
-    using match_fields_t = decltype(&T::MatchFields);
-
-    template <typename T>
-    using match_fields_in_x1_t = decltype(&T::MatchFieldsInX1);
-
-    template <typename T>
-    using match_fields_in_x2_t = decltype(&T::MatchFieldsInX2);
-
-    template <typename T>
-    using match_fields_in_x3_t = decltype(&T::MatchFieldsInX3);
-
-    template <typename T>
-    using match_fields_const_t = decltype(&T::MatchFieldsConst);
-
-    template <typename T>
-    using fix_fields_t = decltype(&T::FixFields);
-
-    template <typename T>
-    using fix_fields_const_t = decltype(&T::FixFieldsConst);
-
-    template <typename T>
-    using perfect_conductor_fields_t = decltype(&T::PerfectConductorFields);
-
-    template <typename T>
-    using perfect_conductor_fields_const_t = decltype(&T::PerfectConductorFieldsConst);
-
-    template <typename T>
-    using perfect_conductor_currents_t = decltype(&T::PerfectConductorCurrents);
-
-    template <typename T>
-    using perfect_conductor_currents_const_t = decltype(&T::PerfectConductorCurrentsConst);
-
-    template <typename T>
-    using custom_fields_t = decltype(&T::CustomFields);
-
-    template <typename T>
-    using custom_poststep_t = decltype(&T::CustomPostStep);
-
-    template <typename T>
-    using custom_field_output_t = decltype(&T::CustomFieldOutput);
-
-    template <typename T>
-    using custom_stat_t = decltype(&T::CustomStat);
-  } // namespace pgen
-
   // for pgen extforce
   template <typename T>
   using species_t = decltype(&T::species);
@@ -179,6 +100,242 @@ namespace traits {
   template <typename T, typename U>
   struct is_pair<std::pair<T, U>> : std::true_type {};
 
+  // c++20
+  namespace engine {
+    template <class E>
+    concept HasRun = requires(E& engine) {
+      { engine.run() } -> std::same_as<void>;
+    };
+  } // namespace engine
+
+  namespace params {
+    template <class P>
+    concept HasToString = requires(const P& params) {
+      { params.to_string() } -> std::convertible_to<std::string>;
+    };
+  } // namespace params
+
+  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 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) {
+      { 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 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,
+                                      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 pgen {
+
+    template <class PG>
+    concept HasInitFlds = requires(const PG& pgen) { pgen.init_flds; };
+
+    template <class PG, class D>
+    concept HasInitPrtls = requires(const PG& pgen, D& domain) {
+      pgen.InitPrtls(domain);
+    };
+
+    template <class PG>
+    concept HasExtForce = requires(const PG& pgen) { pgen.ext_force; };
+
+    template <class PG>
+    concept HasExtCurrent = requires(const PG& pgen) { pgen.ext_current; };
+
+    template <class PG>
+    concept HasAtmFields = requires(const PG& pgen, simtime_t time) {
+      pgen.AtmFields(time);
+    };
+
+    template <class PG>
+    concept HasMatchFields = requires(const PG& pgen, simtime_t time) {
+      pgen.MatchFields(time);
+    };
+
+    template <class PG>
+    concept HasMatchFieldsInX1 = requires(const PG& pgen, simtime_t time) {
+      pgen.MatchFieldsInX1(time);
+    };
+
+    template <class PG>
+    concept HasMatchFieldsInX2 = requires(const PG& pgen, simtime_t time) {
+      pgen.MatchFieldsInX2(time);
+    };
+
+    template <class PG>
+    concept HasMatchFieldsInX3 = requires(const PG& pgen, simtime_t time) {
+      pgen.MatchFieldsInX3(time);
+    };
+
+    template <class PG>
+    concept HasFixFieldsConst = requires(const PG&      pgen,
+                                         const bc_in&   bc,
+                                         const ntt::em& comp) {
+      pgen.FixFieldsConst(bc, comp);
+    };
+
+    template <class PG, class D>
+    concept HasCustomPostStep = requires(const PG&  pgen,
+                                         timestep_t s,
+                                         simtime_t  t,
+                                         D&         domain) {
+      pgen.CustomPostStep(s, t, domain);
+    };
+
+    template <class PG, class D, Dimension Dim>
+    concept HasCustomFieldOutput = requires(const PG&          pgen,
+                                            const std::string& name,
+                                            ndfield_t<Dim, 6>& buff,
+                                            index_t            idx,
+                                            timestep_t         step,
+                                            simtime_t          time,
+                                            const D&           dom) {
+      pgen.CustomFieldOutput(name, buff, idx, step, time, dom);
+    };
+
+    template <class PG, class D>
+    concept HasCustomStatOutput = requires(const PG&          pgen,
+                                           const std::string& name,
+                                           timestep_t         s,
+                                           simtime_t          t,
+                                           const D&           dom) {
+      pgen.CustomStat(name, s, t, dom);
+    };
+
+  } // namespace pgen
+
 } // namespace traits
 
 #endif // GLOBAL_ARCH_TRAITS_H
diff --git a/src/global/utils/formatting.h b/src/global/utils/formatting.h
index 46c25c63..aaed184e 100644
--- a/src/global/utils/formatting.h
+++ b/src/global/utils/formatting.h
@@ -75,8 +75,8 @@ namespace fmt {
     oss << "[";
     if (!vec.empty()) {
       if constexpr (traits::is_pair<T>::value) {
-        if constexpr (
-          traits::has_method<traits::to_string_t, typename T::first_type>::value) {
+        if constexpr (traits::params::HasToString<typename T::first_type> and
+                      traits::params::HasToString<typename T::second_type>) {
           oss << "{" << vec[0].first.to_string() << ", "
               << vec[0].second.to_string() << "}";
           for (size_t i = 1; i < vec.size(); ++i) {
diff --git a/src/global/utils/param_container.cpp b/src/global/utils/param_container.cpp
index 126a64e4..895eec18 100644
--- a/src/global/utils/param_container.cpp
+++ b/src/global/utils/param_container.cpp
@@ -5,7 +5,7 @@
   #include "global.h"
 
   #include <adios2.h>
-  #include <adios2/cxx11/KokkosView.h>
+  #include <adios2/cxx/KokkosView.h>
 
   #include <any>
   #include <functional>
diff --git a/src/global/utils/param_container.h b/src/global/utils/param_container.h
index 679b8055..2346de2b 100644
--- a/src/global/utils/param_container.h
+++ b/src/global/utils/param_container.h
@@ -18,7 +18,7 @@
 
 #if defined(OUTPUT_ENABLED)
   #include <adios2.h>
-  #include <adios2/cxx11/KokkosView.h>
+  #include <adios2/cxx/KokkosView.h>
 #endif
 
 #include <any>
diff --git a/src/kernels/ampere_gr.hpp b/src/kernels/ampere_gr.hpp
index 327ce0cd..d4bf569d 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 dbe9d3db..4cb695f2 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 16ea11f4..ac3959e9 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 3bd20594..24273da6 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 c60be564..0635798e 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 b79afa46..1df53ec6 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 be3e60dd..6d87a5bd 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 575bda87..aa252d90 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 a44a7b30..daeb032e 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 cddb1883..bc55b0ff 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 a2145fa8..8fc8616a 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 4432792b..fb96660d 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 23248145..9f9a48a8 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 678fb15a..9c6a3d8b 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 bdc7633d..0f3e130c 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;
diff --git a/src/metrics/kerr_schild.h b/src/metrics/kerr_schild.h
index 14441f99..4df88cf4 100644
--- a/src/metrics/kerr_schild.h
+++ b/src/metrics/kerr_schild.h
@@ -40,6 +40,7 @@ namespace metric {
 
     const real_t dr, dtheta, dphi;
     const real_t dr_inv, dtheta_inv, dphi_inv;
+    const bool   small_angle;
 
     Inline auto Delta(real_t r) const -> real_t {
       return SQR(r) - TWO * r + SQR(a);
@@ -85,7 +86,8 @@ namespace metric {
       , dphi { (x3_max - x3_min) / nx3 }
       , dr_inv { ONE / dr }
       , dtheta_inv { ONE / dtheta }
-      , dphi_inv { ONE / dphi } {
+      , dphi_inv { ONE / dphi }
+      , small_angle { HALF * dtheta < constant::SMALL_ANGLE } {
       set_dxMin(find_dxMin());
     }
 
@@ -443,12 +445,21 @@ namespace metric {
     /**
      * differential area at the pole (used in axisymmetric solvers)
      * @param x1 radial coordinate along the axis (code units)
+     * @note uses small-angle approximation when the resolution is too high
      */
     Inline auto polar_area(real_t x1) const -> real_t {
-      return dr * (SQR(x1 * dr + x1_min) + SQR(a)) *
-             math::sqrt(ONE + TWO * (x1 * dr + x1_min) /
-                                (SQR(x1 * dr + x1_min) + SQR(a))) *
-             (ONE - math::cos(HALF * dtheta));
+      if (small_angle) {
+        return dr * (SQR(x1 * dr + x1_min) + SQR(a)) *
+              math::sqrt(ONE + TWO * (x1 * dr + x1_min) /
+                                  (SQR(x1 * dr + x1_min) + SQR(a))) *
+               (static_cast<real_t>(48) - SQR(dtheta)) * SQR(dtheta) /
+               static_cast<real_t>(384);
+      } else {
+        return dr * (SQR(x1 * dr + x1_min) + SQR(a)) *
+              math::sqrt(ONE + TWO * (x1 * dr + x1_min) /
+                                  (SQR(x1 * dr + x1_min) + SQR(a))) *
+              (ONE - math::cos(HALF * dtheta));
+      }
     }
 
     /**
diff --git a/src/metrics/metric_base.h b/src/metrics/metric_base.h
index 8b0e5cae..6284cd9e 100644
--- a/src/metrics/metric_base.h
+++ b/src/metrics/metric_base.h
@@ -84,7 +84,6 @@ namespace metric {
    */
   template <Dimension D>
   struct MetricBase {
-    static constexpr bool      is_metric { true };
     static constexpr Dimension Dim { D };
 
     MetricBase(const std::vector<ncells_t>& res, const boundaries_t<real_t>& ext)
diff --git a/src/metrics/qkerr_schild.h b/src/metrics/qkerr_schild.h
index 9b608f35..90f366ef 100644
--- a/src/metrics/qkerr_schild.h
+++ b/src/metrics/qkerr_schild.h
@@ -39,6 +39,7 @@ namespace metric {
     const real_t chi_min, eta_min, phi_min;
     const real_t dchi, deta, dphi;
     const real_t dchi_inv, deta_inv, dphi_inv;
+    const bool   small_angle;
 
     Inline auto Delta(real_t r) const -> real_t {
       return SQR(r) - TWO * r + SQR(a);
@@ -89,7 +90,8 @@ namespace metric {
       , dphi { (x3_max - phi_min) / nx3 }
       , dchi_inv { ONE / dchi }
       , deta_inv { ONE / deta }
-      , dphi_inv { ONE / dphi } {
+      , dphi_inv { ONE / dphi } 
+      , small_angle { eta2theta(HALF * deta) < constant::SMALL_ANGLE } {
       set_dxMin(find_dxMin());
     }
 
@@ -503,15 +505,27 @@ namespace metric {
      * differential area at the pole (used in axisymmetric solvers)
      * @note approximate solution for the polar area
      * @param x1 radial coordinate along the axis (code units)
+     * @note uses small-angle approximation when the resolution is too high
      */
     Inline auto polar_area(real_t x1) const -> real_t {
       if constexpr (D != Dim::_1D) {
-        return dchi * math::exp(x1 * dchi + chi_min) *
-               (SQR(r0 + math::exp(x1 * dchi + chi_min)) + SQR(a)) *
-               math::sqrt(
-                 ONE + TWO * (r0 + math::exp(x1 * dchi + chi_min)) /
-                         (SQR(r0 + math::exp(x1 * dchi + chi_min)) + SQR(a))) *
-               (ONE - math::cos(eta2theta(HALF * deta)));
+        if (small_angle) { 
+          const real_t dtheta = eta2theta(HALF * deta);
+          return dchi * math::exp(x1 * dchi + chi_min) *
+                (SQR(r0 + math::exp(x1 * dchi + chi_min)) + SQR(a)) *
+                math::sqrt(
+                  ONE + TWO * (r0 + math::exp(x1 * dchi + chi_min)) /
+                          (SQR(r0 + math::exp(x1 * dchi + chi_min)) + SQR(a))) *
+                (static_cast<real_t>(48) - SQR(dtheta)) * SQR(dtheta) /
+                 static_cast<real_t>(384);
+        } else {
+          return dchi * math::exp(x1 * dchi + chi_min) *
+                (SQR(r0 + math::exp(x1 * dchi + chi_min)) + SQR(a)) *
+                math::sqrt(
+                  ONE + TWO * (r0 + math::exp(x1 * dchi + chi_min)) /
+                          (SQR(r0 + math::exp(x1 * dchi + chi_min)) + SQR(a))) *
+                (ONE - math::cos(eta2theta(HALF * deta)));
+        }
       }
     }
 
diff --git a/src/output/tests/writer-mpi.cpp b/src/output/tests/writer-mpi.cpp
index dad981a4..382f880c 100644
--- a/src/output/tests/writer-mpi.cpp
+++ b/src/output/tests/writer-mpi.cpp
@@ -7,7 +7,7 @@
 
 #include <Kokkos_Core.hpp>
 #include <adios2.h>
-#include <adios2/cxx11/KokkosView.h>
+#include <adios2/cxx/KokkosView.h>
 #include <mpi.h>
 
 #include <filesystem>
diff --git a/src/output/tests/writer-nompi.cpp b/src/output/tests/writer-nompi.cpp
index a185b177..b32b1058 100644
--- a/src/output/tests/writer-nompi.cpp
+++ b/src/output/tests/writer-nompi.cpp
@@ -7,7 +7,7 @@
 
 #include <Kokkos_Core.hpp>
 #include <adios2.h>
-#include <adios2/cxx11/KokkosView.h>
+#include <adios2/cxx/KokkosView.h>
 
 #include <filesystem>
 #include <iostream>
diff --git a/src/output/utils/attr_writer.h b/src/output/utils/attr_writer.h
index 41b64fb7..4e0c0c87 100644
--- a/src/output/utils/attr_writer.h
+++ b/src/output/utils/attr_writer.h
@@ -12,7 +12,7 @@
 #define OUTPUT_UTILS_ATTR_WRITER_H
 
 #include <adios2.h>
-#include <adios2/cxx11/KokkosView.h>
+#include <adios2/cxx/KokkosView.h>
 
 #include <any>
 #include <functional>
diff --git a/src/output/writer.cpp b/src/output/writer.cpp
index 64d0f081..59c560f9 100644
--- a/src/output/writer.cpp
+++ b/src/output/writer.cpp
@@ -10,7 +10,7 @@
 
 #include <Kokkos_Core.hpp>
 #include <adios2.h>
-#include <adios2/cxx11/KokkosView.h>
+#include <adios2/cxx/KokkosView.h>
 
 #if defined(MPI_ENABLED)
   #include "arch/mpi_aliases.h"
diff --git a/src/output/writer.h b/src/output/writer.h
index 4fe0c194..190cbb8e 100644
--- a/src/output/writer.h
+++ b/src/output/writer.h
@@ -17,7 +17,7 @@
 #include "output/spectra.h"
 
 #include <adios2.h>
-#include <adios2/cxx11/KokkosView.h>
+#include <adios2/cxx/KokkosView.h>
 
 #if defined(MPI_ENABLED)
   #include <mpi.h>