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

[WeiqunZhang]

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

[zingale]

this is worried about finest_level changing during the hscf solve -- does it change? are we creating new levels here?

[maximumcats]

My intuition here is that this is absolutely a thing that could happen, so I do think there should be a PR to address this (once I've had time to review the proposed patch).