diff --git a/input.example.toml b/input.example.toml
index 4a671fab..26a5d1c9 100644
--- a/input.example.toml
+++ b/input.example.toml
@@ -453,6 +453,10 @@
     #   @default: ["B^2", "E^2", "ExB", "Rho", "T00"]
     #   @note: Same notation as for `output.fields.quantities`
     quantities = ""
+    # Custom (user-defined) stats:
+    #   @type: array of strings
+    #   @default: []
+    custom = ""
 
 [checkpoint]
   # Number of timesteps between checkpoints:
diff --git a/pgens/reconnection/reconnection.toml b/pgens/reconnection/reconnection.toml
index 62c1ffe8..76acbf08 100644
--- a/pgens/reconnection/reconnection.toml
+++ b/pgens/reconnection/reconnection.toml
@@ -3,29 +3,35 @@
   engine  = "srpic"
   runtime = 10.0
 
+  [simulation.domain]
+    decomposition = [-1, 2]
+
 [grid]
-  resolution = [1024, 512]
-  extent     = [[-1.0, 1.0], [-0.5, 0.5]]
+  resolution = [512, 512]
+  extent     = [[-1.0, 1.0], [-1.0, 1.0]]
 
   [grid.metric]
     metric = "minkowski"
 
   [grid.boundaries]
-    fields    = [["PERIODIC"], ["MATCH", "MATCH"], ["PERIODIC"]]
-    particles = [["PERIODIC"], ["ABSORB", "ABSORB"], ["PERIODIC"]]
+    fields    = [["PERIODIC"], ["MATCH", "MATCH"]]
+    particles = [["PERIODIC"], ["ABSORB", "ABSORB"]]
+
+    [grid.boundaries.match]
+      ds = [[0.04], [0.1]]
 
 [scales]
   larmor0    = 2e-4
   skindepth0 = 2e-3
 
 [algorithms]
-  current_filters = 4
+  current_filters = 8
 
   [algorithms.timestep]
     CFL = 0.5
 
 [particles]
-  ppc0 = 2.0
+  ppc0 = 8.0
 
   [[particles.species]]
     label    = "e-"
@@ -40,10 +46,12 @@
     maxnpart = 1e7
 
 [setup]
-  Bmag        = 1.0
-  width       = 0.01
-  bg_temp     = 1e-4
-  overdensity = 3.0
+  bg_B           = 1.0
+  bg_Bguide      = 0.0
+  bg_temperature = 1e-4
+  inj_ypad       = 0.25
+  cs_width       = 0.05
+  cs_overdensity = 3.0
 
 [output]
   format        = "hdf5"
@@ -51,3 +59,13 @@
 
   [output.fields]
     quantities = ["N_1", "N_2", "E", "B", "J"]
+
+  [output.particles]
+    enable = false
+
+  [output.spectra]
+    enable = false
+
+[diagnostics]
+  colored_stdout = true
+  interval       = 10
diff --git a/src/engines/engine.hpp b/src/engines/engine.hpp
index 9525874a..5956ea5d 100644
--- a/src/engines/engine.hpp
+++ b/src/engines/engine.hpp
@@ -23,7 +23,6 @@
 
 #include "arch/traits.h"
 #include "utils/error.h"
-#include "utils/progressbar.h"
 #include "utils/timer.h"
 #include "utils/toml.h"
 
diff --git a/src/engines/engine_run.cpp b/src/engines/engine_run.cpp
index 472acae4..65c87e93 100644
--- a/src/engines/engine_run.cpp
+++ b/src/engines/engine_run.cpp
@@ -72,23 +72,41 @@ namespace ntt {
           auto lambda_custom_field_output = [&](const std::string&    name,
                                                 ndfield_t<M::Dim, 6>& buff,
                                                 index_t               idx,
+                                                timestep_t            step,
+                                                simtime_t             time,
                                                 const Domain<S, M>&   dom) {
-            m_pgen.CustomFieldOutput(name, buff, idx, dom);
+            m_pgen.CustomFieldOutput(name, buff, idx, step, time, dom);
           };
-          print_output = m_metadomain.Write(m_params,
-                                            step,
-                                            step - 1,
-                                            time,
-                                            time - dt,
-                                            lambda_custom_field_output);
+          print_output &= m_metadomain.Write(m_params,
+                                             step,
+                                             step - 1,
+                                             time,
+                                             time - dt,
+                                             lambda_custom_field_output);
         } else {
-          print_output = m_metadomain.Write(m_params, step, step - 1, time, time - dt);
+          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) {
+          auto lambda_custom_stat = [&](const std::string&  name,
+                                        timestep_t          step,
+                                        simtime_t           time,
+                                        const Domain<S, M>& dom) -> real_t {
+            return m_pgen.CustomStat(name, step, time, dom);
+          };
+          print_output &= m_metadomain.WriteStats(m_params,
+                                                  step,
+                                                  step - 1,
+                                                  time,
+                                                  time - dt,
+                                                  lambda_custom_stat);
+        } else {
+          print_output &= m_metadomain.WriteStats(m_params,
+                                                  step,
+                                                  step - 1,
+                                                  time,
+                                                  time - dt);
         }
-        print_output &= m_metadomain.WriteStats(m_params,
-                                                step,
-                                                step - 1,
-                                                time,
-                                                time - dt);
         timers.stop("Output");
 
         timers.start("Checkpoint");
@@ -133,4 +151,5 @@ namespace ntt {
   template void Engine<SimEngine::GRPIC, metric::KerrSchild<Dim::_2D>>::run();
   template void Engine<SimEngine::GRPIC, metric::KerrSchild0<Dim::_2D>>::run();
   template void Engine<SimEngine::GRPIC, metric::QKerrSchild<Dim::_2D>>::run();
+
 } // namespace ntt
diff --git a/src/framework/domain/checkpoint.cpp b/src/framework/domain/checkpoint.cpp
index 656ff57d..a39f18e6 100644
--- a/src/framework/domain/checkpoint.cpp
+++ b/src/framework/domain/checkpoint.cpp
@@ -474,13 +474,26 @@ namespace ntt {
       HERE);
   }
 
-  template struct Metadomain<SimEngine::SRPIC, metric::Minkowski<Dim::_1D>>;
-  template struct Metadomain<SimEngine::SRPIC, metric::Minkowski<Dim::_2D>>;
-  template struct Metadomain<SimEngine::SRPIC, metric::Minkowski<Dim::_3D>>;
-  template struct Metadomain<SimEngine::SRPIC, metric::Spherical<Dim::_2D>>;
-  template struct Metadomain<SimEngine::SRPIC, metric::QSpherical<Dim::_2D>>;
-  template struct Metadomain<SimEngine::GRPIC, metric::KerrSchild<Dim::_2D>>;
-  template struct Metadomain<SimEngine::GRPIC, metric::QKerrSchild<Dim::_2D>>;
-  template struct Metadomain<SimEngine::GRPIC, metric::KerrSchild0<Dim::_2D>>;
+#define METADOMAIN_CHECKPOINTS(S, M)                                             \
+  template void Metadomain<S, M>::InitCheckpointWriter(adios2::ADIOS*,           \
+                                                       const SimulationParams&); \
+  template auto Metadomain<S, M>::WriteCheckpoint(const SimulationParams&,       \
+                                                  timestep_t,                    \
+                                                  timestep_t,                    \
+                                                  simtime_t,                     \
+                                                  simtime_t) -> bool;            \
+  template void Metadomain<S, M>::ContinueFromCheckpoint(adios2::ADIOS*,         \
+                                                         const SimulationParams&);
+
+  METADOMAIN_CHECKPOINTS(SimEngine::SRPIC, metric::Minkowski<Dim::_1D>)
+  METADOMAIN_CHECKPOINTS(SimEngine::SRPIC, metric::Minkowski<Dim::_2D>)
+  METADOMAIN_CHECKPOINTS(SimEngine::SRPIC, metric::Minkowski<Dim::_3D>)
+  METADOMAIN_CHECKPOINTS(SimEngine::SRPIC, metric::Spherical<Dim::_2D>)
+  METADOMAIN_CHECKPOINTS(SimEngine::SRPIC, metric::QSpherical<Dim::_2D>)
+  METADOMAIN_CHECKPOINTS(SimEngine::GRPIC, metric::KerrSchild<Dim::_2D>)
+  METADOMAIN_CHECKPOINTS(SimEngine::GRPIC, metric::QKerrSchild<Dim::_2D>)
+  METADOMAIN_CHECKPOINTS(SimEngine::GRPIC, metric::KerrSchild0<Dim::_2D>)
+
+#undef METADOMAIN_CHECKPOINTS
 
 } // namespace ntt
