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 }