From 561b9163c17f351e2f06ab8371d3e660fef1474d Mon Sep 17 00:00:00 2001
From: liangwang0734 <liang.wang.phys@gmail.com>
Date: Sun, 25 Jun 2023 23:11:05 -0400
Subject: [PATCH 01/30] tfield: fix obsolete tfield_average_length

---
 src/include/OutputFieldsDefault.h | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/include/OutputFieldsDefault.h b/src/include/OutputFieldsDefault.h
index d3ff048462..f754c7d448 100644
--- a/src/include/OutputFieldsDefault.h
+++ b/src/include/OutputFieldsDefault.h
@@ -70,7 +70,6 @@ struct OutputFieldsItemParams
   int pfield_first = 0;
   int tfield_interval = 0;
   int tfield_first = 0;
-  int tfield_average_length = 1000000;
   int tfield_average_every = 1;
   Int3 rn = {};
   Int3 rx = {10000000, 10000000, 10000000};
@@ -121,7 +120,7 @@ class OutputFieldsItem : public OutputFieldsItemParams
     bool do_tfield = tfield_interval > 0 && timestep >= tfield_next_;
     bool doaccum_tfield =
       tfield_interval > 0 &&
-      (((timestep >= (tfield_next_ - tfield_average_length + 1)) &&
+      (((timestep >= (tfield_next_ - tfield_interval + 1)) &&
         timestep % tfield_average_every == 0) ||
        timestep == 0);
 

From 747f02894b24372bad25d9b4a4380732b9cf23a2 Mon Sep 17 00:00:00 2001
From: liangwang0734 <liang.wang.phys@gmail.com>
Date: Mon, 26 Jun 2023 00:36:30 -0400
Subject: [PATCH 02/30] fix tfield flag logic

---
 src/include/OutputFieldsDefault.h | 12 +++++++-----
 1 file changed, 7 insertions(+), 5 deletions(-)

diff --git a/src/include/OutputFieldsDefault.h b/src/include/OutputFieldsDefault.h
index f754c7d448..907acf2e56 100644
--- a/src/include/OutputFieldsDefault.h
+++ b/src/include/OutputFieldsDefault.h
@@ -113,16 +113,18 @@ class OutputFieldsItem : public OutputFieldsItemParams
       if (timestep != 0) {
         pfield_next_ = timestep + pfield_interval;
         tfield_next_ = timestep + tfield_interval;
+
+        tfield_next_ = timestep / tfield_interval * tfield_interval;
+        if (timestep % tfield_interval > 0)
+          tfield_next_ += tfield_interval;
       }
     }
 
     bool do_pfield = pfield_interval > 0 && timestep >= pfield_next_;
     bool do_tfield = tfield_interval > 0 && timestep >= tfield_next_;
-    bool doaccum_tfield =
-      tfield_interval > 0 &&
-      (((timestep >= (tfield_next_ - tfield_interval + 1)) &&
-        timestep % tfield_average_every == 0) ||
-       timestep == 0);
+    bool doaccum_tfield = timestep >= (tfield_next_ - tfield_average_every + 1);
+    doaccum_tfield = doaccum_tfield && (tfield_interval > 0);
+    doaccum_tfield = doaccum_tfield || (timestep == 0);
 
     if (do_pfield || doaccum_tfield) {
       prof_start(pr_eval);

From 3a0667086d273ba18530575da260626ddda472ea Mon Sep 17 00:00:00 2001
From: liangwang0734 <liang.wang.phys@gmail.com>
Date: Tue, 27 Jun 2023 16:35:50 -0400
Subject: [PATCH 03/30] Move writer_checkpoint function to the Checkpointing
 class

to reuse the IOAdios2 object

seems to get bp5 checkpointing to work in some cases
---
 src/include/checkpoint.hxx | 74 +++++++++++++++++++-------------------
 1 file changed, 37 insertions(+), 37 deletions(-)

diff --git a/src/include/checkpoint.hxx b/src/include/checkpoint.hxx
index e6a7d285ca..715af7d8f1 100644
--- a/src/include/checkpoint.hxx
+++ b/src/include/checkpoint.hxx
@@ -7,43 +7,6 @@
 #include "particles_simple.inl"
 #include <kg/io.h>
 
-// ----------------------------------------------------------------------
-// write_checkpoint
-//
-
-template <typename Mparticles, typename MfieldsState>
-void write_checkpoint(const Grid_t& grid, Mparticles& mprts,
-                      MfieldsState& mflds)
-{
-  static int pr, pr_A;
-  if (!pr) {
-    pr = prof_register("cp_write", 1., 0, 0);
-    pr_A = prof_register("cp_write_barier", 1., 0, 0);
-  }
-
-  prof_start(pr);
-  mpi_printf(grid.comm(), "**** Writing checkpoint...\n");
-#if defined(PSC_HAVE_ADIOS2) && !defined(VPIC)
-  prof_start(pr_A);
-  MPI_Barrier(grid.comm()); // not really necessary
-  prof_stop(pr_A);
-
-  std::string filename =
-    "checkpoint_" + std::to_string(grid.timestep()) + ".bp";
-
-  auto io = kg::io::IOAdios2{};
-  auto writer =
-    io.open(filename, kg::io::Mode::Write, grid.comm(), "checkpoint");
-  writer.put("grid", grid);
-  writer.put("mprts", mprts);
-  writer.put("mflds", mflds);
-  writer.close();
-#else
-  std::cerr << "write_checkpoint not available without adios2" << std::endl;
-  std::abort();
-#endif
-  prof_stop(pr);
-}
 
 // ----------------------------------------------------------------------
 // read_checkpoint
@@ -120,7 +83,44 @@ public:
     write_checkpoint(grid, mprts, mflds);
   }
 
+  template <typename Mparticles, typename MfieldsState>
+  void write_checkpoint(const Grid_t& grid, Mparticles& mprts,
+                        MfieldsState& mflds)
+  {
+    static int pr, pr_A;
+    if (!pr) {
+      pr = prof_register("cp_write", 1., 0, 0);
+      pr_A = prof_register("cp_write_barier", 1., 0, 0);
+    }
+
+    prof_start(pr);
+    mpi_printf(grid.comm(), "**** Writing checkpoint...\n");
+#if defined(PSC_HAVE_ADIOS2) && !defined(VPIC)
+    prof_start(pr_A);
+    MPI_Barrier(grid.comm()); // not really necessary
+    prof_stop(pr_A);
+
+    std::string filename =
+      "checkpoint_" + std::to_string(grid.timestep()) + ".bp";
+
+    mpi_printf(MPI_COMM_WORLD, "??? checkpoint IOAdios2 %p\n", &io);
+    writer = io.open(filename, kg::io::Mode::Write, grid.comm(), "checkpoint");
+    writer.put("grid", grid);
+    writer.put("mprts", mprts);
+    writer.put("mflds", mflds);
+    writer.performPuts();
+    writer.close();
+#else
+    std::cerr << "write_checkpoint not available without adios2" << std::endl;
+    std::abort();
+#endif
+    prof_stop(pr);
+  }
+
+
 private:
   int interval_; // write checkpoint every so many steps
   bool first_time_ = true;
+  kg::io::IOAdios2 io;
+  kg::io::Engine writer;
 };

From 343097f00f2178bb90dca2b6a09a05ae71d27f35 Mon Sep 17 00:00:00 2001
From: liangwang0734 <liang.wang.phys@gmail.com>
Date: Wed, 28 Jun 2023 01:00:34 -0400
Subject: [PATCH 04/30] checkpointing: add beginStep and endStep to support bp5

checkpoint: remove performPuts as it was suggested against for bp5 in the adios2 doc
---
 src/include/checkpoint.hxx | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/src/include/checkpoint.hxx b/src/include/checkpoint.hxx
