scritps

miniapps/plasma/pic/electrostatic-pic.cpp

@@ -93,6 +93,21 @@ struct PICContext
    bool reproduce = true;      ///< Enable reproducible results.
 } ctx;
 
+/** Accumulating timers for the main phases of the PIC loop:
+    (1) Calculate forces  - interpolate E field to particle positions
+    (2) Push particles    - leap-frog momentum/position update
+    (3) Redistribute      - migrate particles across MPI ranks
+    (4) Update fields     - charge deposit, Poisson solve, E = -∇φ
+    (5) FindPoints        - locate particles in the mesh (all call sites) */
+struct PICTimers
+{
+   mfem::StopWatch forces;        ///< E-field interpolation to particles.
+   mfem::StopWatch push;          ///< Particle push (leap-frog update only).
+   mfem::StopWatch redistribute;  ///< Particle redistribution across ranks.
+   mfem::StopWatch fields;        ///< Field update (deposit + Poisson + grad).
+   mfem::StopWatch findpts;       ///< FindPointsGSLIB particle location.
+} timers;
+
 /** This class implements explicit time integration for charged particles
     in an electric field using ParticleSet. */
 class ParticleMover
@@ -116,13 +131,16 @@ protected:
    /// ParticleSet of charged particles
    std::unique_ptr<ParticleSet> charged_particles;
 
+   /// Particle pusher on cuda or not
+   bool use_device = false;
+
    /// Temporary vectors for particle computation
    mutable Vector pm_, pp_;
 
 public:
    ParticleMover(MPI_Comm comm, ParGridFunction* E_gf_,
                  FindPointsGSLIB& E_finder_, int num_particles,
-                 Ordering::Type pdata_ordering);
+                 Ordering::Type pdata_ordering, bool use_device_);
 
    /// Initialize charged particles with given parameters
    void InitializeChargedParticles(const real_t& k, const real_t& alpha,
@@ -132,9 +150,12 @@ public:
    /// Find Particles in mesh corresponding to E and field
    void FindParticles();
 
-   /// Advance particles one time step using Boris algorithm
+   /// Advance particles one time step using Boris algorithm on cpu.
    void Step(real_t& t, real_t dt, real_t L, bool first_step = false);
 
+   /// Advance particles one time step using MFEM on GPU.
+   void StepDevice(real_t& t, real_t dt, real_t L, bool first_step = false);
+
    /// Redistribute particles across processors
    void Redistribute();
 
@@ -157,8 +178,13 @@ private:
    real_t neutralizing_const;
    ParLinearForm* precomputed_neutralizing_lf = nullptr;
    bool precompute_neutralizing_const = false;
-   // Diffusion matrix
-   HypreParMatrix* diffusion_matrix;
+
+   HypreParMatrix* diffusion_matrix = nullptr;
+   ParBilinearForm* diffusion_form = nullptr;
+
+   bool use_full_assembly = false;
+   bool use_partial_assembly = false;
+
    // Gradient operator for computing E = -∇φ
    ParDiscreteLinearOperator* grad_interpolator;
    FindPointsGSLIB& E_finder;
@@ -178,7 +204,9 @@ protected:
 public:
    FieldSolver(ParFiniteElementSpace* phi_fes, ParFiniteElementSpace* E_fes,
                FindPointsGSLIB& E_finder_,
-               bool precompute_neutralizing_const_ = false);
+               bool precompute_neutralizing_const_ = false,
+               bool use_full_assembly_ = false,
+               bool use_partial_assembly_ = false);
 
    ~FieldSolver();
 
@@ -207,7 +235,20 @@ int main(int argc, char* argv[])
 
    if (Mpi::Root()) { display_banner(cout); }
 
+   const char *device_config = "cpu";
+   bool fa = false;
+   bool pa = false;
+
    OptionsParser args(argc, argv);
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, e.g. cpu, cuda, ceed-cuda.");
+   args.AddOption(&fa, "-fa", "--full-assembly", "-no-fa",
+                  "--no-full-assembly",
+                  "Use MFEM full assembly path with the selected device backend. "
+                  "This keeps the existing Hypre matrix solver.");
+   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly",
+                  "Use MFEM partial assembly path with a matrix-free CG solve.");
    args.AddOption(&ctx.dim, "-dim", "--dimension",
                   "Spatial dimension (2 or 3)");
    args.AddOption(&ctx.order, "-O", "--order",
@@ -248,6 +289,13 @@ int main(int argc, char* argv[])
    }
    if (Mpi::Root()) { args.PrintOptions(cout); }
 
+   MFEM_VERIFY(!(pa && fa),
+               "Cannot use both partial and full assembly.");
+
+   Device device(device_config);
+   if (Mpi::Root()) { device.Print(); }
+   const bool use_device = (std::string(device_config) != "cpu") && Device::IsEnabled();
+
    // Assert that dimension is 2 or 3
    MFEM_VERIFY(ctx.dim == 2 || ctx.dim == 3,
                "Dimension must be 2 or 3, got " << ctx.dim);
@@ -299,9 +347,11 @@ int main(int argc, char* argv[])
    ParGridFunction E_gf(&E_fespace);
    phi_gf = 0.0;  // Initialize phi_gf to zero
    E_gf = 0.0;    // Initialize E_gf to zero
+   phi_gf.UseDevice(use_device);
+   E_gf.UseDevice(use_device);
 
    // 6. Construct the field solver
-   FieldSolver field_solver(&phi_fespace, &E_fespace, E_finder, true);
+   FieldSolver field_solver(&phi_fespace, &E_fespace, E_finder, true, fa, pa);
 
    // 7. Initialize ParticleMover
    Ordering::Type ordering_type =
@@ -309,7 +359,7 @@ int main(int argc, char* argv[])
    int num_particles =
       ctx.npt / num_ranks + (rank < (ctx.npt % num_ranks) ? 1 : 0);
    ParticleMover particle_mover(MPI_COMM_WORLD, &E_gf, E_finder, num_particles,
-                                ordering_type);
+                                ordering_type, use_device);
    particle_mover.InitializeChargedParticles(ctx.k, ctx.alpha, ctx.m, ctx.q,
                                              ctx.L, ctx.reproduce);
 
@@ -326,14 +376,19 @@ int main(int argc, char* argv[])
           (step % ctx.redist_interval == 0 || step == 1) &&
           particle_mover.GetParticles().GetGlobalNParticles() > 0)
       {
-         // Redistribute
+         // (3) Redistribute particles across MPI ranks
+         //     (timed inside Redistribute(); the FindPoints call it makes
+         //      is accumulated into timers.findpts)
          particle_mover.Redistribute();
 
+         // (4) Update fields: deposit charge, solve Poisson, E = -∇φ
+         timers.fields.Start();
          // Update phi_gf from particles
          field_solver.UpdatePhiGridFunction(particle_mover.GetParticles(),
                                             phi_gf);
          // Update E_gf from phi_gf
          field_solver.UpdateEGridFunction(phi_gf, E_gf);
+         timers.fields.Stop();
 
          // Visualize fields if requested
          if (ctx.visualization)
@@ -347,7 +402,16 @@ int main(int argc, char* argv[])
       }
 
       // Step the ParticleMover
-      particle_mover.Step(t, dt, ctx.L, step == 1);
+      // (1) Calculate forces and (2) push particles are timed inside
+      //     Step() / StepDevice() via the global 'timers' object.
+      if (use_device)
+      {
+         particle_mover.StepDevice(t, dt, ctx.L, step == 1);
+      }
+      else
+      {
+         particle_mover.Step(t, dt, ctx.L, step == 1);
+      }
       if (Mpi::Root())
       {
          mfem::out << "Step: " << step << " | Time: " << t;
@@ -395,12 +459,58 @@ int main(int argc, char* argv[])
          }
       }
    }