diff --git a/src/framework/domain/communications.cpp b/src/framework/domain/communications.cpp
index 841c9a7d..a538f853 100644
--- a/src/framework/domain/communications.cpp
+++ b/src/framework/domain/communications.cpp
@@ -36,10 +36,10 @@ namespace ntt {
   using comm_params_t = std::pair<address_t, std::vector<range_tuple_t>>;
 
   template <SimEngine::type S, class M>
-  auto GetSendRecvRanks(
-    Metadomain<S, M>*        metadomain,
-    Domain<S, M>&            domain,
-    dir::direction_t<M::Dim> direction) -> std::pair<address_t, address_t> {
+  auto GetSendRecvRanks(Metadomain<S, M>*        metadomain,
+                        Domain<S, M>&            domain,
+                        dir::direction_t<M::Dim> direction)
+    -> std::pair<address_t, address_t> {
     Domain<S, M>* send_to_nghbr_ptr   = nullptr;
     Domain<S, M>* recv_from_nghbr_ptr = nullptr;
     // set pointers to the correct send/recv domains
@@ -119,11 +119,11 @@ namespace ntt {
   }
 
   template <SimEngine::type S, class M>
-  auto GetSendRecvParams(
-    Metadomain<S, M>*        metadomain,
-    Domain<S, M>&            domain,
-    dir::direction_t<M::Dim> direction,
-    bool synchronize) -> std::pair<comm_params_t, comm_params_t> {
+  auto GetSendRecvParams(Metadomain<S, M>*        metadomain,
+                         Domain<S, M>&            domain,
+                         dir::direction_t<M::Dim> direction,
+                         bool                     synchronize)
+    -> std::pair<comm_params_t, comm_params_t> {
     const auto [send_indrank,
                 recv_indrank] = GetSendRecvRanks(metadomain, domain, direction);
     const auto [send_ind, send_rank] = send_indrank;
@@ -641,13 +641,23 @@ namespace ntt {
     }
   }
 
-  template struct Metadomain<SimEngine::SRPIC, metric::Minkowski<Dim::_1D>>;
-  template struct Metadomain<SimEngine::SRPIC, metric::Minkowski<Dim::_2D>>;
-  template struct Metadomain<SimEngine::SRPIC, metric::Minkowski<Dim::_3D>>;
-  template struct Metadomain<SimEngine::SRPIC, metric::Spherical<Dim::_2D>>;
-  template struct Metadomain<SimEngine::SRPIC, metric::QSpherical<Dim::_2D>>;
-  template struct Metadomain<SimEngine::GRPIC, metric::KerrSchild<Dim::_2D>>;
-  template struct Metadomain<SimEngine::GRPIC, metric::QKerrSchild<Dim::_2D>>;
-  template struct Metadomain<SimEngine::GRPIC, metric::KerrSchild0<Dim::_2D>>;
+#define METADOMAIN_COMM(S, M)                                                  \
+  template void Metadomain<S, M>::CommunicateFields(Domain<S, M>&, CommTags);  \
+  template void Metadomain<S, M>::SynchronizeFields(Domain<S, M>&,             \
+                                                    CommTags,                  \
+                                                    const range_tuple_t&);     \
+  template void Metadomain<S, M>::CommunicateParticles(Domain<S, M>&);         \
+  template void Metadomain<S, M>::RemoveDeadParticles(Domain<S, M>&);
+
+  METADOMAIN_COMM(SimEngine::SRPIC, metric::Minkowski<Dim::_1D>)
+  METADOMAIN_COMM(SimEngine::SRPIC, metric::Minkowski<Dim::_2D>)
+  METADOMAIN_COMM(SimEngine::SRPIC, metric::Minkowski<Dim::_3D>)
+  METADOMAIN_COMM(SimEngine::SRPIC, metric::Spherical<Dim::_2D>)
+  METADOMAIN_COMM(SimEngine::SRPIC, metric::QSpherical<Dim::_2D>)
+  METADOMAIN_COMM(SimEngine::GRPIC, metric::KerrSchild<Dim::_2D>)
+  METADOMAIN_COMM(SimEngine::GRPIC, metric::QKerrSchild<Dim::_2D>)
+  METADOMAIN_COMM(SimEngine::GRPIC, metric::KerrSchild0<Dim::_2D>)
+
+#undef METADOMAIN_COMM
 
 } // namespace ntt
diff --git a/src/framework/domain/grid.cpp b/src/framework/domain/grid.cpp
index ddbfc38f..f087df6f 100644
--- a/src/framework/domain/grid.cpp
+++ b/src/framework/domain/grid.cpp
@@ -86,8 +86,8 @@ namespace ntt {
   }
 
   template <Dimension D>
-  auto Grid<D>::rangeCellsOnHost(const box_region_t<D>& region) const
-    -> range_h_t<D> {
+  auto Grid<D>::rangeCellsOnHost(
+    const box_region_t<D>& region) const -> range_h_t<D> {
     tuple_t<ncells_t, D> imin, imax;
     for (auto i { 0u }; i < D; i++) {
       switch (region[i]) {
@@ -163,8 +163,8 @@ namespace ntt {
   }
 
   template <Dimension D>
-  auto Grid<D>::rangeCells(const tuple_t<list_t<int, 2>, D>& ranges) const
-    -> range_t<D> {
+  auto Grid<D>::rangeCells(
+    const tuple_t<list_t<int, 2>, D>& ranges) const -> range_t<D> {
     tuple_t<ncells_t, D> imin, imax;
     for (auto i { 0u }; i < D; i++) {
       raise::ErrorIf((ranges[i][0] < -(int)N_GHOSTS) ||
diff --git a/src/framework/domain/grid.h b/src/framework/domain/grid.h
index af5b6c8d..87b21b1f 100644
--- a/src/framework/domain/grid.h
+++ b/src/framework/domain/grid.h
@@ -130,6 +130,15 @@ namespace ntt {
       return m_resolution;
     }
 
+    [[nodiscard]]
+    auto num_active() const -> ncells_t {
+      ncells_t total_active = 1u;
+      for (const auto& res : m_resolution) {
+        total_active *= res;
+      }
+      return total_active;
+    }
+
     [[nodiscard]]
     auto n_all(in i) const -> ncells_t {
       switch (i) {
@@ -154,6 +163,15 @@ namespace ntt {
       return nall;
     }
 
+    [[nodiscard]]
+    auto num_all() const -> ncells_t {
+      ncells_t total_all = 1u;
+      for (const auto& res : n_all()) {
+        total_all *= res;
+      }
+      return total_all;
+    }
+
     /* Ranges in the device execution space --------------------------------- */
     /**
      * @brief Loop over all active cells (disregard ghost cells)
diff --git a/src/framework/domain/metadomain.cpp b/src/framework/domain/metadomain.cpp
index 4b9057e2..a0326913 100644
--- a/src/framework/domain/metadomain.cpp
+++ b/src/framework/domain/metadomain.cpp
@@ -23,7 +23,6 @@
 
 #include <limits>
 #include <map>
-#include <memory>
 #include <string>
 #include <vector>
 
diff --git a/src/framework/domain/metadomain.h b/src/framework/domain/metadomain.h
index 1fb6ce00..0cb67096 100644
--- a/src/framework/domain/metadomain.h
+++ b/src/framework/domain/metadomain.h
@@ -19,7 +19,6 @@
 #include "global.h"
 
 #include "arch/kokkos_aliases.h"
-#include "utils/timer.h"
 
 #include "framework/containers/species.h"
 #include "framework/domain/domain.h"
@@ -125,8 +124,10 @@ namespace ntt {
                simtime_t,
                std::function<void(const std::string&,
                                   ndfield_t<M::Dim, 6>&,
-                                  std::size_t,
-                                  const Domain<S, M>&)> = {}) -> bool;
+                                  index_t,
+                                  timestep_t,
+                                  simtime_t,
+                                  const Domain<S, M>&)> = nullptr) -> bool;
     void InitCheckpointWriter(adios2::ADIOS*, const SimulationParams&);
     auto WriteCheckpoint(const SimulationParams&,
                          timestep_t,
@@ -138,8 +139,14 @@ namespace ntt {
 #endif
 
     void InitStatsWriter(const SimulationParams&, bool);
-    auto WriteStats(const SimulationParams&, timestep_t, timestep_t, simtime_t, simtime_t)
-      -> bool;
+    auto WriteStats(
+      const SimulationParams&,
+      timestep_t,
+      timestep_t,
+      simtime_t,
+      simtime_t,
+      std::function<real_t(const std::string&, timestep_t, simtime_t, const Domain<S, M>&)> =
+        nullptr) -> bool;
 
     /* setters -------------------------------------------------------------- */
     void setFldsBC(const bc_in&, const FldsBC&);
diff --git a/src/framework/domain/output.cpp b/src/framework/domain/output.cpp
index 6903c719..ce5a1316 100644
--- a/src/framework/domain/output.cpp
+++ b/src/framework/domain/output.cpp
@@ -188,14 +188,17 @@ namespace ntt {
 
   template <SimEngine::type S, class M>
   auto Metadomain<S, M>::Write(
-    const SimulationParams& params,
-    timestep_t              current_step,
-    timestep_t              finished_step,
-    simtime_t               current_time,
-    simtime_t               finished_time,
-    std::function<
-      void(const std::string&, ndfield_t<M::Dim, 6>&, std::size_t, const Domain<S, M>&)>
-      CustomFieldOutput) -> bool {
+    const SimulationParams&                  params,
+    timestep_t                               current_step,
+    timestep_t                               finished_step,
+    simtime_t                                current_time,
+    simtime_t                                finished_time,
+    std::function<void(const std::string&,
+                       ndfield_t<M::Dim, 6>&,
+                       index_t,
+                       timestep_t,
+                       simtime_t,
+                       const Domain<S, M>&)> CustomFieldOutput) -> bool {
     raise::ErrorIf(
       l_subdomain_indices().size() != 1,
       "Output for now is only supported for one subdomain per rank",
@@ -382,6 +385,8 @@ namespace ntt {
               CustomFieldOutput(fld.name().substr(1),
                                 local_domain->fields.bckp,
                                 addresses.back(),
+                                finished_step,
+                                finished_time,
                                 *local_domain);
             } else {
               raise::Error("Custom output requested but no function provided",
@@ -708,13 +713,32 @@ namespace ntt {
     return true;
   }
 
-  template struct Metadomain<SimEngine::SRPIC, metric::Minkowski<Dim::_1D>>;
-  template struct Metadomain<SimEngine::SRPIC, metric::Minkowski<Dim::_2D>>;
-  template struct Metadomain<SimEngine::SRPIC, metric::Minkowski<Dim::_3D>>;
-  template struct Metadomain<SimEngine::SRPIC, metric::Spherical<Dim::_2D>>;
-  template struct Metadomain<SimEngine::SRPIC, metric::QSpherical<Dim::_2D>>;
-  template struct Metadomain<SimEngine::GRPIC, metric::KerrSchild<Dim::_2D>>;
-  template struct Metadomain<SimEngine::GRPIC, metric::QKerrSchild<Dim::_2D>>;
-  template struct Metadomain<SimEngine::GRPIC, metric::KerrSchild0<Dim::_2D>>;
+#define METADOMAIN_OUTPUT(S, M)                                                \
+  template void Metadomain<S, M>::InitWriter(adios2::ADIOS*,                   \
+                                             const SimulationParams&,          \
+                                             bool);                            \
+  template auto Metadomain<S, M>::Write(                                       \
+    const SimulationParams&,                                                   \
+    timestep_t,                                                                \
+    timestep_t,                                                                \
+    simtime_t,                                                                 \
+    simtime_t,                                                                 \
+    std::function<void(const std::string&,                                     \
+                       ndfield_t<M::Dim, 6>&,                                  \
+                       index_t,                                                \
+                       timestep_t,                                             \
+                       simtime_t,                                              \
+                       const Domain<S, M>&)>) -> bool;
+
+  METADOMAIN_OUTPUT(SimEngine::SRPIC, metric::Minkowski<Dim::_1D>)
+  METADOMAIN_OUTPUT(SimEngine::SRPIC, metric::Minkowski<Dim::_2D>)
+  METADOMAIN_OUTPUT(SimEngine::SRPIC, metric::Minkowski<Dim::_3D>)
+  METADOMAIN_OUTPUT(SimEngine::SRPIC, metric::Spherical<Dim::_2D>)
+  METADOMAIN_OUTPUT(SimEngine::SRPIC, metric::QSpherical<Dim::_2D>)
+  METADOMAIN_OUTPUT(SimEngine::GRPIC, metric::KerrSchild<Dim::_2D>)
+  METADOMAIN_OUTPUT(SimEngine::GRPIC, metric::QKerrSchild<Dim::_2D>)
+  METADOMAIN_OUTPUT(SimEngine::GRPIC, metric::KerrSchild0<Dim::_2D>)
+
+#undef METADOMAIN_OUTPUT
 
 } // namespace ntt
diff --git a/src/framework/domain/stats.cpp b/src/framework/domain/stats.cpp
index 60acd64a..a252e975 100644
--- a/src/framework/domain/stats.cpp
+++ b/src/framework/domain/stats.cpp
@@ -1,7 +1,7 @@
 #include "enums.h"
 #include "global.h"
 
-#include "arch/kokkos_aliases.h"
+#include "utils/comparators.h"
 #include "utils/error.h"
 #include "utils/log.h"
 #include "utils/numeric.h"
@@ -18,6 +18,8 @@
 #include "framework/domain/metadomain.h"
 #include "framework/parameters.h"
 
+#include "kernels/reduced_stats.hpp"
+
 #include <Kokkos_Core.hpp>
 #include <Kokkos_ScatterView.hpp>
 #include <Kokkos_StdAlgorithms.hpp>
@@ -51,23 +53,26 @@ namespace ntt {
     }
     const auto stats_to_write = params.template get<std::vector<std::string>>(
       "output.stats.quantities");
+    const auto custom_stats_to_write = params.template get<std::vector<std::string>>(
+      "output.stats.custom");
     g_stats_writer.init(
       params.template get<timestep_t>("output.stats.interval"),
       params.template get<simtime_t>("output.stats.interval_time"));
     g_stats_writer.defineStatsFilename(filename);
-    g_stats_writer.defineStatsOutputs(stats_to_write);
+    g_stats_writer.defineStatsOutputs(stats_to_write, false);
+    g_stats_writer.defineStatsOutputs(custom_stats_to_write, true);
 
     if (not std::filesystem::exists(filename)) {
       g_stats_writer.writeHeader();
     }
   }
 
-  template <SimEngine::type S, class M, typename T, StatsID::type F>
+  template <SimEngine::type S, class M, StatsID::type P>
   auto ComputeMoments(const SimulationParams& params,
                       const Mesh<M>&          mesh,
                       const std::vector<Particles<M::Dim, M::CoordType>>& prtl_species,
                       const std::vector<spidx_t>&        species,
-                      const std::vector<unsigned short>& components) -> T {
+                      const std::vector<unsigned short>& components) -> real_t {
     std::vector<spidx_t> specs = species;
     if (specs.size() == 0) {
       // if no species specified, take all massive species
@@ -86,47 +91,103 @@ namespace ntt {
     const auto use_weights = params.template get<bool>("particles.use_weights");
     const auto inv_n0      = ONE / params.template get<real_t>("scales.n0");
 
-    T buffer = static_cast<T>(0);
+    real_t buffer = static_cast<real_t>(0);
     for (const auto& sp : specs) {
       auto& prtl_spec = prtl_species[sp - 1];
-      // Kokkos::parallel_reduce(
-      //   "ComputeMoments",
-      //   prtl_spec.rangeActiveParticles(),
-      //   // clang-format off
-      //   kernel::ReducedParticleMoments_kernel<S, M, T, F>(components,
-      //                                                     prtl_spec.i1, prtl_spec.i2, prtl_spec.i3,
-      //                                                     prtl_spec.dx1, prtl_spec.dx2, prtl_spec.dx3,
-      //                                                     prtl_spec.ux1, prtl_spec.ux2, prtl_spec.ux3,
-      //                                                     prtl_spec.phi, prtl_spec.weight, prtl_spec.tag,
-      //                                                     prtl_spec.mass(), prtl_spec.charge(),
-      //                                                     use_weights, mesh.metric, mesh.flds_bc(), inv_n0),
-      //   // clang-format on
-      //   buffer);
+      if (P == StatsID::Charge and cmp::AlmostZero_host(prtl_spec.charge())) {
+        continue;
+      }
+      if (P == StatsID::Rho and cmp::AlmostZero_host(prtl_spec.mass())) {
+        continue;
+      }
+      Kokkos::parallel_reduce(
+        "ComputeMoments",
+        prtl_spec.rangeActiveParticles(),
+        // clang-format off
+        kernel::ReducedParticleMoments_kernel<S, M, P>(components,
+                                                       prtl_spec.i1, prtl_spec.i2, prtl_spec.i3,
+                                                       prtl_spec.dx1, prtl_spec.dx2, prtl_spec.dx3,
+                                                       prtl_spec.ux1, prtl_spec.ux2, prtl_spec.ux3,
+                                                       prtl_spec.phi, prtl_spec.weight, prtl_spec.tag,
+                                                       prtl_spec.mass(), prtl_spec.charge(),
+                                                       use_weights, mesh.metric),
+        // clang-format on
+        buffer);
     }
     return buffer;
   }
 
   template <SimEngine::type S, class M, StatsID::type F>
-  auto ComputeFields(Domain<S, M>*                      domain,
-                     const std::vector<unsigned short>& components) -> real_t {
+  auto ReduceFields(Domain<S, M>*                      domain,
+                    const M&                           global_metric,
+                    const std::vector<unsigned short>& components) -> real_t {
     auto buffer { ZERO };
-    // Kokkos::parallel_reduce(
-    //   "ComputeMoments",
-    //   prtl_spec.rangeActiveParticles(),
-    //   kernel::ReducedFields_kernel<S, M, F>(components,
-    //                                         domain->fields.em,
-    //                                         domain->fields.cur,
-    //                                         domain->mesh.metric),
-    //   buffer);
-    return buffer;
+    if constexpr (F == StatsID::JdotE) {
+      if (components.size() == 0) {
+        Kokkos::parallel_reduce(
+          "ReduceFields",
+          domain->mesh.rangeActiveCells(),
+          kernel::ReducedFields_kernel<S, M, F, 0>(domain->fields.em,
+                                                   domain->fields.cur,
+                                                   domain->mesh.metric),
+          buffer);
+      } else {
+        raise::Error("Components not supported for JdotE", HERE);
+      }
+    } else if constexpr (
+      (S == SimEngine::SRPIC) and
+      (F == StatsID::B2 or F == StatsID::E2 or F == StatsID::ExB)) {
+      raise::ErrorIf(components.size() != 1,
+                     "Components must be of size 1 for B2, E2 or ExB stats",
+                     HERE);
+      const auto comp = components[0];
+      if (comp == 1) {
+        Kokkos::parallel_reduce(
+          "ReduceFields",
+          domain->mesh.rangeActiveCells(),
+          kernel::ReducedFields_kernel<S, M, F, 1>(domain->fields.em,
+                                                   domain->fields.cur,
+                                                   domain->mesh.metric),
+          buffer);
+      } else if (comp == 2) {
+        Kokkos::parallel_reduce(
+          "ReduceFields",
+          domain->mesh.rangeActiveCells(),
+          kernel::ReducedFields_kernel<S, M, F, 2>(domain->fields.em,
+                                                   domain->fields.cur,
+                                                   domain->mesh.metric),
+          buffer);
+      } else if (comp == 3) {
+        Kokkos::parallel_reduce(
+          "ReduceFields",
+          domain->mesh.rangeActiveCells(),
+          kernel::ReducedFields_kernel<S, M, F, 3>(domain->fields.em,
+                                                   domain->fields.cur,
+                                                   domain->mesh.metric),
+          buffer);
+      } else {
+        raise::Error(
+          "Invalid component for B2, E2 or ExB stats: " + std::to_string(comp),
+          HERE);
+      }
+    } else {
+      raise::Error("ReduceFields not implemented for this stats ID + SimEngine "
+                   "combination",
+                   HERE);
+    }
+
+    return buffer / global_metric.totVolume();
   }
 
   template <SimEngine::type S, class M>
-  auto Metadomain<S, M>::WriteStats(const SimulationParams& params,
-                                    timestep_t              current_step,
-                                    timestep_t              finished_step,
-                                    simtime_t               current_time,
-                                    simtime_t finished_time) -> bool {
+  auto Metadomain<S, M>::WriteStats(
+    const SimulationParams& params,
+    timestep_t              current_step,
+    timestep_t              finished_step,
+    simtime_t               current_time,
+    simtime_t               finished_time,
+    std::function<real_t(const std::string&, timestep_t, simtime_t, const Domain<S, M>&)> CustomStat)
+    -> bool {
     if (not(params.template get<bool>("output.stats.enable") and
             g_stats_writer.shouldWrite(finished_step, finished_time))) {
       return false;
@@ -136,69 +197,101 @@ namespace ntt {
     g_stats_writer.write(current_step);
     g_stats_writer.write(current_time);
     for (const auto& stat : g_stats_writer.statsWriters()) {
-      if (stat.id() == StatsID::N) {
-        g_stats_writer.write(
-          ComputeMoments<S, M, real_t, StatsID::N>(params,
-                                                   local_domain->mesh,
-                                                   local_domain->species,
-                                                   stat.species,
-                                                   {}));
+      if (stat.id() == StatsID::Custom) {
+        if (CustomStat != nullptr) {
+          g_stats_writer.write(
+            CustomStat(stat.name(), finished_step, finished_time, *local_domain));
+        } else {
+          raise::Error("Custom output requested but no function provided", HERE);
+        }
+      } else if (stat.id() == StatsID::N) {
+        g_stats_writer.write(ComputeMoments<S, M, StatsID::N>(params,
+                                                              local_domain->mesh,
+                                                              local_domain->species,
+                                                              stat.species,
+                                                              {}));
       } else if (stat.id() == StatsID::Npart) {
         g_stats_writer.write(
-          ComputeMoments<S, M, npart_t, StatsID::Npart>(params,
-                                                        local_domain->mesh,
-                                                        local_domain->species,
-                                                        stat.species,
-                                                        {}));
+          ComputeMoments<S, M, StatsID::Npart>(params,
+                                               local_domain->mesh,
+                                               local_domain->species,
+                                               stat.species,
+                                               {}));
       } else if (stat.id() == StatsID::Rho) {
         g_stats_writer.write(
-          ComputeMoments<S, M, real_t, StatsID::Rho>(params,
-                                                     local_domain->mesh,
-                                                     local_domain->species,
-                                                     stat.species,
-                                                     {}));
+          ComputeMoments<S, M, StatsID::Rho>(params,
+                                             local_domain->mesh,
+                                             local_domain->species,
+                                             stat.species,
+                                             {}));
       } else if (stat.id() == StatsID::Charge) {
         g_stats_writer.write(
-          ComputeMoments<S, M, real_t, StatsID::Charge>(params,
-                                                        local_domain->mesh,
-                                                        local_domain->species,
-                                                        stat.species,
-                                                        {}));
+          ComputeMoments<S, M, StatsID::Charge>(params,
+                                                local_domain->mesh,
+                                                local_domain->species,
+                                                stat.species,
+                                                {}));
       } else if (stat.id() == StatsID::T) {
         for (const auto& comp : stat.comp) {
           g_stats_writer.write(
-            ComputeMoments<S, M, real_t, StatsID::T>(params,
-                                                     local_domain->mesh,
-                                                     local_domain->species,
-                                                     stat.species,
-                                                     comp));
+            ComputeMoments<S, M, StatsID::T>(params,
+                                             local_domain->mesh,
+                                             local_domain->species,
+                                             stat.species,
+                                             comp));
         }
       } else if (stat.id() == StatsID::JdotE) {
-        g_stats_writer.write(ComputeFields<S, M, StatsID::JdotE>(local_domain, {}));
-      } else if (stat.id() == StatsID::E2) {
-        g_stats_writer.write(ComputeFields<S, M, StatsID::E2>(local_domain, {}));
-      } else if (stat.id() == StatsID::B2) {
-        g_stats_writer.write(ComputeFields<S, M, StatsID::B2>(local_domain, {}));
-      } else if (stat.id() == StatsID::ExB) {
-        for (const auto& comp : stat.comp) {
-          g_stats_writer.write(
-            ComputeFields<S, M, StatsID::ExB>(local_domain, comp));
+        g_stats_writer.write(
+          ReduceFields<S, M, StatsID::JdotE>(local_domain, g_mesh.metric, {}));
+      } else if (S == SimEngine::SRPIC) {
+        if (stat.id() == StatsID::E2) {
+          for (const auto& comp : stat.comp) {
+            g_stats_writer.write(
+              ReduceFields<S, M, StatsID::E2>(local_domain, g_mesh.metric, comp));
+          }
+        } else if (stat.id() == StatsID::B2) {
+          for (const auto& comp : stat.comp) {
+            g_stats_writer.write(
+              ReduceFields<S, M, StatsID::B2>(local_domain, g_mesh.metric, comp));
+          }
+        } else if (stat.id() == StatsID::ExB) {
+          for (const auto& comp : stat.comp) {
+            g_stats_writer.write(
+              ReduceFields<S, M, StatsID::ExB>(local_domain, g_mesh.metric, comp));
+          }
+        } else {
+          raise::Error("Unrecognized stats ID " + stat.name(), HERE);
         }
       } else {
-        raise::Error("Unrecognized stats ID " + stat.name(), HERE);
+        raise::Error("StatsID not implemented for particular SimEngine: " +
+                       std::to_string(static_cast<int>(S)),
+                     HERE);
       }
     }
     g_stats_writer.endWriting();
     return true;
   }
 
-  template struct Metadomain<SimEngine::SRPIC, metric::Minkowski<Dim::_1D>>;
-  template struct Metadomain<SimEngine::SRPIC, metric::Minkowski<Dim::_2D>>;
-  template struct Metadomain<SimEngine::SRPIC, metric::Minkowski<Dim::_3D>>;
-  template struct Metadomain<SimEngine::SRPIC, metric::Spherical<Dim::_2D>>;
-  template struct Metadomain<SimEngine::SRPIC, metric::QSpherical<Dim::_2D>>;
-  template struct Metadomain<SimEngine::GRPIC, metric::KerrSchild<Dim::_2D>>;
-  template struct Metadomain<SimEngine::GRPIC, metric::QKerrSchild<Dim::_2D>>;
-  template struct Metadomain<SimEngine::GRPIC, metric::KerrSchild0<Dim::_2D>>;
+#define METADOMAIN_STATS(S, M)                                                             \
+  template void Metadomain<S, M>::InitStatsWriter(const SimulationParams&, bool);          \
+  template auto Metadomain<S, M>::WriteStats(                                              \
+    const SimulationParams&,                                                               \
+    timestep_t,                                                                            \
+    timestep_t,                                                                            \
+    simtime_t,                                                                             \
+    simtime_t,                                                                             \
+    std::function<real_t(const std::string&, timestep_t, simtime_t, const Domain<S, M>&)>) \
+    -> bool;
+
+  METADOMAIN_STATS(SimEngine::SRPIC, metric::Minkowski<Dim::_1D>)
+  METADOMAIN_STATS(SimEngine::SRPIC, metric::Minkowski<Dim::_2D>)
+  METADOMAIN_STATS(SimEngine::SRPIC, metric::Minkowski<Dim::_3D>)
+  METADOMAIN_STATS(SimEngine::SRPIC, metric::Spherical<Dim::_2D>)
+  METADOMAIN_STATS(SimEngine::SRPIC, metric::QSpherical<Dim::_2D>)
+  METADOMAIN_STATS(SimEngine::GRPIC, metric::KerrSchild<Dim::_2D>)
+  METADOMAIN_STATS(SimEngine::GRPIC, metric::QKerrSchild<Dim::_2D>)
+  METADOMAIN_STATS(SimEngine::GRPIC, metric::KerrSchild0<Dim::_2D>)
+
+#undef METADOMAIN_STATS
 
 } // namespace ntt
diff --git a/src/framework/parameters.cpp b/src/framework/parameters.cpp
index 168640e1..40779d73 100644
--- a/src/framework/parameters.cpp
+++ b/src/framework/parameters.cpp
@@ -31,10 +31,10 @@
 namespace ntt {
 
   template <typename M>
-  auto get_dx0_V0(
-    const std::vector<ncells_t>&         resolution,
-    const boundaries_t<real_t>&          extent,
-    const std::map<std::string, real_t>& params) -> std::pair<real_t, real_t> {
+  auto get_dx0_V0(const std::vector<ncells_t>&         resolution,
+                  const boundaries_t<real_t>&          extent,
+                  const std::map<std::string, real_t>& params)
+    -> std::pair<real_t, real_t> {
     const auto      metric = M(resolution, extent, params);
     const auto      dx0    = metric.dxMin();
     coord_t<M::Dim> x_corner { ZERO };
@@ -549,6 +549,12 @@ namespace ntt {
                       "stats",
                       "quantities",
                       defaults::output::stats_quantities));
+    set("output.stats.custom",
+        toml::find_or(toml_data,
+                      "output",
+                      "stats",
+                      "custom",
+                      std::vector<std::string> {}));
 
     // intervals
     for (const auto& type : { "fields", "particles", "spectra", "stats" }) {
@@ -873,6 +879,11 @@ namespace ntt {
                         "gamma_rad",
                         defaults::synchrotron::gamma_rad));
     }
+
+    // @TODO: disabling stats for non-Cartesian
+    if (coord_enum != Coord::Cart) {
+      set("output.stats.enable", false);
+    }
   }
 
   void SimulationParams::setSetupParams(const toml::value& toml_data) {
diff --git a/src/global/arch/traits.h b/src/global/arch/traits.h
index 6d6c51f7..d9e5e931 100644
--- a/src/global/arch/traits.h
+++ b/src/global/arch/traits.h
@@ -99,7 +99,7 @@ namespace traits {
     using ext_force_t = decltype(&T::ext_force);
 
     template <typename T>
-    using ext_current_t = decltype(&T::ext_current);	    
+    using ext_current_t = decltype(&T::ext_current);
 
     template <typename T>
     using atm_fields_t = decltype(&T::AtmFields);
@@ -145,6 +145,9 @@ namespace traits {
 
     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
diff --git a/src/global/enums.h b/src/global/enums.h
index 7e06d4a8..08130a2c 100644
--- a/src/global/enums.h
+++ b/src/global/enums.h
@@ -14,6 +14,8 @@
  *   - enum ntt::Cooling           // synchrotron, none
  *   - enum ntt::FldsID            // e, dive, d, divd, b, h, j,
  *                                    a, t, rho, charge, n, nppc, v, custom
+ *   - enum ntt::StatsID           // b^2, e^2, exb, j.e, t, rho,
+ *                                    charge, n, npart
  * @namespaces:
  *   - ntt::
  * @note Enums of the same type can be compared with each other and with strings
@@ -321,15 +323,16 @@ namespace ntt {
       Charge  = 7,
       N       = 8,
       Npart   = 9,
+      Custom  = 10,
     };
 
     constexpr StatsID(uint8_t c) : enums_hidden::BaseEnum<StatsID> { c } {}
 
     static constexpr type        variants[] = { B2,  E2,     ExB, JdotE, T,
-                                                Rho, Charge, N,   Npart };
-    static constexpr const char* lookup[]   = { "b^2",    "e^2", "exb",
-                                                "j.e",    "t",   "rho",
-                                                "charge", "n",   "npart" };
+                                                Rho, Charge, N,   Npart, Custom };
+    static constexpr const char* lookup[] = { "b^2",   "e^2",   "exb",    "j.e",
+                                              "t",     "rho",   "charge", "n",
+                                              "npart", "custom" };
     static constexpr std::size_t total = sizeof(variants) / sizeof(variants[0]);
   };
 
diff --git a/src/global/tests/enums.cpp b/src/global/tests/enums.cpp
index 1190417e..f653f472 100644
--- a/src/global/tests/enums.cpp
+++ b/src/global/tests/enums.cpp
@@ -71,8 +71,8 @@ auto main() -> int {
                               "h",      "j",    "a",    "t",    "rho",
                               "charge", "n",    "nppc", "v",    "custom" };
 
-  enum_str_t all_out_stats = { "b^2", "e^2",    "exb", "j.e",  "t",
-                               "rho", "charge", "n",   "npart" };
+  enum_str_t all_out_stats = { "b^2", "e^2",    "exb", "j.e",   "t",
+                               "rho", "charge", "n",   "npart", "custom" };
 
   checkEnum<Coord>(all_coords);
   checkEnum<Metric>(all_metrics);
diff --git a/src/global/utils/error.h b/src/global/utils/error.h
index 9d5afed2..5f378bf1 100644
--- a/src/global/utils/error.h
+++ b/src/global/utils/error.h
@@ -34,6 +34,7 @@
 namespace raise {
   using namespace files;
 
+  [[noreturn]]
   inline void Error(const std::string& msg,
                     const std::string& file,
                     const std::string& func,
diff --git a/src/kernels/particle_moments.hpp b/src/kernels/particle_moments.hpp
index 8be68f5d..9c66aed0 100644
--- a/src/kernels/particle_moments.hpp
+++ b/src/kernels/particle_moments.hpp
@@ -59,7 +59,6 @@ namespace kernel {
     const bool               use_weights;
     const M                  metric;
     const int                ni2;
-    const real_t             inv_n0;
     const unsigned short     window;
 
     const real_t contrib;
@@ -113,11 +112,10 @@ namespace kernel {
       , use_weights { use_weights }
       , metric { metric }
       , ni2 { static_cast<int>(ni2) }
-      , inv_n0 { inv_n0 }
       , window { window }
       , contrib { get_contrib<F>(mass, charge) }
-      , smooth { ONE / (real_t)(math::pow(TWO * (real_t)window + ONE,
-                                          static_cast<int>(D))) } {
+      , smooth { inv_n0 / (real_t)(math::pow(TWO * (real_t)window + ONE,
+                                             static_cast<int>(D))) } {
       raise::ErrorIf(buff_idx >= N, "Invalid buffer index", HERE);
       raise::ErrorIf(window > N_GHOSTS, "Window size too large", HERE);
       raise::ErrorIf(((F == FldsID::Rho) || (F == FldsID::Charge)) && (mass == ZERO),
@@ -296,19 +294,21 @@ namespace kernel {
         // for nppc calculation ...
         // ... do not take volume, weights or smoothing into account
         if constexpr (D == Dim::_1D) {
-          coeff *= inv_n0 /
+          coeff *= smooth /
                    metric.sqrt_det_h({ static_cast<real_t>(i1(p)) + HALF });
         } else if constexpr (D == Dim::_2D) {
-          coeff *= inv_n0 /
+          coeff *= smooth /
                    metric.sqrt_det_h({ static_cast<real_t>(i1(p)) + HALF,
                                        static_cast<real_t>(i2(p)) + HALF });
         } else if constexpr (D == Dim::_3D) {
-          coeff *= inv_n0 /
+          coeff *= smooth /
                    metric.sqrt_det_h({ static_cast<real_t>(i1(p)) + HALF,
                                        static_cast<real_t>(i2(p)) + HALF,
                                        static_cast<real_t>(i3(p)) + HALF });
         }
-        coeff *= weight(p) * smooth;
+        if (use_weights) {
+          coeff *= weight(p);
+        }
       }
       auto buff_access = Buff.access();
       if constexpr (D == Dim::_1D) {
diff --git a/src/kernels/reduced_stats.hpp b/src/kernels/reduced_stats.hpp
new file mode 100644
index 00000000..68e6aa97
--- /dev/null
+++ b/src/kernels/reduced_stats.hpp
@@ -0,0 +1,551 @@
+/**
+ * @file kernels/reduced_stats.hpp
+ * @brief Compute reduced field/moment quantities for stats output
+ * @implements
+ *   - kernel::PrtlToPhys_kernel<>
+ * @namespaces:
+ *   - kernel::
+ */
+
+#ifndef KERNELS_REDUCED_STATS_HPP
+#define KERNELS_REDUCED_STATS_HPP
+
+#include "enums.h"
+#include "global.h"
+
+#include "arch/kokkos_aliases.h"
+#include "utils/numeric.h"
+
+namespace kernel {
+  using namespace ntt;
+
+  template <SimEngine::type S, class M, StatsID::type F, unsigned short I = 0>
+  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;
+    ndfield_t<D, 3> J;
+    const M         metric;
+
+  public:
+    ReducedFields_kernel(const ndfield_t<D, 6>& EM,
+                         const ndfield_t<D, 3>& J,
+                         const M&               metric)
+      : EM { EM }
+      , J { J }
+      , metric { metric } {}
+
+    Inline void operator()(index_t i1, real_t& buff) const {
+      const auto i1_ = COORD(i1);
+      if constexpr (F == StatsID::B2) {
+        if constexpr (I == 1) {
+          const auto b1_u = EM(i1, em::bx1);
+          const auto b1_d = metric.template transform<1, Idx::U, Idx::D>({ i1_ },
+                                                                         b1_u);
+          buff += b1_u * b1_d * metric.sqrt_det_h({ i1_ });
+        } else if constexpr (I == 2) {
+          const auto b2_u = EM(i1, em::bx2);
+          const auto b2_d = metric.template transform<2, Idx::U, Idx::D>(
+            { i1_ + HALF },
+            b2_u);
+          buff += b2_u * b2_d * metric.sqrt_det_h({ i1_ + HALF });
+        } else {
+          const auto b3_u = EM(i1, em::bx3);
+          const auto b3_d = metric.template transform<3, Idx::U, Idx::D>(
+            { i1_ + HALF },
+            b3_u);
+          buff += b3_u * b3_d * metric.sqrt_det_h({ i1_ + HALF });
+        }
+      } else if constexpr (F == StatsID::E2) {
+        if constexpr (I == 1) {
+          const auto e1_u = EM(i1, em::ex1);
+          const auto e1_d = metric.template transform<1, Idx::U, Idx::D>(
+            { i1_ + HALF },
+            e1_u);
+          buff += e1_u * e1_d * metric.sqrt_det_h({ i1_ + HALF });
+        } else if constexpr (I == 2) {
+          const auto e2_u = EM(i1, em::ex2);
+          const auto e2_d = metric.template transform<2, Idx::U, Idx::D>({ i1_ },
+                                                                         e2_u);
+          buff += e2_u * e2_d * metric.sqrt_det_h({ i1_ });
+        } else {
+          const auto e3_u = EM(i1, em::ex3);
+          const auto e3_d = metric.template transform<3, Idx::U, Idx::D>({ i1_ },
+                                                                         e3_u);
+          buff += e3_u * e3_d * metric.sqrt_det_h({ i1_ });
+        }
+      } else if constexpr (F == StatsID::ExB) {
+        if constexpr (I == 1) {
+          const auto e2_t = metric.template transform<2, Idx::U, Idx::T>(
+            { i1_ + HALF },
+            HALF * (EM(i1, em::ex2) + EM(i1 + 1, em::ex2)));
+          const auto e3_t = metric.template transform<3, Idx::U, Idx::T>(
+            { i1_ + HALF },
+            HALF * (EM(i1, em::ex3) + EM(i1 + 1, em::ex3)));
+          const auto b2_t = metric.template transform<2, Idx::U, Idx::T>(
+            { i1_ + HALF },
+            EM(i1, em::bx2));
+          const auto b3_t = metric.template transform<3, Idx::U, Idx::T>(
+            { i1_ + HALF },
+            EM(i1, em::bx3));
+          buff += (e2_t * b3_t - e3_t * b2_t) * metric.sqrt_det_h({ i1_ + HALF });
+        } else if constexpr (I == 2) {
+          const auto e1_t = metric.template transform<1, Idx::U, Idx::T>(
+            { i1_ + HALF },
+            EM(i1, em::ex1));
+          const auto e3_t = metric.template transform<3, Idx::U, Idx::T>(
+            { i1_ + HALF },
+            HALF * (EM(i1, em::ex3) + EM(i1 + 1, em::ex3)));
+          const auto b1_t = metric.template transform<1, Idx::U, Idx::T>(
+            { i1_ + HALF },
+            HALF * (EM(i1, em::bx1) + EM(i1 + 1, em::bx1)));
+          const auto b3_t = metric.template transform<3, Idx::U, Idx::T>(
+            { i1_ + HALF },
+            EM(i1, em::bx3));
+          buff += (e3_t * b1_t - e1_t * b3_t) * metric.sqrt_det_h({ i1_ + HALF });
+        } else {
+          const auto e1_t = metric.template transform<1, Idx::U, Idx::T>(
+            { i1_ + HALF },
+            EM(i1, em::ex1));
+          const auto e2_t = metric.template transform<2, Idx::U, Idx::T>(
+            { i1_ + HALF },
+            HALF * (EM(i1, em::ex2) + EM(i1 + 1, em::ex2)));
+          const auto b1_t = metric.template transform<1, Idx::U, Idx::T>(
+            { i1_ + HALF },
+            HALF * (EM(i1, em::bx1) + EM(i1 + 1, em::bx1)));
+          const auto b2_t = metric.template transform<2, Idx::U, Idx::T>(
+            { i1_ + HALF },
+            EM(i1, em::bx2));
+          buff += (e1_t * b2_t - e2_t * b1_t) * metric.sqrt_det_h({ i1_ + HALF });
+        }
+      } else if constexpr (F == StatsID::JdotE) {
+        vec_t<Dim::_3D> e_t { ZERO };
+        vec_t<Dim::_3D> j_t { ZERO };
+        metric.template transform<Idx::U, Idx::T>(
+          { i1_ + HALF },
+          { EM(i1, em::ex1),
+            HALF * (EM(i1, em::ex2) + EM(i1 + 1, em::ex2)),
+            HALF * (EM(i1, em::ex3) + EM(i1 + 1, em::ex3)) },
+          e_t);
+        metric.template transform<Idx::U, Idx::T>(
+          { i1_ + HALF },
+          { J(i1, cur::jx1),
+            HALF * (J(i1, cur::jx2) + J(i1 + 1, cur::jx2)),
+            HALF * (J(i1, cur::jx3) + J(i1 + 1, cur::jx3)) },
+          j_t);
+        buff += (e_t[0] * j_t[0] + e_t[1] * j_t[1] + e_t[2] * j_t[2]) *
+                metric.sqrt_det_h({ i1_ + HALF });
+      }
+    }
+
+    Inline void operator()(index_t i1, index_t i2, real_t& buff) const {
+      const auto i1_ = COORD(i1);
+      const auto i2_ = COORD(i2);
+      if constexpr (F == StatsID::B2) {
+        if constexpr (I == 1) {
+          const auto b1_u = EM(i1, i2, em::bx1);
+          const auto b1_d = metric.template transform<1, Idx::U, Idx::D>(
+            { i1_, i2_ + HALF },
+            b1_u);
+          buff += b1_u * b1_d * metric.sqrt_det_h({ i1_, i2_ + HALF });
+        } else if constexpr (I == 2) {
+          const auto b2_u = EM(i1, i2, em::bx2);
+          const auto b2_d = metric.template transform<2, Idx::U, Idx::D>(
+            { i1_ + HALF, i2_ },
+            b2_u);
+          buff += b2_u * b2_d * metric.sqrt_det_h({ i1_ + HALF, i2_ });
+        } else {
+          const auto b3_u = EM(i1, i2, em::bx3);
+          const auto b3_d = metric.template transform<3, Idx::U, Idx::D>(
+            { i1_ + HALF, i2_ + HALF },
+            b3_u);
+          buff += b3_u * b3_d * metric.sqrt_det_h({ i1_ + HALF, i2_ + HALF });
+        }
+      } else if constexpr (F == StatsID::E2) {
+        if constexpr (I == 1) {
+          const auto e1_u = EM(i1, i2, em::ex1);
+          const auto e1_d = metric.template transform<1, Idx::U, Idx::D>(
+            { i1_ + HALF, i2_ },
+            e1_u);
+          buff += e1_u * e1_d * metric.sqrt_det_h({ i1_ + HALF, i2_ });
+        } else if constexpr (I == 2) {
+          const auto e2_u = EM(i1, i2, em::ex2);
+          const auto e2_d = metric.template transform<2, Idx::U, Idx::D>(
+            { i1_, i2_ + HALF },
+            e2_u);
+          buff += e2_u * e2_d * metric.sqrt_det_h({ i1_, i2_ + HALF });
+        } else {
+          const auto e3_u = EM(i1, i2, em::ex3);
+          const auto e3_d = metric.template transform<3, Idx::U, Idx::D>(
+            { i1_, i2_ },
+            e3_u);
+          buff += e3_u * e3_d * metric.sqrt_det_h({ i1_, i2_ });
+        }
+      } else if constexpr (F == StatsID::ExB) {
+        if constexpr (I == 1) {
+          const auto e2_t = metric.template transform<2, Idx::U, Idx::T>(
+            { i1_ + HALF, i2_ + HALF },
+            HALF * (EM(i1, i2, em::ex2) + EM(i1 + 1, i2, em::ex2)));
+          const auto e3_t = metric.template transform<3, Idx::U, Idx::T>(
+            { i1_ + HALF, i2_ + HALF },
+            INV_4 * (EM(i1, i2, em::ex3) + EM(i1 + 1, i2, em::ex3) +
+                     EM(i1, i2 + 1, em::ex3) + EM(i1 + 1, i2 + 1, em::ex3)));
+          const auto b2_t = metric.template transform<2, Idx::U, Idx::T>(
+            { i1_ + HALF, i2_ + HALF },
+            HALF * (EM(i1, i2, em::bx2) + EM(i1, i2 + 1, em::bx2)));
+          const auto b3_t = metric.template transform<3, Idx::U, Idx::T>(
+            { i1_ + HALF, i2_ + HALF },
+            EM(i1, i2, em::bx3));
+          buff += (e2_t * b3_t - e3_t * b2_t) *
+                  metric.sqrt_det_h({ i1_ + HALF, i2_ + HALF });
+        } else if constexpr (I == 2) {
+          const auto e1_t = metric.template transform<1, Idx::U, Idx::T>(
+            { i1_ + HALF, i2_ + HALF },
+            HALF * (EM(i1, i2, em::ex1) + EM(i1, i2 + 1, em::ex1)));
+          const auto e3_t = metric.template transform<3, Idx::U, Idx::T>(
+            { i1_ + HALF, i2_ + HALF },
+            INV_4 * (EM(i1, i2, em::ex3) + EM(i1 + 1, i2, em::ex3) +
+                     EM(i1, i2 + 1, em::ex3) + EM(i1 + 1, i2 + 1, em::ex3)));
+          const auto b1_t = metric.template transform<1, Idx::U, Idx::T>(
+            { i1_ + HALF, i2_ + HALF },
+            HALF * (EM(i1, i2, em::bx1) + EM(i1 + 1, i2, em::bx1)));
+          const auto b3_t = metric.template transform<3, Idx::U, Idx::T>(
+            { i1_ + HALF, i2_ + HALF },
+            EM(i1, i2, em::bx3));
+          buff += (e3_t * b1_t - e1_t * b3_t) *
+                  metric.sqrt_det_h({ i1_ + HALF, i2_ + HALF });
+        } else {
+          const auto e1_t = metric.template transform<1, Idx::U, Idx::T>(
+            { i1_ + HALF, i2_ + HALF },
+            HALF * (EM(i1, i2, em::ex1) + EM(i1, i2 + 1, em::ex1)));
+          const auto e2_t = metric.template transform<2, Idx::U, Idx::T>(
+            { i1_ + HALF, i2_ + HALF },
+            HALF * (EM(i1, i2, em::ex2) + EM(i1 + 1, i2, em::ex2)));
+          const auto b1_t = metric.template transform<1, Idx::U, Idx::T>(
+            { i1_ + HALF, i2_ + HALF },
+            HALF * (EM(i1, i2, em::bx1) + EM(i1 + 1, i2, em::bx1)));
+          const auto b2_t = metric.template transform<2, Idx::U, Idx::T>(
+            { i1_ + HALF, i2_ + HALF },
+            HALF * (EM(i1, i2, em::bx2) + EM(i1, i2 + 1, em::bx2)));
+          buff += (e1_t * b2_t - e2_t * b1_t) *
+                  metric.sqrt_det_h({ i1_ + HALF, i2_ + HALF });
+        }
+      } else if constexpr (F == StatsID::JdotE) {
+        vec_t<Dim::_3D> e_t { ZERO };
+        vec_t<Dim::_3D> j_t { ZERO };
+        metric.template transform<Idx::U, Idx::T>(
+          { i1_ + HALF, i2_ + HALF },
+          { HALF * (EM(i1, i2, em::ex1) + EM(i1, i2 + 1, em::ex1)),
+            HALF * (EM(i1, i2, em::ex2) + EM(i1 + 1, i2, em::ex2)),
+            INV_4 * (EM(i1, i2, em::ex3) + EM(i1 + 1, i2, em::ex3) +
+                     EM(i1, i2 + 1, em::ex3) + EM(i1 + 1, i2 + 1, em::ex3)) },
+          e_t);
+        metric.template transform<Idx::U, Idx::T>(
+          { i1_ + HALF, i2_ + HALF },
+          { HALF * (J(i1, i2, cur::jx1) + J(i1, i2 + 1, cur::jx1)),
+            HALF * (J(i1, i2, cur::jx2) + J(i1 + 1, i2, cur::jx2)),
+            INV_4 * (J(i1, i2, cur::jx3) + J(i1 + 1, i2, cur::jx3) +
+                     J(i1, i2 + 1, cur::jx3) + J(i1 + 1, i2 + 1, cur::jx3)) },
+          j_t);
+        buff += (e_t[0] * j_t[0] + e_t[1] * j_t[1] + e_t[2] * j_t[2]) *
+                metric.sqrt_det_h({ i1_ + HALF, i2_ + HALF });
+      }
+    }
+
+    Inline void operator()(index_t i1, index_t i2, index_t i3, real_t& buff) const {
+      const auto i1_ = COORD(i1);
+      const auto i2_ = COORD(i2);
+      const auto i3_ = COORD(i3);
+      if constexpr (F == StatsID::B2) {
+        if constexpr (I == 1) {
+          const auto b1_u = EM(i1, i2, i3, em::bx1);
+          const auto b1_d = metric.template transform<1, Idx::U, Idx::D>(
+            { i1_, i2_ + HALF, i3_ + HALF },
+            b1_u);
+          buff += b1_u * b1_d * metric.sqrt_det_h({ i1_, i2_ + HALF, i3_ + HALF });
+        } else if constexpr (I == 2) {
+          const auto b2_u = EM(i1, i2, i3, em::bx2);
+          const auto b2_d = metric.template transform<2, Idx::U, Idx::D>(
+            { i1_ + HALF, i2_, i3_ + HALF },
+            b2_u);
+          buff += b2_u * b2_d * metric.sqrt_det_h({ i1_ + HALF, i2_, i3_ + HALF });
+        } else {
+          const auto b3_u = EM(i1, i2, i3, em::bx3);
+          const auto b3_d = metric.template transform<3, Idx::U, Idx::D>(
+            { i1_ + HALF, i2_ + HALF, i3_ },
+            b3_u);
+          buff += b3_u * b3_d * metric.sqrt_det_h({ i1_ + HALF, i2_ + HALF, i3_ });
+        }
+      } else if constexpr (F == StatsID::E2) {
+        if constexpr (I == 1) {
+          const auto e1_u = EM(i1, i2, i3, em::ex1);
+          const auto e1_d = metric.template transform<1, Idx::U, Idx::D>(
+            { i1_ + HALF, i2_, i3_ },
+            e1_u);
+          buff += e1_u * e1_d * metric.sqrt_det_h({ i1_ + HALF, i2_, i3_ });
+        } else if constexpr (I == 2) {
+          const auto e2_u = EM(i1, i2, i3, em::ex2);
+          const auto e2_d = metric.template transform<2, Idx::U, Idx::D>(
+            { i1_, i2_ + HALF, i3_ },
+            e2_u);
+          buff += e2_u * e2_d * metric.sqrt_det_h({ i1_, i2_ + HALF, i3_ });
+        } else {
+          const auto e3_u = EM(i1, i2, i3, em::ex3);
+          const auto e3_d = metric.template transform<3, Idx::U, Idx::D>(
+            { i1_, i2_, i3_ + HALF },
+            e3_u);
+          buff += e3_u * e3_d * metric.sqrt_det_h({ i1_, i2_, i3_ + HALF });
+        }
+      } else if constexpr (F == StatsID::ExB) {
+        if constexpr (I == 1) {
+          const auto e2_t = metric.template transform<2, Idx::U, Idx::T>(
+            { i1_ + HALF, i2_ + HALF, i3_ + HALF },
+            INV_4 *
+              (EM(i1, i2, i3, em::ex2) + EM(i1 + 1, i2, i3, em::ex2) +
+               EM(i1, i2, i3 + 1, em::ex2) + EM(i1 + 1, i2, i3 + 1, em::ex2)));
+          const auto e3_t = metric.template transform<3, Idx::U, Idx::T>(
+            { i1_ + HALF, i2_ + HALF, i3_ + HALF },
+            INV_4 *
+              (EM(i1, i2, i3, em::ex3) + EM(i1 + 1, i2, i3, em::ex3) +
+               EM(i1, i2 + 1, i3, em::ex3) + EM(i1 + 1, i2 + 1, i3, em::ex3)));
+          const auto b2_t = metric.template transform<2, Idx::U, Idx::T>(
+            { i1_ + HALF, i2_ + HALF, i3_ + HALF },
+            HALF * (EM(i1, i2, i3, em::bx2) + EM(i1, i2 + 1, i3, em::bx2)));
+          const auto b3_t = metric.template transform<3, Idx::U, Idx::T>(
+            { i1_ + HALF, i2_ + HALF, i3_ + HALF },
+            HALF * (EM(i1, i2, i3, em::bx3) + EM(i1, i2, i3 + 1, em::bx3)));
+          buff += (e2_t * b3_t - e3_t * b2_t) *
+                  metric.sqrt_det_h({ i1_ + HALF, i2_ + HALF, i3_ + HALF });
+        } else if constexpr (I == 2) {
+          const auto e1_t = metric.template transform<1, Idx::U, Idx::T>(
+            { i1_ + HALF, i2_ + HALF, i3_ + HALF },
+            INV_4 *
+              (EM(i1, i2, i3, em::ex1) + EM(i1, i2 + 1, i3, em::ex1) +
+               EM(i1, i2, i3 + 1, em::ex1) + EM(i1, i2 + 1, i3 + 1, em::ex1)));
+          const auto e3_t = metric.template transform<3, Idx::U, Idx::T>(
+            { i1_ + HALF, i2_ + HALF, i3_ + HALF },
+            INV_4 *
+              (EM(i1, i2, i3, em::ex3) + EM(i1 + 1, i2, i3, em::ex3) +
+               EM(i1, i2 + 1, i3, em::ex3) + EM(i1 + 1, i2 + 1, i3, em::ex3)));
+          const auto b1_t = metric.template transform<1, Idx::U, Idx::T>(
+            { i1_ + HALF, i2_ + HALF, i3_ + HALF },
+            HALF * (EM(i1, i2, i3, em::bx1) + EM(i1 + 1, i2, i3, em::bx1)));
+          const auto b3_t = metric.template transform<3, Idx::U, Idx::T>(
+            { i1_ + HALF, i2_ + HALF, i3_ + HALF },
+            HALF * (EM(i1, i2, i3, em::bx3) + EM(i1, i2, i3 + 1, em::bx3)));
+          buff += (e3_t * b1_t - e1_t * b3_t) *
+                  metric.sqrt_det_h({ i1_ + HALF, i2_ + HALF, i3_ + HALF });
+        } else {
+          const auto e1_t = metric.template transform<1, Idx::U, Idx::T>(
+            { i1_ + HALF, i2_ + HALF, i3_ + HALF },
+            INV_4 *
+              (EM(i1, i2, i3, em::ex1) + EM(i1, i2 + 1, i3, em::ex1) +
+               EM(i1, i2, i3 + 1, em::ex1) + EM(i1, i2 + 1, i3 + 1, em::ex1)));
+          const auto e2_t = metric.template transform<2, Idx::U, Idx::T>(
+            { i1_ + HALF, i2_ + HALF, i3_ + HALF },
+            INV_4 *
+              (EM(i1, i2, i3, em::ex2) + EM(i1 + 1, i2, i3, em::ex2) +
+               EM(i1, i2, i3 + 1, em::ex2) + EM(i1 + 1, i2, i3 + 1, em::ex2)));
+          const auto b1_t = metric.template transform<1, Idx::U, Idx::T>(
+            { i1_ + HALF, i2_ + HALF, i3_ + HALF },
+            HALF * (EM(i1, i2, i3, em::bx1) + EM(i1 + 1, i2, i3, em::bx1)));
+          const auto b2_t = metric.template transform<2, Idx::U, Idx::T>(
+            { i1_ + HALF, i2_ + HALF, i3_ + HALF },
+            HALF * (EM(i1, i2, i3, em::bx2) + EM(i1, i2 + 1, i3, em::bx2)));
+          buff += (e1_t * b2_t - e2_t * b1_t) *
+                  metric.sqrt_det_h({ i1_ + HALF, i2_ + HALF, i3_ + HALF });
+        }
+      } else if constexpr (F == StatsID::JdotE) {
+        vec_t<Dim::_3D> e_t { ZERO };
+        vec_t<Dim::_3D> j_t { ZERO };
+        metric.template transform<Idx::U, Idx::T>(
+          { i1_ + HALF, i2_ + HALF, i3_ + HALF },
+          { INV_4 * (EM(i1, i2, i3, em::ex1) + EM(i1, i2 + 1, i3, em::ex1) +
+                     EM(i1, i2, i3 + 1, em::ex1) + EM(i1, i2 + 1, i3 + 1, em::ex1)),
+            INV_4 * (EM(i1, i2, i3, em::ex2) + EM(i1 + 1, i2, i3, em::ex2) +
+                     EM(i1, i2, i3 + 1, em::ex2) + EM(i1 + 1, i2, i3 + 1, em::ex2)),
+            INV_4 * (EM(i1, i2, i3, em::ex3) + EM(i1 + 1, i2, i3, em::ex3) +
+                     EM(i1, i2 + 1, i3, em::ex3) + EM(i1 + 1, i2 + 1, i3, em::ex3)) },
+          e_t);
+        metric.template transform<Idx::U, Idx::T>(
+          { i1_ + HALF, i2_ + HALF, i3_ + HALF },
+          { INV_4 * (J(i1, i2, i3, cur::jx1) + J(i1, i2 + 1, i3, cur::jx1) +
+                     J(i1, i2, i3 + 1, cur::jx1) + J(i1, i2 + 1, i3 + 1, cur::jx1)),
+            INV_4 * (J(i1, i2, i3, cur::jx2) + J(i1 + 1, i2, i3, cur::jx2) +
+                     J(i1, i2, i3 + 1, cur::jx2) + J(i1 + 1, i2, i3 + 1, cur::jx2)),
+            INV_4 * (J(i1, i2, i3, cur::jx3) + J(i1 + 1, i2, i3, cur::jx3) +
+                     J(i1, i2 + 1, i3, cur::jx3) + J(i1 + 1, i2 + 1, i3, cur::jx3)) },
+          j_t);
+        buff += (e_t[0] * j_t[0] + e_t[1] * j_t[1] + e_t[2] * j_t[2]) *
+                metric.sqrt_det_h({ i1_ + HALF, i2_ + HALF, i3_ + HALF });
+      }
+    }
+  };
+
+  template <StatsID::type P>
+  auto get_contrib(float mass, float charge) -> real_t {
+    if constexpr (P == StatsID::Rho) {
+      return mass;
+    } else if constexpr (P == StatsID::Charge) {
+      return charge;
+    } else {
+      return ONE;
+    }
+  }
+
+  template <SimEngine::type S, class M, StatsID::type P>
+  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;
+    const array_t<real_t*>   ux1, ux2, ux3;
+    const array_t<real_t*>   phi;
+    const array_t<real_t*>   weight;
+    const array_t<short*>    tag;
+    const float              mass;
+    const float              charge;
+    const bool               use_weights;
+    const M                  metric;
+
+    const real_t contrib;
+
+  public:
+    ReducedParticleMoments_kernel(const std::vector<unsigned short>& components,
+                                  const array_t<int*>&               i1,
+                                  const array_t<int*>&               i2,
+                                  const array_t<int*>&               i3,
+                                  const array_t<prtldx_t*>&          dx1,
+                                  const array_t<prtldx_t*>&          dx2,
+                                  const array_t<prtldx_t*>&          dx3,
+                                  const array_t<real_t*>&            ux1,
+                                  const array_t<real_t*>&            ux2,
+                                  const array_t<real_t*>&            ux3,
+                                  const array_t<real_t*>&            phi,
+                                  const array_t<real_t*>&            weight,
+                                  const array_t<short*>&             tag,
+                                  float                              mass,
+                                  float                              charge,
+                                  bool     use_weights,
+                                  const M& metric)
+      : c1 { (components.size() > 0) ? components[0]
+                                     : static_cast<unsigned short>(0) }
+      , c2 { (components.size() == 2) ? components[1]
+                                      : static_cast<unsigned short>(0) }
+      , i1 { i1 }
+      , i2 { i2 }
+      , i3 { i3 }
+      , dx1 { dx1 }
+      , dx2 { dx2 }
+      , dx3 { dx3 }
+      , ux1 { ux1 }
+      , ux2 { ux2 }
+      , ux3 { ux3 }
+      , phi { phi }
+      , weight { weight }
+      , tag { tag }
+      , mass { mass }
+      , charge { charge }
+      , use_weights { use_weights }
+      , metric { metric }
+      , contrib { get_contrib<P>(mass, charge) } {
+      raise::ErrorIf(((P == StatsID::Rho) || (P == StatsID::Charge)) &&
+                       (mass == ZERO),
+                     "Rho & Charge for massless particles not defined",
+                     HERE);
+    }
+
+    Inline void operator()(index_t p, real_t& buff) const {
+      if (tag(p) == ParticleTag::dead) {
+        return;
+      }
+      if constexpr (P == StatsID::Npart) {
+        buff += ONE;
+        return;
+      } else if constexpr (P == StatsID::N or P == StatsID::Rho or
+                           P == StatsID::Charge) {
+        buff += use_weights ? weight(p) : contrib;
+        return;
+      } else {
+        // for stress-energy tensor
+        real_t          energy { ZERO };
+        vec_t<Dim::_3D> u_Phys { ZERO };
+        if constexpr (S == SimEngine::SRPIC) {
+          // SR
+          // stress-energy tensor for SR is computed in the tetrad (hatted) basis
+          if constexpr (M::CoordType == Coord::Cart) {
+            u_Phys[0] = ux1(p);
+            u_Phys[1] = ux2(p);
+            u_Phys[2] = ux3(p);
+          } else {
+            static_assert(D != Dim::_1D, "non-Cartesian SRPIC 1D");
+            coord_t<M::PrtlDim> x_Code { ZERO };
+            x_Code[0] = static_cast<real_t>(i1(p)) + static_cast<real_t>(dx1(p));
+            x_Code[1] = static_cast<real_t>(i2(p)) + static_cast<real_t>(dx2(p));
+            if constexpr (D == Dim::_3D) {
+              x_Code[2] = static_cast<real_t>(i3(p)) + static_cast<real_t>(dx3(p));
+            } else {
+              x_Code[2] = phi(p);
+            }
+            metric.template transform_xyz<Idx::XYZ, Idx::T>(
+              x_Code,
+              { ux1(p), ux2(p), ux3(p) },
+              u_Phys);
+          }
+          if (mass == ZERO) {
+            energy = NORM(u_Phys[0], u_Phys[1], u_Phys[2]);
+          } else {
+            energy = mass *
+                     math::sqrt(ONE + NORM_SQR(u_Phys[0], u_Phys[1], u_Phys[2]));
+          }
+        } else {
+          // GR
+          // stress-energy tensor for GR is computed in contravariant basis
+          static_assert(D != Dim::_1D, "GRPIC 1D");
+          coord_t<D> x_Code { ZERO };
+          x_Code[0] = static_cast<real_t>(i1(p)) + static_cast<real_t>(dx1(p));
+          x_Code[1] = static_cast<real_t>(i2(p)) + static_cast<real_t>(dx2(p));
+          if constexpr (D == Dim::_3D) {
+            x_Code[2] = static_cast<real_t>(i3(p)) + static_cast<real_t>(dx3(p));
+          }
+          vec_t<Dim::_3D> u_Cntrv { ZERO };
+          // compute u_i u^i for energy
+          metric.template transform<Idx::D, Idx::U>(x_Code,
+                                                    { ux1(p), ux2(p), ux3(p) },
+                                                    u_Cntrv);
+          energy = u_Cntrv[0] * ux1(p) + u_Cntrv[1] * ux2(p) + u_Cntrv[2] * ux3(p);
+          if (mass == ZERO) {
+            energy = math::sqrt(energy);
+          } else {
+            energy = mass * math::sqrt(ONE + energy);
+          }
+          metric.template transform<Idx::U, Idx::PU>(x_Code, u_Cntrv, u_Phys);
+        }
+        // compute the corresponding moment
+        real_t coeff = ONE;
+#pragma unroll
+        for (const auto& c : { c1, c2 }) {
+          if (c == 0) {
+            coeff *= energy;
+          } else {
+            coeff *= u_Phys[c - 1];
+          }
+        }
+        buff += coeff / energy;
+      }
+    }
+  };
+
+} // namespace kernel
+
+#endif
diff --git a/src/kernels/tests/CMakeLists.txt b/src/kernels/tests/CMakeLists.txt
index 551f9012..2702e052 100644
--- a/src/kernels/tests/CMakeLists.txt
+++ b/src/kernels/tests/CMakeLists.txt
@@ -35,3 +35,4 @@ gen_test(prtl_bc)
 gen_test(flds_bc)
 gen_test(pusher)
 gen_test(ext_force)
+gen_test(reduced_stats)
diff --git a/src/kernels/tests/particle_moments.cpp b/src/kernels/tests/particle_moments.cpp
index 25ed7d4d..ca3c2a7a 100644
--- a/src/kernels/tests/particle_moments.cpp
+++ b/src/kernels/tests/particle_moments.cpp
@@ -12,6 +12,8 @@
 #include "metrics/qspherical.h"
 #include "metrics/spherical.h"
 
+#include "kernels/reduced_stats.hpp"
+
 #include <Kokkos_Core.hpp>
 #include <Kokkos_ScatterView.hpp>
 
@@ -20,7 +22,6 @@
 #include <map>
 #include <stdexcept>
 #include <string>
-#include <utility>
 #include <vector>
 
 using namespace ntt;
@@ -56,7 +57,7 @@ void testParticleMoments(const std::vector<std::size_t>&      res,
   } else {
     extent = {
       ext[0],
-      {ZERO, constant::PI}
+      { ZERO, constant::PI }
     };
   }
 
@@ -106,8 +107,8 @@ void testParticleMoments(const std::vector<std::size_t>&      res,
   auto boundaries = boundaries_t<FldsBC> {};
   if constexpr (M::CoordType != Coord::Cart) {
     boundaries = {
-      {FldsBC::CUSTOM, FldsBC::CUSTOM},
-      {  FldsBC::AXIS,   FldsBC::AXIS}
+      { FldsBC::CUSTOM, FldsBC::CUSTOM },
+      {   FldsBC::AXIS,   FldsBC::AXIS }
     };
   }
 
@@ -117,84 +118,85 @@ void testParticleMoments(const std::vector<std::size_t>&      res,
   const unsigned short              window = 1;
 
   auto scatter_buff = Kokkos::Experimental::create_scatter_view(buff);
+  // clang-format off
   Kokkos::parallel_for(
-    "ParticleMoments",
-    10,
-    kernel::ParticleMoments_kernel<S, M, FldsID::T, 3>(comp1,
-                                                       scatter_buff,
-                                                       0,
-                                                       i1,
-                                                       i2,
-                                                       i3,
-                                                       dx1,
-                                                       dx2,
-                                                       dx3,
-                                                       ux1,
-                                                       ux2,
-                                                       ux3,
-                                                       phi,
-                                                       weight,
-                                                       tag,
-                                                       mass,
-                                                       charge,
+    "ParticleMoments", 10,
+    kernel::ParticleMoments_kernel<S, M, FldsID::T, 3>(comp1, scatter_buff, 0,
+                                                       i1, i2, i3,
+                                                       dx1, dx2, dx3,
+                                                       ux1, ux2, ux3,
+                                                       phi, weight, tag,
+                                                       mass, charge,
                                                        use_weights,
                                                        metric,
-                                                       boundaries,
-                                                       nx2,
-                                                       inv_n0,
-                                                       window));
+                                                       boundaries, nx2, inv_n0, window));
+
   Kokkos::parallel_for(
-    "ParticleMoments",
-    10,
-    kernel::ParticleMoments_kernel<S, M, FldsID::T, 3>(comp2,
-                                                       scatter_buff,
-                                                       1,
-                                                       i1,
-                                                       i2,
-                                                       i3,
-                                                       dx1,
-                                                       dx2,
-                                                       dx3,
-                                                       ux1,
-                                                       ux2,
-                                                       ux3,
-                                                       phi,
-                                                       weight,
-                                                       tag,
-                                                       mass,
-                                                       charge,
+    "ParticleMoments", 10,
+    kernel::ParticleMoments_kernel<S, M, FldsID::T, 3>(comp2, scatter_buff, 1,
+                                                       i1, i2, i3,
+                                                       dx1, dx2, dx3,
+                                                       ux1, ux2, ux3,
+                                                       phi, weight, tag,
+                                                       mass, charge,
                                                        use_weights,
                                                        metric,
-                                                       boundaries,
-                                                       nx2,
-                                                       inv_n0,
-                                                       window));
+                                                       boundaries, nx2, inv_n0, window));
   Kokkos::parallel_for(
-    "ParticleMoments",
-    10,
-    kernel::ParticleMoments_kernel<S, M, FldsID::T, 3>(comp3,
-                                                       scatter_buff,
-                                                       2,
-                                                       i1,
-                                                       i2,
-                                                       i3,
-                                                       dx1,
-                                                       dx2,
-                                                       dx3,
-                                                       ux1,
-                                                       ux2,
-                                                       ux3,
-                                                       phi,
-                                                       weight,
-                                                       tag,
-                                                       mass,
-                                                       charge,
+    "ParticleMoments", 10,
+    kernel::ParticleMoments_kernel<S, M, FldsID::T, 3>(comp3, scatter_buff, 2,
+                                                       i1, i2, i3,
+                                                       dx1, dx2, dx3,
+                                                       ux1, ux2, ux3,
+                                                       phi, weight, tag,
+                                                       mass, charge,
                                                        use_weights,
                                                        metric,
-                                                       boundaries,
-                                                       nx2,
-                                                       inv_n0,
-                                                       window));
+                                                       boundaries, nx2, inv_n0, window));
+
+  real_t n = ZERO, npart = ZERO, rho = ZERO, t00 = ZERO;
+  Kokkos::parallel_reduce(
+    "ReducedParticleMoments", 10,
+    kernel::ReducedParticleMoments_kernel<S, M, StatsID::N>({}, 
+                                                            i1, i2, i3,
+                                                            dx1, dx2, dx3,
+                                                            ux1, ux2, ux3,
+                                                            phi, weight, tag,
+                                                            mass, charge,
+                                                            use_weights,
+                                                            metric), n);
+
+  Kokkos::parallel_reduce(
+    "ReducedParticleMoments", 10,
+    kernel::ReducedParticleMoments_kernel<S, M, StatsID::Npart>({}, 
+                                                                i1, i2, i3,
+                                                                dx1, dx2, dx3,
+                                                                ux1, ux2, ux3,
+                                                                phi, weight, tag,
+                                                                mass, charge,
+                                                                use_weights,
+                                                                metric), npart);
+  Kokkos::parallel_reduce(
+    "ReducedParticleMoments", 10,
+    kernel::ReducedParticleMoments_kernel<S, M, StatsID::Rho>({}, 
+                                                              i1, i2, i3,
+                                                              dx1, dx2, dx3,
+                                                              ux1, ux2, ux3,
+                                                              phi, weight, tag,
+                                                              mass, charge,
+                                                              use_weights,
+                                                              metric), rho);
+  Kokkos::parallel_reduce(
+    "ReducedParticleMoments", 10,
+    kernel::ReducedParticleMoments_kernel<S, M, StatsID::T>({0u, 0u}, 
+                                                            i1, i2, i3,
+                                                            dx1, dx2, dx3,
+                                                            ux1, ux2, ux3,
+                                                            phi, weight, tag,
+                                                            mass, charge,
+                                                            use_weights,
+                                                            metric), t00);
+  // clang-format on
   Kokkos::Experimental::contribute(buff, scatter_buff);
 
   auto i1_h = Kokkos::create_mirror_view(i1);
@@ -231,6 +233,9 @@ void testParticleMoments(const std::vector<std::size_t>&      res,
     const real_t gammaSQR_1_expect = 15.0;
     const real_t gammaSQR_2_expect = 15.0;
 
+    const real_t n_expect   = 2.0;
+    const real_t t00_expect = 2.0 * math::sqrt(15.0);
+
     errorIf(not cmp::AlmostEqual_host(gammaSQR_1, gammaSQR_1_expect, epsilon * acc),
             fmt::format("wrong gamma_1 %.12e %.12e for %dD %s",
                         gammaSQR_1,
@@ -243,6 +248,31 @@ void testParticleMoments(const std::vector<std::size_t>&      res,
                         gammaSQR_2_expect,
                         metric.Dim,
                         metric.Label));
+
+    errorIf(not cmp::AlmostEqual_host(n, n_expect, epsilon * acc),
+            fmt::format("wrong n reduction %.12e %.12e for %dD %s",
+                        n,
+                        n_expect,
+                        metric.Dim,
+                        metric.Label));
+    errorIf(not cmp::AlmostEqual_host(npart, n_expect, epsilon * acc),
+            fmt::format("wrong npart reduction %.12e %.12e for %dD %s",
+                        npart,
+                        n_expect,
+                        metric.Dim,
+                        metric.Label));
+    errorIf(not cmp::AlmostEqual_host(rho, n_expect, epsilon * acc),
+            fmt::format("wrong rho reduction %.12e %.12e for %dD %s",
+                        rho,
+                        n_expect,
+                        metric.Dim,
+                        metric.Label));
+    errorIf(not cmp::AlmostEqual_host(t00, t00_expect, epsilon * acc),
+            fmt::format("wrong t00 reduction %.12e %.12e for %dD %s",
+                        t00,
+                        t00_expect,
+                        metric.Dim,
+                        metric.Label));
   }
 }
 
@@ -286,4 +316,4 @@ auto main(int argc, char* argv[]) -> int {
   }
   Kokkos::finalize();
   return 0;
-}
\ No newline at end of file
+}
diff --git a/src/kernels/tests/reduced_stats.cpp b/src/kernels/tests/reduced_stats.cpp
new file mode 100644
index 00000000..caae51cc
--- /dev/null
+++ b/src/kernels/tests/reduced_stats.cpp
@@ -0,0 +1,322 @@
+#include "kernels/reduced_stats.hpp"
+
+#include "enums.h"
+#include "global.h"
+
+#include "arch/kokkos_aliases.h"
+#include "utils/comparators.h"
+#include "utils/error.h"
+
+#include "metrics/minkowski.h"
+
+#include <iostream>
+#include <limits>
+#include <utility>
+
+using namespace ntt;
+using namespace metric;
+
+template <Dimension D, unsigned short N>
+class Fill_kernel {
+  ndfield_t<D, N>& arr;
+  real_t           v;
+  unsigned short   c;
+
+public:
+  Fill_kernel(ndfield_t<D, N>& arr_, real_t v_, unsigned short c_)
+    : arr { arr_ }
+    , v { v_ }
+    , c { c_ } {}
+
+  Inline void operator()(index_t i1) const {
+    arr(i1, c) = v;
+  }
+
+  Inline void operator()(index_t i1, index_t i2) const {
+    arr(i1, i2, c) = v;
+  }
+
+  Inline void operator()(index_t i1, index_t i2, index_t i3) const {
+    arr(i1, i2, i3, c) = v;
+  }
+};
+
+template <Dimension D, unsigned short N>
+void put_value(ndfield_t<D, N>& arr, real_t v, unsigned short c) {
+  range_t<D> range;
+  if constexpr (D == Dim::_1D) {
+    range = {
+      { 0u, arr.extent(0) }
+    };
+  } else if constexpr (D == Dim::_2D) {
+    range = {
+      {            0u,            0u },
+      { arr.extent(0), arr.extent(1) }
+    };
+  } else {
+    range = {
+      {            0u,            0u,            0u },
+      { arr.extent(0), arr.extent(1), arr.extent(2) }
+    };
+  }
+  Kokkos::parallel_for("Fill", range, Fill_kernel<D, N>(arr, v, c));
+}
+
+template <SimEngine::type S, class M, StatsID::type F, unsigned short I = 0>
+auto compute_field_stat(const M&                    metric,
+                        const ndfield_t<M::Dim, 6>& em,
+                        const ndfield_t<M::Dim, 3>& j,
+                        const range_t<M::Dim>&      range,
+                        ncells_t                    num_cells) -> real_t {
+  real_t buff = ZERO;
+  Kokkos::parallel_reduce("ReduceFields",
+                          range,
+                          kernel::ReducedFields_kernel<S, M, F, I>(em, j, metric),
+                          buff);
+  return buff / metric.totVolume();
+}
+
+auto almost_equal(real_t a, real_t b, real_t acc) -> bool {
+  return (math::fabs(a - b) < acc * math::max(math::fabs(a), math::fabs(b))) +
+         (real_t)1e-10;
+}
+
+template <SimEngine::type S, class M>
+void testReducedStats(const std::vector<std::size_t>& res,
+                      const boundaries_t<real_t>&     ext,
+                      const real_t                    acc) {
+  raise::ErrorIf(res.size() != M::Dim, "Invalid resolution size", HERE);
+
+  M metric { res, ext, {} };
+
+  std::vector<ncells_t> x_s, y_s, z_s;
+
+  coord_t<M::Dim>     dummy { ZERO };
+  std::vector<real_t> values;
+  values.push_back(metric.template transform<1, Idx::T, Idx::U>(dummy, ONE));
+  values.push_back(metric.template transform<2, Idx::T, Idx::U>(dummy, TWO));
+  values.push_back(metric.template transform<3, Idx::T, Idx::U>(dummy, THREE));
+
+  values.push_back(metric.template transform<1, Idx::T, Idx::U>(dummy, FOUR * ONE));
+  values.push_back(metric.template transform<2, Idx::T, Idx::U>(dummy, FOUR * TWO));
+  values.push_back(
+    metric.template transform<3, Idx::T, Idx::U>(dummy, FOUR * THREE));
+
+  values.push_back(metric.template transform<1, Idx::T, Idx::U>(dummy, -ONE));
+  values.push_back(metric.template transform<2, Idx::T, Idx::U>(dummy, -TWO));
+  values.push_back(metric.template transform<3, Idx::T, Idx::U>(dummy, THREE));
+
+  values.push_back(metric.template transform<1, Idx::T, Idx::U>(dummy, FOUR));
+  values.push_back(metric.template transform<2, Idx::T, Idx::U>(dummy, TWO));
+  values.push_back(metric.template transform<3, Idx::T, Idx::U>(dummy, ONE));
+
+  ndfield_t<M::Dim, 6> EM;
+  ndfield_t<M::Dim, 3> J;
+  range_t<M::Dim>      cell_range;
+  ncells_t             num_cells;
+
+  if constexpr (M::Dim == Dim::_1D) {
+    EM         = ndfield_t<M::Dim, 6> { "EM", res[0] + 2 * N_GHOSTS };
+    J          = ndfield_t<M::Dim, 3> { "J", res[0] + 2 * N_GHOSTS };
+    cell_range = { N_GHOSTS, res[0] + N_GHOSTS };
+    num_cells  = res[0];
+
+    put_value<M::Dim, 6>(EM, values[0], em::ex1);
+    put_value<M::Dim, 6>(EM, values[1], em::ex2);
+    put_value<M::Dim, 6>(EM, values[2], em::ex3);
+
+    put_value<M::Dim, 6>(EM, values[6], em::bx1);
+    put_value<M::Dim, 6>(EM, values[7], em::bx2);
+    put_value<M::Dim, 6>(EM, values[8], em::bx3);
+
+    put_value<M::Dim, 3>(J, values[9], cur::jx1);
+    put_value<M::Dim, 3>(J, values[10], cur::jx2);
+    put_value<M::Dim, 3>(J, values[11], cur::jx3);
+  } else if constexpr (M::Dim == Dim::_2D) {
+    EM = ndfield_t<M::Dim, 6> { "EM", res[0] + 2 * N_GHOSTS, res[1] + 2 * N_GHOSTS };
+    J = ndfield_t<M::Dim, 3> { "J", res[0] + 2 * N_GHOSTS, res[1] + 2 * N_GHOSTS };
+
+    cell_range = {
+      {          N_GHOSTS,          N_GHOSTS },
+      { res[0] + N_GHOSTS, res[1] + N_GHOSTS }
+    };
+    num_cells = res[0] * res[1];
+
+    put_value<M::Dim, 6>(EM, values[0], em::ex1);
+    put_value<M::Dim, 6>(EM, values[1], em::ex2);
+    put_value<M::Dim, 6>(EM, values[2], em::ex3);
+
+    put_value<M::Dim, 6>(EM, values[6], em::bx1);
+    put_value<M::Dim, 6>(EM, values[7], em::bx2);
+    put_value<M::Dim, 6>(EM, values[8], em::bx3);
+
+    put_value<M::Dim, 3>(J, values[9], cur::jx1);
+    put_value<M::Dim, 3>(J, values[10], cur::jx2);
+    put_value<M::Dim, 3>(J, values[11], cur::jx3);
+  } else {
+    EM = ndfield_t<M::Dim, 6> { "EM",
+                                res[0] + 2 * N_GHOSTS,
+                                res[1] + 2 * N_GHOSTS,
+                                res[2] + 2 * N_GHOSTS };
+    J  = ndfield_t<M::Dim, 3> { "J",
+                                res[0] + 2 * N_GHOSTS,
+                                res[1] + 2 * N_GHOSTS,
+                                res[2] + 2 * N_GHOSTS };
+
+    cell_range = {
+      {          N_GHOSTS,          N_GHOSTS,          N_GHOSTS },
+      { res[0] + N_GHOSTS, res[1] + N_GHOSTS, res[2] + N_GHOSTS }
+    };
+    num_cells = res[0] * res[1] * res[2];
+
+    put_value<M::Dim, 6>(EM, values[0], em::ex1);
+    put_value<M::Dim, 6>(EM, values[1], em::ex2);
+    put_value<M::Dim, 6>(EM, values[2], em::ex3);
+
+    put_value<M::Dim, 6>(EM, values[6], em::bx1);
+    put_value<M::Dim, 6>(EM, values[7], em::bx2);
+    put_value<M::Dim, 6>(EM, values[8], em::bx3);
+
+    put_value<M::Dim, 3>(J, values[9], cur::jx1);
+    put_value<M::Dim, 3>(J, values[10], cur::jx2);
+    put_value<M::Dim, 3>(J, values[11], cur::jx3);
+  }
+
+  {
+    const auto Ex_Sq = compute_field_stat<S, M, StatsID::E2, 1>(metric,
+                                                                EM,
+                                                                J,
+                                                                cell_range,
+                                                                num_cells);
+    raise::ErrorIf(not almost_equal(Ex_Sq, (real_t)(1), acc),
+                   "Ex_Sq does not match expected value",
+                   HERE);
+  }
+
+  {
+    const auto Ey_Sq = compute_field_stat<S, M, StatsID::E2, 2>(metric,
+                                                                EM,
+                                                                J,
+                                                                cell_range,
+                                                                num_cells);
+    raise::ErrorIf(not almost_equal(Ey_Sq, (real_t)(4), acc),
+                   "Ey_Sq does not match expected value",
+                   HERE);
+  }
+
+  {
+    const auto Ez_Sq = compute_field_stat<S, M, StatsID::E2, 3>(metric,
+                                                                EM,
+                                                                J,
+                                                                cell_range,
+                                                                num_cells);
+    raise::ErrorIf(not almost_equal(Ez_Sq, (real_t)(9), acc),
+                   "Ez_Sq does not match expected value",
+                   HERE);
+  }
+
+  {
+    const auto Bx_Sq = compute_field_stat<S, M, StatsID::B2, 1>(metric,
+                                                                EM,
+                                                                J,
+                                                                cell_range,
+                                                                num_cells);
+    raise::ErrorIf(not almost_equal(Bx_Sq, (real_t)(1), acc),
+                   "Bx_Sq does not match expected value",
+                   HERE);
+  }
+
+  {
+    const auto By_Sq = compute_field_stat<S, M, StatsID::B2, 2>(metric,
+                                                                EM,
+                                                                J,
+                                                                cell_range,
+                                                                num_cells);
+    raise::ErrorIf(not almost_equal(By_Sq, (real_t)(4), acc),
+                   "By_Sq does not match expected value",
+                   HERE);
+  }
+
+  {
+    const auto Bz_Sq = compute_field_stat<S, M, StatsID::B2, 3>(metric,
+                                                                EM,
+                                                                J,
+                                                                cell_range,
+                                                                num_cells);
+    raise::ErrorIf(not almost_equal(Bz_Sq, (real_t)(9), acc),
+                   "Bz_Sq does not match expected value",
+                   HERE);
+  }
+
+  {
+    const auto ExB_x = compute_field_stat<S, M, StatsID::ExB, 1>(metric,
+                                                                 EM,
+                                                                 J,
+                                                                 cell_range,
+                                                                 num_cells);
+    raise::ErrorIf(not almost_equal(ExB_x, (real_t)(12), acc),
+                   "ExB_x does not match expected value",
+                   HERE);
+  }
+
+  {
+    const auto ExB_y = compute_field_stat<S, M, StatsID::ExB, 2>(metric,
+                                                                 EM,
+                                                                 J,
+                                                                 cell_range,
+                                                                 num_cells);
+    raise::ErrorIf(not almost_equal(ExB_y, (real_t)(-6), acc),
+                   "ExB_y does not match expected value",
+                   HERE);
+  }
+
+  {
+    const auto ExB_z = compute_field_stat<S, M, StatsID::ExB, 3>(metric,
+                                                                 EM,
+                                                                 J,
+                                                                 cell_range,
+                                                                 num_cells);
+    raise::ErrorIf(not almost_equal(ExB_z, (real_t)(0), acc),
+                   "ExB_z does not match expected value",
+                   HERE);
+  }
+
+  {
+    const auto JdotE = compute_field_stat<S, M, StatsID::JdotE>(metric,
+                                                                EM,
+                                                                J,
+                                                                cell_range,
+                                                                num_cells);
+    raise::ErrorIf(not almost_equal(JdotE, (real_t)(11), acc),
+                   "JdotE does not match expected value",
+                   HERE);
+  }
+}
+
+auto main(int argc, char* argv[]) -> int {
+  Kokkos::initialize(argc, argv);
+
+  try {
+    using namespace ntt;
+
+    const ncells_t            nx = 100, ny = 123, nz = 52;
+    std::pair<real_t, real_t> x_ext { -2.0, 2.0 };
+    std::pair<real_t, real_t> y_ext { 0.0, 4.92 };
+    std::pair<real_t, real_t> z_ext { 0.0, 2.08 };
+
+    testReducedStats<SimEngine::SRPIC, Minkowski<Dim::_1D>>({ nx }, { x_ext }, 1e-6);
+    testReducedStats<SimEngine::SRPIC, Minkowski<Dim::_2D>>({ nx, ny },
+                                                            { x_ext, y_ext },
+                                                            1e-6);
+    testReducedStats<SimEngine::SRPIC, Minkowski<Dim::_3D>>({ nx, ny, nz },
+                                                            { x_ext, y_ext, z_ext },
+                                                            1e-6);
+
+  } catch (std::exception& e) {
+    std::cerr << e.what() << std::endl;
+    Kokkos::finalize();
+    return 1;
+  }
+  Kokkos::finalize();
+  return 0;
+}
diff --git a/src/metrics/kerr_schild.h b/src/metrics/kerr_schild.h
index 13294dd7..5edd43f2 100644
--- a/src/metrics/kerr_schild.h
+++ b/src/metrics/kerr_schild.h
@@ -127,6 +127,15 @@ namespace metric {
       return min_dx;
     }
 
+    /**
+     * total volume of the region described by the metric (in physical units)
+     */
+    [[nodiscard]]
+    auto totVolume() const -> real_t override {
+      // @TODO: Ask Alisa
+      return ZERO;
+    }
+
     /**
      * metric component with lower indices: h_ij
      * @param x coordinate array in code units
diff --git a/src/metrics/kerr_schild_0.h b/src/metrics/kerr_schild_0.h
index fdd2fd84..6e4f86c1 100644
--- a/src/metrics/kerr_schild_0.h
+++ b/src/metrics/kerr_schild_0.h
@@ -104,6 +104,20 @@ namespace metric {
       return min_dx;
     }
 
+    /**
+     * total volume of the region described by the metric (in physical units)
+     */
+    [[nodiscard]]
+    auto totVolume() const -> real_t override {
+      if constexpr (D == Dim::_1D) {
+        raise::Error("1D spherical metric not applicable", HERE);
+      } else if constexpr (D == Dim::_2D) {
+        return (SQR(x1_max) - SQR(x1_min)) * (x2_max - x2_min);
+      } else {
+        return (SQR(x1_max) - SQR(x1_min)) * (x2_max - x2_min) * (x3_max - x3_min);
+      }
+    }
+
     /**
      * metric component with lower indices: h_ij
      * @param x coordinate array in code units
diff --git a/src/metrics/metric_base.h b/src/metrics/metric_base.h
index 66cc2891..b386d46a 100644
--- a/src/metrics/metric_base.h
+++ b/src/metrics/metric_base.h
@@ -77,6 +77,9 @@ namespace metric {
     [[nodiscard]]
     virtual auto find_dxMin() const -> real_t = 0;
 
+    [[nodiscard]]
+    virtual auto totVolume() const -> real_t = 0;
+
     [[nodiscard]]
     auto dxMin() const -> real_t {
       return dx_min;
diff --git a/src/metrics/minkowski.h b/src/metrics/minkowski.h
index 244ce9f2..c1eda49d 100644
--- a/src/metrics/minkowski.h
+++ b/src/metrics/minkowski.h
@@ -55,16 +55,18 @@ namespace metric {
       , dx_inv { ONE / dx } {
       set_dxMin(find_dxMin());
       const auto epsilon = std::numeric_limits<real_t>::epsilon() *
-                         static_cast<real_t>(100.0);
+                           static_cast<real_t>(100.0);
       if constexpr (D != Dim::_1D) {
-        raise::ErrorIf(not cmp::AlmostEqual((x2_max - x2_min) / (real_t)(nx2), dx, epsilon),
-                       "dx2 must be equal to dx1 in 2D",
-                       HERE);
+        raise::ErrorIf(
+          not cmp::AlmostEqual((x2_max - x2_min) / (real_t)(nx2), dx, epsilon),
+          "dx2 must be equal to dx1 in 2D",
+          HERE);
       }
       if constexpr (D == Dim::_3D) {
-        raise::ErrorIf(not cmp::AlmostEqual((x3_max - x3_min) / (real_t)(nx3), dx, epsilon),
-                       "dx3 must be equal to dx1 in 3D",
-                       HERE);
+        raise::ErrorIf(
+          not cmp::AlmostEqual((x3_max - x3_min) / (real_t)(nx3), dx, epsilon),
+          "dx3 must be equal to dx1 in 3D",
+          HERE);
       }
     }
 
@@ -78,6 +80,20 @@ namespace metric {
       return dx / math::sqrt(static_cast<real_t>(D));
     }
 
+    /**
+     * total volume of the region described by the metric (in physical units)
+     */
+    [[nodiscard]]
+    auto totVolume() const -> real_t override {
+      if constexpr (D == Dim::_1D) {
+        return x1_max - x1_min;
+      } else if constexpr (D == Dim::_2D) {
+        return (x1_max - x1_min) * (x2_max - x2_min);
+      } else {
+        return (x1_max - x1_min) * (x2_max - x2_min) * (x3_max - x3_min);
+      }
+    }
+
     /**
      * metric component with lower indices: h_ij
      * @param x coordinate array in code units
@@ -242,7 +258,8 @@ namespace metric {
      * @note tetrad/cart <-> cntrv <-> cov
      */
     template <idx_t i, Idx in, Idx out>
-    Inline auto transform(const coord_t<D>& xi, const real_t& v_in) const -> real_t {
+    Inline auto transform(const coord_t<D>& xi, const real_t& v_in) const
+      -> real_t {
       static_assert(i > 0 && i <= 3, "Invalid index i");
       static_assert(in != out, "Invalid vector transformation");
       if constexpr (i > static_cast<idx_t>(D)) {
diff --git a/src/metrics/qkerr_schild.h b/src/metrics/qkerr_schild.h
index d64fe05d..7d327b89 100644
--- a/src/metrics/qkerr_schild.h
+++ b/src/metrics/qkerr_schild.h
@@ -132,6 +132,15 @@ namespace metric {
       return min_dx;
     }
 
+    /**
+     * total volume of the region described by the metric (in physical units)
+     */
+    [[nodiscard]]
+    auto totVolume() const -> real_t override {
+      // @TODO: Ask Alisa
+      return ZERO;
+    }
+
     /**
      * metric component with lower indices: h_ij
      * @param x coordinate array in code units
diff --git a/src/metrics/qspherical.h b/src/metrics/qspherical.h
index 4df5db86..2af1f02a 100644
--- a/src/metrics/qspherical.h
+++ b/src/metrics/qspherical.h
@@ -97,6 +97,20 @@ namespace metric {
       return min_dx;
     }
 
+    /**
+     * total volume of the region described by the metric (in physical units)
+     */
+    [[nodiscard]]
+    auto totVolume() const -> real_t override {
+      if constexpr (D == Dim::_1D) {
+        raise::Error("1D spherical metric not applicable", HERE);
+      } else if constexpr (D == Dim::_2D) {
+        return (SQR(x1_max) - SQR(x1_min)) * (x2_max - x2_min);
+      } else {
+        return (SQR(x1_max) - SQR(x1_min)) * (x2_max - x2_min) * (x3_max - x3_min);
+      }
+    }
+
     /**
      * metric component with lower indices: h_ij
      * @param x coordinate array in code units
@@ -284,7 +298,8 @@ namespace metric {
      * @note tetrad/sph <-> cntrv <-> cov
      */
     template <idx_t i, Idx in, Idx out>
-    Inline auto transform(const coord_t<D>& xi, const real_t& v_in) const -> real_t {
+    Inline auto transform(const coord_t<D>& xi, const real_t& v_in) const
+      -> real_t {
       static_assert(i > 0 && i <= 3, "Invalid index i");
       static_assert(in != out, "Invalid vector transformation");
       if constexpr ((in == Idx::T && out == Idx::Sph) ||
diff --git a/src/metrics/spherical.h b/src/metrics/spherical.h
index 54e12367..bf4c3453 100644
--- a/src/metrics/spherical.h
+++ b/src/metrics/spherical.h
@@ -76,6 +76,20 @@ namespace metric {
       return ONE / math::sqrt(ONE / SQR(dx1) + ONE / SQR(dx2));
     }
 
+    /**
+     * total volume of the region described by the metric (in physical units)
+     */
+    [[nodiscard]]
+    auto totVolume() const -> real_t override {
+      if constexpr (D == Dim::_1D) {
+        raise::Error("1D spherical metric not applicable", HERE);
+      } else if constexpr (D == Dim::_2D) {
+        return (SQR(x1_max) - SQR(x1_min)) * (x2_max - x2_min);
+      } else {
+        return (SQR(x1_max) - SQR(x1_min)) * (x2_max - x2_min) * (x3_max - x3_min);
+      }
+    }
+
     /**
      * metric component with lower indices: h_ij
      * @param x coordinate array in code units
@@ -252,7 +266,8 @@ namespace metric {
      * @note tetrad/sph <-> cntrv <-> cov
      */
     template <idx_t i, Idx in, Idx out>
-    Inline auto transform(const coord_t<D>& xi, const real_t& v_in) const -> real_t {
+    Inline auto transform(const coord_t<D>& xi, const real_t& v_in) const
+      -> real_t {
       static_assert(i > 0 && i <= 3, "Invalid index i");
       static_assert(in != out, "Invalid vector transformation");
       if constexpr ((in == Idx::T && out == Idx::Sph) ||
diff --git a/src/output/stats.cpp b/src/output/stats.cpp
index ba3e85c4..6aa65067 100644
--- a/src/output/stats.cpp
+++ b/src/output/stats.cpp
@@ -11,6 +11,7 @@
 
 #include <Kokkos_Core.hpp>
 
+#include <iomanip>
 #include <string>
 #include <vector>
 
@@ -19,7 +20,14 @@ using namespace out;
 
 namespace stats {
 
-  OutputStats::OutputStats(const std::string& name) : m_name { name } {
+  OutputStats::OutputStats(const std::string& name, bool is_custom)
+    : m_name { name } {
+    if (is_custom) {
+      m_id    = StatsID::Custom;
+      comp    = {};
+      species = {};
+      return;
+    }
     // determine the stats ID
     const auto pos = name.find("_");
     auto name_raw  = (pos == std::string::npos) ? name : name.substr(0, pos);
@@ -29,7 +37,7 @@ namespace stats {
     if (StatsID::contains(fmt::toLower(name_raw).c_str())) {
       m_id = StatsID::pick(fmt::toLower(name_raw).c_str());
     } else {
-      raise::Error("Unrecognized stats ID " + fmt::toLower(name_raw), HERE);
+      raise::Error("Unrecognized stats name: " + name, HERE);
     }
     // determine the species and components to output
     if (is_moment()) {
@@ -38,13 +46,13 @@ namespace stats {
       species = {};
     }
     if (is_vector()) {
-      // always write all the ExB and V components
+      // always write all the E^2, B^2, ExB components
       comp = { { 1 }, { 2 }, { 3 } };
     } else if (id() == StatsID::T) {
       // energy-momentum tensor
       comp = InterpretComponents({ name.substr(1, 1), name.substr(2, 1) });
     } else {
-      // scalar (e.g., Rho, E^2, etc.)
+      // scalar (e.g., Rho, Npart, etc.)
       comp = {};
     }
   }
@@ -61,9 +69,10 @@ namespace stats {
     m_fname = filename;
   }
 
-  void Writer::defineStatsOutputs(const std::vector<std::string>& stats_to_write) {
+  void Writer::defineStatsOutputs(const std::vector<std::string>& stats_to_write,
+                                  bool is_custom) {
     for (const auto& stat : stats_to_write) {
-      m_stat_writers.emplace_back(stat);
+      m_stat_writers.emplace_back(stat, is_custom);
     }
   }
 
@@ -71,14 +80,15 @@ namespace stats {
     CallOnce(
       [](auto& fname, auto& stat_writers) {
         std::fstream StatsOut(fname, std::fstream::out | std::fstream::app);
-        StatsOut << "step,time,";
+        StatsOut << std::setw(14) << "step" << "," << std::setw(14) << "time"
+                 << ",";
         for (const auto& stat : stat_writers) {
           if (stat.is_vector()) {
             for (auto i { 0u }; i < stat.comp.size(); ++i) {
-              StatsOut << stat.name(i) << ",";
+              StatsOut << std::setw(14) << stat.name(i) << ",";
             }
           } else {
-            StatsOut << stat.name() << ",";
+            StatsOut << std::setw(14) << stat.name() << ",";
           }
         }
         StatsOut << std::endl;
diff --git a/src/output/stats.h b/src/output/stats.h
index c81b9b7b..fc5bbf3a 100644
--- a/src/output/stats.h
+++ b/src/output/stats.h
@@ -3,10 +3,11 @@
  * @brief Class defining the metadata necessary to prepare the stats for output
  * @implements
  *   - out::OutputStats
+ *   - out::Writer
  * @cpp:
  *   - stats.cpp
  * @namespaces:
- *   - out::
+ *   - stats::
  */
 
 #ifndef OUTPUT_STATS_H
@@ -16,10 +17,17 @@
 #include "global.h"
 
 #include "utils/error.h"
+#include "utils/formatting.h"
 #include "utils/tools.h"
 
+#if defined(MPI_ENABLED)
+  #include "arch/mpi_aliases.h"
+
+  #include <mpi.h>
+#endif
+
 #include <fstream>
-#include <iostream>
+#include <iomanip>
 #include <string>
 #include <vector>
 
@@ -35,7 +43,7 @@ namespace stats {
     std::vector<std::vector<unsigned short>> comp {};
     std::vector<spidx_t>                     species {};
 
-    OutputStats(const std::string&);
+    OutputStats(const std::string&, bool);
 
     ~OutputStats() = default;
 
@@ -47,11 +55,19 @@ namespace stats {
 
     [[nodiscard]]
     auto is_vector() const -> bool {
-      return id() == StatsID::ExB;
+      return id() == StatsID::ExB || id() == StatsID::E2 || id() == StatsID::B2;
+    }
+
+    [[nodiscard]]
+    auto is_custom() const -> bool {
+      return id() == StatsID::Custom;
     }
 
     [[nodiscard]]
     inline auto name() const -> std::string {
+      if (id() == StatsID::Custom) {
+        return m_name;
+      }
       // generate the name
       std::string tmp = std::string(id().to_string());
       if (tmp == "exb") {
@@ -65,7 +81,11 @@ namespace stats {
       if (id() == StatsID::T) {
         tmp += m_name.substr(1, 2);
       } else if (is_vector()) {
-        tmp += "i";
+        if (id() == StatsID::E2 || id() == StatsID::B2) {
+          tmp = fmt::format("%ci^2", tmp[0]);
+        } else {
+          tmp += "i";
+        }
       }
       if (species.size() > 0) {
         tmp += "_";
@@ -101,16 +121,20 @@ namespace stats {
         // capitalize the first letter
         tmp[0] = std::toupper(tmp[0]);
       }
-      for (auto& c : comp[ci]) {
-        tmp += std::to_string(c);
-      }
-      if (species.size() > 0) {
-        tmp += "_";
-        for (auto& s : species) {
-          tmp += std::to_string(s);
+      if (tmp == "E^2" or tmp == "B^2") {
+        tmp = fmt::format("%c%d^2", tmp[0], comp[ci][0]);
+      } else {
+        for (auto& c : comp[ci]) {
+          tmp += std::to_string(c);
+        }
+        if (species.size() > 0) {
           tmp += "_";
+          for (auto& s : species) {
+            tmp += std::to_string(s);
+            tmp += "_";
+          }
+          tmp.pop_back();
         }
-        tmp.pop_back();
       }
       return tmp;
     }
@@ -135,7 +159,7 @@ namespace stats {
 
     void init(timestep_t, simtime_t);
     void defineStatsFilename(const std::string&);
-    void defineStatsOutputs(const std::vector<std::string>&);
+    void defineStatsOutputs(const std::vector<std::string>&, bool);
 
     void writeHeader();
 
@@ -144,17 +168,26 @@ namespace stats {
 
     template <typename T>
     inline void write(const T& value) const {
+      auto tot_value { static_cast<T>(0) };
 #if defined(MPI_ENABLED)
-        // @TODO: reduce
+      MPI_Reduce(&value,
+                 &tot_value,
+                 1,
+                 mpi::get_type<T>(),
+                 MPI_SUM,
+                 MPI_ROOT_RANK,
+                 MPI_COMM_WORLD);
+#else
+      tot_value = value;
 #endif
       CallOnce(
-        [](auto& fname, auto& value) {
+        [](auto&& fname, auto&& value) {
           std::fstream StatsOut(fname, std::fstream::out | std::fstream::app);
-          StatsOut << value << ",";
+          StatsOut << std::setw(14) << value << ",";
           StatsOut.close();
         },
         m_fname,
-        value);
+        tot_value);
     }
 
     void endWriting();
diff --git a/src/output/tests/stats.cpp b/src/output/tests/stats.cpp
index db6730a8..11e9961e 100644
--- a/src/output/tests/stats.cpp
+++ b/src/output/tests/stats.cpp
@@ -13,17 +13,17 @@ auto main() -> int {
   using namespace ntt;
   try {
     {
-      const auto e = OutputStats("E^2");
-      raise::ErrorIf(e.is_vector(), "E^2 should not be a vector quantity", HERE);
+      const auto e = OutputStats("E^2", false);
+      raise::ErrorIf(not e.is_vector(), "E^2 should be a vector quantity", HERE);
       raise::ErrorIf(e.is_moment(), "E^2 should not be a moment", HERE);
       raise::ErrorIf(e.id() != StatsID::E2, "E^2 should have ID StatsID::E2", HERE);
       raise::ErrorIf(e.species.size() != 0, "E^2 should have no species", HERE);
-      raise::ErrorIf(e.comp.size() != 0, "E^2 should have no components", HERE);
-      raise::ErrorIf(e.name() != "E^2", "E^2 should have name `E^2`", HERE);
+      raise::ErrorIf(e.comp.size() != 3, "E^2 should have 3 components", HERE);
+      raise::ErrorIf(e.name() != "Ei^2", "E^2 should have name `Ei^2`", HERE);
     }
 
     {
-      const auto e = OutputStats("ExB");
+      const auto e = OutputStats("ExB", false);
       raise::ErrorIf(not e.is_vector(), "ExB should be a vector quantity", HERE);
       raise::ErrorIf(e.is_moment(), "ExB should not be a moment", HERE);
       raise::ErrorIf(e.id() != StatsID::ExB, "ExB should have ID StatsID::ExB", HERE);
@@ -33,7 +33,7 @@ auto main() -> int {
     }
 
     {
-      const auto e = OutputStats("J.E");
+      const auto e = OutputStats("J.E", false);
       raise::ErrorIf(e.is_vector(), "J.E should not be a vector quantity", HERE);
       raise::ErrorIf(e.is_moment(), "J.E should not be a moment", HERE);
       raise::ErrorIf(e.id() != StatsID::JdotE,
@@ -45,7 +45,7 @@ auto main() -> int {
     }
 
     {
-      const auto rho = OutputStats("Rho_1_3");
+      const auto rho = OutputStats("Rho_1_3", false);
       raise::ErrorIf(not rho.is_moment(), "Rho should be a moment", HERE);
       raise::ErrorIf(rho.id() != StatsID::Rho,
                      "Rho should have ID StatsID::Rho",
@@ -58,7 +58,7 @@ auto main() -> int {
     }
 
     {
-      const auto t = OutputStats("Tti_2_3");
+      const auto t = OutputStats("Tti_2_3", false);
       raise::ErrorIf(not t.is_moment(), "T should be a moment", HERE);
       raise::ErrorIf(t.is_vector(), "T should not be a vector quantity", HERE);
       raise::ErrorIf(t.id() != StatsID::T, "T should have ID StatsID::T", HERE);
@@ -91,9 +91,21 @@ auto main() -> int {
     }
 
     {
-      const auto t = OutputStats("Tij");
+      const auto t = OutputStats("Tij", false);
       raise::ErrorIf(t.comp.size() != 6, "T should have 6 component", HERE);
     }
+
+    {
+      const auto custom = OutputStats("C2x_12", true);
+      raise::ErrorIf(custom.name() != "C2x_12",
+                     "Custom should have name `C2x_12`",
+                     HERE);
+      raise::ErrorIf(not custom.is_custom(),
+                     "Custom should be... well... a custom",
+                     HERE);
+      raise::ErrorIf(custom.is_moment(), "Custom should not be a moment", HERE);
+      raise::ErrorIf(custom.is_vector(), "Custom should not be a vector", HERE);
+    }
   } catch (const std::exception& e) {
     std::cerr << e.what() << std::endl;
     return 1;