index 715af7d8f1..d2b140d45b 100644
--- a/src/include/checkpoint.hxx
+++ b/src/include/checkpoint.hxx
@@ -22,6 +22,7 @@ inline void read_checkpoint(const std::string& filename, Grid_t& grid,
 #ifdef PSC_HAVE_ADIOS2
   auto io = kg::io::IOAdios2{};
   auto reader = io.open(filename, kg::io::Mode::Read);
+  reader.beginStep(kg::io::StepMode::Read);
   reader.get("grid", grid);
   mprts.~Mparticles();
   mflds.~MfieldsState();
@@ -29,6 +30,7 @@ inline void read_checkpoint(const std::string& filename, Grid_t& grid,
   new (&mflds) MfieldsState(grid);
   reader.get("mprts", mprts);
   reader.get("mflds", mflds);
+  reader.endStep();
   reader.close();
 
   // FIXME, when we read back a rebalanced grid, other existing fields will
@@ -105,10 +107,11 @@ public:
 
     mpi_printf(MPI_COMM_WORLD, "??? checkpoint IOAdios2 %p\n", &io);
     writer = io.open(filename, kg::io::Mode::Write, grid.comm(), "checkpoint");
+    writer.beginStep(kg::io::StepMode::Append);
     writer.put("grid", grid);
     writer.put("mprts", mprts);
     writer.put("mflds", mflds);
-    writer.performPuts();
+    writer.endStep();
     writer.close();
 #else
     std::cerr << "write_checkpoint not available without adios2" << std::endl;

From af53ca26ec0b4527d2e8a809637b569a4b69f9bb Mon Sep 17 00:00:00 2001
From: liangwang0734 <liang.wang.phys@gmail.com>
Date: Fri, 30 Jun 2023 17:48:51 -0400
Subject: [PATCH 05/30] TestIOAdios2: add .bp suffix; seems needed to pass bp4
 tests

---
 src/kg/testing/io/TestIOAdios2.cxx | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/kg/testing/io/TestIOAdios2.cxx b/src/kg/testing/io/TestIOAdios2.cxx
index 4c3c616247..65da036c1d 100644
--- a/src/kg/testing/io/TestIOAdios2.cxx
+++ b/src/kg/testing/io/TestIOAdios2.cxx
@@ -22,7 +22,7 @@ TEST(IOAdios2, OpenWriteThenRead)
 {
   auto io = kg::io::IOAdios2{};
   {
-    auto file = io.openFile("test2", kg::io::Mode::Write);
+    auto file = io.openFile("test2.bp", kg::io::Mode::Write);
   }
   {
     auto file = io.openFile("test2.bp", kg::io::Mode::Read);
@@ -33,7 +33,7 @@ TEST(IOAdios2, FilePutGetVariable)
 {
   auto io = kg::io::IOAdios2{};
   {
-    auto file = io.openFile("test3", kg::io::Mode::Write);
+    auto file = io.openFile("test3.bp", kg::io::Mode::Write);
     auto dbl = std::vector<double>{1., 2., 3., 4., 5.};
     file.putVariable("dbl", dbl.data(), kg::io::Mode::NonBlocking, {5},
                      {{0}, {5}}, {});
@@ -68,7 +68,7 @@ TEST(IOAdios2, FilePutGetAttribute)
 {
   auto io = kg::io::IOAdios2{};
   {
-    auto file = io.openFile("test4", kg::io::Mode::Write);
+    auto file = io.openFile("test4.bp", kg::io::Mode::Write);
     auto dbl = std::vector<double>{1., 2., 3., 4., 5.};
     file.putAttribute("attr_dbl", dbl.data(), dbl.size());
   }

From 617741aff25c7996cec10b8c4324e180942fddc1 Mon Sep 17 00:00:00 2001
From: James <jm1667@wildcats.unh.edu>
Date: Thu, 10 Nov 2022 10:24:46 -0500
Subject: [PATCH 06/30] setup_particles: consolidate init_npt funcs

may incur a slight slowdown during initialization phase,
but should make things more flexible
---
 src/include/psc.hxx             | 42 +++++++-----------
 src/include/setup_particles.hxx | 77 ++++++++++++++++++++-------------
 src/psc_bgk.cxx                 |  7 ++-
 3 files changed, 65 insertions(+), 61 deletions(-)

diff --git a/src/include/psc.hxx b/src/include/psc.hxx
index d05b87aefb..cbb1ca0d4e 100644
--- a/src/include/psc.hxx
+++ b/src/include/psc.hxx
@@ -781,39 +781,28 @@ Psc<PscConfig, Diagnostics, InjectParticles> makePscIntegrator(
 
 // ======================================================================
 // partitionAndSetupParticles
-// convenience/backward compatibility method to allow simpler function signature
-// init_npt signature: (int kind, Double3 pos, npt) -> void
 
-template <typename SetupParticles, typename Balance, typename Mparticles,
-          typename FUNC>
+template <typename SetupParticles, typename Balance, typename Mparticles>
 void partitionAndSetupParticles(SetupParticles& setup_particles,
                                 Balance& balance, Grid_t*& grid_ptr,
-                                Mparticles& mprts, FUNC&& init_npt)
+                                Mparticles& mprts, InitNptFunc init_npt)
 {
-  auto init_npt_general = [&](int kind, Double3 pos, int, Int3,
-                              psc_particle_npt& npt) {
-    init_npt(kind, pos, npt);
-  };
-  partitionParticlesGeneralInit(setup_particles, balance, grid_ptr, mprts,
-                                init_npt_general);
-  setupParticlesGeneralInit(setup_particles, mprts, init_npt_general);
+  partitionParticles(setup_particles, balance, grid_ptr, mprts, init_npt);
+  setupParticles(setup_particles, mprts, init_npt);
 }
 
 // ----------------------------------------------------------------------
-// partitionParticlesGeneralInit
-// init_npt signature: (int kind, Double3 pos, int p, Int3 index, npt) -> void
-
-template <typename SetupParticles, typename Balance, typename Mparticles,
-          typename FUNC>
-void partitionParticlesGeneralInit(SetupParticles& setup_particles,
-                                   Balance& balance, Grid_t*& grid_ptr,
-                                   Mparticles& mprts, FUNC&& init_npt)
+// partitionParticles
+
+template <typename SetupParticles, typename Balance, typename Mparticles>
+void partitionParticles(SetupParticles& setup_particles, Balance& balance,
+                        Grid_t*& grid_ptr, Mparticles& mprts,
+                        InitNptFunc init_npt)
 {
   auto comm = grid_ptr->comm();
   mpi_printf(comm, "**** Partitioning...\n");
 
-  auto n_prts_by_patch =
-    setup_particles.partition_general_init(*grid_ptr, init_npt);
+  auto n_prts_by_patch = setup_particles.partition(*grid_ptr, init_npt);
 
   balance.initial(grid_ptr, n_prts_by_patch);
   // !!! FIXME! grid is now invalid
@@ -824,12 +813,11 @@ void partitionParticlesGeneralInit(SetupParticles& setup_particles,
 }
 
 // ----------------------------------------------------------------------
-// setupParticlesGeneralInit
-// init_npt signature: (int kind, Double3 pos, int p, Int3 index, npt) -> void
+// setupParticles
 
-template <typename SetupParticles, typename Mparticles, typename FUNC>
-void setupParticlesGeneralInit(SetupParticles& setup_particles,
-                               Mparticles& mprts, FUNC&& init_npt)
+template <typename SetupParticles, typename Mparticles>
+void setupParticles(SetupParticles& setup_particles, Mparticles& mprts,
+                    InitNptFunc init_npt)
 {
   mpi_printf(MPI_COMM_WORLD, "**** Setting up particles...\n");
   setup_particles.setupParticles(mprts, init_npt);
diff --git a/src/include/setup_particles.hxx b/src/include/setup_particles.hxx
index 2f3e8580cc..27f3267899 100644
--- a/src/include/setup_particles.hxx
+++ b/src/include/setup_particles.hxx
@@ -1,6 +1,8 @@
 
 #pragma once
 
+#include <functional>
+#include <type_traits>
 #include "centering.hxx"
 
 struct psc_particle_npt
@@ -20,6 +22,43 @@ namespace
 const Centering::Centerer defaultCenterer(Centering::CC);
 }
 
+struct InitNptFunc
+{
+  // Initialize particles according to a Maxwellian.
+  // (int kind, Double3 pos, psc_particle_npt& npt) -> void
+  template <
+    typename INIT_NPT,
+    std::enable_if_t<
+      std::is_convertible<
+        INIT_NPT, std::function<void(int, Double3, psc_particle_npt&)>>::value,
+      bool> = true>
+  InitNptFunc(INIT_NPT init_npt_simple)
+    : func{[&](int kind, Double3 pos, int p, Int3 idx, psc_particle_npt& npt) {
+        init_npt_simple(kind, pos, npt);
+      }}
+  {}
+
+  // Initialize particles according to a Maxwellian.
+  // (int kind, Double3 pos, int patch, Int3 idx, psc_particle_npt& npt) -> void
+  template <
+    typename INIT_NPT,
+    std::enable_if_t<std::is_convertible<
+                       INIT_NPT, std::function<void(int, Double3, int, Int3,
+                                                    psc_particle_npt&)>>::value,
+                     bool> = true>
+  InitNptFunc(INIT_NPT init_npt) : func{init_npt}
+  {}
+
+  void operator()(int kind, Double3 pos, int patch, Int3 index,
+                  psc_particle_npt& npt)
+  {
+    func(kind, pos, patch, index, npt);
+  }
+
+private:
+  std::function<void(int, Double3, int, Int3, psc_particle_npt&)> func;
+};
+
 template <typename MP>
 struct SetupParticles
 {
@@ -62,11 +101,10 @@ struct SetupParticles
   // ----------------------------------------------------------------------
   // op_cellwise
   // Performs a given operation in each cell.
-  // init_npt signature: (int kind, Double3 pos, int p, Int3 index, npt) -> void
   // op signature: (int n_in_cell, npt, Double3 pos) -> void
 
-  template <typename InitFunc, typename OpFunc>
-  void op_cellwise(const Grid_t& grid, int patch, InitFunc&& init_npt,
+  template <typename OpFunc>
+  void op_cellwise(const Grid_t& grid, int patch, InitNptFunc init_npt,
                    OpFunc&& op)
   {
     Int3 ilo = Int3{}, ihi = grid.ldims;
@@ -157,22 +195,16 @@ struct SetupParticles
 
   // ----------------------------------------------------------------------
   // setup_particles
-  // init_npt signature: (int kind, Double3 pos, npt) -> void
 
-  template <typename FUNC>
-  void operator()(Mparticles& mprts, FUNC&& init_npt)
+  void operator()(Mparticles& mprts, InitNptFunc init_npt)
   {
-    setupParticles(mprts,
-                   [&](int kind, Double3 pos, int p, Int3 idx,
-                       psc_particle_npt& npt) { init_npt(kind, pos, npt); });
+    setupParticles(mprts, init_npt);
   }
 
   // ----------------------------------------------------------------------
   // setupParticles
-  // init_npt signature: (int kind, Double3 pos, int p, Int3 index, npt) -> void
 
-  template <typename FUNC>
-  void setupParticles(Mparticles& mprts, FUNC&& init_npt)
+  void setupParticles(Mparticles& mprts, InitNptFunc init_npt)
   {
     static int pr;
     if (!pr) {
@@ -189,7 +221,7 @@ struct SetupParticles
     for (int p = 0; p < mprts.n_patches(); ++p) {
       auto injector = inj[p];
 
-      op_cellwise(grid, p, std::forward<FUNC>(init_npt),
+      op_cellwise(grid, p, init_npt,
                   [&](int n_in_cell, psc_particle_npt& npt, Double3& pos) {
                     for (int cnt = 0; cnt < n_in_cell; cnt++) {
                       real_t wni;
@@ -209,28 +241,13 @@ struct SetupParticles
 
   // ----------------------------------------------------------------------
   // partition
-  // init_npt signature: (int kind, Double3 pos, npt) -> void
-
-  template <typename FUNC>
-  std::vector<uint> partition(const Grid_t& grid, FUNC&& init_npt)
-  {
-    return partition_general_init(
-      grid, [&](int kind, Double3 pos, int, Int3, psc_particle_npt& npt) {
-        init_npt(kind, pos, npt);
-      });
-  }
-
-  // ----------------------------------------------------------------------
-  // partition_general_init
-  // init_npt signature: (int kind, Double3 pos, int p, Int3 index, npt) -> void
 
-  template <typename FUNC>
-  std::vector<uint> partition_general_init(const Grid_t& grid, FUNC&& init_npt)
+  std::vector<uint> partition(const Grid_t& grid, InitNptFunc init_npt)
   {
     std::vector<uint> n_prts_by_patch(grid.n_patches());
 
     for (int p = 0; p < grid.n_patches(); ++p) {
-      op_cellwise(grid, p, std::forward<FUNC>(init_npt),
+      op_cellwise(grid, p, init_npt,
                   [&](int n_in_cell, psc_particle_npt&, Double3&) {
                     n_prts_by_patch[p] += n_in_cell;
                   });
diff --git a/src/psc_bgk.cxx b/src/psc_bgk.cxx
index 30469c4261..ba7da4a295 100644
--- a/src/psc_bgk.cxx
+++ b/src/psc_bgk.cxx
@@ -220,7 +220,7 @@ void initializeParticles(Balance& balance, Grid_t*& grid_ptr, Mparticles& mprts,
 
   auto&& qDensity = -psc::mflds::interior(divGradPhi.grid(), divGradPhi.gt());
 
-  auto npt_init = [&](int kind, double crd[3], int p, Int3 idx,
+  auto npt_init = [&](int kind, Double3 crd, int p, Int3 idx,
                       psc_particle_npt& npt) {
     double y = getCoord(crd[1]);
     double z = getCoord(crd[2]);
@@ -255,9 +255,8 @@ void initializeParticles(Balance& balance, Grid_t*& grid_ptr, Mparticles& mprts,
     }
   };
 
-  partitionParticlesGeneralInit(setup_particles, balance, grid_ptr, mprts,
-                                npt_init);
-  setupParticlesGeneralInit(setup_particles, mprts, npt_init);
+  partitionAndSetupParticles(setup_particles, balance, grid_ptr, mprts,
+                             npt_init);
 }
 
 // ======================================================================

From 2cdf867fccecb49e50755fa1d0d3dabe026d1409 Mon Sep 17 00:00:00 2001
From: James <jm1667@wildcats.unh.edu>
Date: Sun, 13 Nov 2022 11:54:15 -0500
Subject: [PATCH 07/30] setup_particles/distribution: move rng to latter

---
 src/include/distribution.hxx    | 30 ++++++++++++++++++++++++++++++
 src/include/setup_particles.hxx | 31 +++++++++++--------------------
 2 files changed, 41 insertions(+), 20 deletions(-)
 create mode 100644 src/include/distribution.hxx

diff --git a/src/include/distribution.hxx b/src/include/distribution.hxx
new file mode 100644
index 0000000000..8ef89eef98
--- /dev/null
+++ b/src/include/distribution.hxx
@@ -0,0 +1,30 @@
+#include <random>
+
+namespace distribution
+{
+
+// ======================================================================
+// Distribution
+// Base class for distributions
+
+template <typename Real>
+struct Distribution
+{
+  virtual Real get() = 0;
+};
+
+// ======================================================================
+// Uniform
+
+template <typename Real>
+struct Uniform : public Distribution<Real>
+{
+  Uniform(Real min, Real max) : dist{min, max} {}
+  Real get() override { return dist(generator); }
+
+private:
+  std::default_random_engine generator;
+  std::uniform_real_distribution<Real> dist;
+};
+
+} // namespace distribution
\ No newline at end of file
diff --git a/src/include/setup_particles.hxx b/src/include/setup_particles.hxx
index 27f3267899..fd48c9d970 100644
--- a/src/include/setup_particles.hxx
+++ b/src/include/setup_particles.hxx
@@ -4,6 +4,7 @@
 #include <functional>
 #include <type_traits>
 #include "centering.hxx"
+#include "distribution.hxx"
 
 struct psc_particle_npt
 {
@@ -154,31 +155,21 @@ struct SetupParticles
   psc::particle::Inject setupParticle(const psc_particle_npt& npt, Double3 pos,
                                       double wni)
   {
+    static distribution::Uniform<float> dist(0, 1);
     double beta = norm_.beta;
 
     assert(npt.kind >= 0 && npt.kind < kinds_.size());
     double m = kinds_[npt.kind].m;
 
-    float ran1, ran2, ran3, ran4, ran5, ran6;
-    do {
-      ran1 = random() / ((float)RAND_MAX + 1);
-      ran2 = random() / ((float)RAND_MAX + 1);
-      ran3 = random() / ((float)RAND_MAX + 1);
-      ran4 = random() / ((float)RAND_MAX + 1);
-      ran5 = random() / ((float)RAND_MAX + 1);
-      ran6 = random() / ((float)RAND_MAX + 1);
-    } while (ran1 >= 1.f || ran2 >= 1.f || ran3 >= 1.f || ran4 >= 1.f ||
-             ran5 >= 1.f || ran6 >= 1.f);
-
-    double pxi =
-      npt.p[0] + sqrtf(-2.f * npt.T[0] / m * sqr(beta) * logf(1.0 - ran1)) *
-                   cosf(2.f * M_PI * ran2);
-    double pyi =
-      npt.p[1] + sqrtf(-2.f * npt.T[1] / m * sqr(beta) * logf(1.0 - ran3)) *
-                   cosf(2.f * M_PI * ran4);
-    double pzi =
-      npt.p[2] + sqrtf(-2.f * npt.T[2] / m * sqr(beta) * logf(1.0 - ran5)) *
-                   cosf(2.f * M_PI * ran6);
+    double pxi = npt.p[0] + sqrtf(-2.f * npt.T[0] / m * sqr(beta) *
+                                  logf(1.0 - dist.get())) *
+                              cosf(2.f * M_PI * dist.get());
+    double pyi = npt.p[1] + sqrtf(-2.f * npt.T[1] / m * sqr(beta) *
+                                  logf(1.0 - dist.get())) *
+                              cosf(2.f * M_PI * dist.get());
+    double pzi = npt.p[2] + sqrtf(-2.f * npt.T[2] / m * sqr(beta) *
+                                  logf(1.0 - dist.get())) *
+                              cosf(2.f * M_PI * dist.get());
 
     if (initial_momentum_gamma_correction) {
       double gam;

From 8335e873de907f6c1c583688dd8d8f3b47a9b8ca Mon Sep 17 00:00:00 2001
From: James <jm1667@wildcats.unh.edu>
Date: Sun, 13 Nov 2022 12:42:44 -0500
Subject: [PATCH 08/30] setup_particles/distribution: impl normal dist

---
 src/include/distribution.hxx    | 33 +++++++++++++++++++++++++++++++++
 src/include/setup_particles.hxx | 14 ++++----------
 2 files changed, 37 insertions(+), 10 deletions(-)

diff --git a/src/include/distribution.hxx b/src/include/distribution.hxx
index 8ef89eef98..d0cfa7e63c 100644
--- a/src/include/distribution.hxx
+++ b/src/include/distribution.hxx
@@ -1,4 +1,5 @@
 #include <random>
+#include <cmath>
 
 namespace distribution
 {
@@ -27,4 +28,36 @@ private:
   std::uniform_real_distribution<Real> dist;
 };
 
+// ======================================================================
+// Normal
+
+template <typename Real>
+struct Normal : public Distribution<Real>
+{
+  Normal(Real mean, Real stdev) : mean{mean}, stdev{stdev} {}
+  Normal() : Normal(0, 1) {}
+
+  Real get() override { return get(mean, stdev); }
+  Real get(Real mean, Real stdev)
+  {
+    if (has_next_value) {
+      has_next_value = false;
+      return next_value;
+    }
+
+    // Box-Muller transform to generate 2 independent values at once
+    Real r = stdev * std::sqrt(-2 * std::log(1 - uniform.get()));
+    Real theta = uniform.get() * 2 * M_PI;
+
+    next_value = r * std::cos(theta) + mean;
+    return r * std::sin(theta) + mean;
+  }
+
+private:
+  Uniform<Real> uniform{0, 1};
+  Real mean, stdev;
+  Real next_value;
+  bool has_next_value;
+};
+
 } // namespace distribution
\ No newline at end of file
diff --git a/src/include/setup_particles.hxx b/src/include/setup_particles.hxx
index fd48c9d970..cd0f82040e 100644
--- a/src/include/setup_particles.hxx
+++ b/src/include/setup_particles.hxx
@@ -155,21 +155,15 @@ struct SetupParticles
   psc::particle::Inject setupParticle(const psc_particle_npt& npt, Double3 pos,
                                       double wni)
   {
-    static distribution::Uniform<float> dist(0, 1);
+    static distribution::Normal<double> dist;
     double beta = norm_.beta;
 
     assert(npt.kind >= 0 && npt.kind < kinds_.size());
     double m = kinds_[npt.kind].m;
 
-    double pxi = npt.p[0] + sqrtf(-2.f * npt.T[0] / m * sqr(beta) *
-                                  logf(1.0 - dist.get())) *
-                              cosf(2.f * M_PI * dist.get());
-    double pyi = npt.p[1] + sqrtf(-2.f * npt.T[1] / m * sqr(beta) *
-                                  logf(1.0 - dist.get())) *
-                              cosf(2.f * M_PI * dist.get());
-    double pzi = npt.p[2] + sqrtf(-2.f * npt.T[2] / m * sqr(beta) *
-                                  logf(1.0 - dist.get())) *
-                              cosf(2.f * M_PI * dist.get());
+    double pxi = dist.get(npt.p[0], beta * std::sqrt(npt.T[0] / m));
+    double pyi = dist.get(npt.p[1], beta * std::sqrt(npt.T[1] / m));
+    double pzi = dist.get(npt.p[2], beta * std::sqrt(npt.T[2] / m));
 
     if (initial_momentum_gamma_correction) {
       double gam;

From 3d3c9e3367d79135ddee3ce032d622cfd595384d Mon Sep 17 00:00:00 2001
From: James <jm1667@wildcats.unh.edu>
Date: Sun, 13 Nov 2022 18:13:48 -0500
Subject: [PATCH 09/30] setup_particles: cleanup get_n_in_cell

---
 src/include/setup_particles.hxx | 9 ++-------
 1 file changed, 2 insertions(+), 7 deletions(-)

diff --git a/src/include/setup_particles.hxx b/src/include/setup_particles.hxx
index cd0f82040e..c0074929a1 100644
--- a/src/include/setup_particles.hxx
+++ b/src/include/setup_particles.hxx
@@ -84,17 +84,12 @@ struct SetupParticles
 
   int get_n_in_cell(const psc_particle_npt& npt)
   {
+    static distribution::Uniform<float> dist{0, 1};
     if (npt.n == 0) {
       return 0;
     }
     if (fractional_n_particles_per_cell) {
-      int n_prts = npt.n / norm_.cori;
-      float rmndr = npt.n / norm_.cori - n_prts;
-      float ran = random() / ((float)RAND_MAX + 1);
-      if (ran < rmndr) {
-        n_prts++;
-      }
-      return n_prts;
+      return npt.n / norm_.cori + dist.get();
     }
     return npt.n / norm_.cori + .5;
   }

From d265c469359c6656e03e57c6a093382f0c4b2650 Mon Sep 17 00:00:00 2001
From: James <jm1667@wildcats.unh.edu>
Date: Sun, 13 Nov 2022 18:18:12 -0500
Subject: [PATCH 10/30] setup_particles: cleanup setupParticle

---
 src/include/setup_particles.hxx | 18 ++++++++----------
 1 file changed, 8 insertions(+), 10 deletions(-)

diff --git a/src/include/setup_particles.hxx b/src/include/setup_particles.hxx
index c0074929a1..40944c1b6e 100644
--- a/src/include/setup_particles.hxx
+++ b/src/include/setup_particles.hxx
@@ -156,21 +156,19 @@ struct SetupParticles
     assert(npt.kind >= 0 && npt.kind < kinds_.size());
     double m = kinds_[npt.kind].m;
 
-    double pxi = dist.get(npt.p[0], beta * std::sqrt(npt.T[0] / m));
-    double pyi = dist.get(npt.p[1], beta * std::sqrt(npt.T[1] / m));
-    double pzi = dist.get(npt.p[2], beta * std::sqrt(npt.T[2] / m));
+    Double3 p;
+    for (int i = 0; i < 3; i++)
+      p[i] = dist.get(npt.p[i], beta * std::sqrt(npt.T[i] / m));
 
     if (initial_momentum_gamma_correction) {
-      double gam;
-      if (sqr(pxi) + sqr(pyi) + sqr(pzi) < 1.) {
-        gam = 1. / sqrt(1. - sqr(pxi) - sqr(pyi) - sqr(pzi));
-        pxi *= gam;
-        pyi *= gam;
-        pzi *= gam;
+      double p_squared = sqr(p[0]) + sqr(p[1]) + sqr(p[2]);
+      if (p_squared < 1.) {
+        double gamma = 1. / sqrt(1. - p_squared);
+        p *= gamma;
       }
     }
 
-    return psc::particle::Inject{pos, {pxi, pyi, pzi}, wni, npt.kind, npt.tag};
+    return psc::particle::Inject{pos, p, wni, npt.kind, npt.tag};
   }
 
   // ----------------------------------------------------------------------

From b6670ed6a1e6d7089b25b07440e27560faa112eb Mon Sep 17 00:00:00 2001
From: James <jm1667@wildcats.unh.edu>
Date: Sun, 13 Nov 2022 19:04:04 -0500
Subject: [PATCH 11/30] setup_particles: mv p gen to psc_particle_np

---
 src/include/setup_particles.hxx | 55 +++++++++++++++++++++++----------
 1 file changed, 38 insertions(+), 17 deletions(-)

diff --git a/src/include/setup_particles.hxx b/src/include/setup_particles.hxx
index 40944c1b6e..e81b07f582 100644
--- a/src/include/setup_particles.hxx
+++ b/src/include/setup_particles.hxx
@@ -15,6 +15,36 @@ struct psc_particle_npt
   psc::particle::Tag tag;
 };
 
+struct psc_particle_np
+{
+  int kind;                   ///< particle kind
+  double n;                   ///< density
+  std::function<Double3()> p; ///< returns a random momentum
+  psc::particle::Tag tag;
+
+  psc_particle_np(const psc_particle_npt& npt, double m, double beta,
+                  bool initial_momentum_gamma_correction)
+    : kind{npt.kind}, n{npt.n}, tag{npt.tag}
+  {
+    p = [=]() {
+      static distribution::Normal<double> dist;
+
+      Double3 p;
+      for (int i = 0; i < 3; i++)
+        p[i] = dist.get(npt.p[i], beta * std::sqrt(npt.T[i] / m));
+
+      if (initial_momentum_gamma_correction) {
+        double p_squared = sqr(p[0]) + sqr(p[1]) + sqr(p[2]);
+        if (p_squared < 1.) {
+          double gamma = 1. / sqrt(1. - p_squared);
+          p *= gamma;
+        }
+      }
+      return p;
+    };
+  };
+};
+
 // ======================================================================
 // SetupParticles
 
@@ -150,25 +180,16 @@ struct SetupParticles
   psc::particle::Inject setupParticle(const psc_particle_npt& npt, Double3 pos,
                                       double wni)
   {
-    static distribution::Normal<double> dist;
-    double beta = norm_.beta;
-
     assert(npt.kind >= 0 && npt.kind < kinds_.size());
-    double m = kinds_[npt.kind].m;
-
-    Double3 p;
-    for (int i = 0; i < 3; i++)
-      p[i] = dist.get(npt.p[i], beta * std::sqrt(npt.T[i] / m));
-
-    if (initial_momentum_gamma_correction) {
-      double p_squared = sqr(p[0]) + sqr(p[1]) + sqr(p[2]);
-      if (p_squared < 1.) {
-        double gamma = 1. / sqrt(1. - p_squared);
-        p *= gamma;
-      }
-    }
+    return setupParticle(psc_particle_np(npt, kinds_[npt.kind].m, norm_.beta,
+                                         initial_momentum_gamma_correction),
+                         pos, wni);
+  }
 
-    return psc::particle::Inject{pos, p, wni, npt.kind, npt.tag};
+  psc::particle::Inject setupParticle(const psc_particle_np& np, Double3 pos,
+                                      double wni)
+  {
+    return psc::particle::Inject{pos, np.p(), wni, np.kind, np.tag};
   }
 
   // ----------------------------------------------------------------------

From 36f4a1474b9bb18de718328f0590d651af50f818 Mon Sep 17 00:00:00 2001
From: James <jm1667@wildcats.unh.edu>
Date: Sun, 13 Nov 2022 19:07:16 -0500
Subject: [PATCH 12/30] setup_particles: add helper npt -> np conversion

---
 src/include/setup_particles.hxx | 15 +++++++++++----
 1 file changed, 11 insertions(+), 4 deletions(-)

diff --git a/src/include/setup_particles.hxx b/src/include/setup_particles.hxx
index e81b07f582..2051e802d7 100644
--- a/src/include/setup_particles.hxx
+++ b/src/include/setup_particles.hxx
@@ -180,10 +180,7 @@ struct SetupParticles
   psc::particle::Inject setupParticle(const psc_particle_npt& npt, Double3 pos,
                                       double wni)
   {
-    assert(npt.kind >= 0 && npt.kind < kinds_.size());
-    return setupParticle(psc_particle_np(npt, kinds_[npt.kind].m, norm_.beta,
-                                         initial_momentum_gamma_correction),
-                         pos, wni);
+    return setupParticle(npt_to_np(npt), pos, wni);
   }
 
   psc::particle::Inject setupParticle(const psc_particle_np& np, Double3 pos,
@@ -192,6 +189,16 @@ struct SetupParticles
     return psc::particle::Inject{pos, np.p(), wni, np.kind, np.tag};
   }
 
+  // ----------------------------------------------------------------------
+  // npt_to_np
+
+  psc_particle_np npt_to_np(const psc_particle_npt& npt)
+  {
+    assert(npt.kind >= 0 && npt.kind < kinds_.size());
+    return psc_particle_np(npt, kinds_[npt.kind].m, norm_.beta,
+                           initial_momentum_gamma_correction);
+  }
+
   // ----------------------------------------------------------------------
   // setup_particles
 

From 44cfc84039de9d842cc1d2a814f0f559d8bb4c8d Mon Sep 17 00:00:00 2001
From: James <jm1667@wildcats.unh.edu>
Date: Sun, 13 Nov 2022 19:09:25 -0500
Subject: [PATCH 13/30] setup_particles: rm redundant setupParticle

---
 src/include/setup_particles.hxx | 8 +-------
 1 file changed, 1 insertion(+), 7 deletions(-)

diff --git a/src/include/setup_particles.hxx b/src/include/setup_particles.hxx
index 2051e802d7..ff2b9b5bc8 100644
--- a/src/include/setup_particles.hxx
+++ b/src/include/setup_particles.hxx
@@ -177,12 +177,6 @@ struct SetupParticles
   // ----------------------------------------------------------------------
   // setupParticle
 
-  psc::particle::Inject setupParticle(const psc_particle_npt& npt, Double3 pos,
-                                      double wni)
-  {
-    return setupParticle(npt_to_np(npt), pos, wni);
-  }
-
   psc::particle::Inject setupParticle(const psc_particle_np& np, Double3 pos,
                                       double wni)
   {
@@ -236,7 +230,7 @@ struct SetupParticles
                       } else {
                         wni = npt.n / (n_in_cell * norm_.cori);
                       }
-                      auto prt = setupParticle(npt, pos, wni);
+                      auto prt = setupParticle(npt_to_np(npt), pos, wni);
                       injector(prt);
                     }
                   });

From 4cd0ee8a7a96c1bda7c7f6060830ee42a3273e46 Mon Sep 17 00:00:00 2001
From: James <jm1667@wildcats.unh.edu>
Date: Sun, 13 Nov 2022 19:19:14 -0500
Subject: [PATCH 14/30] setup_particles: make & use np in op_cellwise

---
 src/include/setup_particles.hxx | 27 ++++++++++++++-------------
 1 file changed, 14 insertions(+), 13 deletions(-)

diff --git a/src/include/setup_particles.hxx b/src/include/setup_particles.hxx
index ff2b9b5bc8..e810bb134e 100644
--- a/src/include/setup_particles.hxx
+++ b/src/include/setup_particles.hxx
@@ -112,22 +112,22 @@ struct SetupParticles
   //
   // helper function for partition / particle setup
 
-  int get_n_in_cell(const psc_particle_npt& npt)
+  int get_n_in_cell(const psc_particle_np& np)
   {
     static distribution::Uniform<float> dist{0, 1};
-    if (npt.n == 0) {
+    if (np.n == 0) {
       return 0;
     }
     if (fractional_n_particles_per_cell) {
-      return npt.n / norm_.cori + dist.get();
+      return np.n / norm_.cori + dist.get();
     }
-    return npt.n / norm_.cori + .5;
+    return np.n / norm_.cori + .5;
   }
 
   // ----------------------------------------------------------------------
   // op_cellwise
   // Performs a given operation in each cell.
-  // op signature: (int n_in_cell, npt, Double3 pos) -> void
+  // op signature: (int n_in_cell, np, Double3 pos) -> void
 
   template <typename OpFunc>
   void op_cellwise(const Grid_t& grid, int patch, InitNptFunc init_npt,
@@ -156,18 +156,19 @@ struct SetupParticles
               npt.kind = pop;
             }
             init_npt(pop, pos, patch, index, npt);
+            auto np = npt_to_np(npt);
 
             int n_in_cell;
             if (pop != neutralizing_population) {
-              n_in_cell = get_n_in_cell(npt);
-              n_q_in_cell += kinds_[npt.kind].q * n_in_cell;
+              n_in_cell = get_n_in_cell(np);
+              n_q_in_cell += kinds_[np.kind].q * n_in_cell;
             } else {
               // FIXME, should handle the case where not the last population
               // is neutralizing
               assert(neutralizing_population == n_populations_ - 1);
-              n_in_cell = -n_q_in_cell / kinds_[npt.kind].q;
+              n_in_cell = -n_q_in_cell / kinds_[np.kind].q;
             }
-            op(n_in_cell, npt, pos);
+            op(n_in_cell, np, pos);
           }
         }
       }
@@ -222,15 +223,15 @@ struct SetupParticles
       auto injector = inj[p];
 
       op_cellwise(grid, p, init_npt,
-                  [&](int n_in_cell, psc_particle_npt& npt, Double3& pos) {
+                  [&](int n_in_cell, psc_particle_np& np, Double3& pos) {
                     for (int cnt = 0; cnt < n_in_cell; cnt++) {
                       real_t wni;
                       if (fractional_n_particles_per_cell) {
                         wni = 1.;
                       } else {
-                        wni = npt.n / (n_in_cell * norm_.cori);
+                        wni = np.n / (n_in_cell * norm_.cori);
                       }
-                      auto prt = setupParticle(npt_to_np(npt), pos, wni);
+                      auto prt = setupParticle(np, pos, wni);
                       injector(prt);
                     }
                   });
@@ -248,7 +249,7 @@ struct SetupParticles
 
     for (int p = 0; p < grid.n_patches(); ++p) {
       op_cellwise(grid, p, init_npt,
-                  [&](int n_in_cell, psc_particle_npt&, Double3&) {
+                  [&](int n_in_cell, psc_particle_np&, Double3&) {
                     n_prts_by_patch[p] += n_in_cell;
                   });
     }

From 36417b426e7bc918b3a4250b5985cc3384288c94 Mon Sep 17 00:00:00 2001
From: James <jm1667@wildcats.unh.edu>
Date: Mon, 14 Nov 2022 09:43:43 -0500
Subject: [PATCH 15/30] setup_particles: mv maxwellian creation new func

---
 src/include/setup_particles.hxx | 51 +++++++++++++++++----------------
 1 file changed, 27 insertions(+), 24 deletions(-)

diff --git a/src/include/setup_particles.hxx b/src/include/setup_particles.hxx
index e810bb134e..660c1d861e 100644
--- a/src/include/setup_particles.hxx
+++ b/src/include/setup_particles.hxx
@@ -21,28 +21,6 @@ struct psc_particle_np
   double n;                   ///< density
   std::function<Double3()> p; ///< returns a random momentum
   psc::particle::Tag tag;
-
-  psc_particle_np(const psc_particle_npt& npt, double m, double beta,
-                  bool initial_momentum_gamma_correction)
-    : kind{npt.kind}, n{npt.n}, tag{npt.tag}
-  {
-    p = [=]() {
-      static distribution::Normal<double> dist;
-
-      Double3 p;
-      for (int i = 0; i < 3; i++)
-        p[i] = dist.get(npt.p[i], beta * std::sqrt(npt.T[i] / m));
-
-      if (initial_momentum_gamma_correction) {
-        double p_squared = sqr(p[0]) + sqr(p[1]) + sqr(p[2]);
-        if (p_squared < 1.) {
-          double gamma = 1. / sqrt(1. - p_squared);
-          p *= gamma;
-        }
-      }
-      return p;
-    };
-  };
 };
 
 // ======================================================================
@@ -188,10 +166,35 @@ struct SetupParticles
   // npt_to_np
 
   psc_particle_np npt_to_np(const psc_particle_npt& npt)
+  {
+    return psc_particle_np{npt.kind, npt.n, createMaxwellian(npt), npt.tag};
+  }
+
+  // ----------------------------------------------------------------------
+  // createMaxwellian
+
+  std::function<Double3()> createMaxwellian(const psc_particle_npt& npt)
   {
     assert(npt.kind >= 0 && npt.kind < kinds_.size());
-    return psc_particle_np(npt, kinds_[npt.kind].m, norm_.beta,
-                           initial_momentum_gamma_correction);
+    double beta = norm_.beta;
+    double m = kinds_[npt.kind].m;
+
+    return [=]() {
+      static distribution::Normal<double> dist;
+
+      Double3 p;
+      for (int i = 0; i < 3; i++)
+        p[i] = dist.get(npt.p[i], beta * std::sqrt(npt.T[i] / m));
+
+      if (initial_momentum_gamma_correction) {
+        double p_squared = sqr(p[0]) + sqr(p[1]) + sqr(p[2]);
+        if (p_squared < 1.) {
+          double gamma = 1. / sqrt(1. - p_squared);
+          p *= gamma;
+        }
+      }
+      return p;
+    };
   }
 
   // ----------------------------------------------------------------------

From 56bb0ae6280b6a9b75f1cfd532be001d4eb43932 Mon Sep 17 00:00:00 2001
From: James <jm1667@wildcats.unh.edu>
Date: Mon, 14 Nov 2022 09:45:33 -0500
Subject: [PATCH 16/30] setup_particles: rm npt_to_np

---
 src/include/setup_particles.hxx | 10 +---------
 1 file changed, 1 insertion(+), 9 deletions(-)

diff --git a/src/include/setup_particles.hxx b/src/include/setup_particles.hxx
index 660c1d861e..cc2626b2ad 100644
--- a/src/include/setup_particles.hxx
+++ b/src/include/setup_particles.hxx
@@ -134,7 +134,7 @@ struct SetupParticles
               npt.kind = pop;
             }
             init_npt(pop, pos, patch, index, npt);
-            auto np = npt_to_np(npt);
+            psc_particle_np np{npt.kind, npt.n, createMaxwellian(npt), npt.tag};
 
             int n_in_cell;
             if (pop != neutralizing_population) {
@@ -162,14 +162,6 @@ struct SetupParticles
     return psc::particle::Inject{pos, np.p(), wni, np.kind, np.tag};
   }
 
-  // ----------------------------------------------------------------------
-  // npt_to_np
-
-  psc_particle_np npt_to_np(const psc_particle_npt& npt)
-  {
-    return psc_particle_np{npt.kind, npt.n, createMaxwellian(npt), npt.tag};
-  }
-
   // ----------------------------------------------------------------------
   // createMaxwellian
 

From 52fb45438b32c674c2ab4109b27f58332e2dd74a Mon Sep 17 00:00:00 2001
From: James <jm1667@wildcats.unh.edu>
Date: Mon, 14 Nov 2022 09:59:05 -0500
Subject: [PATCH 17/30] setup_particles: introduce InitNpFunc

---
 src/include/setup_particles.hxx | 52 ++++++++++++++++++++++++++++-----
 1 file changed, 45 insertions(+), 7 deletions(-)

diff --git a/src/include/setup_particles.hxx b/src/include/setup_particles.hxx
index cc2626b2ad..141dde4584 100644
--- a/src/include/setup_particles.hxx
+++ b/src/include/setup_particles.hxx
@@ -68,6 +68,29 @@ private:
   std::function<void(int, Double3, int, Int3, psc_particle_npt&)> func;
 };
 
+struct InitNpFunc
+{
+  // Initialize particles according to an arbitrary momentum disribution.
+  // (int kind, Double3 pos, int patch, Int3 idx, psc_particle_np& np) -> void
+  template <
+    typename INIT_NP,
+    std::enable_if_t<std::is_convertible<
+                       INIT_NP, std::function<void(int, Double3, int, Int3,
+                                                   psc_particle_np&)>>::value,
+                     bool> = true>
+  InitNpFunc(INIT_NP init_np) : func{init_np}
+  {}
+
+  void operator()(int kind, Double3 pos, int patch, Int3 index,
+                  psc_particle_np& np)
+  {
+    func(kind, pos, patch, index, np);
+  }
+
+private:
+  std::function<void(int, Double3, int, Int3, psc_particle_np&)> func;
+};
+
 template <typename MP>
 struct SetupParticles
 {
@@ -108,7 +131,7 @@ struct SetupParticles
   // op signature: (int n_in_cell, np, Double3 pos) -> void
 
   template <typename OpFunc>
-  void op_cellwise(const Grid_t& grid, int patch, InitNptFunc init_npt,
+  void op_cellwise(const Grid_t& grid, int patch, InitNpFunc init_np,
                    OpFunc&& op)
   {
     Int3 ilo = Int3{}, ihi = grid.ldims;
@@ -129,12 +152,11 @@ struct SetupParticles
 
           int n_q_in_cell = 0;
           for (int pop = 0; pop < n_populations_; pop++) {
-            psc_particle_npt npt{};
+            psc_particle_np np{};
             if (pop < kinds_.size()) {
-              npt.kind = pop;
+              np.kind = pop;
             }
-            init_npt(pop, pos, patch, index, npt);
-            psc_particle_np np{npt.kind, npt.n, createMaxwellian(npt), npt.tag};
+            init_np(pop, pos, patch, index, np);
 
             int n_in_cell;
             if (pop != neutralizing_population) {
@@ -197,6 +219,22 @@ struct SetupParticles
     setupParticles(mprts, init_npt);
   }
 
+  // ----------------------------------------------------------------------
+  // initNpt_to_initNp
+
+  InitNpFunc initNpt_to_initNp(InitNptFunc& init_npt)
+  {
+    return InitNpFunc(
+      [&](int kind, Double3 pos, int patch, Int3 idx, psc_particle_np& np) {
+        psc_particle_npt npt{};
+        npt.kind = np.kind;
+        init_npt(kind, pos, patch, idx, npt);
+        np.n = npt.n;
+        np.p = createMaxwellian(npt);
+        np.tag = npt.tag;
+      });
+  }
+
   // ----------------------------------------------------------------------
   // setupParticles
 
@@ -217,7 +255,7 @@ struct SetupParticles
     for (int p = 0; p < mprts.n_patches(); ++p) {
       auto injector = inj[p];
 
-      op_cellwise(grid, p, init_npt,
+      op_cellwise(grid, p, initNpt_to_initNp(init_npt),
                   [&](int n_in_cell, psc_particle_np& np, Double3& pos) {
                     for (int cnt = 0; cnt < n_in_cell; cnt++) {
                       real_t wni;
@@ -243,7 +281,7 @@ struct SetupParticles
     std::vector<uint> n_prts_by_patch(grid.n_patches());
 
     for (int p = 0; p < grid.n_patches(); ++p) {
-      op_cellwise(grid, p, init_npt,
+      op_cellwise(grid, p, initNpt_to_initNp(init_npt),
                   [&](int n_in_cell, psc_particle_np&, Double3&) {
                     n_prts_by_patch[p] += n_in_cell;
                   });

From 971095f16fe56e5f9e39a7640fecbfa3a8f32c2a Mon Sep 17 00:00:00 2001
From: James <jm1667@wildcats.unh.edu>
Date: Mon, 14 Nov 2022 10:16:34 -0500
Subject: [PATCH 18/30] setup_particles: add methods accepting InitNpFuncs

---
 src/include/psc.hxx             | 23 ++++++++++++-----------
 src/include/setup_particles.hxx | 19 +++++++++++++++++--
 2 files changed, 29 insertions(+), 13 deletions(-)

diff --git a/src/include/psc.hxx b/src/include/psc.hxx
index cbb1ca0d4e..8b2c3b812b 100644
--- a/src/include/psc.hxx
+++ b/src/include/psc.hxx
@@ -782,27 +782,28 @@ Psc<PscConfig, Diagnostics, InjectParticles> makePscIntegrator(
 // ======================================================================
 // partitionAndSetupParticles
 
-template <typename SetupParticles, typename Balance, typename Mparticles>
+template <typename SetupParticles, typename Balance, typename Mparticles,
+          typename InitFunc>
 void partitionAndSetupParticles(SetupParticles& setup_particles,
                                 Balance& balance, Grid_t*& grid_ptr,
-                                Mparticles& mprts, InitNptFunc init_npt)
+                                Mparticles& mprts, InitFunc init_np)
 {
-  partitionParticles(setup_particles, balance, grid_ptr, mprts, init_npt);
-  setupParticles(setup_particles, mprts, init_npt);
+  partitionParticles(setup_particles, balance, grid_ptr, mprts, init_np);
+  setupParticles(setup_particles, mprts, init_np);
 }
 
 // ----------------------------------------------------------------------
 // partitionParticles
 
-template <typename SetupParticles, typename Balance, typename Mparticles>
+template <typename SetupParticles, typename Balance, typename Mparticles,
+          typename InitFunc>
 void partitionParticles(SetupParticles& setup_particles, Balance& balance,
-                        Grid_t*& grid_ptr, Mparticles& mprts,
-                        InitNptFunc init_npt)
+                        Grid_t*& grid_ptr, Mparticles& mprts, InitFunc init_np)
 {
   auto comm = grid_ptr->comm();
   mpi_printf(comm, "**** Partitioning...\n");
 
-  auto n_prts_by_patch = setup_particles.partition(*grid_ptr, init_npt);
+  auto n_prts_by_patch = setup_particles.partition(*grid_ptr, init_np);
 
   balance.initial(grid_ptr, n_prts_by_patch);
   // !!! FIXME! grid is now invalid
@@ -815,12 +816,12 @@ void partitionParticles(SetupParticles& setup_particles, Balance& balance,
 // ----------------------------------------------------------------------
 // setupParticles
 
-template <typename SetupParticles, typename Mparticles>
+template <typename SetupParticles, typename Mparticles, typename InitFunc>
 void setupParticles(SetupParticles& setup_particles, Mparticles& mprts,
-                    InitNptFunc init_npt)
+                    InitFunc init_np)
 {
   mpi_printf(MPI_COMM_WORLD, "**** Setting up particles...\n");
-  setup_particles.setupParticles(mprts, init_npt);
+  setup_particles.setupParticles(mprts, init_np);
 }
 
 // ======================================================================
diff --git a/src/include/setup_particles.hxx b/src/include/setup_particles.hxx
index 141dde4584..343cba9ef2 100644
--- a/src/include/setup_particles.hxx
+++ b/src/include/setup_particles.hxx
@@ -214,6 +214,11 @@ struct SetupParticles
   // ----------------------------------------------------------------------
   // setup_particles
 
+  void operator()(Mparticles& mprts, InitNpFunc init_np)
+  {
+    setupParticles(mprts, init_np);
+  }
+
   void operator()(Mparticles& mprts, InitNptFunc init_npt)
   {
     setupParticles(mprts, init_npt);
@@ -239,6 +244,11 @@ struct SetupParticles
   // setupParticles
 
   void setupParticles(Mparticles& mprts, InitNptFunc init_npt)
+  {
+    setupParticles(mprts, initNpt_to_initNp(init_npt));
+  }
+
+  void setupParticles(Mparticles& mprts, InitNpFunc init_np)
   {
     static int pr;
     if (!pr) {
@@ -255,7 +265,7 @@ struct SetupParticles
     for (int p = 0; p < mprts.n_patches(); ++p) {
       auto injector = inj[p];
 
-      op_cellwise(grid, p, initNpt_to_initNp(init_npt),
+      op_cellwise(grid, p, init_np,
                   [&](int n_in_cell, psc_particle_np& np, Double3& pos) {
                     for (int cnt = 0; cnt < n_in_cell; cnt++) {
                       real_t wni;
@@ -277,11 +287,16 @@ struct SetupParticles
   // partition
 
   std::vector<uint> partition(const Grid_t& grid, InitNptFunc init_npt)
+  {
+    return partition(grid, initNpt_to_initNp(init_npt));
+  }
+
+  std::vector<uint> partition(const Grid_t& grid, InitNpFunc init_np)
   {
     std::vector<uint> n_prts_by_patch(grid.n_patches());
 
     for (int p = 0; p < grid.n_patches(); ++p) {
-      op_cellwise(grid, p, initNpt_to_initNp(init_npt),
+      op_cellwise(grid, p, init_np,
                   [&](int n_in_cell, psc_particle_np&, Double3&) {
                     n_prts_by_patch[p] += n_in_cell;
                   });

From 3b1c990423502bddd2cf95032f30065dc7ac16d4 Mon Sep 17 00:00:00 2001
From: James <jm1667@wildcats.unh.edu>
Date: Wed, 16 Nov 2022 11:04:45 -0500
Subject: [PATCH 19/30] distribution: impl cdf inverter

---
 src/include/distribution.hxx | 50 ++++++++++++++++++++++++++++++++++++
 1 file changed, 50 insertions(+)

diff --git a/src/include/distribution.hxx b/src/include/distribution.hxx
index d0cfa7e63c..3f0897d59a 100644
--- a/src/include/distribution.hxx
+++ b/src/include/distribution.hxx
@@ -1,4 +1,5 @@
 #include <random>
+#include <vector>
 #include <cmath>
 
 namespace distribution
@@ -60,4 +61,53 @@ private:
   bool has_next_value;
 };
 
+// ======================================================================
+// InvertedCdf
+
+template <typename Real>
+struct InvertedCdf : public Distribution<Real>
+{
+  InvertedCdf(std::function<Real(Real)> cdf, int n_points = 100,
+              Real eps = .0001)
+  {
+    Real xmin = -1, xmax = 1;
+
+    while (cdf(xmin) < eps)
+      xmin += .5;
+    while (cdf(xmin) > eps)
+      xmin -= 1 / 16.;
+    while (cdf(xmax) > 1 - eps)
+      xmax -= .5;
+    while (cdf(xmax) < 1 - eps)
+      xmax += 1 / 16.;
+
+    x.push_back(xmin);
+    this->cdf.push_back(0);
+    Real dx = (xmax - xmin) / (n_points - 1);
+    for (int i = 1; i < n_points - 1; i++) {
+      x.push_back(xmin + i * dx);
+      this->cdf.push_back(cdf(xmin + i * dx));
+    }
+    x.push_back(xmax);
+    this->cdf.push_back(1);
+  }
+
+  Real get()
+  {
+    Real u = uniform.get();
+    int guess_idx = cdf.size() * u;
+
+    while (cdf[guess_idx] > u)
+      --guess_idx;
+    while (cdf[guess_idx] < u)
+      ++guess_idx;
+
+    return x[guess_idx];
+  }
+
+private:
+  std::vector<Real> cdf, x;
+  Uniform<Real> uniform{0, 1};
+};
+
 } // namespace distribution
\ No newline at end of file

From 651885fd34d915d3762ef6da6d88f087ae941394 Mon Sep 17 00:00:00 2001
From: James <jm1667@wildcats.unh.edu>
Date: Wed, 16 Nov 2022 11:05:29 -0500
Subject: [PATCH 20/30] distribution: interpolate in get()

---
 src/include/distribution.hxx | 22 +++++++++++++++++-----
 1 file changed, 17 insertions(+), 5 deletions(-)

diff --git a/src/include/distribution.hxx b/src/include/distribution.hxx
index 3f0897d59a..84a5699725 100644
--- a/src/include/distribution.hxx
+++ b/src/include/distribution.hxx
@@ -97,12 +97,24 @@ struct InvertedCdf : public Distribution<Real>
     Real u = uniform.get();
     int guess_idx = cdf.size() * u;
 
-    while (cdf[guess_idx] > u)
-      --guess_idx;
-    while (cdf[guess_idx] < u)
-      ++guess_idx;
+    while (u < cdf[guess_idx]) {
+      --guess_idx; // decrease cdf
+    }
+    while (u > cdf[guess_idx]) {
+      ++guess_idx; // increase cdf
+    }
+
+    Real d_cdf = cdf[guess_idx] - cdf[guess_idx - 1];
+    if (d_cdf == 0.)
+      return x[guess_idx];
+
+    // now: cdf[guess_idx-1] <= u < cdf[guess_idx]
+    // so linearly interpolate between points
+
+    Real w1 = (cdf[guess_idx] - u) / d_cdf;
+    Real w2 = (u - cdf[guess_idx - 1]) / d_cdf;
 
-    return x[guess_idx];
+    return w1 * x[guess_idx - 1] + w2 * x[guess_idx];
   }
 
 private:

From 67a7a127686da1721daf628d9e8ab8334632dbeb Mon Sep 17 00:00:00 2001
From: James <jm1667@wildcats.unh.edu>
Date: Wed, 16 Nov 2022 11:06:35 -0500
Subject: [PATCH 21/30] distribution: limit n_iter

---
 src/include/distribution.hxx | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/include/distribution.hxx b/src/include/distribution.hxx
index 84a5699725..cbe06f8bfb 100644
--- a/src/include/distribution.hxx
+++ b/src/include/distribution.hxx
@@ -68,17 +68,17 @@ template <typename Real>
 struct InvertedCdf : public Distribution<Real>
 {
   InvertedCdf(std::function<Real(Real)> cdf, int n_points = 100,
-              Real eps = .0001)
+              Real eps = .0001, int max_n_iter = 50)
   {
     Real xmin = -1, xmax = 1;
 
-    while (cdf(xmin) < eps)
+    for (int n_iter = 0; cdf(xmin) < eps && n_iter < max_n_iter; ++n_iter)
       xmin += .5;
-    while (cdf(xmin) > eps)
+    for (int n_iter = 0; cdf(xmin) > eps && n_iter < max_n_iter; ++n_iter)
       xmin -= 1 / 16.;
-    while (cdf(xmax) > 1 - eps)
+    for (int n_iter = 0; cdf(xmax) > 1 - eps && n_iter < max_n_iter; ++n_iter)
       xmax -= .5;
-    while (cdf(xmax) < 1 - eps)
+    for (int n_iter = 0; cdf(xmax) < 1 - eps && n_iter < max_n_iter; ++n_iter)
       xmax += 1 / 16.;
 
     x.push_back(xmin);

From 5773df34558f43d51af18c6a6a2d343bd3fd6d7d Mon Sep 17 00:00:00 2001
From: James <jm1667@wildcats.unh.edu>
Date: Wed, 16 Nov 2022 12:24:02 -0500
Subject: [PATCH 22/30] distribution: fix seeds

---
 src/include/distribution.hxx | 41 ++++++++++++++++++------------------
 1 file changed, 20 insertions(+), 21 deletions(-)

diff --git a/src/include/distribution.hxx b/src/include/distribution.hxx
index cbe06f8bfb..afc75a896e 100644
--- a/src/include/distribution.hxx
+++ b/src/include/distribution.hxx
@@ -1,6 +1,8 @@
-#include <random>
 #include <vector>
 #include <cmath>
+#include "stdlib.h"
+#include <random>
+#include <chrono>
 
 namespace distribution
 {
@@ -21,11 +23,15 @@ struct Distribution
 template <typename Real>
 struct Uniform : public Distribution<Real>
 {
-  Uniform(Real min, Real max) : dist{min, max} {}
-  Real get() override { return dist(generator); }
+  Uniform(Real min, Real max)
+    : dist(min, max),
+      gen(std::chrono::system_clock::now().time_since_epoch().count())
+  {}
+
+  Real get() override { return dist(gen); }
 
 private:
-  std::default_random_engine generator;
+  std::default_random_engine gen;
   std::uniform_real_distribution<Real> dist;
 };
 
@@ -35,30 +41,23 @@ private:
 template <typename Real>
 struct Normal : public Distribution<Real>
 {
-  Normal(Real mean, Real stdev) : mean{mean}, stdev{stdev} {}
+  Normal(Real mean, Real stdev)
+    : dist(mean, stdev),
+      gen(std::chrono::system_clock::now().time_since_epoch().count())
+  {}
   Normal() : Normal(0, 1) {}
 
-  Real get() override { return get(mean, stdev); }
+  Real get() override { return dist(gen); }
+  // FIXME remove me, or make standalone func
   Real get(Real mean, Real stdev)
   {
-    if (has_next_value) {
-      has_next_value = false;
-      return next_value;
-    }
-
-    // Box-Muller transform to generate 2 independent values at once
-    Real r = stdev * std::sqrt(-2 * std::log(1 - uniform.get()));
-    Real theta = uniform.get() * 2 * M_PI;
-
-    next_value = r * std::cos(theta) + mean;
-    return r * std::sin(theta) + mean;
+    // should be possible to pass params to existing dist
+    return std::normal_distribution<Real>(mean, stdev)(gen);
   }
 
 private:
-  Uniform<Real> uniform{0, 1};
-  Real mean, stdev;
-  Real next_value;
-  bool has_next_value;
+  std::default_random_engine gen;
+  std::normal_distribution<Real> dist;
 };
 
 // ======================================================================

From c6299309e71b9eb311416167590d3426c46e581a Mon Sep 17 00:00:00 2001
From: James <jm1667@wildcats.unh.edu>
Date: Mon, 19 Dec 2022 14:32:23 -0500
Subject: [PATCH 23/30] distribution: fix finding pdf support

this version is a disaster of repeated code, but hopefully works
---
 src/include/distribution.hxx | 82 +++++++++++++++++++++++++++++++-----
 1 file changed, 71 insertions(+), 11 deletions(-)

diff --git a/src/include/distribution.hxx b/src/include/distribution.hxx
index afc75a896e..c59ee0ad5b 100644
--- a/src/include/distribution.hxx
+++ b/src/include/distribution.hxx
@@ -63,22 +63,82 @@ private:
 // ======================================================================
 // InvertedCdf
 
+template <typename Real>
+void findPdfSupport(std::function<Real(Real)>& cdf, Real& xmin, Real& xmax,
+                    Real tolerance, int max_n_iter)
+{
+  Real xmin_ub = 0, xmax_lb = 0;
+  Real sf;
+  int n_iter;
+
+  // find where cfd > tolerance
+  n_iter = 0;
+  sf = 1;
+  while (n_iter < max_n_iter && cdf(xmin_ub) < tolerance) {
+    xmin_ub += (sf *= 1.5);
+    n_iter++;
+  }
+
+  // find where cfd < 1-tolerance
+  n_iter = 0;
+  sf = 1;
+  while (n_iter < max_n_iter && cdf(xmax_lb) > 1 - tolerance) {
+    xmax_lb -= (sf *= 1.5);
+    n_iter++;
+  }
+
+  Real xmin_lb = xmin_ub, xmax_ub = xmax_lb;
+
+  // find where cdf < tolerance
+  n_iter = 0;
+  sf = 1;
+  while (n_iter < max_n_iter && cdf(xmin_lb) > tolerance) {
+    xmin_lb -= (sf *= 1.5);
+    n_iter++;
+  }
+
+  // find where cdf > 1-tolerance
+  n_iter = 0;
+  sf = 1;
+  while (n_iter < max_n_iter && cdf(xmax_ub) < 1 - tolerance) {
+    xmax_ub += (sf *= 1.5);
+    n_iter++;
+  }
+
+  // bisect search for where cdf = tolerance
+  n_iter = 0;
+  xmin = (xmin_lb + xmin_ub) / 2.;
+  xmax = (xmax_lb + xmax_ub) / 2.;
+  // logDiff measures the precision of the range relative to that of each bound
+  Real logDiff =
+    std::log2((xmax - xmin) / std::max(xmax_ub - xmax_lb, xmin_ub - xmin_lb));
+  // if logDiff were recalculated, it would go up by ~1 with each iteration, so
+  // logDiff(t) ~= n_iter(t) + logDiff(0)
+  while (n_iter < max_n_iter && n_iter + logDiff < 10) {
+    if (cdf(xmin) < tolerance) {
+      xmin_lb = xmin;
+    } else {
+      xmin_ub = xmin;
+    }
+    if (cdf(xmax) > 1 - tolerance) {
+      xmax_ub = xmax;
+    } else {
+      xmax_lb = xmax;
+    }
+    xmin = (xmin_lb + xmin_ub) / 2.;
+    xmax = (xmax_lb + xmax_ub) / 2.;
+    n_iter++;
+  }
+}
+
 template <typename Real>
 struct InvertedCdf : public Distribution<Real>
 {
   InvertedCdf(std::function<Real(Real)> cdf, int n_points = 100,
-              Real eps = .0001, int max_n_iter = 50)
+              Real tolerance = .0001, int max_n_iter = 50)
   {
-    Real xmin = -1, xmax = 1;
-
-    for (int n_iter = 0; cdf(xmin) < eps && n_iter < max_n_iter; ++n_iter)
-      xmin += .5;
-    for (int n_iter = 0; cdf(xmin) > eps && n_iter < max_n_iter; ++n_iter)
-      xmin -= 1 / 16.;
-    for (int n_iter = 0; cdf(xmax) > 1 - eps && n_iter < max_n_iter; ++n_iter)
-      xmax -= .5;
-    for (int n_iter = 0; cdf(xmax) < 1 - eps && n_iter < max_n_iter; ++n_iter)
-      xmax += 1 / 16.;
+    Real xmin, xmax;
+    findPdfSupport(cdf, xmin, xmax, tolerance, max_n_iter);
 
     x.push_back(xmin);
     this->cdf.push_back(0);

From 88445e85f4463b487359de3e0d746b09935da269 Mon Sep 17 00:00:00 2001
From: James <jm1667@wildcats.unh.edu>
Date: Wed, 24 May 2023 19:01:32 -0400
Subject: [PATCH 24/30] distribution: add default Uniform

---
 src/include/distribution.hxx | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/include/distribution.hxx b/src/include/distribution.hxx
index c59ee0ad5b..737f10b02a 100644
--- a/src/include/distribution.hxx
+++ b/src/include/distribution.hxx
@@ -28,6 +28,8 @@ struct Uniform : public Distribution<Real>
       gen(std::chrono::system_clock::now().time_since_epoch().count())
   {}
 
+  Uniform() : Uniform(0, 1) {}
+
   Real get() override { return dist(gen); }
 
 private:

From fe58c77caaa23561b8c5beeeecbda12802c9f5b3 Mon Sep 17 00:00:00 2001
From: James <jm1667@wildcats.unh.edu>
Date: Wed, 24 May 2023 19:01:45 -0400
Subject: [PATCH 25/30] distribution: fix includes

---
 src/include/distribution.hxx | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/src/include/distribution.hxx b/src/include/distribution.hxx
index 737f10b02a..aa462ea2ca 100644
--- a/src/include/distribution.hxx
+++ b/src/include/distribution.hxx
@@ -1,8 +1,11 @@
+#pragma once
+
 #include <vector>
 #include <cmath>
 #include "stdlib.h"
 #include <random>
 #include <chrono>
+#include <functional>
 
 namespace distribution
 {

From 0545dea932579604508c969c371c624056ecca51 Mon Sep 17 00:00:00 2001
From: James <jm1667@wildcats.unh.edu>
Date: Wed, 24 May 2023 19:13:51 -0400
Subject: [PATCH 26/30] distribution: rm base class

---
 src/include/distribution.hxx | 20 +++++---------------
 1 file changed, 5 insertions(+), 15 deletions(-)

diff --git a/src/include/distribution.hxx b/src/include/distribution.hxx
index aa462ea2ca..6bd685ab91 100644
--- a/src/include/distribution.hxx
+++ b/src/include/distribution.hxx
@@ -10,21 +10,11 @@
 namespace distribution
 {
 
-// ======================================================================
-// Distribution
-// Base class for distributions
-
-template <typename Real>
-struct Distribution
-{
-  virtual Real get() = 0;
-};
-
 // ======================================================================
 // Uniform
 
 template <typename Real>
-struct Uniform : public Distribution<Real>
+struct Uniform
 {
   Uniform(Real min, Real max)
     : dist(min, max),
@@ -33,7 +23,7 @@ struct Uniform : public Distribution<Real>
 
   Uniform() : Uniform(0, 1) {}
 
-  Real get() override { return dist(gen); }
+  Real get() { return dist(gen); }
 
 private:
   std::default_random_engine gen;
@@ -44,7 +34,7 @@ private:
 // Normal
 
 template <typename Real>
-struct Normal : public Distribution<Real>
+struct Normal
 {
   Normal(Real mean, Real stdev)
     : dist(mean, stdev),
@@ -52,7 +42,7 @@ struct Normal : public Distribution<Real>
   {}
   Normal() : Normal(0, 1) {}
 
-  Real get() override { return dist(gen); }
+  Real get() { return dist(gen); }
   // FIXME remove me, or make standalone func
   Real get(Real mean, Real stdev)
   {
@@ -137,7 +127,7 @@ void findPdfSupport(std::function<Real(Real)>& cdf, Real& xmin, Real& xmax,
 }
 
 template <typename Real>
-struct InvertedCdf : public Distribution<Real>
+struct InvertedCdf
 {
   InvertedCdf(std::function<Real(Real)> cdf, int n_points = 100,
               Real tolerance = .0001, int max_n_iter = 50)

From fcffa24247af1a49fbc024663d2ee59d81a0862d Mon Sep 17 00:00:00 2001
From: James <jm1667@wildcats.unh.edu>
Date: Wed, 24 May 2023 19:20:44 -0400
Subject: [PATCH 27/30] distribution, *: rename to rng

---
 src/include/{distribution.hxx => rng.hxx} | 4 ++--
 src/include/setup_particles.hxx           | 6 +++---
 2 files changed, 5 insertions(+), 5 deletions(-)
 rename src/include/{distribution.hxx => rng.hxx} (98%)

diff --git a/src/include/distribution.hxx b/src/include/rng.hxx
similarity index 98%
rename from src/include/distribution.hxx
rename to src/include/rng.hxx
index 6bd685ab91..13e24b22e1 100644
--- a/src/include/distribution.hxx
+++ b/src/include/rng.hxx
@@ -7,7 +7,7 @@
 #include <chrono>
 #include <functional>
 
-namespace distribution
+namespace rng
 {
 
 // ======================================================================
@@ -176,4 +176,4 @@ private:
   Uniform<Real> uniform{0, 1};
 };
 
-} // namespace distribution
\ No newline at end of file
+} // namespace rng
\ No newline at end of file
diff --git a/src/include/setup_particles.hxx b/src/include/setup_particles.hxx
index 343cba9ef2..03fd233e1b 100644
--- a/src/include/setup_particles.hxx
+++ b/src/include/setup_particles.hxx
@@ -4,7 +4,7 @@
 #include <functional>
 #include <type_traits>
 #include "centering.hxx"
-#include "distribution.hxx"
+#include "rng.hxx"
 
 struct psc_particle_npt
 {
@@ -115,7 +115,7 @@ struct SetupParticles
 
   int get_n_in_cell(const psc_particle_np& np)
   {
-    static distribution::Uniform<float> dist{0, 1};
+    static rng::Uniform<float> dist{0, 1};
     if (np.n == 0) {
       return 0;
     }
@@ -194,7 +194,7 @@ struct SetupParticles
     double m = kinds_[npt.kind].m;
 
     return [=]() {
-      static distribution::Normal<double> dist;
+      static rng::Normal<double> dist;
 
       Double3 p;
       for (int i = 0; i < 3; i++)

From 01804a43940dd181dd79ba108007044c0c447221 Mon Sep 17 00:00:00 2001
From: liangwang0734 <liang.wang.phys@gmail.com>
Date: Fri, 30 Jun 2023 22:46:59 -0400
Subject: [PATCH 28/30] fix formatting

---
 src/include/checkpoint.hxx | 2 --
 1 file changed, 2 deletions(-)

diff --git a/src/include/checkpoint.hxx b/src/include/checkpoint.hxx
index d2b140d45b..cb44ecaeb5 100644
--- a/src/include/checkpoint.hxx
+++ b/src/include/checkpoint.hxx
@@ -7,7 +7,6 @@
 #include "particles_simple.inl"
 #include <kg/io.h>
 
-
 // ----------------------------------------------------------------------
 // read_checkpoint
 //
@@ -120,7 +119,6 @@ public:
     prof_stop(pr);
   }
 
-
 private:
   int interval_; // write checkpoint every so many steps
   bool first_time_ = true;

From 97e84ada7e6c81545c8ef62511717a60ab64f3dc Mon Sep 17 00:00:00 2001
From: liangwang0734 <liang.wang.phys@gmail.com>
Date: Fri, 30 Jun 2023 23:29:49 -0400
Subject: [PATCH 29/30] TestIOAdios2: add beginStep/endStep pairs to support
 bp5

---
 src/kg/testing/io/TestIO.cxx       | 12 ++++++++++++
 src/kg/testing/io/TestIOAdios2.cxx |  8 ++++++++
 2 files changed, 20 insertions(+)

diff --git a/src/kg/testing/io/TestIO.cxx b/src/kg/testing/io/TestIO.cxx
index 785be49c24..2e6fb14e73 100644
--- a/src/kg/testing/io/TestIO.cxx
+++ b/src/kg/testing/io/TestIO.cxx
@@ -9,14 +9,18 @@ TEST(IO, WriteRead)
 
   {
     auto writer = io.open("test.bp", kg::io::Mode::Write);
+    writer.beginStep(kg::io::StepMode::Append);
     writer.putLocal("var_double", 99.);
+    writer.endStep();
     writer.close();
   }
 
   {
     auto reader = io.open("test.bp", kg::io::Mode::Read);
+    reader.beginStep(kg::io::StepMode::Read);
     double dbl;
     reader.getLocal("var_double", dbl);
+    reader.endStep();
     reader.close();
     EXPECT_EQ(dbl, 99.);
   }
@@ -28,14 +32,18 @@ TEST(IO, WriteReadAttr)
 
   {
     auto writer = io.open("test.bp", kg::io::Mode::Write);
+    writer.beginStep(kg::io::StepMode::Append);
     writer.put("attr_double", 99.);
+    writer.endStep();
     writer.close();
   }
 
   {
     auto reader = io.open("test.bp", kg::io::Mode::Read);
+    reader.beginStep(kg::io::StepMode::Read);
     double dbl;
     reader.get("attr_double", dbl);
+    reader.endStep();
     reader.close();
     EXPECT_EQ(dbl, 99.);
   }
@@ -72,15 +80,19 @@ TEST(IO, WriteReadCustom)
 
   {
     auto writer = io.open("test.bp", kg::io::Mode::Write);
+    writer.beginStep(kg::io::StepMode::Append);
     auto c = Custom{3, 99.};
     writer.put("var_custom", c);
+    writer.endStep();
     writer.close();
   }
 
   {
     auto reader = io.open("test.bp", kg::io::Mode::Read);
+    reader.beginStep(kg::io::StepMode::Read);
     auto c = Custom{};
     reader.get("var_custom", c);
+    reader.endStep();
     reader.close();
     EXPECT_EQ(c.i, 3);
     EXPECT_EQ(c.d, 99.);
diff --git a/src/kg/testing/io/TestIOAdios2.cxx b/src/kg/testing/io/TestIOAdios2.cxx
index 65da036c1d..927e06e7d2 100644
--- a/src/kg/testing/io/TestIOAdios2.cxx
+++ b/src/kg/testing/io/TestIOAdios2.cxx
@@ -34,6 +34,7 @@ TEST(IOAdios2, FilePutGetVariable)
   auto io = kg::io::IOAdios2{};
   {
     auto file = io.openFile("test3.bp", kg::io::Mode::Write);
+    file.beginStep(kg::io::StepMode::Append);
     auto dbl = std::vector<double>{1., 2., 3., 4., 5.};
     file.putVariable("dbl", dbl.data(), kg::io::Mode::NonBlocking, {5},
                      {{0}, {5}}, {});
@@ -42,9 +43,11 @@ TEST(IOAdios2, FilePutGetVariable)
     file.putVariable("dbl2", dbl.data() + 3, kg::io::Mode::NonBlocking, {5},
                      {{3}, {2}}, {});
     file.performPuts();
+    file.endStep();
   }
   {
     auto file = io.openFile("test3.bp", kg::io::Mode::Read);
+    file.beginStep(kg::io::StepMode::Read);
 
     auto shape = file.shapeVariable("dbl");
     EXPECT_EQ(shape, kg::io::Dims{5});
@@ -58,6 +61,7 @@ TEST(IOAdios2, FilePutGetVariable)
     file.getVariable("dbl2", dbl2.data(), kg::io::Mode::NonBlocking, {{0}, {5}},
                      {});
     file.performGets();
+    file.endStep();
 
     EXPECT_EQ(dbl, (std::vector<double>{1., 2., 3., 4., 5.}));
     EXPECT_EQ(dbl2, (std::vector<double>{1., 2., 0., 4., 5.}));
@@ -69,15 +73,19 @@ TEST(IOAdios2, FilePutGetAttribute)
   auto io = kg::io::IOAdios2{};
   {
     auto file = io.openFile("test4.bp", kg::io::Mode::Write);
+    file.beginStep(kg::io::StepMode::Append);
     auto dbl = std::vector<double>{1., 2., 3., 4., 5.};
     file.putAttribute("attr_dbl", dbl.data(), dbl.size());
+    file.endStep();
   }
   {
     auto file = io.openFile("test4.bp", kg::io::Mode::Read);
+    file.beginStep(kg::io::StepMode::Read);
     auto size = file.sizeAttribute("attr_dbl");
     auto dbl = std::vector<double>(size);
     file.getAttribute("attr_dbl", dbl.data());
     EXPECT_EQ(dbl, (std::vector<double>{1., 2., 3., 4., 5.}));
+    file.endStep();
   }
 }
 

From eaa19995be257442fc3f7717a1d3a11b71d3d0f2 Mon Sep 17 00:00:00 2001
From: liangwang0734 <liang.wang.phys@gmail.com>
Date: Sun, 2 Jul 2023 00:09:13 -0400
Subject: [PATCH 30/30] checkpoint: check if adios2 is available

---
 src/include/checkpoint.hxx | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/include/checkpoint.hxx b/src/include/checkpoint.hxx
index cb44ecaeb5..af05a60c0c 100644
--- a/src/include/checkpoint.hxx
+++ b/src/include/checkpoint.hxx
@@ -122,6 +122,8 @@ public:
 private:
   int interval_; // write checkpoint every so many steps
   bool first_time_ = true;
+#ifdef PSC_HAVE_ADIOS2
   kg::io::IOAdios2 io;
   kg::io::Engine writer;
+#endif
 };