From 5e62c2a9b0a91bc3db3612190da868f52807e1ae Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Mon, 14 Sep 2026 06:39:31 -0700 Subject: [PATCH] fix issues 111 112 113 114 115 116 - #111 construct_hydro_source (EB): decide the regular-box fallback on grow(bx,4), since the non-EB kernels read q out to distance 4 and flatten() disagrees with flatten_eb() within 3 cells of a cut cell. Without this the fluxes on a shared face differ between the two sides and the answer depends on the box layout. - #112 main.cpp moving-EB loop: stop calling finalize_EB2() every step, which destroyed the EB2 index space that the live level factories still point into. Build the new space only once time has advanced and erase the previous one after the step, when the forced regrid has rebuilt every level. Abort up front unless amr.regrid_int = 1 and amr.force_regrid_level_zero = 1, which the scheme requires. - #113 clean_state now takes an explicit ghost count instead of using S.nGrow(): with CAMR.state_nghost > 0 it ran computeTemp over ghost cells of the state data that nothing ever fills. Also corrects the stale enforce_min_density comment, which described an S_old argument and a return value the function does not have and claimed it touches ghost zones. - #114 setPlotVariables no longer strips the species from an explicit amr.plot_vars list when CAMR.plot_rhoy was never set. - #115 read_tagging_params: the 24 tagging.* keys were parsed into a struct nothing reads, silently ignoring user thresholds. Refuse them with a message pointing at amr.refinement_indicators. - #116 delete Source/CAMR_initial_redist.cpp, an orphaned stub that is in no build file, declared in no header, called from nowhere, and would not compile. Built and ran Exec/Sod (2D, non-EB), Exec/SodPlusSphere (2D, EB) and Exec/MovingEBCases/MovingPlane_Shock (3D, moving EB). Co-Authored-By: Claude Opus 5 (1M context) --- Source/CAMR.H | 7 +- Source/CAMR.cpp | 30 +++++---- Source/CAMR_advance.cpp | 6 +- Source/CAMR_initial_redist.cpp | 67 -------------------- Source/Hydro/CAMR_construct_hydro_source.cpp | 9 ++- Source/Utils/CAMR_io.cpp | 9 ++- Source/Utils/Tagging.cpp | 42 ++++++------ Source/main.cpp | 29 ++++++++- 8 files changed, 88 insertions(+), 111 deletions(-) delete mode 100644 Source/CAMR_initial_redist.cpp diff --git a/Source/CAMR.H b/Source/CAMR.H index bf0ac39..47b3005 100644 --- a/Source/CAMR.H +++ b/Source/CAMR.H @@ -388,7 +388,12 @@ protected: static void error_setup (); void expand_state(amrex::MultiFab& S, const amrex::Real time, const int ng); - void clean_state(amrex::MultiFab& S); + + // Floor the density and renormalize the species on the valid region, then + // recompute T on the valid region plus ng ghost cells. ng must not exceed + // the number of ghost cells that have actually been filled (0 for state + // data, numGrow() for Sborder). + void clean_state(amrex::MultiFab& S, int ng); #ifdef CAMR_USE_MOVING_EB void ZeroingOutForPlotting(amrex::MultiFab& S); #endif diff --git a/Source/CAMR.cpp b/Source/CAMR.cpp index 7bb6799..c414ec2 100644 --- a/Source/CAMR.cpp +++ b/Source/CAMR.cpp @@ -704,7 +704,7 @@ CAMR::post_timestep(int /*iteration*/) // Clean up any aberrant state data generated by the reflux and average-down, // and then update quantities like temperature to be consistent. amrex::MultiFab& S_new_crse = get_new_data(State_Type); - clean_state(S_new_crse); + clean_state(S_new_crse, 0); #ifdef AMREX_USE_EB // If we redistribute then the ML redistribution algorithm may change @@ -715,7 +715,7 @@ CAMR::post_timestep(int /*iteration*/) if ( level < finest_level && (redistribution_type != "NoRedist") ) { CAMR& fine_level = getLevel(level + 1); amrex::MultiFab& S_new_fine = fine_level.get_new_data(State_Type); - clean_state(S_new_fine); + clean_state(S_new_fine, 0); } #endif @@ -948,15 +948,14 @@ CAMR::enforce_consistent_e(amrex::MultiFab& S) void CAMR::enforce_min_density(amrex::MultiFab& S_new) { - // This routine sets the density in S_new to be larger than the density - // floor. Note that it will operate everywhere on S_new, including ghost - // zones. S_old is present so that, after the hydro call, we know what the old - // density was so that we have a reference for comparison. If you are calling - // it elsewhere and there's no meaningful reference state, just pass in the - // same amrex::MultiFab twice. - // @return The return value is the the negative fractional change in the - // state that has the largest magnitude. If there is no reference state, this - // is meaningless. + // This routine resets any zone of S_new whose density is below the density + // floor. It operates on the valid region only; the ghost cells of S_new are + // left untouched. + // + // A zone below the floor is overwritten with a copy of whichever of its 27 + // neighbors within the same tile has the highest (pre-reset) density, or, if + // no such neighbor reaches the floor, with a state at (small_dens, + // small_temp) and zero velocity. #ifdef AMREX_USE_EB auto const& fact = dynamic_cast(S_new.Factory()); @@ -1276,18 +1275,21 @@ CAMR::expand_state(amrex::MultiFab& S, const amrex::Real time, const int ng) AmrLevel::FillPatch(*this,S,ng,time,State_Type,0,S.nComp()); - clean_state(S); + clean_state(S, ng); } void -CAMR::clean_state(amrex::MultiFab& S) +CAMR::clean_state(amrex::MultiFab& S, int ng) { + AMREX_ALWAYS_ASSERT(S.nGrow() >= ng); + // Enforce a minimum density. enforce_min_density(S); normalize_species(S); - int ng = S.nGrow(); + // Only touch ghost cells the caller has actually filled: the state data + // itself carries state_nghost unfilled ghost cells. computeTemp(S,ng); } diff --git a/Source/CAMR_advance.cpp b/Source/CAMR_advance.cpp index 15f9431..82d3c9d 100644 --- a/Source/CAMR_advance.cpp +++ b/Source/CAMR_advance.cpp @@ -80,7 +80,7 @@ CAMR::CAMR_advance (Real time, // trusted to respect the consistency between certain state variables // (e.g. UEINT and UEDEN) that we demand in every zone. - clean_state(get_old_data(State_Type)); + clean_state(get_old_data(State_Type), 0); MultiFab& S_old = get_old_data(State_Type); amrex::ignore_unused(S_old); @@ -158,14 +158,14 @@ CAMR::CAMR_advance (Real time, } // Sync up state after old sources and hydro source. - clean_state(S_new); + clean_state(S_new, 0); // "new source" is actually the correction to the old source we've already added for (int n = 0; n < src_list.size(); ++n) { construct_new_source(src_list[n], time, dt); MultiFab::Saxpy(S_new, dt, *new_sources[src_list[n]], 0, 0, NVAR, 0); - clean_state(S_new); + clean_state(S_new, 0); } Sborder.clear(); diff --git a/Source/CAMR_initial_redist.cpp b/Source/CAMR_initial_redist.cpp deleted file mode 100644 index 516fcab..0000000 --- a/Source/CAMR_initial_redist.cpp +++ /dev/null @@ -1,67 +0,0 @@ -#include "CAMR.H" -#include "CAMR_hydro.H" -#include "CAMR_utils_K.H" -#include "Godunov.H" -#include "MOL_umeth.H" -#include "CAMR_Constants.H" - -#include -#include - -using namespace amrex; - -#ifdef AMREX_USE_EB -void -CAMR::ApplyInitialMLRedistribution( Box const& bx, - Geometry const& geom, - const EBFArrayBoxFactory* ebfact, - const int* bclo, const int* bchi, - Array4 const& uin_arr, - Array4< Real> const& dsdt_arr, - Array4 const& vf_arr, - Array4 const& flag_arr, - int as_crse, Array4 const& drho_as_crse, - Array4 const& rrflag_as_crse, - int as_fine, Array4 const& dm_as_fine, - Array4 const& lev_mask, - const BCRec* bcs_d_ptr) -{ - BL_PROFILE_VAR("ApplyInitialMLRedist()", ApplyInitialMLRedist); - - const Box& bxg_ii = grow(bxg_i,1); - - Array4 AMREX_D_DECL(fcx, fcy, fcz), AMREX_D_DECL(apx, apy, apz), ccc; - AMREX_D_TERM(fcx = ebfact->getFaceCent()[0]->const_array(mfi);, - fcy = ebfact->getFaceCent()[1]->const_array(mfi);, - fcz = ebfact->getFaceCent()[2]->const_array(mfi);); - AMREX_D_TERM(apx = ebfact->getAreaFrac()[0]->const_array(mfi);, - apy = ebfact->getAreaFrac()[1]->const_array(mfi);, - apz = ebfact->getAreaFrac()[2]->const_array(mfi);); - ccc = ebfact->getCentroid().const_array(mfi); - - int l_ncomp = dsdt_arr.nComp(); - int level_mask_not_covered = CAMRConstants::level_mask_notcovered; - bool use_wts_in_divnc = false; - - // We need to set fac_for_redist to 1/2 for MOL because we will - // compute this twice per time step, so the contribution of - // each needs to be weighted by 1/2 - Real fac_for_redist = (do_mol) ? Real(0.5) : Real(1.0); - if (redistribution_type == "StateRedist") { - ApplyInitialMLRedistribution(bx, l_ncomp, - dsdt_arr, divc_arr, uin_arr, redistwgt_arr, - flag_arr, - AMREX_D_DECL(apx, apy, apz), - vf_arr, - AMREX_D_DECL(fcx, fcy, fcz), - ccc, bcs_d_ptr, geom, dt, - l_redistribution_type, - as_crse, drho_as_crse, rrflag_as_crse, - as_fine, dm_as_fine, lev_mask, - level_mask_not_covered, - fac_for_redist); - } - - BL_PROFILE_VAR_STOP(ApplyInitialMLRedist); -} -#endif diff --git a/Source/Hydro/CAMR_construct_hydro_source.cpp b/Source/Hydro/CAMR_construct_hydro_source.cpp index 9ac05f6..473835d 100644 --- a/Source/Hydro/CAMR_construct_hydro_source.cpp +++ b/Source/Hydro/CAMR_construct_hydro_source.cpp @@ -166,7 +166,14 @@ CAMR::construct_hydro_source (const MultiFab& S, ngrow_bx = 2; } const Box& bxg_i = grow(bx,ngrow_bx); - if (flagfab.getType(bxg_i) != FabType::regular) { + + // The non-EB kernels used in the fallback read q up to 4 cells outside + // bx (flatten: +-3 around the bxg1 ring) and flatten_eb() returns 1 + // within 3 cells of a cut cell, so the two paths give the same fluxes + // on the faces of bx only if grow(bx,4) is entirely regular. numGrow() + // is at least 5 in EB builds, so the flag fab covers this box. + const Box& bx_reg = grow(bx, amrex::max(ngrow_bx,4)); + if (flagfab.getType(bx_reg) != FabType::regular) { EBFluxRegister* fr_as_crse = nullptr; if (do_reflux && level < parent->finestLevel()) { diff --git a/Source/Utils/CAMR_io.cpp b/Source/Utils/CAMR_io.cpp index d939e9a..6044ef6 100644 --- a/Source/Utils/CAMR_io.cpp +++ b/Source/Utils/CAMR_io.cpp @@ -90,13 +90,18 @@ CAMR::setPlotVariables() amrex::ParmParse pp("CAMR"); + // Species are dropped from the default ("all state variables") plot list + // unless plot_rhoy is set. If the user gave an explicit amr.plot_vars list + // and did not say anything about plot_rhoy, honor that list as written + // instead of silently removing the species from it. bool plot_rhoy = false; - pp.query("plot_rhoy", plot_rhoy); + const bool have_plot_rhoy = pp.query("plot_rhoy", plot_rhoy); + const bool have_plot_vars = amrex::ParmParse("amr").contains("plot_vars"); if (plot_rhoy) { for (int i = 0; i < NUM_SPECIES; i++) { amrex::Amr::addStatePlotVar(desc_lst[State_Type].name(UFS + i)); } - } else { + } else if (have_plot_rhoy || !have_plot_vars) { for (int i = 0; i < NUM_SPECIES; i++) { amrex::Amr::deleteStatePlotVar(desc_lst[State_Type].name(UFS + i)); } diff --git a/Source/Utils/Tagging.cpp b/Source/Utils/Tagging.cpp index 2e0330a..43cb776 100644 --- a/Source/Utils/Tagging.cpp +++ b/Source/Utils/Tagging.cpp @@ -6,28 +6,28 @@ void CAMR::read_tagging_params() { + // Nothing in CAMR reads TaggingParm: CAMR::errorEst tags only through the + // amr.refinement_indicators mechanism (CAMR_error.cpp). Rather than query + // these keys -- which hides them from AMReX's unused-inputs report and + // silently ignores them -- refuse them with a pointer to what works. amrex::ParmParse pp("tagging"); - pp.query("denerr", tagging_parm->denerr); - pp.query("max_denerr_lev", tagging_parm->max_denerr_lev); - pp.query("dengrad", tagging_parm->dengrad); - pp.query("max_dengrad_lev", tagging_parm->max_dengrad_lev); + static const char* const unsupported[] = { + "denerr", "max_denerr_lev", "dengrad", "max_dengrad_lev", + "presserr", "max_presserr_lev", "pressgrad", "max_pressgrad_lev", + "velerr", "max_velerr_lev", "velgrad", "max_velgrad_lev", + "vorterr", "max_vorterr_lev", + "temperr", "max_temperr_lev", "tempgrad", "max_tempgrad_lev", + "ftracerr", "max_ftracerr_lev", "ftracgrad", "max_ftracgrad_lev", + "vfracerr", "max_vfracerr_lev"}; - pp.query("presserr", tagging_parm->presserr); - pp.query("max_presserr_lev", tagging_parm->max_presserr_lev); - pp.query("pressgrad", tagging_parm->pressgrad); - pp.query("max_pressgrad_lev", tagging_parm->max_pressgrad_lev); - - pp.query("velerr", tagging_parm->velerr); - pp.query("max_velerr_lev", tagging_parm->max_velerr_lev); - pp.query("velgrad", tagging_parm->velgrad); - pp.query("max_velgrad_lev", tagging_parm->max_velgrad_lev); - - pp.query("temperr", tagging_parm->temperr); - pp.query("max_temperr_lev", tagging_parm->max_temperr_lev); - pp.query("tempgrad", tagging_parm->tempgrad); - pp.query("max_tempgrad_lev", tagging_parm->max_tempgrad_lev); - - pp.query("vfracerr", tagging_parm->vfracerr); - pp.query("max_vfracerr_lev", tagging_parm->max_vfracerr_lev); + for (const char* name : unsupported) { + if (pp.contains(name)) { + amrex::Abort(std::string("CAMR: input tagging.") + name + + " is not supported -- CAMR::errorEst ignores all tagging.*" + " keys. Use amr.refinement_indicators with value_greater /" + " value_less / adjacent_difference_greater / vorticity_greater" + " (see Source/Utils/CAMR_error.cpp)."); + } + } } diff --git a/Source/main.cpp b/Source/main.cpp index 6da39e5..9fa9883 100644 --- a/Source/main.cpp +++ b/Source/main.cpp @@ -116,6 +116,20 @@ main(int argc, char* argv[]) amrptr->init(strt_time, stop_time); +#ifdef CAMR_USE_MOVING_EB + { + // The geometry is rebuilt at the start of every coarse step; the level + // factories only pick it up when every level (level 0 included) is + // regridded at the top of that step. + int force_regrid_level_zero = 0; + amrex::ParmParse ppa("amr"); + ppa.query("force_regrid_level_zero", force_regrid_level_zero); + if (amrptr->regridInt(0) != 1 || !force_regrid_level_zero) { + amrex::Abort("CAMR_USE_MOVING_EB requires amr.regrid_int = 1 and amr.force_regrid_level_zero = 1"); + } + } +#endif + // If we set the regrid_on_restart flag and if we are *not* going to take // a time step then we want to go ahead and regrid here. if ( @@ -131,12 +145,23 @@ main(int argc, char* argv[]) (amrptr->levelSteps(0) < max_step || max_step < 0) && (amrptr->cumTime() < stop_time || stop_time < 0.0)) { #ifdef CAMR_USE_MOVING_EB - initialize_EB2(amrptr->Geom(amrptr->maxLevel()), amrptr->maxLevel(), amrptr->maxLevel(), amrptr->cumTime()); + // Build the geometry for this step's start time. The existing levels + // still point into the previous index space until the regrid inside + // coarseTimeStep replaces them, so it is erased only after the step. + // The initial index space (built above at strt_time) serves the first + // step, which does not regrid. + amrex::EB2::IndexSpace const* old_is = nullptr; + if (amrptr->cumTime() > strt_time) { + old_is = amrex::EB2::TopIndexSpaceIfPresent(); + initialize_EB2(amrptr->Geom(amrptr->maxLevel()), amrptr->maxLevel(), amrptr->maxLevel(), amrptr->cumTime()); + } #endif // Do a timestep amrptr->coarseTimeStep(stop_time); #ifdef CAMR_USE_MOVING_EB - finalize_EB2(); + if (old_is) { + amrex::EB2::IndexSpace::erase(const_cast(old_is)); + } #endif }