Stats output

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:

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,14 +46,26 @@
     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"
   interval_time = 0.1
 
   [output.fields]
     quantities = ["N_1", "N_2", "E", "B", "J"]
+
+  [output.particles]
+    enable = false
+
+  [output.spectra]
+    enable = false
+
+[diagnostics]
+  colored_stdout = true
+  interval       = 10

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"
 

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

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

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

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

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)

src/framework/domain/metadomain.cpp

@@ -23,7 +23,6 @@
 
 #include <limits>
 #include <map>
-#include <memory>
 #include <string>
 #include <vector>
 

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