`do_hscf_solve()` caches `finest_level` once, but calls `RegridOnly()` inside the iteration loop

Summary

do_hscf_solve() stores const int finest_level = parent->finestLevel(); before entering the SCF iteration loop, but each iteration can call parent->RegridOnly(time, do_io). If regridding changes AMR level topology, subsequent loops and the gravity solve still use the stale finest_level.

Location

• Source/scf/scf_relax.cpp:63
• Source/scf/scf_relax.cpp:156
• Source/scf/scf_relax.cpp:600
• Source/scf/scf_relax.cpp:605

Problem Details

Current control flow:

const int finest_level = parent->finestLevel();
...
while (ctr <= scf_max_iterations) {
    ...
    if (finest_level > 0) {
        parent->RegridOnly(time, do_io);
    }
    gravity->multilevel_solve_for_new_phi(0, finest_level);
    for (int lev = 0; lev <= finest_level; ++lev) { ... }
}

This assumes AMR level count never changes during the SCF loop. If RegridOnly() refines/coarsens levels, the SCF iteration continues with outdated bounds.

Impact

• New refined levels can be ignored by subsequent SCF operations.
• Coarsened-away levels can still be referenced via stale loop bounds.
• Gravity updates can be computed with an out-of-date finest level index.

Suggested Patch

Re-evaluate parent->finestLevel() each iteration and after regridding, then use the refreshed bound for all per-level loops and gravity calls.

diff --git a/Source/scf/scf_relax.cpp b/Source/scf/scf_relax.cpp
--- a/Source/scf/scf_relax.cpp
+++ b/Source/scf/scf_relax.cpp
@@
-    const int finest_level = parent->finestLevel();
-    const int n_levs = finest_level + 1;
+    int finest_level = parent->finestLevel();
@@
-    Vector< std::unique_ptr<MultiFab> > psi(n_levs);
-    Vector< std::unique_ptr<MultiFab> > enthalpy(n_levs);
-    Vector< std::unique_ptr<MultiFab> > state_vec(n_levs);
-    Vector< std::unique_ptr<MultiFab> > phi(n_levs);
+    Vector< std::unique_ptr<MultiFab> > psi;
+    Vector< std::unique_ptr<MultiFab> > enthalpy;
+    Vector< std::unique_ptr<MultiFab> > state_vec;
+    Vector< std::unique_ptr<MultiFab> > phi;
@@  
     while (ctr <= scf_max_iterations) {
+        finest_level = parent->finestLevel();
+        const int n_levs = finest_level + 1;
+        psi.resize(n_levs);
+        enthalpy.resize(n_levs);
+        state_vec.resize(n_levs);
+        phi.resize(n_levs);
@@
         if (finest_level > 0) {
             parent->RegridOnly(time, do_io);
         }

-        gravity->multilevel_solve_for_new_phi(0, finest_level);
+        finest_level = parent->finestLevel();
+        gravity->multilevel_solve_for_new_phi(0, finest_level);

Prepared by Codex

Start in Source/scf/scf_relax.cpp at do_hscf_solve(), especially lines 63, 156, 600, and 605. Trace how RegridOnly() changes the level topology and how the per-level containers, gravity solve, and loops use finest_level. Done means subsequent SCF work uses the current level bound after regridding without referencing invalid or omitted levels.