+   sw.Stop();
+
+   // 9. Print timing summary
+   // Reduce with MPI_MAX so the reported times reflect the slowest rank
+   // (the one that determines wall-clock time for each phase).
+   real_t t_local[6] = {timers.forces.RealTime(),
+                        timers.push.RealTime(),
+                        timers.redistribute.RealTime(),
+                        timers.fields.RealTime(),
+                        timers.findpts.RealTime(),
+                        sw.RealTime()
+                       };
+   real_t t_max[6];
+   MPI_Reduce(t_local, t_max, 6, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD);
+
+   if (Mpi::Root())
+   {
+      const real_t total = t_max[5];
+      const real_t accounted =
+         t_max[0] + t_max[1] + t_max[2] + t_max[3] + t_max[4];
+      const real_t other = total - accounted;
+
+      auto print_row = [&](const char* name, real_t time)
+      {
+         cout << "  " << left << setw(28) << name << right << setw(12)
+              << fixed << setprecision(4) << time << " s"
+              << setw(9) << setprecision(1)
+              << (total > 0.0 ? 100.0 * time / total : 0.0) << " %" << endl;
+      };
+
+      cout << "\n=========== PIC Timing Summary (max over ranks) ==========="
+           << endl;
+      print_row("1. Calculate forces (E@p)", t_max[0]);
+      print_row("2. Push particles", t_max[1]);
+      print_row("3. Redistribute particles", t_max[2]);
+      print_row("4. Update fields (Poisson)", t_max[3]);
+      print_row("5. FindPoints (locate p)", t_max[4]);
+      print_row("Other (I/O, vis, energies)", other);
+      cout << "  " << string(51, '-') << endl;
+      print_row("Total time", total);
+      cout << "  " << left << setw(28) << "Avg time per step" << right
+           << setw(12) << fixed << setprecision(4)
+           << (ctx.nt > 0 ? total / ctx.nt : 0.0) << " s" << endl;
+      cout << "============================================================"
+           << endl;
+   }
 }
 
 ParticleMover::ParticleMover(MPI_Comm comm, ParGridFunction* E_gf_,
                              FindPointsGSLIB& E_finder_, int num_particles,
-                             Ordering::Type pdata_ordering)
-   : E_gf(E_gf_), E_finder(E_finder_)
+                             Ordering::Type pdata_ordering, bool use_device_)
+   : E_gf(E_gf_), E_finder(E_finder_), use_device(use_device_)
 {
    MFEM_ASSERT(E_gf, "Must pass an E field to ParticleMover.");
 
@@ -413,7 +523,8 @@ ParticleMover::ParticleMover(MPI_Comm comm, ParGridFunction* E_gf_,
    // 2 vectors of size space dim for momentum and e field
    Array<int> field_vdims({1, 1, dim, dim});
    charged_particles = std::make_unique<ParticleSet>(
-                          comm, num_particles, dim, field_vdims, 1, pdata_ordering);
+                          comm, num_particles, dim, field_vdims, 1,
+                          pdata_ordering, use_device);
 }
 
 void ParticleMover::InitializeChargedParticles(const real_t& k,
@@ -436,6 +547,11 @@ void ParticleMover::InitializeChargedParticles(const real_t& k,
    ParticleVector& M = charged_particles->Field(ParticleMover::MASS);
    ParticleVector& Q = charged_particles->Field(ParticleMover::CHARGE);
 
+   X.HostWrite();
+   P.HostWrite();
+   M.HostWrite();
+   Q.HostWrite();
+
    for (int i = 0; i < charged_particles->GetNParticles(); i++)
    {
       // Initialize momentum
@@ -461,19 +577,38 @@ void ParticleMover::InitializeChargedParticles(const real_t& k,
       M(i) = m;
       Q(i) = q;
    }
+
+   X.Read();
+   P.Read();
+   M.Read();
+   Q.Read();
+
    FindParticles();
 }
 
 void ParticleMover::FindParticles()
 {
+   timers.findpts.Start();
    E_finder.FindPoints(charged_particles->Coords());
+#if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
+   if (Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
+   {
+      MFEM_DEVICE_SYNC;
+   }
+#endif
+   timers.findpts.Stop();
 }
 
 void ParticleMover::Step(real_t& t, real_t dt, real_t L, bool first_step)
 {
-   // Update E field at particles
+   // ---- (1) Calculate forces: interpolate E field at particles ----
+   timers.forces.Start();
    ParticleVector& E = charged_particles->Field(EFIELD);
    E_finder.Interpolate(*E_gf, E, E.GetOrdering());
+   timers.forces.Stop();
+
+   // ---- (2) Push particles: leap-frog update + relocate in mesh ----
+   timers.push.Start();
 
    // Extract particle data
    ParticleVector& X = charged_particles->Coords();
@@ -505,16 +640,96 @@ void ParticleMover::Step(real_t& t, real_t dt, real_t L, bool first_step)
       }
    }
 
-   FindParticles();
+   timers.push.Stop();
+
+   FindParticles();  // timed separately by timers.findpts
 
    // Update time
    t += dt;
 }
 
+
+void ParticleMover::StepDevice(real_t& t, real_t dt, real_t L, bool first_step)
+{
+   // ---- (1) Calculate forces: interpolate E field at particles ----
+   timers.forces.Start();
+   ParticleVector& E = charged_particles->Field(EFIELD);
+   E_finder.Interpolate(*E_gf, E, E.GetOrdering());
+#if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
+   // Only sync if a real GPU backend is active at runtime. The "debug"
+   // backend enables Device::IsEnabled() but runs on the host, and calling
+   // cudaDeviceSynchronize() there triggers a CUDA error.
+   if (Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
+   {
+      MFEM_DEVICE_SYNC;  // ensure accurate timing of async device kernels
+   }
+#endif
+   timers.forces.Stop();
+
+   // ---- (2) Push particles: leap-frog update + relocate in mesh ----
+   timers.push.Start();
+
+   ParticleVector& X = charged_particles->Coords();
+   ParticleVector& P = charged_particles->Field(MOM);
+   ParticleVector& M = charged_particles->Field(MASS);
+   ParticleVector& Q = charged_particles->Field(CHARGE);
+
+   const int npt = charged_particles->GetNParticles();
+   const int dim = X.GetVDim();
+   const real_t accel_dt = first_step ? dt / 2.0 : dt;
+
+   const bool byVDIM_X = (X.GetOrdering() == Ordering::byVDIM);
+   const bool byVDIM_P = (P.GetOrdering() == Ordering::byVDIM);
+   const bool byVDIM_E = (E.GetOrdering() == Ordering::byVDIM);
+
+   auto d_x = X.ReadWrite();
+   auto d_p = P.ReadWrite();
+   auto d_m = M.Read();
+   auto d_q = Q.Read();
+   auto d_e = E.Read();
+
+   MFEM_FORALL(particle, npt,
+   {
+      const real_t m = d_m[particle];
+      const real_t q = d_q[particle];
+
+      for (int d = 0; d < dim; ++d)
+      {
+         const int x_idx = byVDIM_X ? particle * dim + d : particle + d * npt;
+         const int p_idx = byVDIM_P ? particle * dim + d : particle + d * npt;
+         const int e_idx = byVDIM_E ? particle * dim + d : particle + d * npt;
+
+         real_t p_new = d_p[p_idx] + accel_dt * q * d_e[e_idx];
+         real_t x_new = d_x[x_idx] + dt / m * p_new;
+
+         while (x_new >= L) { x_new -= L; }
+         while (x_new < 0.0) { x_new += L; }
+
+         d_p[p_idx] = p_new;
+         d_x[x_idx] = x_new;
+      }
+   });
+
+#if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
+   if (Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
+   {
+      MFEM_DEVICE_SYNC;  // ensure accurate timing of async device kernels
+   }
+#endif
+   timers.push.Stop();
+
+   FindParticles();  // timed separately by timers.findpts
+
+   t += dt;
+}
+
 void ParticleMover::Redistribute()
 {
+   timers.redistribute.Start();
    charged_particles->Redistribute(E_finder.GetProc());
-   FindParticles();
+   timers.redistribute.Stop();
+
+   FindParticles();  // timed separately by timers.findpts
 }
 
 real_t ParticleMover::ComputeKineticEnergy(real_t dt) const
@@ -528,15 +743,49 @@ real_t ParticleMover::ComputeKineticEnergy(real_t dt) const
    // update from Step() is used directly.
 
    real_t kinetic_energy = 0.0;
-   for (int p = 0; p < charged_particles->GetNParticles(); ++p)
+   const int npt = charged_particles->GetNParticles();
+   const int dim = P.GetVDim();
+
+   if (this->use_device)
    {
-      real_t p_square_p = 0.0;
-      for (int d = 0; d < P.GetVDim(); ++d)
+      Vector ke(npt);
+      auto d_ke = ke.Write();
+      auto d_p = P.Read();
+      auto d_m = M.Read();
+      auto d_q = Q.Read();
+      auto d_e = E.Read();
+
+      const bool byVDIM_P = (P.GetOrdering() == Ordering::byVDIM);
+      const bool byVDIM_E = (E.GetOrdering() == Ordering::byVDIM);
+
+      MFEM_FORALL(particle, npt,
       {
-         const real_t P_m = P(p, d) + dt * Q(p) * E(p, d);
-         p_square_p += P_m * P_m;
+         real_t p_square_p = 0.0;
+         const real_t q = d_q[particle];
+         for (int d = 0; d < dim; ++d)
+         {
+            const int p_idx = byVDIM_P ? particle * dim + d : particle + d * npt;
+            const int e_idx = byVDIM_E ? particle * dim + d : particle + d * npt;
+            const real_t p_m = d_p[p_idx] + dt * q * d_e[e_idx];
+            p_square_p += p_m * p_m;
+         }
+         d_ke[particle] = 0.5 * p_square_p / d_m[particle];
+      });
+
+      kinetic_energy = ke.Sum();
+   }
+   else
+   {
+      for (int p = 0; p < npt; ++p)
+      {
+         real_t p_square_p = 0.0;
+         for (int d = 0; d < dim; ++d)
+         {
+            const real_t P_m = P(p, d) + dt * Q(p) * E(p, d);
+            p_square_p += P_m * P_m;
+         }
+         kinetic_energy += 0.5 * p_square_p / M(p);
       }
-      kinetic_energy += 0.5 * p_square_p / M(p);
    }
 
    real_t global_kinetic_energy = 0.0;
@@ -548,10 +797,14 @@ real_t ParticleMover::ComputeKineticEnergy(real_t dt) const
 FieldSolver::FieldSolver(ParFiniteElementSpace* phi_fes,
                          ParFiniteElementSpace* E_fes,
                          FindPointsGSLIB& E_finder_,
-                         bool precompute_neutralizing_const_)
+                         bool precompute_neutralizing_const_,
+                         bool use_full_assembly_,
+                         bool use_partial_assembly_)
    : precompute_neutralizing_const(precompute_neutralizing_const_),
      E_finder(E_finder_),
-     b(phi_fes)
+     b(phi_fes),
+     use_full_assembly(use_full_assembly_),
+     use_partial_assembly(use_partial_assembly_)
 {
    // compute domain volume
    ParMesh* pmesh = phi_fes->GetParMesh();
@@ -564,16 +817,27 @@ FieldSolver::FieldSolver(ParFiniteElementSpace* phi_fes,
                  phi_fes->GetParMesh()->GetComm());
 
    {
-      // Par bilinear form for the gradgrad matrix
-      ParBilinearForm dm(phi_fes);
+      diffusion_form = new ParBilinearForm(phi_fes);
+
+      if (use_partial_assembly)
+      {
+         diffusion_form->SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      }
+      else if (use_full_assembly)
+      {
+         diffusion_form->SetAssemblyLevel(AssemblyLevel::FULL);
+         diffusion_form->EnableSparseMatrixSorting(Device::IsEnabled());
+      }
+
       ConstantCoefficient epsilon(EPSILON);  // ε_0
-      dm.AddDomainIntegrator(
-         new DiffusionIntegrator(epsilon));  // ∫ ∇φ_i · ∇φ_j
+      diffusion_form->AddDomainIntegrator(new DiffusionIntegrator(epsilon));  // ∫ ∇φ_i · ∇φ_j
+      diffusion_form->Assemble();
 
-      dm.Assemble();
-      dm.Finalize();
-
-      diffusion_matrix = dm.ParallelAssemble();  // global gradgrad matrix
+      if (!use_partial_assembly)
+      {
+         diffusion_form->Finalize();
+         diffusion_matrix = diffusion_form->ParallelAssemble();
+      }
    }
 
    {
@@ -587,6 +851,7 @@ FieldSolver::FieldSolver(ParFiniteElementSpace* phi_fes,
 FieldSolver::~FieldSolver()
 {
    delete diffusion_matrix;
+   delete diffusion_form;
    delete precomputed_neutralizing_lf;
    delete grad_interpolator;
 }
@@ -594,9 +859,10 @@ FieldSolver::~FieldSolver()
 const ParLinearForm& FieldSolver::ComputeNeutralizingRHS(
    ParFiniteElementSpace* pfes, const ParticleVector& Q, MPI_Comm comm)
 {
-   int npt = Q.Size();
-   // Get E_finder references
-   const Array<unsigned int>& code = E_finder.GetCode();
+    Q.HostRead();
+    int npt = Q.Size();
+    const Array<unsigned int>& code = E_finder.GetCode();
+    code.HostRead();
 
    if (!precompute_neutralizing_const || precomputed_neutralizing_lf == nullptr)
    {
@@ -639,6 +905,7 @@ const ParLinearForm& FieldSolver::ComputeNeutralizingRHS(
 void FieldSolver::DepositCharge(ParFiniteElementSpace* pfes,
                                 const ParticleVector& Q)
 {
+    Q.HostRead();
    int npt = Q.Size();
    ParMesh* pmesh = pfes->GetParMesh();
    int dim = pmesh->SpaceDimension();
@@ -647,10 +914,15 @@ void FieldSolver::DepositCharge(ParFiniteElementSpace* pfes,
 
    // Get E_finder references
    // 0: inside, 1: boundary, 2: not found
-   const Array<unsigned int>& code = E_finder.GetCode();
-   const Array<unsigned int>& proc = E_finder.GetProc();  // owning MPI rank
-   const Array<unsigned int>& elem = E_finder.GetElem();  // local element id
-   const Vector& rref = E_finder.GetReferencePosition();  // (r,s,t) byVDIM
+    const Array<unsigned int>& code = E_finder.GetCode();
+    const Array<unsigned int>& proc = E_finder.GetProc();
+    const Array<unsigned int>& elem = E_finder.GetElem();
+    const Vector& rref = E_finder.GetReferencePosition();
+
+    code.HostRead();
+    proc.HostRead();
+    elem.HostRead();
+    rref.HostRead();
 
    Array<int> dofs;
 
@@ -716,25 +988,52 @@ void FieldSolver::UpdatePhiGridFunction(ParticleSet& particles,
    // 3) Solve A * phi = B with zero-mean enforcement via OrthoSolver
    // ------------------------------------------------------------------
    phi_gf = 0.0;
-   HypreParVector Phi_true(pfes);
-   Phi_true = 0.0;
 
-   HyprePCG solver(diffusion_matrix->GetComm());
-   solver.SetOperator(*diffusion_matrix);
-   solver.SetTol(1e-12);
-   solver.SetMaxIter(200);
-   solver.SetPrintLevel(0);
+   if (use_partial_assembly)
+   {
+      Array<int> ess_tdof_list;
 
-   HypreBoomerAMG prec(*diffusion_matrix);
-   prec.SetPrintLevel(0);
-   solver.SetPreconditioner(prec);
+      OperatorPtr A;
+      Vector X, B_pa;
+      diffusion_form->FormLinearSystem(ess_tdof_list, phi_gf, b, A, X, B_pa);
 
-   OrthoSolver ortho(comm);
-   ortho.SetSolver(solver);
-   ortho.Mult(B, Phi_true);
+      CGSolver solver(comm);
+      solver.SetRelTol(1e-10);
+      solver.SetAbsTol(0.0);
+      solver.SetMaxIter(10000);
+      solver.SetPrintLevel(0);
+      solver.SetOperator(*A);
 
-   // Map true-dof solution back to the ParGridFunction
-   phi_gf.Distribute(Phi_true);
+      OperatorJacobiSmoother prec(*diffusion_form, ess_tdof_list);
+      solver.SetPreconditioner(prec);
+
+      OrthoSolver ortho(comm);
+      ortho.SetSolver(solver);
+      ortho.Mult(B_pa, X);
+
+      diffusion_form->RecoverFEMSolution(X, b, phi_gf);
+   }
+   else
+   {
+      HypreParVector Phi_true(pfes);
+      Phi_true = 0.0;
+
+      HyprePCG solver(diffusion_matrix->GetComm());
+      solver.SetOperator(*diffusion_matrix);
+      solver.SetTol(1e-12);
+      solver.SetMaxIter(200);
+      solver.SetPrintLevel(0);
+
+      HypreBoomerAMG prec(*diffusion_matrix);
+      prec.SetPrintLevel(0);
+      solver.SetPreconditioner(prec);
+
+      OrthoSolver ortho(comm);
+      ortho.SetSolver(solver);
+      ortho.Mult(B, Phi_true);
+
+      phi_gf.Distribute(Phi_true);
+   }
 }
 
 void FieldSolver::UpdateEGridFunction(ParGridFunction& phi_gf,
@@ -785,4 +1084,4 @@ void display_banner(ostream& os)
          )"
       << endl
       << flush;
-}
+}

miniapps/plasma/pic/logs/cuda_2d.log

@@ -0,0 +1,276 @@
+#############################
+ranks=1  mesh=32^2  npt=409600  q=m=0.0011816406253684498
+#############################
+
+      ██████╗░██╗░█████╗░
+      ██╔══██╗██║██╔══██╗
+      ██████╔╝██║██║░░╚═╝
+      ██╔═══╝░██║██║░░██╗
+      ██║░░░░░██║╚█████╔╝
+      ╚═╝░░░░░╚═╝░╚════╝░
+         
+Options used:
+   --device cpu
+   --full-assembly
+   --no-partial-assembly
+   --dimension 2
+   --order 1
+   --num-x 32
+   --num-y 32
+   --num-z 100
+   --charge 0.00118164
+   --mass 0.00118164
+   --time-step 0.1
+   --num-timesteps 10
+   --num-particles 409600
+   --k 0.285599
+   --alpha 0.05
+   --ordering 1
+   --redist-interval 1
+   --output-csv-interval -1
+   --no-visualization
+   --send-port 19916
+   --reproduce
+Device configuration: cpu
+Memory configuration: host-std
+Total charge: 484, Domain volume: 484, Neutralizing constant: -1
+Further updates will use this precomputed neutralizing constant.
+Step: 1 | Time: 0.1 | Time per step: 0.410409
+Kinetic energy: 484.208	Field energy: 7.2893	Total energy: 491.497
+Step: 2 | Time: 0.2 | Time per step: 0.418639
+Kinetic energy: 484.311	Field energy: 7.18476	Total energy: 491.496
+Step: 3 | Time: 0.3 | Time per step: 0.417743
+Kinetic energy: 484.568	Field energy: 6.92838	Total energy: 491.496
+Step: 4 | Time: 0.4 | Time per step: 0.417233
+Kinetic energy: 484.964	Field energy: 6.53039	Total energy: 491.495
+Step: 5 | Time: 0.5 | Time per step: 0.41692
+Kinetic energy: 485.482	Field energy: 6.01109	Total energy: 491.493
+Step: 6 | Time: 0.6 | Time per step: 0.41677
+Kinetic energy: 486.098	Field energy: 5.39433	Total energy: 491.493
+Step: 7 | Time: 0.7 | Time per step: 0.41671
+Kinetic energy: 486.785	Field energy: 4.70615	Total energy: 491.491
+Step: 8 | Time: 0.8 | Time per step: 0.416604
+Kinetic energy: 487.512	Field energy: 3.97614	Total energy: 491.488
+Step: 9 | Time: 0.9 | Time per step: 0.41669
+Kinetic energy: 488.249	Field energy: 3.2369	Total energy: 491.486
+Step: 10 | Time: 1 | Time per step: 0.416736
+Kinetic energy: 488.965	Field energy: 2.51996	Total energy: 491.485
+
+=========== PIC Timing Summary (max over ranks) ===========
+  1. Calculate forces (E@p)         1.3285 s     31.8 %
+  2. Push particles                 0.1487 s      3.6 %
+  3. Redistribute particles         0.0028 s      0.1 %
+  4. Update fields (Poisson)        0.3524 s      8.4 %
+  5. FindPoints (locate p)          2.3713 s     56.8 %
+  Other (I/O, vis, energies)       -0.0286 s     -0.7 %
+  ---------------------------------------------------
+  Total time                        4.1752 s    100.0 %
+  Avg time per step                 0.4175 s
+============================================================
+#############################
+ranks=4  mesh=64^2  npt=1638400  q=m=0.00029541015634211245
+#############################
+
+      ██████╗░██╗░█████╗░
+      ██╔══██╗██║██╔══██╗
+      ██████╔╝██║██║░░╚═╝
+      ██╔═══╝░██║██║░░██╗
+      ██║░░░░░██║╚█████╔╝
+      ╚═╝░░░░░╚═╝░╚════╝░
+         
+Options used:
+   --device cpu
+   --full-assembly
+   --no-partial-assembly
+   --dimension 2
+   --order 1
+   --num-x 64
+   --num-y 64
+   --num-z 100
+   --charge 0.00029541
+   --mass 0.00029541
+   --time-step 0.1
+   --num-timesteps 10
+   --num-particles 1638400
+   --k 0.285599
+   --alpha 0.05
+   --ordering 1
+   --redist-interval 1
+   --output-csv-interval -1
+   --no-visualization
+   --send-port 19916
+   --reproduce
+Device configuration: cpu
+Memory configuration: host-std
+Total charge: 484, Domain volume: 484, Neutralizing constant: -1
+Further updates will use this precomputed neutralizing constant.
+Step: 1 | Time: 0.1 | Time per step: 0.677767
+Kinetic energy: 483.977	Field energy: 7.82945	Total energy: 491.807
+Step: 2 | Time: 0.2 | Time per step: 0.653161
+Kinetic energy: 484.061	Field energy: 7.74561	Total energy: 491.807
+Step: 3 | Time: 0.3 | Time per step: 0.643092
+Kinetic energy: 484.309	Field energy: 7.49724	Total energy: 491.806
+Step: 4 | Time: 0.4 | Time per step: 0.636938
+Kinetic energy: 484.709	Field energy: 7.09583	Total energy: 491.805
+Step: 5 | Time: 0.5 | Time per step: 0.634511
+Kinetic energy: 485.244	Field energy: 6.55892	Total energy: 491.803
+Step: 6 | Time: 0.6 | Time per step: 0.631674
+Kinetic energy: 485.891	Field energy: 5.91033	Total energy: 491.802
+Step: 7 | Time: 0.7 | Time per step: 0.628884
+Kinetic energy: 486.621	Field energy: 5.1789	Total energy: 491.8
+Step: 8 | Time: 0.8 | Time per step: 0.626798
+Kinetic energy: 487.401	Field energy: 4.39656	Total energy: 491.797
+Step: 9 | Time: 0.9 | Time per step: 0.625584
+Kinetic energy: 488.198	Field energy: 3.59751	Total energy: 491.795
+Step: 10 | Time: 1 | Time per step: 0.624849
+Kinetic energy: 488.978	Field energy: 2.81567	Total energy: 491.794
+
+=========== PIC Timing Summary (max over ranks) ===========
+  1. Calculate forces (E@p)         2.0789 s     33.2 %
+  2. Push particles                 0.2043 s      3.3 %
+  3. Redistribute particles         0.1613 s      2.6 %
+  4. Update fields (Poisson)        0.6066 s      9.7 %
+  5. FindPoints (locate p)          3.5290 s     56.4 %
+  Other (I/O, vis, energies)       -0.3195 s     -5.1 %
+  ---------------------------------------------------
+  Total time                        6.2607 s    100.0 %
+  Avg time per step                 0.6261 s
+============================================================
+#############################
+ranks=16  mesh=128^2  npt=6553600  q=m=7.3852539085528113e-05
+#############################
+
+      ██████╗░██╗░█████╗░
+      ██╔══██╗██║██╔══██╗
+      ██████╔╝██║██║░░╚═╝
+      ██╔═══╝░██║██║░░██╗
+      ██║░░░░░██║╚█████╔╝
+      ╚═╝░░░░░╚═╝░╚════╝░
+         
+Options used:
+   --device cpu
+   --full-assembly
+   --no-partial-assembly
+   --dimension 2
+   --order 1
+   --num-x 128
+   --num-y 128
+   --num-z 100
+   --charge 7.38525e-05
+   --mass 7.38525e-05
+   --time-step 0.1
+   --num-timesteps 10
+   --num-particles 6553600
+   --k 0.285599
+   --alpha 0.05
+   --ordering 1
+   --redist-interval 1
+   --output-csv-interval -1
+   --no-visualization
+   --send-port 19916
+   --reproduce
+Device configuration: cpu
+Memory configuration: host-std
+Total charge: 484, Domain volume: 484, Neutralizing constant: -1
+Further updates will use this precomputed neutralizing constant.
+Step: 1 | Time: 0.1 | Time per step: 1.04666
+Kinetic energy: 484.098	Field energy: 7.45399	Total energy: 491.552
+Step: 2 | Time: 0.2 | Time per step: 0.899693
+Kinetic energy: 484.175	Field energy: 7.37733	Total energy: 491.552
+Step: 3 | Time: 0.3 | Time per step: 0.848083
+Kinetic energy: 484.408	Field energy: 7.1433	Total energy: 491.551
+Step: 4 | Time: 0.4 | Time per step: 0.824958
+Kinetic energy: 484.788	Field energy: 6.76264	Total energy: 491.55
+Step: 5 | Time: 0.5 | Time per step: 0.809805
+Kinetic energy: 485.297	Field energy: 6.25213	Total energy: 491.549
+Step: 6 | Time: 0.6 | Time per step: 0.798508
+Kinetic energy: 485.912	Field energy: 5.63448	Total energy: 491.547
+Step: 7 | Time: 0.7 | Time per step: 0.791095
+Kinetic energy: 486.608	Field energy: 4.93722	Total energy: 491.545
+Step: 8 | Time: 0.8 | Time per step: 0.785146
+Kinetic energy: 487.352	Field energy: 4.19078	Total energy: 491.543
+Step: 9 | Time: 0.9 | Time per step: 0.778998
+Kinetic energy: 488.113	Field energy: 3.42757	Total energy: 491.541
+Step: 10 | Time: 1 | Time per step: 0.773919
+Kinetic energy: 488.859	Field energy: 2.68018	Total energy: 491.539
+
+=========== PIC Timing Summary (max over ranks) ===========
+  1. Calculate forces (E@p)         2.3439 s     30.2 %
+  2. Push particles                 0.2228 s      2.9 %
+  3. Redistribute particles         0.4353 s      5.6 %
+  4. Update fields (Poisson)        0.7318 s      9.4 %
+  5. FindPoints (locate p)          4.6342 s     59.8 %
+  Other (I/O, vis, energies)       -0.6154 s     -7.9 %
+  ---------------------------------------------------
+  Total time                        7.7526 s    100.0 %
+  Avg time per step                 0.7753 s
+============================================================
+#############################
+ranks=64  mesh=256^2  npt=26214400  q=m=1.8463134771382028e-05
+#############################
+
+      ██████╗░██╗░█████╗░
+      ██╔══██╗██║██╔══██╗
+      ██████╔╝██║██║░░╚═╝
+      ██╔═══╝░██║██║░░██╗
+      ██║░░░░░██║╚█████╔╝
+      ╚═╝░░░░░╚═╝░╚════╝░
+         
+Options used:
+   --device cpu
+   --full-assembly
+   --no-partial-assembly
+   --dimension 2
+   --order 1
+   --num-x 256
+   --num-y 256
+   --num-z 100
+   --charge 1.84631e-05
+   --mass 1.84631e-05
+   --time-step 0.1
+   --num-timesteps 10
+   --num-particles 26214400
+   --k 0.285599
+   --alpha 0.05
+   --ordering 1
+   --redist-interval 1
+   --output-csv-interval -1
+   --no-visualization
+   --send-port 19916
+   --reproduce
+Device configuration: cpu
+Memory configuration: host-std
+Total charge: 484, Domain volume: 484, Neutralizing constant: -1
+Further updates will use this precomputed neutralizing constant.
+Step: 1 | Time: 0.1 | Time per step: 1.50239
+Kinetic energy: 484.063	Field energy: 7.39334	Total energy: 491.456
+Step: 2 | Time: 0.2 | Time per step: 1.17161
+Kinetic energy: 484.139	Field energy: 7.31687	Total energy: 491.456
+Step: 3 | Time: 0.3 | Time per step: 1.06237
+Kinetic energy: 484.371	Field energy: 7.08415	Total energy: 491.455
+Step: 4 | Time: 0.4 | Time per step: 1.00664
+Kinetic energy: 484.748	Field energy: 6.70566	Total energy: 491.454
+Step: 5 | Time: 0.5 | Time per step: 0.97318
+Kinetic energy: 485.254	Field energy: 6.19836	Total energy: 491.453
+Step: 6 | Time: 0.6 | Time per step: 0.950163
+Kinetic energy: 485.866	Field energy: 5.58482	Total energy: 491.451
+Step: 7 | Time: 0.7 | Time per step: 0.934015
+Kinetic energy: 486.557	Field energy: 4.89226	Total energy: 491.449
+Step: 8 | Time: 0.8 | Time per step: 0.922009
+Kinetic energy: 487.296	Field energy: 4.15094	Total energy: 491.447
+Step: 9 | Time: 0.9 | Time per step: 0.912949
+Kinetic energy: 488.052	Field energy: 3.39314	Total energy: 491.445
+Step: 10 | Time: 1 | Time per step: 0.905082
+Kinetic energy: 488.791	Field energy: 2.65125	Total energy: 491.443
+
+=========== PIC Timing Summary (max over ranks) ===========
+  1. Calculate forces (E@p)         2.6915 s     29.7 %
+  2. Push particles                 0.2120 s      2.3 %
+  3. Redistribute particles         1.1491 s     12.7 %
+  4. Update fields (Poisson)        0.8290 s      9.1 %
+  5. FindPoints (locate p)          5.5258 s     61.0 %
+  Other (I/O, vis, energies)       -1.3417 s    -14.8 %
+  ---------------------------------------------------
+  Total time                        9.0657 s    100.0 %
+  Avg time per step                 0.9066 s
+============================================================

miniapps/plasma/pic/logs/cuda_3d.log

@@ -0,0 +1,210 @@
+#############################
+ranks=1  mesh=16^3  npt=409600  q=m=0.0048447307312968462
+#############################
+
+      ██████╗░██╗░█████╗░
+      ██╔══██╗██║██╔══██╗
+      ██████╔╝██║██║░░╚═╝
+      ██╔═══╝░██║██║░░██╗
+      ██║░░░░░██║╚█████╔╝
+      ╚═╝░░░░░╚═╝░╚════╝░
+         
+Options used:
+   --device cuda
+   --full-assembly
+   --no-partial-assembly
+   --dimension 3
+   --order 1
+   --num-x 16
+   --num-y 16
+   --num-z 16
+   --charge 0.00484473
+   --mass 0.00484473
+   --time-step 0.02
+   --num-timesteps 10
+   --num-particles 409600
+   --k 0.5
+   --alpha 0.01
+   --ordering 1
+   --redist-interval 1
+   --output-csv-interval -1
+   --no-visualization
+   --send-port 19916
+   --reproduce
+Device configuration: cuda,cpu
+Memory configuration: host-std,cuda
+Use GPU-aware MPI:    no
+Total charge: 1984.4, Domain volume: 1984.4, Neutralizing constant: -1
+Further updates will use this precomputed neutralizing constant.
+Step: 1 | Time: 0.02 | Time per step: 0.35749
+Kinetic energy: 2979.46	Field energy: 1.65511	Total energy: 2981.11
+Step: 2 | Time: 0.04 | Time per step: 0.348279
+Kinetic energy: 2979.46	Field energy: 1.65308	Total energy: 2981.11
+Step: 3 | Time: 0.06 | Time per step: 0.346449
+Kinetic energy: 2979.46	Field energy: 1.65022	Total energy: 2981.11
+Step: 4 | Time: 0.08 | Time per step: 0.343719
+Kinetic energy: 2979.46	Field energy: 1.64676	Total energy: 2981.11
+Step: 5 | Time: 0.1 | Time per step: 0.342003
+Kinetic energy: 2979.47	Field energy: 1.64308	Total energy: 2981.11
+Step: 6 | Time: 0.12 | Time per step: 0.342526
+Kinetic energy: 2979.47	Field energy: 1.63865	Total energy: 2981.11
+Step: 7 | Time: 0.14 | Time per step: 0.341468
+Kinetic energy: 2979.48	Field energy: 1.6331	Total energy: 2981.11
+Step: 8 | Time: 0.16 | Time per step: 0.341208
+Kinetic energy: 2979.49	Field energy: 1.62567	Total energy: 2981.11
+Step: 9 | Time: 0.18 | Time per step: 0.342033
+Kinetic energy: 2979.49	Field energy: 1.61658	Total energy: 2981.11
+Step: 10 | Time: 0.2 | Time per step: 0.341634
+Kinetic energy: 2979.5	Field energy: 1.60714	Total energy: 2981.11
+
+=========== PIC Timing Summary (max over ranks) ===========
+  1. Calculate forces (E@p)         2.4046 s     70.1 %
+  2. Push particles                 0.0022 s      0.1 %
+  3. Redistribute particles         0.0046 s      0.1 %
+  4. Update fields (Poisson)        0.7747 s     22.6 %
+  5. FindPoints (locate p)          0.1194 s      3.5 %
+  Other (I/O, vis, energies)        0.1256 s      3.7 %
+  ---------------------------------------------------
+  Total time                        3.4310 s    100.0 %
+  Avg time per step                 0.3431 s
+============================================================
+#############################
+ranks=8  mesh=32^3  npt=3276800  q=m=0.00060559134141210578
+#############################
+
+      ██████╗░██╗░█████╗░
+      ██╔══██╗██║██╔══██╗
+      ██████╔╝██║██║░░╚═╝
+      ██╔═══╝░██║██║░░██╗
+      ██║░░░░░██║╚█████╔╝
+      ╚═╝░░░░░╚═╝░╚════╝░
+         
+Options used:
+   --device cuda
+   --full-assembly
+   --no-partial-assembly
+   --dimension 3
+   --order 1
+   --num-x 32
+   --num-y 32
+   --num-z 32
+   --charge 0.000605591
+   --mass 0.000605591
+   --time-step 0.02
+   --num-timesteps 10
+   --num-particles 3276800
+   --k 0.5
+   --alpha 0.01
+   --ordering 1
+   --redist-interval 1
+   --output-csv-interval -1
+   --no-visualization
+   --send-port 19916
+   --reproduce
+Device configuration: cuda,cpu
+Memory configuration: host-std,cuda
+Use GPU-aware MPI:    no
+Total charge: 1984.4, Domain volume: 1984.4, Neutralizing constant: -1
+Further updates will use this precomputed neutralizing constant.
+Step: 1 | Time: 0.02 | Time per step: 0.746432
+Kinetic energy: 2978.05	Field energy: 0.799737	Total energy: 2978.85
+Step: 2 | Time: 0.04 | Time per step: 0.690279
+Kinetic energy: 2978.05	Field energy: 0.800322	Total energy: 2978.85
+Step: 3 | Time: 0.06 | Time per step: 0.667386
+Kinetic energy: 2978.05	Field energy: 0.799868	Total energy: 2978.85
+Step: 4 | Time: 0.08 | Time per step: 0.651101
+Kinetic energy: 2978.05	Field energy: 0.798389	Total energy: 2978.85
+Step: 5 | Time: 0.1 | Time per step: 0.64199
+Kinetic energy: 2978.05	Field energy: 0.796093	Total energy: 2978.85
+Step: 6 | Time: 0.12 | Time per step: 0.636094
+Kinetic energy: 2978.06	Field energy: 0.793324	Total energy: 2978.85
+Step: 7 | Time: 0.14 | Time per step: 0.630078
+Kinetic energy: 2978.06	Field energy: 0.789973	Total energy: 2978.85
+Step: 8 | Time: 0.16 | Time per step: 0.625063
+Kinetic energy: 2978.06	Field energy: 0.786052	Total energy: 2978.85
+Step: 9 | Time: 0.18 | Time per step: 0.621014
+Kinetic energy: 2978.07	Field energy: 0.781706	Total energy: 2978.85
+Step: 10 | Time: 0.2 | Time per step: 0.618461
+Kinetic energy: 2978.07	Field energy: 0.776765	Total energy: 2978.85
+
+=========== PIC Timing Summary (max over ranks) ===========
+  1. Calculate forces (E@p)         3.7631 s     60.6 %
+  2. Push particles                 0.0059 s      0.1 %
+  3. Redistribute particles         0.2159 s      3.5 %
+  4. Update fields (Poisson)        1.4331 s     23.1 %
+  5. FindPoints (locate p)          0.9369 s     15.1 %
+  Other (I/O, vis, energies)       -0.1493 s     -2.4 %
+  ---------------------------------------------------
+  Total time                        6.2056 s    100.0 %
+  Avg time per step                 0.6206 s
+============================================================
+#############################
+ranks=64  mesh=64^3  npt=26214400  q=m=7.5698917676513222e-05
+#############################
+
+      ██████╗░██╗░█████╗░
+      ██╔══██╗██║██╔══██╗
+      ██████╔╝██║██║░░╚═╝
+      ██╔═══╝░██║██║░░██╗
+      ██║░░░░░██║╚█████╔╝
+      ╚═╝░░░░░╚═╝░╚════╝░
+         
+Options used:
+   --device cuda
+   --full-assembly
+   --no-partial-assembly
+   --dimension 3
+   --order 1
+   --num-x 64
+   --num-y 64
+   --num-z 64
+   --charge 7.56989e-05
+   --mass 7.56989e-05
+   --time-step 0.02
+   --num-timesteps 10
+   --num-particles 26214400
+   --k 0.5
+   --alpha 0.01
+   --ordering 1
+   --redist-interval 1
+   --output-csv-interval -1
+   --no-visualization
+   --send-port 19916
+   --reproduce
+Device configuration: cuda,cpu
+Memory configuration: host-std,cuda
+Use GPU-aware MPI:    no
+Total charge: 1984.4, Domain volume: 1984.4, Neutralizing constant: -1
+Further updates will use this precomputed neutralizing constant.
+Step: 1 | Time: 0.02 | Time per step: 2.10004
+Kinetic energy: 2977.19	Field energy: 0.636153	Total energy: 2977.83
+Step: 2 | Time: 0.04 | Time per step: 1.74643
+Kinetic energy: 2977.19	Field energy: 0.635865	Total energy: 2977.83
+Step: 3 | Time: 0.06 | Time per step: 1.55546
+Kinetic energy: 2977.19	Field energy: 0.634917	Total energy: 2977.83
+Step: 4 | Time: 0.08 | Time per step: 1.44911
+Kinetic energy: 2977.19	Field energy: 0.63339	Total energy: 2977.83
+Step: 5 | Time: 0.1 | Time per step: 1.37504
+Kinetic energy: 2977.2	Field energy: 0.631344	Total energy: 2977.83
+Step: 6 | Time: 0.12 | Time per step: 1.32539
+Kinetic energy: 2977.2	Field energy: 0.628712	Total energy: 2977.83
+Step: 7 | Time: 0.14 | Time per step: 1.28712
+Kinetic energy: 2977.2	Field energy: 0.625445	Total energy: 2977.83
+Step: 8 | Time: 0.16 | Time per step: 1.25707
+Kinetic energy: 2977.21	Field energy: 0.621599	Total energy: 2977.83
+Step: 9 | Time: 0.18 | Time per step: 1.23524
+Kinetic energy: 2977.21	Field energy: 0.61722	Total energy: 2977.83
+Step: 10 | Time: 0.2 | Time per step: 1.21563
+Kinetic energy: 2977.21	Field energy: 0.612302	Total energy: 2977.83
+
+=========== PIC Timing Summary (max over ranks) ===========
+  1. Calculate forces (E@p)         4.6997 s     38.5 %
+  2. Push particles                 0.0202 s      0.2 %
+  3. Redistribute particles         1.1160 s      9.2 %
+  4. Update fields (Poisson)        3.3711 s     27.6 %
+  5. FindPoints (locate p)          3.8824 s     31.8 %
+  Other (I/O, vis, energies)       -0.8947 s     -7.3 %
+  ---------------------------------------------------
+  Total time                       12.1947 s    100.0 %
+  Avg time per step                 1.2195 s
+============================================================

miniapps/plasma/pic/logs/weak-2d.log

@@ -0,0 +1,276 @@
+#############################
+ranks=1  mesh=32^2  npt=409600  q=m=0.0011816406253684498
+#############################
+
+      ██████╗░██╗░█████╗░
+      ██╔══██╗██║██╔══██╗
+      ██████╔╝██║██║░░╚═╝
+      ██╔═══╝░██║██║░░██╗
+      ██║░░░░░██║╚█████╔╝
+      ╚═╝░░░░░╚═╝░╚════╝░
+         
+Options used:
+   --device cpu
+   --full-assembly
+   --no-partial-assembly
+   --dimension 2
+   --order 1
+   --num-x 32
+   --num-y 32
+   --num-z 100
+   --charge 0.00118164
+   --mass 0.00118164
+   --time-step 0.1
+   --num-timesteps 10
+   --num-particles 409600
+   --k 0.285599
+   --alpha 0.05
+   --ordering 1
+   --redist-interval 1
+   --output-csv-interval -1
+   --no-visualization
+   --send-port 19916
+   --reproduce
+Device configuration: cpu
+Memory configuration: host-std
+Total charge: 484, Domain volume: 484, Neutralizing constant: -1
+Further updates will use this precomputed neutralizing constant.
+Step: 1 | Time: 0.1 | Time per step: 0.704219
+Kinetic energy: 484.208	Field energy: 7.2893	Total energy: 491.497
+Step: 2 | Time: 0.2 | Time per step: 0.708807
+Kinetic energy: 484.311	Field energy: 7.18476	Total energy: 491.496
+Step: 3 | Time: 0.3 | Time per step: 0.711702
+Kinetic energy: 484.568	Field energy: 6.92838	Total energy: 491.496
+Step: 4 | Time: 0.4 | Time per step: 0.712573
+Kinetic energy: 484.964	Field energy: 6.53039	Total energy: 491.495
+Step: 5 | Time: 0.5 | Time per step: 0.71302
+Kinetic energy: 485.482	Field energy: 6.01109	Total energy: 491.493
+Step: 6 | Time: 0.6 | Time per step: 0.713366
+Kinetic energy: 486.098	Field energy: 5.39433	Total energy: 491.493
+Step: 7 | Time: 0.7 | Time per step: 0.713481
+Kinetic energy: 486.785	Field energy: 4.70615	Total energy: 491.491
+Step: 8 | Time: 0.8 | Time per step: 0.713655
+Kinetic energy: 487.512	Field energy: 3.97614	Total energy: 491.488
+Step: 9 | Time: 0.9 | Time per step: 0.713807
+Kinetic energy: 488.249	Field energy: 3.2369	Total energy: 491.486
+Step: 10 | Time: 1 | Time per step: 0.713923
+Kinetic energy: 488.965	Field energy: 2.51996	Total energy: 491.485
+
+=========== PIC Timing Summary (max over ranks) ===========
+  1. Calculate forces (E@p)         2.2907 s     32.0 %
+  2. Push particles                 0.2338 s      3.3 %
+  3. Redistribute particles         0.0047 s      0.1 %
+  4. Update fields (Poisson)        0.5954 s      8.3 %
+  5. FindPoints (locate p)          4.0973 s     57.3 %
+  Other (I/O, vis, energies)       -0.0694 s     -1.0 %
+  ---------------------------------------------------
+  Total time                        7.1526 s    100.0 %
+  Avg time per step                 0.7153 s
+============================================================
+#############################
+ranks=4  mesh=64^2  npt=1638400  q=m=0.00029541015634211245
+#############################
+
+      ██████╗░██╗░█████╗░
+      ██╔══██╗██║██╔══██╗
+      ██████╔╝██║██║░░╚═╝
+      ██╔═══╝░██║██║░░██╗
+      ██║░░░░░██║╚█████╔╝
+      ╚═╝░░░░░╚═╝░╚════╝░
+         
+Options used:
+   --device cpu
+   --full-assembly
+   --no-partial-assembly
+   --dimension 2
+   --order 1
+   --num-x 64
+   --num-y 64
+   --num-z 100
+   --charge 0.00029541
+   --mass 0.00029541
+   --time-step 0.1
+   --num-timesteps 10
+   --num-particles 1638400
+   --k 0.285599
+   --alpha 0.05
+   --ordering 1
+   --redist-interval 1
+   --output-csv-interval -1
+   --no-visualization
+   --send-port 19916
+   --reproduce
+Device configuration: cpu
+Memory configuration: host-std
+Total charge: 484, Domain volume: 484, Neutralizing constant: -1
+Further updates will use this precomputed neutralizing constant.
+Step: 1 | Time: 0.1 | Time per step: 0.905431
+Kinetic energy: 483.977	Field energy: 7.82945	Total energy: 491.807
+Step: 2 | Time: 0.2 | Time per step: 0.870283
+Kinetic energy: 484.061	Field energy: 7.74561	Total energy: 491.807
+Step: 3 | Time: 0.3 | Time per step: 0.857828
+Kinetic energy: 484.309	Field energy: 7.49724	Total energy: 491.806
+Step: 4 | Time: 0.4 | Time per step: 0.850545
+Kinetic energy: 484.709	Field energy: 7.09583	Total energy: 491.805
+Step: 5 | Time: 0.5 | Time per step: 0.846357
+Kinetic energy: 485.244	Field energy: 6.55892	Total energy: 491.803
+Step: 6 | Time: 0.6 | Time per step: 0.840311
+Kinetic energy: 485.891	Field energy: 5.91033	Total energy: 491.802
+Step: 7 | Time: 0.7 | Time per step: 0.834989
+Kinetic energy: 486.621	Field energy: 5.1789	Total energy: 491.8
+Step: 8 | Time: 0.8 | Time per step: 0.831411
+Kinetic energy: 487.401	Field energy: 4.39656	Total energy: 491.797
+Step: 9 | Time: 0.9 | Time per step: 0.831248
+Kinetic energy: 488.198	Field energy: 3.59751	Total energy: 491.795
+Step: 10 | Time: 1 | Time per step: 0.827664
+Kinetic energy: 488.978	Field energy: 2.81567	Total energy: 491.794
+
+=========== PIC Timing Summary (max over ranks) ===========
+  1. Calculate forces (E@p)         2.7633 s     33.3 %
+  2. Push particles                 0.2409 s      2.9 %
+  3. Redistribute particles         0.3018 s      3.6 %
+  4. Update fields (Poisson)        0.6681 s      8.1 %
+  5. FindPoints (locate p)          4.7120 s     56.8 %
+  Other (I/O, vis, energies)       -0.3958 s     -4.8 %
+  ---------------------------------------------------
+  Total time                        8.2903 s    100.0 %
+  Avg time per step                 0.8290 s
+============================================================
+#############################
+ranks=16  mesh=128^2  npt=6553600  q=m=7.3852539085528113e-05
+#############################
+
+      ██████╗░██╗░█████╗░
+      ██╔══██╗██║██╔══██╗
+      ██████╔╝██║██║░░╚═╝
+      ██╔═══╝░██║██║░░██╗
+      ██║░░░░░██║╚█████╔╝
+      ╚═╝░░░░░╚═╝░╚════╝░
+         
+Options used:
+   --device cpu
+   --full-assembly
+   --no-partial-assembly
+   --dimension 2
+   --order 1
+   --num-x 128
+   --num-y 128
+   --num-z 100
+   --charge 7.38525e-05
+   --mass 7.38525e-05
+   --time-step 0.1
+   --num-timesteps 10
+   --num-particles 6553600
+   --k 0.285599
+   --alpha 0.05
+   --ordering 1
+   --redist-interval 1
+   --output-csv-interval -1
+   --no-visualization
+   --send-port 19916
+   --reproduce
+Device configuration: cpu
+Memory configuration: host-std
+Total charge: 484, Domain volume: 484, Neutralizing constant: -1
+Further updates will use this precomputed neutralizing constant.
+Step: 1 | Time: 0.1 | Time per step: 1.17015
+Kinetic energy: 484.098	Field energy: 7.45399	Total energy: 491.552
+Step: 2 | Time: 0.2 | Time per step: 1.05184
+Kinetic energy: 484.175	Field energy: 7.37733	Total energy: 491.552
+Step: 3 | Time: 0.3 | Time per step: 1.01096
+Kinetic energy: 484.408	Field energy: 7.1433	Total energy: 491.551
+Step: 4 | Time: 0.4 | Time per step: 0.990213
+Kinetic energy: 484.788	Field energy: 6.76264	Total energy: 491.55
+Step: 5 | Time: 0.5 | Time per step: 0.975491
+Kinetic energy: 485.297	Field energy: 6.25213	Total energy: 491.549
+Step: 6 | Time: 0.6 | Time per step: 0.963004
+Kinetic energy: 485.912	Field energy: 5.63448	Total energy: 491.547
+Step: 7 | Time: 0.7 | Time per step: 0.955492
+Kinetic energy: 486.608	Field energy: 4.93722	Total energy: 491.545
+Step: 8 | Time: 0.8 | Time per step: 0.947635
+Kinetic energy: 487.352	Field energy: 4.19078	Total energy: 491.543
+Step: 9 | Time: 0.9 | Time per step: 0.94217
+Kinetic energy: 488.113	Field energy: 3.42757	Total energy: 491.541
+Step: 10 | Time: 1 | Time per step: 0.937194
+Kinetic energy: 488.859	Field energy: 2.68018	Total energy: 491.539
+
+=========== PIC Timing Summary (max over ranks) ===========
+  1. Calculate forces (E@p)         2.8975 s     30.9 %
+  2. Push particles                 0.2580 s      2.7 %
+  3. Redistribute particles         0.4210 s      4.5 %
+  4. Update fields (Poisson)        0.8276 s      8.8 %
+  5. FindPoints (locate p)          5.6030 s     59.7 %
+  Other (I/O, vis, energies)       -0.6178 s     -6.6 %
+  ---------------------------------------------------
+  Total time                        9.3892 s    100.0 %
+  Avg time per step                 0.9389 s
+============================================================
+#############################
+ranks=64  mesh=256^2  npt=26214400  q=m=1.8463134771382028e-05
+#############################
+
+      ██████╗░██╗░█████╗░
+      ██╔══██╗██║██╔══██╗
+      ██████╔╝██║██║░░╚═╝
+      ██╔═══╝░██║██║░░██╗
+      ██║░░░░░██║╚█████╔╝
+      ╚═╝░░░░░╚═╝░╚════╝░
+         
+Options used:
+   --device cpu
+   --full-assembly
+   --no-partial-assembly
+   --dimension 2
+   --order 1
+   --num-x 256
+   --num-y 256
+   --num-z 100
+   --charge 1.84631e-05
+   --mass 1.84631e-05
+   --time-step 0.1
+   --num-timesteps 10
+   --num-particles 26214400
+   --k 0.285599
+   --alpha 0.05
+   --ordering 1
+   --redist-interval 1
+   --output-csv-interval -1
+   --no-visualization
+   --send-port 19916
+   --reproduce
+Device configuration: cpu
+Memory configuration: host-std
+Total charge: 484, Domain volume: 484, Neutralizing constant: -1
+Further updates will use this precomputed neutralizing constant.
+Step: 1 | Time: 0.1 | Time per step: 1.53165
+Kinetic energy: 484.063	Field energy: 7.39334	Total energy: 491.456
+Step: 2 | Time: 0.2 | Time per step: 1.31346
+Kinetic energy: 484.139	Field energy: 7.31687	Total energy: 491.456
+Step: 3 | Time: 0.3 | Time per step: 1.23876
+Kinetic energy: 484.371	Field energy: 7.08415	Total energy: 491.455
+Step: 4 | Time: 0.4 | Time per step: 1.20432
+Kinetic energy: 484.748	Field energy: 6.70566	Total energy: 491.454
+Step: 5 | Time: 0.5 | Time per step: 1.17897
+Kinetic energy: 485.254	Field energy: 6.19836	Total energy: 491.453
+Step: 6 | Time: 0.6 | Time per step: 1.16336
+Kinetic energy: 485.866	Field energy: 5.58482	Total energy: 491.451
+Step: 7 | Time: 0.7 | Time per step: 1.15155
+Kinetic energy: 486.557	Field energy: 4.89226	Total energy: 491.449
+Step: 8 | Time: 0.8 | Time per step: 1.14081
+Kinetic energy: 487.296	Field energy: 4.15094	Total energy: 491.447
+Step: 9 | Time: 0.9 | Time per step: 1.13135
+Kinetic energy: 488.052	Field energy: 3.39314	Total energy: 491.445
+Step: 10 | Time: 1 | Time per step: 1.12425
+Kinetic energy: 488.791	Field energy: 2.65125	Total energy: 491.443
+
+=========== PIC Timing Summary (max over ranks) ===========
+  1. Calculate forces (E@p)         3.1641 s     28.1 %
+  2. Push particles                 0.2679 s      2.4 %
+  3. Redistribute particles         0.7077 s      6.3 %
+  4. Update fields (Poisson)        0.9426 s      8.4 %
+  5. FindPoints (locate p)          7.0523 s     62.6 %
+  Other (I/O, vis, energies)       -0.8735 s     -7.8 %
+  ---------------------------------------------------
+  Total time                       11.2611 s    100.0 %
+  Avg time per step                 1.1261 s
+============================================================

miniapps/plasma/pic/logs/weak-3d.log

@@ -0,0 +1,207 @@
+#############################
+ranks=1  mesh=16^3  npt=409600  q=m=0.0048447307312968462
+#############################
+
+      ██████╗░██╗░█████╗░
+      ██╔══██╗██║██╔══██╗
+      ██████╔╝██║██║░░╚═╝
+      ██╔═══╝░██║██║░░██╗
+      ██║░░░░░██║╚█████╔╝
+      ╚═╝░░░░░╚═╝░╚════╝░
+         
+Options used:
+   --device cpu
+   --full-assembly
+   --no-partial-assembly
+   --dimension 3
+   --order 1
+   --num-x 16
+   --num-y 16
+   --num-z 16
+   --charge 0.00484473
+   --mass 0.00484473
+   --time-step 0.02
+   --num-timesteps 10
+   --num-particles 409600
+   --k 0.5
+   --alpha 0.01
+   --ordering 1
+   --redist-interval 1
+   --output-csv-interval -1
+   --no-visualization
+   --send-port 19916
+   --reproduce
+Device configuration: cpu
+Memory configuration: host-std
+Total charge: 1984.4, Domain volume: 1984.4, Neutralizing constant: -1
+Further updates will use this precomputed neutralizing constant.
+Step: 1 | Time: 0.02 | Time per step: 1.36198
+Kinetic energy: 2979.46	Field energy: 1.65511	Total energy: 2981.11
+Step: 2 | Time: 0.04 | Time per step: 1.36538
+Kinetic energy: 2979.46	Field energy: 1.65308	Total energy: 2981.11
+Step: 3 | Time: 0.06 | Time per step: 1.36776
+Kinetic energy: 2979.46	Field energy: 1.65022	Total energy: 2981.11
+Step: 4 | Time: 0.08 | Time per step: 1.36901
+Kinetic energy: 2979.46	Field energy: 1.64676	Total energy: 2981.11
+Step: 5 | Time: 0.1 | Time per step: 1.37007
+Kinetic energy: 2979.47	Field energy: 1.64308	Total energy: 2981.11
+Step: 6 | Time: 0.12 | Time per step: 1.3701
+Kinetic energy: 2979.47	Field energy: 1.63865	Total energy: 2981.11
+Step: 7 | Time: 0.14 | Time per step: 1.36988
+Kinetic energy: 2979.48	Field energy: 1.6331	Total energy: 2981.11
+Step: 8 | Time: 0.16 | Time per step: 1.36996
+Kinetic energy: 2979.49	Field energy: 1.62567	Total energy: 2981.11
+Step: 9 | Time: 0.18 | Time per step: 1.37024
+Kinetic energy: 2979.49	Field energy: 1.61658	Total energy: 2981.11
+Step: 10 | Time: 0.2 | Time per step: 1.37026
+Kinetic energy: 2979.5	Field energy: 1.60714	Total energy: 2981.11
+
+=========== PIC Timing Summary (max over ranks) ===========
+  1. Calculate forces (E@p)         3.6764 s     26.8 %
+  2. Push particles                 0.3683 s      2.7 %
+  3. Redistribute particles         0.0047 s      0.0 %
+  4. Update fields (Poisson)        1.1195 s      8.1 %
+  5. FindPoints (locate p)          8.6048 s     62.6 %
+  Other (I/O, vis, energies)       -0.0324 s     -0.2 %
+  ---------------------------------------------------
+  Total time                       13.7413 s    100.0 %
+  Avg time per step                 1.3741 s
+============================================================
+#############################
+ranks=8  mesh=32^3  npt=3276800  q=m=0.00060559134141210578
+#############################
+
+      ██████╗░██╗░█████╗░
+      ██╔══██╗██║██╔══██╗
+      ██████╔╝██║██║░░╚═╝
+      ██╔═══╝░██║██║░░██╗
+      ██║░░░░░██║╚█████╔╝
+      ╚═╝░░░░░╚═╝░╚════╝░
+         
+Options used:
+   --device cpu
+   --full-assembly
+   --no-partial-assembly
+   --dimension 3
+   --order 1
+   --num-x 32
+   --num-y 32
+   --num-z 32
+   --charge 0.000605591
+   --mass 0.000605591
+   --time-step 0.02
+   --num-timesteps 10
+   --num-particles 3276800
+   --k 0.5
+   --alpha 0.01
+   --ordering 1
+   --redist-interval 1
+   --output-csv-interval -1
+   --no-visualization
+   --send-port 19916
+   --reproduce
+Device configuration: cpu
+Memory configuration: host-std
+Total charge: 1984.4, Domain volume: 1984.4, Neutralizing constant: -1
+Further updates will use this precomputed neutralizing constant.
+Step: 1 | Time: 0.02 | Time per step: 1.67352
+Kinetic energy: 2978.05	Field energy: 0.799737	Total energy: 2978.85
+Step: 2 | Time: 0.04 | Time per step: 1.63226
+Kinetic energy: 2978.05	Field energy: 0.800322	Total energy: 2978.85
+Step: 3 | Time: 0.06 | Time per step: 1.60665
+Kinetic energy: 2978.05	Field energy: 0.799868	Total energy: 2978.85
+Step: 4 | Time: 0.08 | Time per step: 1.59359
+Kinetic energy: 2978.05	Field energy: 0.798389	Total energy: 2978.85
+Step: 5 | Time: 0.1 | Time per step: 1.58883
+Kinetic energy: 2978.05	Field energy: 0.796093	Total energy: 2978.85
+Step: 6 | Time: 0.12 | Time per step: 1.58967
+Kinetic energy: 2978.06	Field energy: 0.793324	Total energy: 2978.85
+Step: 7 | Time: 0.14 | Time per step: 1.58496
+Kinetic energy: 2978.06	Field energy: 0.789973	Total energy: 2978.85
+Step: 8 | Time: 0.16 | Time per step: 1.58093
+Kinetic energy: 2978.06	Field energy: 0.786052	Total energy: 2978.85
+Step: 9 | Time: 0.18 | Time per step: 1.57764
+Kinetic energy: 2978.07	Field energy: 0.781706	Total energy: 2978.85
+Step: 10 | Time: 0.2 | Time per step: 1.57648
+Kinetic energy: 2978.07	Field energy: 0.776765	Total energy: 2978.85
+
+=========== PIC Timing Summary (max over ranks) ===========
+  1. Calculate forces (E@p)         4.3329 s     27.4 %
+  2. Push particles                 0.3601 s      2.3 %
+  3. Redistribute particles         0.3988 s      2.5 %
+  4. Update fields (Poisson)        1.3243 s      8.4 %
+  5. FindPoints (locate p)          9.8843 s     62.5 %
+  Other (I/O, vis, energies)       -0.4894 s     -3.1 %
+  ---------------------------------------------------
+  Total time                       15.8110 s    100.0 %
+  Avg time per step                 1.5811 s
+============================================================
+#############################
+ranks=64  mesh=64^3  npt=26214400  q=m=7.5698917676513222e-05
+#############################
+
+      ██████╗░██╗░█████╗░
+      ██╔══██╗██║██╔══██╗
+      ██████╔╝██║██║░░╚═╝
+      ██╔═══╝░██║██║░░██╗
+      ██║░░░░░██║╚█████╔╝
+      ╚═╝░░░░░╚═╝░╚════╝░
+         
+Options used:
+   --device cpu
+   --full-assembly
+   --no-partial-assembly
+   --dimension 3
+   --order 1
+   --num-x 64
+   --num-y 64
+   --num-z 64
+   --charge 7.56989e-05
+   --mass 7.56989e-05
+   --time-step 0.02
+   --num-timesteps 10
+   --num-particles 26214400
+   --k 0.5
+   --alpha 0.01
+   --ordering 1
+   --redist-interval 1
+   --output-csv-interval -1
+   --no-visualization
+   --send-port 19916
+   --reproduce
+Device configuration: cpu
+Memory configuration: host-std
+Total charge: 1984.4, Domain volume: 1984.4, Neutralizing constant: -1
+Further updates will use this precomputed neutralizing constant.
+Step: 1 | Time: 0.02 | Time per step: 2.28528
+Kinetic energy: 2977.19	Field energy: 0.636153	Total energy: 2977.83
+Step: 2 | Time: 0.04 | Time per step: 2.062
+Kinetic energy: 2977.19	Field energy: 0.635865	Total energy: 2977.83
+Step: 3 | Time: 0.06 | Time per step: 1.99404
+Kinetic energy: 2977.19	Field energy: 0.634917	Total energy: 2977.83
+Step: 4 | Time: 0.08 | Time per step: 1.95795
+Kinetic energy: 2977.19	Field energy: 0.63339	Total energy: 2977.83
+Step: 5 | Time: 0.1 | Time per step: 1.93115
+Kinetic energy: 2977.2	Field energy: 0.631344	Total energy: 2977.83
+Step: 6 | Time: 0.12 | Time per step: 1.91559
+Kinetic energy: 2977.2	Field energy: 0.628712	Total energy: 2977.83
+Step: 7 | Time: 0.14 | Time per step: 1.90556
+Kinetic energy: 2977.2	Field energy: 0.625445	Total energy: 2977.83
+Step: 8 | Time: 0.16 | Time per step: 1.8954
+Kinetic energy: 2977.21	Field energy: 0.621599	Total energy: 2977.83
+Step: 9 | Time: 0.18 | Time per step: 1.88887
+Kinetic energy: 2977.21	Field energy: 0.61722	Total energy: 2977.83
+Step: 10 | Time: 0.2 | Time per step: 1.8827
+Kinetic energy: 2977.21	Field energy: 0.612302	Total energy: 2977.83
+
+=========== PIC Timing Summary (max over ranks) ===========
+  1. Calculate forces (E@p)         5.0533 s     26.8 %
+  2. Push particles                 0.4047 s      2.1 %
+  3. Redistribute particles         0.6334 s      3.4 %
+  4. Update fields (Poisson)        1.8219 s      9.6 %
+  5. FindPoints (locate p)         11.7914 s     62.5 %
+  Other (I/O, vis, energies)       -0.8252 s     -4.4 %
+  ---------------------------------------------------
+  Total time                       18.8794 s    100.0 %
+  Avg time per step                 1.8879 s
+============================================================

miniapps/plasma/pic/scaling.sh

@@ -0,0 +1,31 @@
+nx=32
+ny=32
+ 
+for t in 1 2 4 8 16 24 32 64
+do
+
+# echo $t
+
+nz=$(($t * 4))
+npt=$(($nx * $ny * $nz * 100))
+
+echo "#############################"
+echo $t
+echo "#############################"
+
+srun -n $t ./electrostatic-pic \
+    -no-vis \
+    -rdi 1 \
+    -dim 3 \
+    -npt $npt \
+    -k 0.5 -a 0.01 \
+    -nt 10 \
+    -nx $nx -ny $ny -nz $nz \
+    -O 1 \
+    -q 0.00004844730731 \
+    -m 0.00004844730731 \
+    -oci -1 \
+    -dt 0.02 \
+    -fa -d cpu
+
+done

miniapps/plasma/pic/script.sh

@@ -0,0 +1,14 @@
+srun -p 64 ./electrostatic-pic \
+    -no-vis \
+    -rdi 1 \
+    -dim 3 \
+    -npt 26214400 \
+    -k 0.5 -a 0.01 \
+    -nt 10 \
+    -nx 32 -ny 32 -nz 128  \
+    -O 1 \
+    -q 0.00004844730731 \
+    -m 0.00004844730731 \
+    -oci -1 \
+    -dt 0.02 \
+    -fa -d cuda

miniapps/plasma/pic/weak-scale-2d.sh

@@ -0,0 +1,65 @@
+#!/bin/bash
+#
+# 2D weak-scaling sweep for electrostatic-pic.
+#
+# The domain is always a square of side L = 2*pi/k, so the mesh is refined
+# ISOTROPICALLY and the rank count grows by 4x per step (2x per direction):
+#
+#   ranks:            1        4       16       64
+#   mesh:          32^2     64^2    128^2    256^2
+#   elems/rank:    1024     1024     1024     1024
+#   parts/rank:  409600   409600   409600   409600   (400 particles/cell,
+#                                                      as in the reference
+#                                                      2D Landau run)
+#
+# IMPORTANT: q and m are NOT free constants in this test case. The Landau
+# setup assumes unit charge density (unit plasma frequency), i.e.
+#     q = m = L^2 / npt
+# (check: the reference run has q*npt = 0.001181640625 * 409600 = L^2).
+# Since npt changes with scale, q and m are recomputed for every run below.
+
+k=0.2855993321
+alpha=0.05
+dt=0.1
+
+ranks=(1 4 16 64)
+cells=(32 64 128 256)   # n^2 mesh
+
+ppc=400                 # particles per cell (matches the reference run)
+nt=10                   # long enough that one-time setup doesn't dominate
+
+# Set to your node's core count so all runs pack nodes identically.
+PPN=""
+
+for i in "${!ranks[@]}"
+do
+    t=${ranks[$i]}
+    n=${cells[$i]}
+    npt=$(($n * $n * $ppc))
+
+    # q = m = L^2 / npt, with L = 2*pi/k  (unit charge density)
+    q=$(awk -v k=$k -v npt=$npt \
+        'BEGIN { L = 2*atan2(0,-1)/k; printf "%.17g", L*L/npt }')
+
+    echo "#############################"
+    echo "ranks=$t  mesh=${n}^2  npt=$npt  q=m=$q"
+    echo "#############################"
+
+    srun -n $t ${PPN:+--ntasks-per-node=$PPN} --cpu-bind=cores \
+        ./electrostatic-pic \
+        -no-vis \
+        -rdi 1 \
+        -dim 2 \
+        -npt $npt \
+        -k $k -a $alpha \
+        -nt $nt \
+        -nx $n -ny $n \
+        -O 1 \
+        -q $q \
+        -m $q \
+        -oci -1 \
+        -dt $dt \
+        -fa -d cpu
+
+
+done

miniapps/plasma/pic/weak-scale-3d.sh

@@ -0,0 +1,66 @@
+#!/bin/bash
+#
+# 3D weak-scaling sweep for electrostatic-pic.
+#
+# The domain is always a cube of side L = 2*pi/k, so the mesh is refined
+# ISOTROPICALLY and the rank count grows by 8x per step (2x per direction):
+#
+#   ranks:            1        8       64
+#   mesh:          16^3     32^3     64^3
+#   elems/rank:    4096     4096     4096
+#   parts/rank:  409600   409600   409600      (100 particles per cell)
+#
+# IMPORTANT: q and m are NOT free constants in this test case. The Landau
+# setup assumes unit charge density (unit plasma frequency), i.e.
+#     q = m = L^3 / npt
+# (check: the reference 3D run has q*npt = 0.00004844730731 * 40960000 = L^3).
+# Since npt changes with scale, q and m are recomputed for every run below.
+
+k=0.5
+alpha=0.01
+dt=0.02
+
+ranks=(1 8 64)
+cells=(16 32 64)     # n^3 mesh; shift to (32 64 128) for finer resolution --
+                     # per-rank load then grows 8x but stays constant across
+                     # the sweep, which is all weak scaling requires.
+
+ppc=100              # particles per cell
+nt=10                # long enough that one-time setup doesn't dominate
+
+# Set to your node's core count so all runs pack nodes identically.
+PPN=""
+
+for i in "${!ranks[@]}"
+do
+    t=${ranks[$i]}
+    n=${cells[$i]}
+    npt=$(($n * $n * $n * $ppc))
+
+    # q = m = L^3 / npt, with L = 2*pi/k  (unit charge density)
+    q=$(awk -v k=$k -v npt=$npt \
+        'BEGIN { L = 2*atan2(0,-1)/k; printf "%.17g", L*L*L/npt }')
+
+    echo "#############################"
+    echo "ranks=$t  mesh=${n}^3  npt=$npt  q=m=$q"
+    echo "#############################"
+
+    # NOTE: -n sets the task count; -p selects a PARTITION.
+    srun -n $t ${PPN:+--ntasks-per-node=$PPN} --cpu-bind=cores \
+        ./electrostatic-pic \
+        -no-vis \
+        -rdi 1 \
+        -dim 3 \
+        -npt $npt \
+        -k $k -a $alpha \
+        -nt $nt \
+        -nx $n -ny $n -nz $n \
+        -O 1 \
+        -q $q \
+        -m $q \
+        -oci -1 \
+        -dt $dt \
+        -fa -d cuda
+
+
+done