Fix issues #240-#249

[asalmgren]

Fixes #240, Fixes #241, Fixes #242, Fixes #243, Fixes #244, Fixes #245, Fixes #246, Fixes #247, Fixes #248, Fixes #249.

Medium severity

• ApplyCCProjection still allocates cc_phi/cc_gphi with new and never deletes them: the per-call leak reported as #191 was marked closed by PR #213 but the code was not changed #240 ApplyCCProjection allocated cc_phi/cc_gphi with new and never deleted them, leaking 1 + AMREX_SPACEDIM cell-centred fields per call (device memory on GPU). This was reported as ApplyCCProjection leaks cc_phi/cc_gphi MultiFabs (raw new, never deleted) every call #191 and closed by PR fix issues 171 176 181 188 189 191 #213 without the fix landing. They are now held in Vector<MultiFab> and handed to Hydro with GetVecOfPtrs, matching how inv_rho/m_fluxes are already done in the same function.

• DiffusionScalarOp: the single EB scalar op is shared by the tracer and temperature passes, so an EB Dirichlet BC set by one pass (tracer_eb or temperature_eb) is silently applied to the other, with the other field's values and coefficient #241 DiffusionScalarOp owned one EB solve op and one EB apply op, shared by the tracer and temperature passes. MLEBABecLap has no way to revert an EB Dirichlet BC back to homogeneous Neumann, so with only one of eb_flow.tracer/eb_flow.temperature given, the other field silently solved with a Dirichlet EB pinned to the first field's value and coefficient. The temperature now gets its own m_eb_tem_solve_op/m_eb_tem_apply_op, selected by an is_temperature flag passed from diffuse_temperature/compute_laps_T. Pure plumbing — the per-component logic is untouched.

• init_burggraf evaluates the exact solution at x-0.5, y-0.5 (regression from PR #160) while the probtype-16 forcing, lid profile and DiffFromExact use [0,1]: benchmark.burggraf starts half a cavity out of place #242 init_burggraf evaluated the exact solution at x-0.5, y-0.5 while the probtype-16 forcing, the lid profile and DiffFromExact all use [0,1]^2, so the IC was the exact solution translated by (-1/2,-1/2) — an O(1) slip velocity on no-slip walls. A regression from PR wrap constants with Real() #160 (a cosmetic Real()-wrapping pass); the shift is removed.

• incflo_PCInit 2D cylinder cull: ALWAYS_ASSERT rejects cylinder.direction=2 (the only value its disk formula fits) and culls a disk for directions 0/1 where AdvectWithFlow and EB2::CylinderIF use slabs (regression from PR #213) #243 The 2D cylinder cull in incflo_PCInit was self-contradictory: the assert admitted only direction 0 or 1, but the formula was the out-of-plane (direction 2) disk distance and ignored the direction it had just validated. Meanwhile AdvectWithFlow and EB2::CylinderIF treat 2 as a disk and 0/1 as slabs. So direction = 2 aborted outright, and 0/1 seeded a disk while every later step confined particles to a slab. Init now uses the same convention as AdvectWithFlow.

• Fresh start seeds tracer particles on level 0 and never redistributes: initial ErrorEst counts zero particles on level >= 1 (no level 2 from particle tagging at t=0) and the first advection uses level-0 MAC velocities for every particle #244 A fresh start seeds tracer particles on level 0 only and nothing on that path ever calls Redistribute() (the restart path gets it from ParticleContainer::Restart). So initial ErrorEst counted zero particles on level >= 1 — with refine_particles as the criterion the initial hierarchy stopped at level 1 regardless of amr.max_level — and the first advection moved every particle under a fine patch with level-0 MAC velocities. Now redistributes after each level is created from scratch.

Low severity

• ReadParameters accepts steady_state=1 but SteadyStateReached() is an Abort stub: run takes one step, then dies #245 steady_state = 1 was read and accepted, but SteadyStateReached() is an amrex::Abort stub with the real body #if 0'd. Evolve() masks its entry test with && !m_steady_state, so the run paid for full initialization plus one complete time step and then died on what reads as an internal TODO (with the stop_time cap disabled for that step, and steady_state_tol silently ignored). Now refused at read time, exactly as PR fix issue 177 178 179 #212 did for amr.KE_int.

• InitData never records m_last_chk (and restart sets m_last_plt=0): Evolve's final-output block rewrites the step's checkpoint/plotfile, renaming the restart checkpoint to *.old.* #246 InitData wrote chk00000 without setting m_last_chk, and the restart branch left m_last_chk = -1 while recording m_last_plt = 0 instead of m_nstep. When the loop then takes no steps, Evolve's final-output block rewrote the step's checkpoint and plotfile — UtilCreateCleanDirectory renames the existing directory to *.old.* first, so the very checkpoint the run restarted from was renamed and replaced by an identical copy. Both counters now record the step of each write.

• ReadCheckpointFile takes finest_level from the header without checking amr.max_level: restarting with a smaller max_level writes grids[]/m_leveldata[] out of bounds #247 ReadCheckpointFile took finest_level from the header and looped MakeNewLevelFromScratch over it, but grids/dmap/m_leveldata/m_t_new are all sized max_level + 1 from the current inputs. Restarting with a smaller amr.max_level wrote past the end of those vectors — silent heap corruption in a release build. Now aborts with a diagnostic naming both values.

• update_temperature corrector passes the tracer laps (get_laps_new) to compute_laps_T instead of get_laps_tem_new (latent: MOL+temperature is refused at read time) #248 The StepType::Corrector branch of update_temperature passed get_laps_new() (the tracer Laplacian, m_ntrac components) to compute_laps_T instead of get_laps_tem_new(). Unreachable today — ReadParameters refuses use_temperature with advection_type = MOL — but it would clobber the tracer's laps and alias out-of-bounds components for ntrac > 1 the moment that guard is lifted. One-token fix.

• update_density corrector refreshes no ghost cells of density/density_nph, but ApplyCCProjection averages density_nph to faces through its ghost cells: stale predictor values on inter-grid faces (MOL + use_cc_proj + variable density) #249 update_density used ng = 0 in the corrector, so the ghost cells of density/density_nph still held predictor values. ApplyCCProjection averages density_nph to faces via amrex_avg_cc_to_fc, which reads the ghost cell on the low face of every grid, so the dt/rho face coefficient was assembled from two different time levels on faces coinciding with grid boundaries — making MOL + use_cc_proj + variable-density runs depend on amr.max_grid_size. Now one ghost cell in both passes, so the corrector's fillpatch_density refreshes them.

Testing

Built clean with no warnings in four configurations: test_2d (2D EB), test_no_eb_2d (2D no-EB), test_3d (3D EB), test_no_eb_3d (3D no-EB).

Note that #242 (burggraf IC) and #249 (corrector density ghost cells) change numerical output, so the affected benchmarks will need regolding. benchmark.burggraf in particular: the deck runs to stop_time = 5 at Re = 100, so the Norm0 printed by test_convergence_burggraf was previously measuring the decay of the spurious transient rather than the discretisation error of the steady solution.

🤖 Generated with Claude Code

[asalmgren]

Fixes #240, fixes #241, fixes #242, fixes #243, fixes #244, fixes #245, fixes #246, fixes #247, fixes#248, fixes#249