From 5fde90c970bba1059e0834a7b5dce6a46b47a33b Mon Sep 17 00:00:00 2001
From: Liwei Ji <jiliwei.phys@gmail.com>
Date: Fri, 13 Mar 2026 18:56:09 -0400
Subject: [PATCH] TestNuPcsArdBH: iterate all AMR levels when randomizing
 particle momenta

InitRandom distributes particles across all AMR levels via Where(), but
the momentum randomization loop only touched level 0, leaving fine-level
particles with zero momenta. Wrap Phase 2 in a level loop over
finestLevel() and add a defensive Redistribute() after initialization.
Also remove unused lev and dx variables.
---
 TestNuPcsArdBH/src/testnupcsardbh.cxx | 45 ++++++++++++++-------------
 1 file changed, 24 insertions(+), 21 deletions(-)

diff --git a/TestNuPcsArdBH/src/testnupcsardbh.cxx b/TestNuPcsArdBH/src/testnupcsardbh.cxx
index f40972f..1587e26 100644
--- a/TestNuPcsArdBH/src/testnupcsardbh.cxx
+++ b/TestNuPcsArdBH/src/testnupcsardbh.cxx
@@ -23,7 +23,6 @@ extern "C" void TestNuPcsArdBH_InitParticles(CCTK_ARGUMENTS) {
 
   for (int patch = 0; patch < ghext->num_patches(); ++patch) {
     auto &pc = g_nupcs.at(patch);
-    const int lev = 0;
 
     // Phase 1: Use AMReX InitRandom for positions + zero momenta
     // For SoA particles, InitRandom randomly generates positions into SoA
@@ -34,7 +33,6 @@ extern "C" void TestNuPcsArdBH_InitParticles(CCTK_ARGUMENTS) {
         {}, {}, {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0}, {0, 0}};
 
     // Constrain positions to the target cell
-    const auto dx = pc->Geom(lev).CellSizeArray();
     RealBox containing_bx(
         {containing_bx_xmin, containing_bx_ymin, containing_bx_zmin},
         {containing_bx_xmax, containing_bx_ymax, containing_bx_zmax});
@@ -42,25 +40,30 @@ extern "C" void TestNuPcsArdBH_InitParticles(CCTK_ARGUMENTS) {
     pc->InitRandom(num_particles, random_seed, pdata, true, containing_bx);
 
     // Phase 2: Randomize momenta with isotropic sampling
-    for (ParIterType pti(*pc, lev); pti.isValid(); ++pti) {
-      auto &particle_tile = pti.GetParticleTile();
-      auto ptd = particle_tile.getParticleTileData();
-      const int np = pti.numParticles();
-
-      amrex::ParallelForRNG(
-          np, [=] AMREX_GPU_DEVICE(int i,
-                                   amrex::RandomEngine const &engine) noexcept {
-            const Real nu = 1.0;
-            Real costh = Random(engine) * 2 - 1;
-            Real ph = Random(engine) * (2 * M_PI);
-            Real sinth = std::sqrt(amrex::max(Real(0), 1 - costh * costh));
-            ptd.rdata(PIdx::px)[i] = nu * sinth * std::cos(ph);
-            ptd.rdata(PIdx::py)[i] = nu * sinth * std::sin(ph);
-            ptd.rdata(PIdx::pz)[i] = nu * costh;
-            ptd.idata(PIdxInt::species)[i] = 0;
-            ptd.idata(PIdxInt::cell_id)[i] = 0;
-          });
-    }
+    for (int lev = 0; lev <= pc->finestLevel(); ++lev) {
+      for (ParIterType pti(*pc, lev); pti.isValid(); ++pti) {
+        auto &particle_tile = pti.GetParticleTile();
+        auto ptd = particle_tile.getParticleTileData();
+        const int np = pti.numParticles();
+
+        amrex::ParallelForRNG(
+            np, [=] AMREX_GPU_DEVICE(int i,
+                                     amrex::RandomEngine const &engine) noexcept {
+              const Real nu = 1.0;
+              Real costh = Random(engine) * 2 - 1;
+              Real ph = Random(engine) * (2 * M_PI);
+              Real sinth = std::sqrt(amrex::max(Real(0), 1 - costh * costh));
+              ptd.rdata(PIdx::px)[i] = nu * sinth * std::cos(ph);
+              ptd.rdata(PIdx::py)[i] = nu * sinth * std::sin(ph);
+              ptd.rdata(PIdx::pz)[i] = nu * costh;
+              ptd.idata(PIdxInt::species)[i] = 0;
+              ptd.idata(PIdxInt::cell_id)[i] = 0;
+            });
+      }
+    } // for lev
+
+    // Defensive redistribute after initialization
+    pc->Redistribute();
   } // for patch
 
   // IO