AMReX-Astro / AMReX-Astro/Castro
SCF diagnostics use stale `phi` scratch data after `multilevel_solve_for_new_phi()`
Nobody has claimed this yet.
- Dominant language
- C++
- Stars
- 340
- Forks
- 105
- Avg merge
- 3d 8h
- Merged PRs (30d)
- 8
Description
Summary
do_hscf_solve() copies PhiGrav_Type into local scratch phi[lev] early in each iteration, then later calls gravity->multilevel_solve_for_new_phi(...), but the diagnostics still read from the old scratch phi[lev] instead of refreshed level data.
Location
Source/scf/scf_relax.cpp:203Source/scf/scf_relax.cpp:605Source/scf/scf_relax.cpp:639
Problem Details
The sequence in one SCF iteration is:
- Copy
get_new_data(PhiGrav_Type)intophi[lev]. - Modify density and solve gravity with
multilevel_solve_for_new_phi. - Compute
pot_engfromphi_arr = (*phi[lev])[mfi].array();.
So the printed potential energy and virial diagnostic are based on pre-solve phi, not the just-updated gravitational potential.
Impact
- Reported potential energy can lag the actual solved state.
- Virial error output can be misleading during relaxation monitoring.
- SCF convergence diagnostics become harder to interpret.
Suggested Patch
Use current level data in the diagnostic reduce loop (or refresh phi[lev] from PhiGrav_Type immediately after the gravity solve).
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
@@
gravity->multilevel_solve_for_new_phi(0, finest_level);
+ // Refresh local phi scratch with the newly solved potential.
+ for (int lev = 0; lev <= finest_level; ++lev) {
+ MultiFab::Copy((*phi[lev]), getLevel(lev).get_new_data(PhiGrav_Type), 0, 0, 1, 0);
+ if (lev < finest_level) {
+ const MultiFab& mask = getLevel(lev+1).build_fine_mask();
+ MultiFab::Multiply((*phi[lev]), mask, 0, 0, 1, 0);
+ }
+ }
+
// Update diagnostic quantities.
Prepared by Codex
Contributor guide
No contributing guide indexed for this repository
First steps
- Read the whole issue, then the project's contributing guide.
- Comment on the issue to say you are picking it up — it saves two people doing the same work.
- Fork the repository and make your change on a branch.
- Open a pull request that references the issue number.
Research direction
Start in Source/scf/scf_relax.cpp, especially do_hscf_solve() around lines 203, 605, and 639, and trace when phi[lev] is copied versus when multilevel_solve_for_new_phi() runs. Confirm the diagnostic reduce loop uses the refreshed potential, then validate the potential-energy and virial outputs with an SCF run.
Written by the indexing model from the issue text.
Assessment
- Tech stack
- cpp
- Domain
- hpc
- Issue type
- Bug
- Difficulty
- 2/5
- Estimated time
- 1-3 hours
- Activity status
- Stale
- Clarity
- Clearly specified
- Newbie friendliness
- 55/100