From 3651b04a9a6a9583079629efc901b2dc3b7e7ba0 Mon Sep 17 00:00:00 2001 From: Marc Day Date: Thu, 27 Aug 2026 22:12:07 +0200 Subject: [PATCH 1/3] NSCBC 1/7: characteristic ghost-cell boundary conditions, 1-D-normal LODI limit The base characteristic treatment as a ghost-cell fill: one wave decomposition in an outward-normal frame, two boundary types. Outgoing invariants extrapolate from the interior on a minmod-limited slope; the incoming acoustic is modelled as a Poinsot-Lele relaxation written as a gradient increment (a mesh-independent rate, so literature coefficients transfer), sigma at outflows, relax_u/relax_t at inflows, with a stateless feed-forward slot (Target::dudt) for problems that inject signals. Supersonic faces are exact. Flow reversal at an outflow runs the same closure with the material slopes upwinded off across a narrow Mach band -- one code path on both sides of u_out = 0. The fill is a pure function of pre-launch interior data: idempotent, order-independent, bit-reproducible across restarts and decompositions, GPU-resident (POD captures, counted fallbacks instead of device aborts). Every fallback increments a counter and the counters report by default -- a silent fallback is indistinguishable from a healthy boundary. Problems opt in per boundary POINT through a bcnormal_nscbc() hook returning a Target; a face may mix inflow, outflow and wall. EB requires eb_zero_body_state (fatal otherwise, the fill detects covered cells by non-positive density); an AMR fine face touching a characteristic boundary warns. This PR is the 1-D-normal LODI limit on purpose: transverse terms and the reacting-boundary closures arrive as the next two PRs in this stack, each an additive, default-off correction to the modelled incoming wave (some comments here forward-reference them). Verification driver and regression cases follow in the stack; headline number from the case PR: an acoustic pulse leaves through this boundary with 0.8% reflection against 97% for a hard-pressure fill. Also removes the dead nscbc_adv/nscbc_diff keys of the deleted Fortran implementation (aborting, not ignoring, when an inputs file still sets them) and retires the 1 K default for eb_boundary_T in favour of 300 K. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01XU8M23nucKsu1do2WxXFeq --- .codespell-ignore-words | 4 + .gitignore | 1 + Exec/RegTests/MMS/mms-4.inp | 9 +- Source/BCfill.cpp | 325 +++++++- Source/Diffusion.cpp | 50 +- Source/Hydro.cpp | 56 +- Source/NSCBC.H | 738 ++++++++++++++++++ Source/Params/_cpp_parameters | 63 +- Source/Params/param_includes/pelec_defaults.H | 10 +- Source/Params/param_includes/pelec_params.H | 8 +- Source/Params/param_includes/pelec_queries.H | 8 +- Source/PeleC.H | 27 + Source/PeleC.cpp | 119 +++ Source/ProblemSpecificFunctions.H | 32 + Source/SumIQ.cpp | 2 + 15 files changed, 1356 insertions(+), 96 deletions(-) create mode 100644 Source/NSCBC.H diff --git a/.codespell-ignore-words b/.codespell-ignore-words index 49d7e99a6..c90a7c1a9 100644 --- a/.codespell-ignore-words +++ b/.codespell-ignore-words @@ -4,3 +4,7 @@ blocs renewl frop fpr +ue +statics +Rin +anly diff --git a/.gitignore b/.gitignore index 9dea5cb1e..fbbf487af 100644 --- a/.gitignore +++ b/.gitignore @@ -38,3 +38,4 @@ Docs/sphinx/doctrees/ Docs/sphinx_doc/ __pycache__ +build*/ diff --git a/Exec/RegTests/MMS/mms-4.inp b/Exec/RegTests/MMS/mms-4.inp index 8982912ae..ec02f873a 100644 --- a/Exec/RegTests/MMS/mms-4.inp +++ b/Exec/RegTests/MMS/mms-4.inp @@ -16,9 +16,12 @@ amr.n_cell = 16 16 16 pelec.lo_bc = "FOExtrap" "NoSlipWall" "Hard" pelec.hi_bc = "Hard" "NoSlipWall" "FOExtrap" -# We should not have NSCBC turned on for this test case -pelec.nscbc_adv = 0 -pelec.nscbc_diff = 0 +# NOTE: this case must keep its "Hard" faces as exact Dirichlet faces -- the +# manufactured solution is imposed there and any characteristic relaxation of +# the boundary state would invalidate the convergence test. It therefore must +# NOT enable the characteristic (NSCBC) boundary treatment. (This used to be +# spelled pelec.nscbc_adv = 0 / pelec.nscbc_diff = 0, which had no effect; both +# keys have been removed.) # WHICH PHYSICS pelec.do_hydro = 1 diff --git a/Source/BCfill.cpp b/Source/BCfill.cpp index 6d08a815c..51b905345 100644 --- a/Source/BCfill.cpp +++ b/Source/BCfill.cpp @@ -1,22 +1,205 @@ +#include + #include #include +#include #include +#include #include "PeleC.H" #include "prob.H" +#include "NSCBC.H" + +namespace { +// Device-resident fallback counters for the characteristic boundary +// treatment. Allocated on first use, reduced and reported by +// PeleC::nscbc_report_diagnostics(). +// +// Held as a heap pointer released through amrex::ExecOnFinalize rather than as +// a function-local static object. A static Gpu::DeviceVector destructs at +// program exit, which is AFTER amrex::Finalize() has torn down the arena it +// allocated from -- harmless on a CPU build and a use-after-free of the device +// allocator on a GPU one. +amrex::Gpu::DeviceVector* nscbc_diag_p = nullptr; + +amrex::Gpu::DeviceVector& +nscbc_diag() +{ + if (nscbc_diag_p == nullptr) { + nscbc_diag_p = + new amrex::Gpu::DeviceVector(pc_nscbc::Diag::count, 0); + amrex::ExecOnFinalize([]() { + delete nscbc_diag_p; + nscbc_diag_p = nullptr; + }); + } + return *nscbc_diag_p; +} + +} // namespace struct PCHypFillExtDir { ProbParmDevice const* lprobparm; bool m_do_turb_inflow{false}; + bool m_nscbc{false}; + amrex::GpuArray m_nscbc_prm; + amrex::Long* m_nscbc_diag{nullptr}; AMREX_GPU_HOST - constexpr explicit PCHypFillExtDir( - const ProbParmDevice* d_prob_parm, const bool do_turb_inflow) - : lprobparm(d_prob_parm), m_do_turb_inflow(do_turb_inflow) + explicit PCHypFillExtDir( + const ProbParmDevice* d_prob_parm, + const bool do_turb_inflow, + const bool nscbc, + const amrex::GpuArray& nscbc_prm, + amrex::Long* nscbc_diag) + : lprobparm(d_prob_parm), + m_do_turb_inflow(do_turb_inflow), + m_nscbc(nscbc), + m_nscbc_prm(nscbc_prm), + m_nscbc_diag(nscbc_diag) { } + // ------------------------------------------------------------------------- + // Characteristic (NSCBC) fill for one ghost cell. + // + // Returns true if this ghost cell was filled here, in which case the + // ordinary bcnormal() path below is skipped for it entirely. + // + // Corner ownership. A ghost cell may lie outside the domain in more than + // one direction. Such a cell is owned by the LOWEST idir in which it is + // outside on an ext_dir face whose problem hook returns a live target; + // the state is written exactly once. Combined with the clamped stencil + // below this makes the fill a pure function of valid interior data, hence + // independent of the order in which ghost cells are visited and identical + // on CPU and GPU. + // + // (Note that the legacy bcnormal() path further down does NOT have this + // property: at a corner it reads dest() at a tangential index that is + // itself a ghost cell, which another thread in the same launch may be + // writing. That is pre-existing behaviour and is left untouched here so + // that no existing result moves.) + // ------------------------------------------------------------------------- + AMREX_GPU_DEVICE + AMREX_FORCE_INLINE + bool nscbc_fill( + const amrex::IntVect& iv, + amrex::Array4 const& dest, + amrex::GeometryData const& geom, + const amrex::Real time, + const int* bc) const + { + const int* domlo = geom.Domain().loVect(); + const int* domhi = geom.Domain().hiVect(); + const amrex::Real* prob_lo = geom.ProbLo(); + const amrex::Real* dx = geom.CellSize(); + + for (int idir = 0; idir < AMREX_SPACEDIM; ++idir) { + int sgn = 0; + if ((bc[idir] == amrex::BCType::ext_dir) && (iv[idir] < domlo[idir])) { + sgn = +1; + } else if ( + (bc[idir + AMREX_SPACEDIM] == amrex::BCType::ext_dir) && + (iv[idir] > domhi[idir])) { + sgn = -1; + } else { + continue; + } + + const int N_pos = (sgn > 0) ? domlo[idir] : domhi[idir]; + const int layer = sgn * (N_pos - iv[idir]); // 1 = nearest the domain + + const amrex::Dim3 lo3 = amrex::lbound(dest); + const amrex::Dim3 hi3 = amrex::ubound(dest); + const int fab_lo[3] = {lo3.x, lo3.y, lo3.z}; + const int fab_hi[3] = {hi3.x, hi3.y, hi3.z}; + + // Tangential indices are clamped into the domain (and the FAB), in + // every tangential direction: only valid interior cells are read, and + // the fill is a pure function of pre-launch data. + // + // In a PERIODIC tangential direction this leaves a small, measured + // seam residual, and that is a deliberate trade. amrex's corner + // protocol (StateDataPhysBCFunct) recomputes corner ghosts on a strip + // FAB holding only their image band, so a seam-adjacent cell is filled + // twice on two different FABs; exact agreement requires restricting + // the tangential stencil to the band the strip can see. Both + // alternatives were built and measured: a wrap through the resident + // images leaves the array aperiodic at the 2e-2 level, and the + // strip-consistent band is bitwise-periodic and measured equally + // stable (NSCBC-FlameOutflow-DRM, beta = 0.5) -- but costs ~100 lines + // of decomposition-sensitive index machinery to remove a residual that + // measures 2e-4 (inert) to 1.6e-3 (flame on the seam corner). The + // clamp was kept for simplicity; nscbc_check_periodic_wrap() reports + // the residual and aborts above 1e-2, which a broken stencil fails. + amrex::IntVect base(AMREX_D_DECL(iv[0], iv[1], iv[2])); + for (int d = 0; d < AMREX_SPACEDIM; ++d) { + if (d != idir) { + base[d] = amrex::min( + amrex::max(iv[d], amrex::max(domlo[d], fab_lo[d])), + amrex::min(domhi[d], fab_hi[d])); + } + } + + // How deep a normal stencil is available, in the domain and in the FAB. + const int depth_domain = domhi[idir] - domlo[idir] + 1; + const int depth_fab = + (sgn > 0) ? (fab_hi[idir] - N_pos + 1) : (N_pos - fab_lo[idir] + 1); + const int n_stencil = + amrex::min(3, amrex::min(depth_domain, depth_fab)); + if (n_stencil < 1) { + continue; + } + + auto stencil_iv = [&](const int step) { + amrex::IntVect r = base; + r[idir] = N_pos + sgn * amrex::min(step, n_stencil - 1); + return r; + }; + amrex::Real s_N[NVAR], s_Nm1[NVAR], s_Nm2[NVAR]; + const amrex::IntVect ivN = stencil_iv(0); + const amrex::IntVect ivNm1 = stencil_iv(1); + const amrex::IntVect ivNm2 = stencil_iv(2); + for (int n = 0; n < NVAR; n++) { + s_N[n] = dest(ivN, n); + s_Nm1[n] = dest(ivNm1, n); + s_Nm2[n] = dest(ivNm2, n); + } + + // Query the problem for this boundary POINT. x is the location on the + // boundary plane, not the ghost cell centre: the target is a property + // of the boundary point and must not vary with the ghost layer, or the + // relaxation would be applied to a different target in each layer. + amrex::Real x[AMREX_SPACEDIM]; + for (int d = 0; d < AMREX_SPACEDIM; ++d) { + x[d] = prob_lo[d] + (static_cast(base[d]) + 0.5) * dx[d]; + } + x[idir] = prob_lo[idir] + static_cast( + (sgn > 0) ? domlo[idir] : (domhi[idir] + 1)) * + dx[idir]; + + pc_nscbc::Target tgt = ProblemSpecificFunctions::bcnormal_nscbc( + x, s_N, idir, sgn, time, geom, *lprobparm); + if (tgt.type == pc_nscbc::Type::off) { + continue; // this face is not characteristic here; try the next + } + + // Boundary-register composition (level 0 only; the registers change + // once per advance and are read frozen here). The kernel below sees + // only the composed Target -- it stays a pure function. + amrex::Real s_ghost[NVAR]; + pc_nscbc::apply( + s_N, s_Nm1, s_Nm2, n_stencil, dx[idir], idir, sgn, layer, tgt, + m_nscbc_prm[idir], s_ghost, m_nscbc_diag); + for (int n = 0; n < NVAR; n++) { + dest(iv, n) = s_ghost[n]; + } + return true; + } + return false; + } + AMREX_GPU_DEVICE void operator()( const amrex::IntVect& iv, @@ -43,6 +226,11 @@ struct PCHypFillExtDir const int* bc = bcr->data(); + // Characteristic boundary treatment, where the problem asks for it. + if (m_nscbc && nscbc_fill(iv, dest, geom, time, bc)) { + return; + } + amrex::Real s_int[NVAR] = {0.0}; amrex::Real s_ext[NVAR] = {0.0}; amrex::GpuArray turb_fluc{0.0}; @@ -196,8 +384,19 @@ pc_bcfill_hyp( } const ProbParmDevice* lprobparm = PeleC::d_prob_parm_device; + + // Capture the NSCBC parameters HOST-side (a device kernel must + // never touch ParmParse or a class static). + const bool nscbc = PeleC::nscbc_active(); + amrex::GpuArray nscbc_prm; + for (int d = 0; d < AMREX_SPACEDIM; ++d) { + nscbc_prm[d] = PeleC::nscbc_params(d); + } + amrex::Long* diag = nscbc ? nscbc_diag().data() : nullptr; + amrex::GpuBndryFuncFab hyp_bndry_func( - PCHypFillExtDir{lprobparm, PeleC::turb_inflow.is_initialized()}); + PCHypFillExtDir{ + lprobparm, PeleC::turb_inflow.is_initialized(), nscbc, nscbc_prm, diag}); hyp_bndry_func(bx, data, dcomp, numcomp, geom, time, bcr, bcomp, scomp); } @@ -231,3 +430,121 @@ pc_nullfill( const int /*scomp*/) { } + +void +PeleC::nscbc_check_fine_faces() const +{ + if (!bc_nscbc || level == 0) { + return; + } + const amrex::Box& dom = geom.Domain(); + for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { + for (int side = 0; side < 2; ++side) { + const int t = (side == 0) ? phys_bc.lo(dir) : phys_bc.hi(dir); + if ((t != PCPhysBCType::inflow) && (t != PCPhysBCType::user_bc)) { + continue; // not a face the characteristic treatment can reach + } + const int face = (side == 0) ? dom.smallEnd(dir) : dom.bigEnd(dir); + bool touches = false; + for (int i = 0; i < grids.size(); ++i) { + const amrex::Box& b = grids[i]; + if ( + ((side == 0) && (b.smallEnd(dir) == face)) || + ((side == 1) && (b.bigEnd(dir) == face))) { + touches = true; + break; + } + } + // Once per face per run: regridding can recur every few steps, and a + // warning repeated 1600 times is a warning nobody reads. + static bool warned[AMREX_SPACEDIM][2] = {}; + if ( + touches && !warned[dir][side] && + amrex::ParallelDescriptor::IOProcessor()) { + warned[dir][side] = true; + // A warning rather than an abort: the problem hook decides per + // boundary POINT whether a face is characteristic, and the host + // cannot know what it will return. But if any point of this face + // is characteristic, the fill's extrapolation stencil is + // level-local, so the fine patch imposes a DIFFERENT boundary + // condition than the coarse level does on the same face -- a + // level-dependent artefact that refining cannot remove. + amrex::Warning( + "NSCBC: level " + std::to_string(level) + " grids touch the domain " + + (side == 0 ? "lo" : "hi") + " face in direction " + + std::to_string(dir) + + ", which is a Hard/UserBC face with pelec.bc_nscbc = 1. The " + "characteristic fill's stencil is level-local, so a refined patch " + "on a characteristic face makes the boundary condition " + "level-dependent. Keep refinement away from characteristic faces " + "(see the BCs chapter). (This warning is printed once per face.)"); + } + } + } +} + +pc_nscbc::Params +PeleC::nscbc_params(const int idir) +{ + pc_nscbc::Params p; + p.sigma = bc_nscbc_sigma; + p.relax_u = bc_nscbc_relax_u; + p.relax_t = bc_nscbc_relax_t; + p.order = bc_nscbc_order; + p.pin_farfield = bc_nscbc_pin_farfield; + // Only the ratio sigma/L_ref is physical. L_ref is fixed to the domain + // extent along the boundary normal so that sigma keeps the meaning it has + // in the literature, rather than being one of two dials for one degree of + // freedom. Uses probhi - problo, not probhi: the legacy Fortran used + // probhi(idir) and was therefore silently wrong for any domain not + // anchored at the origin. + const auto& geom = amrex::DefaultGeometry(); + p.L_ref = geom.ProbHi(idir) - geom.ProbLo(idir); + return p; +} + +void +PeleC::nscbc_report_diagnostics() +{ + if (!bc_nscbc) { + return; + } + std::vector h(pc_nscbc::Diag::count, 0); + amrex::Gpu::copy( + amrex::Gpu::deviceToHost, nscbc_diag().begin(), nscbc_diag().end(), + h.begin()); + amrex::ParallelDescriptor::ReduceLongSum(h.data(), pc_nscbc::Diag::count); + + // The supersonic path is exact, not a degradation, so it is reported but is + // not a warning. The others mean the boundary is being asked for something + // it cannot cleanly provide. + // transverse_drop and source_drop belong here as much as the rest: a + // beta or beta_s that is silently not being applied looks exactly like a + // beta or beta_s that does nothing, and the only way to tell the two apart + // is to count it. + const amrex::Long total = + h[pc_nscbc::Diag::reversed] + h[pc_nscbc::Diag::body_state] + + h[pc_nscbc::Diag::eos_failure] + h[pc_nscbc::Diag::floored] + + h[pc_nscbc::Diag::transverse_drop] + h[pc_nscbc::Diag::source_drop] + + h[pc_nscbc::Diag::target_incomplete]; + if (amrex::ParallelDescriptor::IOProcessor() && (total > 0 || verbose > 1)) { + amrex::Print() << " NSCBC fallbacks since last report:" << " supersonic " + << h[pc_nscbc::Diag::supersonic] << ", flow reversal " + << h[pc_nscbc::Diag::reversed] << ", EB body state " + << h[pc_nscbc::Diag::body_state] << ", EOS failure " + << h[pc_nscbc::Diag::eos_failure] << ", floored " + << h[pc_nscbc::Diag::floored] << ", transverse dropped " + << h[pc_nscbc::Diag::transverse_drop] << ", source dropped " + << h[pc_nscbc::Diag::source_drop] << ", target incomplete " + << h[pc_nscbc::Diag::target_incomplete] << "\n"; + } + // Settle any counter atomics still in flight on other streams before the + // reset; the blocking Gpu::copy above synchronised only its own stream. + amrex::Gpu::Device::synchronize(); + nscbc_diag().assign(pc_nscbc::Diag::count, 0); + // The reset itself is a device fill on the current stream; the next + // advance's fills bump these counters from MFIter's rotating streams, + // which are not ordered against it. Settle the reset before returning + // so a zero can never land on top of a fresh count. + amrex::Gpu::streamSynchronize(); +} diff --git a/Source/Diffusion.cpp b/Source/Diffusion.cpp index 624042039..aa9665c83 100644 --- a/Source/Diffusion.cpp +++ b/Source/Diffusion.cpp @@ -143,39 +143,23 @@ PeleC::getMOLSrcTerm( } }); } - // TODO deal with NSCBC - /* - for (int dir = 0; dir < AMREX_SPACEDIM ; dir++) { - const amrex::Box& bxtmp = amrex::surroundingNodes(vbox,dir); - amrex::Box TestBox(bxtmp); - for(int d=0; d 0. The +// outward wave speeds are u_out - c (always incoming subsonically), u_out +// (the lambda_0 family: entropy, tangential velocity, species, passives -- +// outgoing at an outflow, incoming at an inflow) and u_out + c (always +// outgoing). Hence a subsonic outflow specifies one quantity (p); a +// subsonic inflow specifies u, T and Y. +// +// We work with linearised invariants, the impedance rho*c FROZEN at the +// boundary cell N (a per-cell impedance would make R differences of two +// different variables): +// +// R_+/- = u_out +/- p/(rho c), S = rho - p/c^2 . +// +// Outgoing quantities extrapolate outward on a minmod-limited slope; the +// incoming one takes the Poinsot-Lele relaxation as a gradient increment, +// +// L_in = K (phi - phi_target), K = coeff (1 - M^2) c / L_ref [1/s] +// R_-,ghost = R_-,N + layer dx L_in / ((c - u_out) rho c) , +// +// so K is a mesh-independent RATE and literature coefficients transfer +// directly. The gradient form gives L_in an additive slot for the +// transverse terms (beta) and the reaction source (beta_s); a hard pin of +// R_- has none. The fill is a pure algebraic function of the current +// interior state: idempotent under repeated FillPatch calls, +// bit-reproducible across restarts and decompositions. +// +// SIGNS. Every user-facing coefficient (sigma, relax_u, relax_t) is +// positive; all sign handling is internal and Params::validate() rejects +// negatives. +// +// GPU. Finiteness guards use amrex::Math::isfinite, never std::isfinite +// (SYCL fast-math makes the latter constant-false, which would silently +// degrade every fill to a zero-gradient copy). No allocation, no virtual +// dispatch, no printf/Abort on the device path; fallbacks are counted via +// atomics (see Diag). Stack footprint is order 10*NVAR doubles per thread +// and will spill past the register budget for large mechanisms; tolerable +// because the launch covers only a thin ghost region. +// +// EOS. Only density-carrying entry points are used (RTY2P, RTY2Cs, RYP2T, +// RTY2E, REY2T, PYT2RE), which keeps the kernel real-gas capable: the +// shortcut forms error out under SRK. +// ============================================================================ + +#include +#include +#include +#include +#include + +#include "IndexDefines.H" +#include "Constants.H" +#include "PelePhysics.H" + +namespace pc_nscbc { + +// --------------------------------------------------------------------------- +// Boundary type, decided per boundary POINT by the problem's bcnormal_nscbc(). +// A single face may carry a mixture (a jet inlet surrounded by outflow). +// --------------------------------------------------------------------------- +enum struct Type : int { off = 0, outflow = 1, inflow = 2 }; + +// --------------------------------------------------------------------------- +// Numerical parameters. POD, captured host-side into the boundary functor; +// never read ParmParse or a class static from inside a device kernel. +// --------------------------------------------------------------------------- +struct Params +{ + // Outflow pressure relaxation, classical Poinsot-Lele sigma: + // K = sigma * (1 - M^2) * c / L_ref [1/s] + // 0 is perfectly non-reflecting and leaves the mean pressure unanchored; + // the literature band is 0.15-0.3 (Rudy & Strikwerda 1980 give ~0.27 as + // optimal for a 1-D pulse). + amrex::Real sigma = 0.25; + + // Inflow normal-velocity relaxation. Larger imposes the target velocity + // more strongly and reflects more; >~10 is a hard Dirichlet in disguise. + amrex::Real relax_u = 2.0; + + // Inflow temperature and tangential-velocity relaxation. + amrex::Real relax_t = 0.2; + + // Boundary-normal reference length, probhi[idir] - problo[idir]. Filled per + // face by the caller. Only the ratio sigma/L_ref is physical, so this is + // not a user knob; it exists so that sigma keeps its literature meaning. + amrex::Real L_ref = 1.0; + + // Order of the outgoing-invariant extrapolation: + // 1 = zeroth-order copy (constant across ghost layers) + // 2 = minmod-limited linear in the ghost layer index (default) + // Exposed so that a bit-comparison against the standalone 1-D driver can be + // done without recompiling. Verification knob, not a physics knob. + int order = 2; + + // Replace the relaxation by a hard pin of R_- to the far-field value + // (u_out = 0, p = p_target): non-reflecting AND anchored, but it anchors to + // p_target + rho c u_out, and as a value constraint its effective rate is + // c/dx -- it does not converge under refinement. For a boundary onto a + // quiescent reservoir; wrong for a through-flow duct. sigma is ignored. + bool pin_farfield = false; + + AMREX_GPU_HOST bool validate(std::string& why) const + { + if (sigma < 0.0) { + why = "bc_nscbc_sigma must be >= 0 (all NSCBC relaxation coefficients " + "are positive in PeleC; internal signs are handled by the kernel)"; + return false; + } + if (relax_u < 0.0) { + why = "bc_nscbc_relax_u must be >= 0"; + return false; + } + if (relax_t < 0.0) { + why = "bc_nscbc_relax_t must be >= 0"; + return false; + } + if ((order != 1) && (order != 2)) { + why = "bc_nscbc_order must be 1 or 2"; + return false; + } + if (L_ref <= 0.0) { + why = "NSCBC reference length must be > 0"; + return false; + } + return true; + } +}; + +// --------------------------------------------------------------------------- +// Target state returned per boundary point by the problem hook. Only the +// entries relevant to `type` are read: outflow -> p; inflow -> u[], T, Y[]. +// Deliberately no relaxation for composition: at an inflow every species +// characteristic is incoming (hard imposition is well posed), at an outflow +// all are outgoing (imposing would over-specify). +// --------------------------------------------------------------------------- +struct Target +{ + Type type = Type::off; + amrex::Real p = 0.0; // outflow: target pressure [dyn/cm^2] + amrex::Real u[3] = {0.0, 0.0, 0.0}; // inflow: target velocity [cm/s] + amrex::Real T = 0.0; // inflow: target temperature [K] + amrex::Real Y[NUM_SPECIES] = {0.0}; // inflow: target mass fractions + + // Inflow acoustic-injection feed-forward: the time derivative of the + // NORMAL velocity target, supplied analytically by the problem hook when + // it is forcing the inlet. The value-relaxation alone has no amplitude + // slot (design doc I.3e) -- an injected signal rides the same relaxation + // that holds the mean, with the measured deterioration of driver t3 -- + // and a register-based NRI transplant does not create one (measured: + // subtracting the outgoing wave from the deviation only exposes the bare + // injection bandwidth limit). This term is that slot, stateless: the + // incoming amplitude gains 2 rho c du^t/dt, the GC transplant of the + // flux form's acoustic forcing, and the relaxation keeps its one job of + // holding the mean. Zero (default) recovers the classical inlet + // identically. + amrex::Real dudt = 0.0; // d(u_normal target)/dt, LAB frame [cm/s^2] +}; + +// --------------------------------------------------------------------------- +// Diagnostic counters. Every fallback is counted; a silent fallback is a +// bug that will not be found. Storage is amrex::Long, not int: with the +// default sum_interval = -1 the counters are never reset, and the supersonic +// counter on a large 3-D run walks past 2^31 in ~1e4 steps. +// --------------------------------------------------------------------------- +namespace Diag { +enum : int { + supersonic = 0, // |M| >= 1: zero-gradient copy is exact, not a failure + reversed, // flow ran the wrong way through the face + body_state, // an EB-covered cell was found in the stencil + eos_failure, // an EOS call returned a non-finite or non-positive value + floored, // a positivity floor was hit + transverse_drop, // a transverse derivative was unavailable or unusable + source_drop, // a modelled source term was unavailable or unusable + target_incomplete, // a supersonic inflow whose Target carries no pressure: + // every characteristic is incoming, so the FULL state is + // required; the interior pressure was substituted, which + // makes the ghost depend on a domain it should be + // causally upstream of. Set Target.p in bcnormal_nscbc. + structure, // ADVISORY, not a fallback: material structure sat in an + // outflow boundary cell (|dS| > 5% of rho). A front at + // the boundary wants the flame closures -- + // extrap_temperature and beta_s = 0; see + // NSCBC-FlameOutflow/README.md. + count +}; +} + +// =========================================================================== +// Small helpers +// =========================================================================== + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real +finite_or(const amrex::Real x, const amrex::Real fallback) noexcept +{ + return amrex::Math::isfinite(x) ? x : fallback; +} + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real +minmod(const amrex::Real a, const amrex::Real b) noexcept +{ + if (a * b <= 0.0) { + return 0.0; + } + return (std::abs(a) < std::abs(b)) ? a : b; +} + +// True if a state looks like EB body state rather than fluid, detected by +// non-positive density. This is why read_params makes eb_zero_body_state = 1 +// mandatory under EB: PeleC's default body state is a sampled fluid state the +// stencil cannot distinguish from interior data. +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool +is_body_state(const amrex::Real s[NVAR]) noexcept +{ + return !(s[URHO] > constants::very_small_num()) || + !amrex::Math::isfinite(s[URHO]) || !amrex::Math::isfinite(s[UEDEN]); +} + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void +copy_state(const amrex::Real src[NVAR], amrex::Real dst[NVAR]) noexcept +{ + for (int n = 0; n < NVAR; n++) { + dst[n] = src[n]; + } +} + +// --------------------------------------------------------------------------- +// Primitive decomposition of one stencil cell. T is taken from the state +// (a maintained component; computeTemp runs after every update), keeping +// every EOS call on this path non-iterative -- which matters under SRK -- +// with the iterative recovery as fallback. +// --------------------------------------------------------------------------- +struct CellPrim +{ + amrex::Real rho, u[3], T, p, c, e; + amrex::Real Y[NUM_SPECIES]; + bool ok; +}; + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE CellPrim +cell_primitives(const amrex::Real s[NVAR]) noexcept +{ + auto eos = pele::physics::PhysicsType::eos(); + CellPrim q; + q.ok = true; + + q.rho = s[URHO]; + const amrex::Real rhoinv = 1.0 / q.rho; + q.u[0] = s[UMX] * rhoinv; + q.u[1] = s[UMY] * rhoinv; + q.u[2] = s[UMZ] * rhoinv; + + amrex::Real Ysum = 0.0; + for (int n = 0; n < NUM_SPECIES; n++) { + q.Y[n] = amrex::max(s[UFS + n] * rhoinv, 0.0); + Ysum += q.Y[n]; + } + if (Ysum > constants::very_small_num()) { + const amrex::Real inv = 1.0 / Ysum; + for (amrex::Real& yv : q.Y) { + yv *= inv; + } + } else { + q.ok = false; + } + + const amrex::Real ke = + 0.5 * (q.u[0] * q.u[0] + q.u[1] * q.u[1] + q.u[2] * q.u[2]); + q.e = s[UEDEN] * rhoinv - ke; + + q.T = s[UTEMP]; + if (!(q.T > 0.0) || !amrex::Math::isfinite(q.T)) { + q.T = 300.0; // initial guess for the iterative fallback + eos.REY2T(q.rho, q.e, q.Y, q.T); + } + + eos.RTY2P(q.rho, q.T, q.Y, q.p); + eos.RTY2Cs(q.rho, q.T, q.Y, q.c); + + if ( + !amrex::Math::isfinite(q.p) || !(q.p > 0.0) || + !amrex::Math::isfinite(q.c) || !(q.c > 0.0) || + !amrex::Math::isfinite(q.T) || !(q.T > 0.0)) { + q.ok = false; + } + return q; +} + +// --------------------------------------------------------------------------- +// Pack a set of ghost primitives into a conserved state. +// +// Every NVAR slot is written, and the algebraic identities that the rest of +// the code assumes are enforced exactly rather than approximately: +// sum_k (rho Y)_k == rho (species renormalised before packing) +// UEDEN == UEINT + 0.5 rho |u|^2 +// The passive families (advected, auxiliary, linear/soot) ride the lambda_0 +// characteristic and are handled by the caller through `pass`. +// --------------------------------------------------------------------------- +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void +pack_ghost( + const amrex::Real rho, + const amrex::Real u[3], + const amrex::Real e, + const amrex::Real T, + const amrex::Real Y[NUM_SPECIES], + const amrex::Real s_N[NVAR], + amrex::Real s_ghost[NVAR]) noexcept +{ + // s_N feeds only the passive families below; in a build with no passive + // fields every #if block is empty and the parameter would warn as unused. + amrex::ignore_unused(s_N); + s_ghost[URHO] = rho; + s_ghost[UMX] = rho * u[0]; + s_ghost[UMY] = rho * u[1]; + s_ghost[UMZ] = rho * u[2]; + s_ghost[UEINT] = rho * e; + s_ghost[UEDEN] = + rho * e + 0.5 * rho * (u[0] * u[0] + u[1] * u[1] + u[2] * u[2]); + s_ghost[UTEMP] = T; + for (int n = 0; n < NUM_SPECIES; n++) { + s_ghost[UFS + n] = rho * Y[n]; + } + // Passive families default to a mass-scaled zero-gradient; the caller + // overwrites them where a target applies. +#if NUM_ADV > 0 + const amrex::Real ratio = + (s_N[URHO] > constants::very_small_num()) ? rho / s_N[URHO] : 1.0; + for (int n = 0; n < NUM_ADV; n++) { + s_ghost[UFA + n] = s_N[UFA + n] * ratio; // per unit mass in ctoprim + } +#endif +#if NUM_AUX > 0 + for (int n = 0; n < NUM_AUX; n++) { + s_ghost[UFX + n] = s_N[UFX + n]; // NOT per unit mass in ctoprim + } +#endif +#if NUM_LIN > 0 + for (int n = 0; n < NUM_LIN; n++) { + s_ghost[ULIN + n] = s_N[ULIN + n]; // NOT per unit mass in ctoprim + } +#endif +} + +// =========================================================================== +// apply() -- fill one ghost cell. +// s_N, s_Nm1, s_Nm2: interior stencil walking inward; tangential indexing +// is the caller's (BCfill.cpp), keeping the fill a deterministic pure +// function of data written before the launch. +// n_stencil: usable cells (3/2/1), reduced further on EB body state. +// idir, sgn: boundary normal and side (+1 lo, -1 hi, as bcnormal). +// layer: 1 = nearest the domain. diag: optional Long[Diag::count]. +// =========================================================================== +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void +apply( + const amrex::Real s_N[NVAR], + const amrex::Real s_Nm1[NVAR], + const amrex::Real s_Nm2[NVAR], + int n_stencil, + const amrex::Real dx_normal, + const int idir, + const int sgn, + const int layer, + Target const& tgt, + Params const& prm, + amrex::Real s_ghost[NVAR], + amrex::Long* diag = nullptr) noexcept +{ + auto eos = pele::physics::PhysicsType::eos(); + + auto bump = [=](const int which) { + if (diag != nullptr) { + amrex::HostDevice::Atomic::Add(&diag[which], amrex::Long(1)); + } + }; + + // ---- 0. Stencil hygiene ------------------------------------------------ + if (is_body_state(s_N)) { + // Nothing usable to build a characteristic state from. Leave whatever the + // caller pre-populated (the reflected NoSlipWall state) in place. + bump(Diag::body_state); + copy_state(s_N, s_ghost); + return; + } + if ((n_stencil > 2) && is_body_state(s_Nm2)) { + n_stencil = 2; + bump(Diag::body_state); + } + if ((n_stencil > 1) && is_body_state(s_Nm1)) { + n_stencil = 1; + bump(Diag::body_state); + } + // Extrapolation order can never exceed what the stencil supports. A + // minmod-limited slope needs three cells; with two it would be unlimited. + const int order = + amrex::min(prm.order, amrex::max(n_stencil - 1, 1)); + const bool linear = (order >= 2) && (n_stencil >= 3); + + const CellPrim qN = cell_primitives(s_N); + if (!qN.ok) { + bump(Diag::eos_failure); + copy_state(s_N, s_ghost); + return; + } + + // ---- 1. Outward-normal frame ------------------------------------------ + const auto n_sgn = static_cast(-sgn); + const amrex::Real u_out_N = n_sgn * qN.u[idir]; + const amrex::Real c_N = qN.c; + const amrex::Real mach = u_out_N / c_N; + + // Frozen acoustic impedance (see the header note: this must be a single + // value shared by every stencil cell). + const amrex::Real rho_c = qN.rho * c_N; + + // ---- 2. Regimes with an exact answer ---------------------------------- + if (std::abs(mach) >= 1.0) { + if (mach >= 1.0) { + // Supersonic outflow: zero-gradient is exact, not an approximation. + bump(Diag::supersonic); + copy_state(s_N, s_ghost); + return; + } + // Supersonic inflow: every characteristic is incoming, so the full state + // is specified. Only meaningful with an inflow target. A Target + // without a pressure under-specifies it: the interior pressure is + // substituted and COUNTED (target_incomplete), because the ghost then + // depends on a domain it should be causally upstream of. + if (tgt.type == Type::inflow) { + if (!(tgt.p > 0.0)) { + bump(Diag::target_incomplete); + } + amrex::Real rho_g = 0.0, e_g = 0.0; + eos.PYT2RE(tgt.p > 0.0 ? tgt.p : qN.p, tgt.Y, tgt.T, rho_g, e_g); + if ( + amrex::Math::isfinite(rho_g) && (rho_g > 0.0) && + amrex::Math::isfinite(e_g)) { + amrex::Real u_g[3] = {tgt.u[0], tgt.u[1], tgt.u[2]}; + pack_ghost(rho_g, u_g, e_g, tgt.T, tgt.Y, s_N, s_ghost); + bump(Diag::supersonic); + return; + } + bump(Diag::eos_failure); + } + bump(Diag::supersonic); + copy_state(s_N, s_ghost); + return; + } + + // ---- 3. Direction sanity ---------------------------------------------- + const bool outflow_face = (tgt.type == Type::outflow); + if (outflow_face && (u_out_N < 0.0)) { + // Transient reversal at an outflow: COUNTED, but no longer a separate + // closure. The acoustic side of the standard outflow path is well + // defined and restoring for u_out < 0 (lambda_in = c - u_out grows, + // 1 - M^2 is even in M, K (p - p_tgt) pulls the ghost pressure toward + // the target AND pushes the ghost normal velocity outward), so one + // closure runs on both sides of u_out = 0; only the lambda_0 material + // SLOPES are upwinded off during reversal (see w_mat below), leaving the + // ghost the interior's material VALUES. Both dedicated branches tried + // here were defects the chamber found: a hard pin of tgt.p is a step + // discontinuity that NaN'd under vent breathing (backflow -144 -> -2628 + // cm/s in 45 steps, a 138 K cell), and the soft pressure-only relaxation + // that replaced it dropped dR+, S_p and T_in relative to the forward + // branch and froze the ghost velocity -- an O(1) discontinuity at + // u_out = 0 that a flame finishing its transit turned into a growing + // dither and a spurious 0.3 atm chamber spike + // (Docs/NSCBC-reversal-branch-defect.md; gates: driver C13 a/b/c). + // Sustained recirculation as a local inflow is Params::backflow_material. + // + // The COUNTER carries a roundoff deadband (the closure itself has no + // branch left to care): a face holding a quiescent charge sits at + // exactly u_out = 0 and the solver's noise dithers the sign, which read + // as 314k "reversals" in the first 100 steps of a chamber whose ignition + // wave was still a centimetre from the boundary. Counting only + // |u_out| > 1e-9 c keeps "reversal count zero" meaning NO PHYSICAL + // BACKFLOW rather than "no roundoff", which is the semantics the + // README's advice ("sustained nonzero means a misplaced outflow") + // depends on. + if (u_out_N < -1.0e-9 * c_N) { + bump(Diag::reversed); + } + } + if (!outflow_face && (u_out_N > 0.0)) { + // Outflow through a configured inflow would impose D+NUM_SPECIES+1 + // conditions where one is admitted; zero-gradient and count it. + bump(Diag::reversed); + copy_state(s_N, s_ghost); + return; + } + + // ---- 4. Outgoing acoustic invariant, R_+ ------------------------------ + // R_+ leaves the domain at both subsonic inflow and subsonic outflow + // (lambda_+ = u_out + c > 0 in both cases), so it is always extrapolated. + const amrex::Real Rp_N = u_out_N + qN.p / rho_c; + amrex::Real dRp = 0.0; + CellPrim qNm1, qNm2; + if (linear) { + qNm1 = cell_primitives(s_Nm1); + qNm2 = cell_primitives(s_Nm2); + if (qNm1.ok && qNm2.ok) { + const amrex::Real Rp_Nm1 = n_sgn * qNm1.u[idir] + qNm1.p / rho_c; + const amrex::Real Rp_Nm2 = n_sgn * qNm2.u[idir] + qNm2.p / rho_c; + dRp = minmod(Rp_N - Rp_Nm1, Rp_Nm1 - Rp_Nm2); + } else { + bump(Diag::eos_failure); + } + } + const amrex::Real Rp_g = Rp_N + static_cast(layer) * dRp; + + // ---- 5. Incoming acoustic invariant, R_- ------------------------------ + // lambda_- < 0 always (subsonic), so R_- is always modelled. Sign checks: + // outflow, p_N > p_target => dR_-/dn > 0 and (dp/dR_- = -rho c/2) the + // ghost pressure falls toward target; inflow, u above target => smaller + // R_-,ghost and the ghost normal velocity falls toward target. + const amrex::Real Rm_N = u_out_N - qN.p / rho_c; + const amrex::Real lambda_in = c_N - u_out_N; // > 0, magnitude of the speed + const amrex::Real one_m_M2 = + amrex::max(1.0 - mach * mach, constants::very_small_num()); + + // Upwinding of the lambda_0 (material) family. Its characteristics leave + // the domain only for u_out > 0; during a reversal an EXTRAPOLATED + // material ghost is advected back in and feeds on itself. That loop is + // the chamber's cold runaway: the outward-cooling T ramp refrigerates the + // boundary cell, the steepened ramp extrapolates colder still, and + // 1/(rho c) amplifies the relaxation as T falls -- 385 -> 241 -> 89 K in + // 42 us, then NaN (Docs/NSCBC-reversal-branch-defect.md). So the material + // SLOPES are used only where their upwind justification holds: full for + // u_out >= 0 (every forward-flow result is untouched to the bit), ramped + // off across a narrow band below zero (1e-3 c -- a continuity + // regularisation of the fill, not a physics knob; C13 gates it), and zero + // under firm reversal, where the ghost keeps the interior VALUES -- which + // is what the surviving closure always did. + const amrex::Real w_mat = + (u_out_N >= 0.0) + ? 1.0 + : amrex::max(0.0, 1.0 + u_out_N / (1.0e-3 * c_N)); + + const amrex::Real dRm = 0.0; + + amrex::Real Rm_g; + if (outflow_face && prm.pin_farfield) { + // Hard far-field pin: R_- takes its far-field value (u_out = 0, + // p = p_target). Value constraint, not a rate, and there is nowhere in it + // to put a transverse term -- so pin_farfield ignores beta. + Rm_g = -tgt.p / rho_c; + } else { + amrex::Real L_in; + if (outflow_face) { + const amrex::Real K = prm.sigma * one_m_M2 * c_N / prm.L_ref; + L_in = K * (qN.p - tgt.p); + } else { + const amrex::Real u_out_tgt = n_sgn * tgt.u[idir]; + const amrex::Real K = prm.relax_u * one_m_M2 * c_N / prm.L_ref; + L_in = -K * qN.rho * c_N * (u_out_N - u_out_tgt); + // Acoustic-injection feed-forward (Target::dudt): the incoming wave + // carries the target's rate directly, so a forced signal does not + // have to fight its way through the relaxation. n_sgn converts the + // lab-frame rate to the outward frame; the sign follows the + // relaxation term's convention (positive L_in raises u_out through + // the Rm inversion), and the factor 2 is u = (R+ + R-)/2: imposing a + // velocity rate through the incoming invariant alone needs twice it. + L_in += 2.0 * rho_c * (n_sgn * tgt.dudt); + } + // There is deliberately NO diffusive source term on L_in. The + // ghost-cell form does not have the flux-form's viscous-condition gap: + // the diffusion operator reads these ghost cells, so a correct ghost + // closure (extrap_temperature, gated by C8 and C12) already carries the + // diffusive physics, and an amplitude-side term counts it twice + // (measured: it moves the C12 driver error 104 -> -911). Re-verify + // against C12 before revisiting. + Rm_g = + Rm_N + static_cast(layer) * dRm + + static_cast(layer) * dx_normal * L_in / (lambda_in * rho_c); + } + + // ---- 6. Invert the acoustic pair -------------------------------------- + amrex::Real u_out_g = 0.5 * (Rp_g + Rm_g); + amrex::Real p_g = 0.5 * rho_c * (Rp_g - Rm_g); + + const amrex::Real p_floor = amrex::max( + 1.0e-3 * (tgt.p > 0.0 ? tgt.p : qN.p), constants::very_small_num()); + if (!amrex::Math::isfinite(p_g) || (p_g < p_floor)) { + p_g = p_floor; + bump(Diag::floored); + } + if (!amrex::Math::isfinite(u_out_g)) { + u_out_g = u_out_N; + bump(Diag::floored); + } + + amrex::Real u_g[3] = {qN.u[0], qN.u[1], qN.u[2]}; + u_g[idir] = n_sgn * u_out_g; + + // ---- 7. The lambda_0 family ------------------------------------------- + // Entropy, tangential velocities, species and passive scalars all convect at + // u_out. At an OUTFLOW they leave the domain and are extrapolated; at an + // INFLOW they enter and are set from the target -- composition hard (all + // NUM_SPECIES characteristics are incoming, so it is well posed), and + // temperature and tangential velocity through a relaxation. + amrex::Real rho_g = 0.0, T_g = 0.0, e_g = 0.0; + amrex::Real Y_g[NUM_SPECIES]; + + if (outflow_face) { + // Entropy invariant extrapolated like R_+; zero slope reduces to the + // linearised isentrope rho_g = rho_N + (p_g - p_N)/c^2. + const amrex::Real inv_c2 = 1.0 / (c_N * c_N); + const amrex::Real S_N = qN.rho - qN.p * inv_c2; + amrex::Real dS = 0.0; + amrex::Real dY[NUM_SPECIES] = {0.0}; + amrex::Real dut[3] = {0.0, 0.0, 0.0}; + if (linear && qNm1.ok && qNm2.ok) { + // All lambda_0 slopes carry the material upwinding factor: extrapolated + // material content is only valid where it advects OUT (see w_mat). + const amrex::Real S_Nm1 = qNm1.rho - qNm1.p * inv_c2; + const amrex::Real S_Nm2 = qNm2.rho - qNm2.p * inv_c2; + dS = w_mat * minmod(S_N - S_Nm1, S_Nm1 - S_Nm2); + for (int n = 0; n < NUM_SPECIES; n++) { + dY[n] = w_mat * minmod(qN.Y[n] - qNm1.Y[n], qNm1.Y[n] - qNm2.Y[n]); + } + for (int d = 0; d < 3; d++) { + if (d != idir) { + dut[d] = w_mat * minmod(qN.u[d] - qNm1.u[d], qNm1.u[d] - qNm2.u[d]); + } + } + } + const auto fl = static_cast(layer); + const amrex::Real S_g = S_N + fl * dS; + rho_g = S_g + p_g * inv_c2; + + amrex::Real Ysum = 0.0; + for (int n = 0; n < NUM_SPECIES; n++) { + Y_g[n] = amrex::max(qN.Y[n] + fl * dY[n], 0.0); + Ysum += Y_g[n]; + } + if (Ysum > constants::very_small_num()) { + const amrex::Real inv = 1.0 / Ysum; + for (amrex::Real& yg : Y_g) { + yg *= inv; + } + } else { + for (int n = 0; n < NUM_SPECIES; n++) { + Y_g[n] = qN.Y[n]; + } + bump(Diag::floored); + } + + for (int d = 0; d < 3; d++) { + if (d != idir) { + u_g[d] = qN.u[d] + fl * dut[d]; + } + } + + const amrex::Real rho_floor = + amrex::max(1.0e-6 * qN.rho, constants::very_small_num()); + if (!amrex::Math::isfinite(rho_g) || (rho_g < rho_floor)) { + rho_g = rho_floor; + bump(Diag::floored); + } + eos.RYP2T(rho_g, Y_g, p_g, T_g); + eos.RTY2E(rho_g, T_g, Y_g, e_g); + + } else { + // Inflow. One dimensionless nudge factor shared by temperature and the + // tangential velocities -- both ride lambda_0, magnitude |u_out|: + // nu = layer dx relax_t c / (L_ref |u_out|), clamped to [0, 1], + // applied as an explicit move toward the target so no sign convention + // can go wrong. + const amrex::Real u_mag = + amrex::max(std::abs(u_out_N), constants::smallu()); + const amrex::Real nu = amrex::min( + static_cast(layer) * dx_normal * prm.relax_t * c_N / + (prm.L_ref * u_mag), + 1.0); + + T_g = qN.T - nu * (qN.T - tgt.T); + for (int d = 0; d < 3; d++) { + if (d != idir) { + u_g[d] = qN.u[d] - nu * (qN.u[d] - tgt.u[d]); + } + } + // Composition: hard. Every species characteristic is incoming. + amrex::Real Ysum = 0.0; + for (int n = 0; n < NUM_SPECIES; n++) { + Y_g[n] = amrex::max(tgt.Y[n], 0.0); + Ysum += Y_g[n]; + } + if (Ysum > constants::very_small_num()) { + const amrex::Real inv = 1.0 / Ysum; + for (amrex::Real& yg : Y_g) { + yg *= inv; + } + } else { + for (int n = 0; n < NUM_SPECIES; n++) { + Y_g[n] = qN.Y[n]; + } + bump(Diag::floored); + } + if (!amrex::Math::isfinite(T_g) || !(T_g > 0.0)) { + T_g = qN.T; + bump(Diag::floored); + } + // Close on (p, T, Y): p from the acoustic pair, T and Y from the target. + eos.PYT2RE(p_g, Y_g, T_g, rho_g, e_g); + } + + if ( + !amrex::Math::isfinite(rho_g) || !(rho_g > 0.0) || + !amrex::Math::isfinite(e_g) || !amrex::Math::isfinite(T_g) || + !(T_g > 0.0)) { + bump(Diag::eos_failure); + copy_state(s_N, s_ghost); + return; + } + + // ---- 8. Pack and sanitise --------------------------------------------- + pack_ghost(rho_g, u_g, e_g, T_g, Y_g, s_N, s_ghost); + + for (int n = 0; n < NVAR; n++) { + if (!amrex::Math::isfinite(s_ghost[n])) { + s_ghost[n] = s_N[n]; + bump(Diag::floored); + } + } +} + +} // namespace pc_nscbc + +#endif diff --git a/Source/Params/_cpp_parameters b/Source/Params/_cpp_parameters index 732bf0b3e..78558dfba 100644 --- a/Source/Params/_cpp_parameters +++ b/Source/Params/_cpp_parameters @@ -67,10 +67,55 @@ use_hybrid_weno bool false # WENO scheme type in PPM method weno_scheme int 1 -# permits Ghost-Cells Navier-Stokes Boundary Conditions to be turned on and off -# for advective terms (adv) and for diffusion terms (diff) -nscbc_adv bool true -nscbc_diff bool false +# --------------------------------------------------------------------------- +# Characteristic (NSCBC) boundary treatment. See Source/NSCBC.H for the +# formulation and Docs/sphinx/BoundaryConditions.rst for guidance. +# +# Master switch. Even when true, the treatment is applied only at boundary +# POINTS for which the problem's ProblemSpecificFunctions::bcnormal_nscbc() +# returns a type other than `off`, so a problem that does not provide that +# hook is unaffected. +bc_nscbc bool false + +# Outflow pressure relaxation, the classical Poinsot-Lele sigma. The +# relaxation rate is K = sigma * (1 - M^2) * c / L, with L the domain extent +# along the boundary normal, so the boundary pressure relaxes toward the +# target over tau = 1/K, which is 1/(sigma (1 - M^2)) acoustic transit times. +# 0 is perfectly non-reflecting and leaves the mean pressure unanchored. +# Literature band is 0.15-0.3. Must be >= 0. +bc_nscbc_sigma Real 0.25 + +# Inflow normal-velocity relaxation coefficient. Larger imposes the target +# velocity more strongly and reflects more. Must be >= 0. +bc_nscbc_relax_u Real 2.0 + +# Inflow temperature and tangential-velocity relaxation coefficient. Must be +# >= 0: unlike the legacy Fortran, every user-facing NSCBC coefficient in +# PeleC is positive and the internal signs are handled by the kernel. +bc_nscbc_relax_t Real 0.2 + +# Order of the outgoing-invariant extrapolation, 1 or 2. Verification and +# debugging knob, not a physics knob; leave at 2. +bc_nscbc_order int 2 + +# Pin the incoming acoustic invariant to its far-field value instead of +# relaxing it toward the target at rate K. Simultaneously non-reflecting and +# anchored, but it anchors to p_target + rho c u rather than to p_target and +# its effective relaxation rate is c/dx, so it does not converge under mesh +# refinement. Appropriate for an open boundary onto a large quiescent +# reservoir; inappropriate for a duct exhausting into a plenum whose true mean +# pressure is not the target. bc_nscbc_sigma is ignored when this is set. +bc_nscbc_pin_farfield bool false + + +# NOTE: pelec.nscbc_adv and pelec.nscbc_diff were removed here. They were +# declared and queried but never read by any live code path -- the +# Ghost-Cell NSCBC implementation they controlled was deleted with the rest +# of the Fortran in 2d3a6f6. nscbc_adv defaulted to true, so leaving them in +# place meant that reconnecting a live NSCBC would silently change the +# behaviour of every Hard/UserBC case that did not explicitly opt out. +# PeleC::read_params aborts with a migration message if either key is present +# in an inputs file. See Docs/sphinx/BoundaryConditions.rst. # if true, define an external source term add_ext_src bool false @@ -219,8 +264,14 @@ PrT Real 1.0 # category: EB #----------------------------------------------------------------------------- -# set the EB boundary temperature for isothermal walls -eb_boundary_T Real 1.0 +# set the EB boundary temperature for isothermal walls. The default is +# ambient-like on purpose: this used to default to 1.0 K, and because +# eb_isothermal also defaults true, any diffusive EB case that did not set +# both parameters was silently conducting heat into one-kelvin walls -- the +# NSCBC-Chamber box variants refrigerated their cut cells to 77 K and NaN'd +# at the corners before this was found (see that case's README). A wall +# temperature is physics; set it deliberately. +eb_boundary_T Real 300.0 # flag for isothermal EB boundary eb_isothermal bool true diff --git a/Source/Params/param_includes/pelec_defaults.H b/Source/Params/param_includes/pelec_defaults.H index d24b70f75..2aeaffe88 100644 --- a/Source/Params/param_includes/pelec_defaults.H +++ b/Source/Params/param_includes/pelec_defaults.H @@ -18,8 +18,12 @@ bool PeleC::do_hydro = true; bool PeleC::do_mol = false; bool PeleC::use_hybrid_weno = false; int PeleC::weno_scheme = 1; -bool PeleC::nscbc_adv = true; -bool PeleC::nscbc_diff = false; +bool PeleC::bc_nscbc = false; +amrex::Real PeleC::bc_nscbc_sigma = 0.25; +amrex::Real PeleC::bc_nscbc_relax_u = 2.0; +amrex::Real PeleC::bc_nscbc_relax_t = 0.2; +int PeleC::bc_nscbc_order = 2; +bool PeleC::bc_nscbc_pin_farfield = false; bool PeleC::add_ext_src = false; amrex::GpuArray PeleC::external_forcing = {0.0}; bool PeleC::add_forcing_src = false; @@ -66,7 +70,7 @@ amrex::Real PeleC::Cs = 0.0; amrex::Real PeleC::Cw = 0.0; amrex::Real PeleC::CI = 0.0; amrex::Real PeleC::PrT = 1.0; -amrex::Real PeleC::eb_boundary_T = 1.0; +amrex::Real PeleC::eb_boundary_T = 300.0; bool PeleC::eb_isothermal = true; bool PeleC::eb_noslip = true; std::string PeleC::redistribution_type = "StateRedist"; diff --git a/Source/Params/param_includes/pelec_params.H b/Source/Params/param_includes/pelec_params.H index 1b45f8318..970779cd0 100644 --- a/Source/Params/param_includes/pelec_params.H +++ b/Source/Params/param_includes/pelec_params.H @@ -18,8 +18,12 @@ static bool do_hydro; static bool do_mol; static bool use_hybrid_weno; static int weno_scheme; -static bool nscbc_adv; -static bool nscbc_diff; +static bool bc_nscbc; +static amrex::Real bc_nscbc_sigma; +static amrex::Real bc_nscbc_relax_u; +static amrex::Real bc_nscbc_relax_t; +static int bc_nscbc_order; +static bool bc_nscbc_pin_farfield; static bool add_ext_src; static amrex::GpuArray external_forcing; static bool add_forcing_src; diff --git a/Source/Params/param_includes/pelec_queries.H b/Source/Params/param_includes/pelec_queries.H index 8aa47a63b..9f850362c 100644 --- a/Source/Params/param_includes/pelec_queries.H +++ b/Source/Params/param_includes/pelec_queries.H @@ -18,8 +18,12 @@ pp.query("do_hydro", do_hydro); pp.query("do_mol", do_mol); pp.query("use_hybrid_weno", use_hybrid_weno); pp.query("weno_scheme", weno_scheme); -pp.query("nscbc_adv", nscbc_adv); -pp.query("nscbc_diff", nscbc_diff); +pp.query("bc_nscbc", bc_nscbc); +pp.query("bc_nscbc_sigma", bc_nscbc_sigma); +pp.query("bc_nscbc_relax_u", bc_nscbc_relax_u); +pp.query("bc_nscbc_relax_t", bc_nscbc_relax_t); +pp.query("bc_nscbc_order", bc_nscbc_order); +pp.query("bc_nscbc_pin_farfield", bc_nscbc_pin_farfield); pp.query("add_ext_src", add_ext_src); { amrex::Vector tmp(AMREX_SPACEDIM, 0.0); diff --git a/Source/PeleC.H b/Source/PeleC.H index 11f36ddd3..df7da62b5 100644 --- a/Source/PeleC.H +++ b/Source/PeleC.H @@ -18,6 +18,7 @@ #include "PeleCUtilities.H" #include "Tagging.H" #include "IndexDefines.H" +#include "NSCBC.H" #include "prob_parm.H" #include "PelePhysics.H" #include "ReactorBase.H" @@ -93,6 +94,32 @@ public: ~PeleC() override; + // NSCBC accessors. pc_bcfill_hyp is a free function and the runtime + // parameters are protected statics, so these exist to let the boundary + // filler capture its parameters HOST-side into a POD functor member. A + // device kernel must never read a class static or ParmParse directly. + static bool nscbc_active() { return bc_nscbc; } + static pc_nscbc::Params nscbc_params(int idir); + // Report and reset the NSCBC fallback counters. A silent fallback is a bug + // that will not be found; this is how it surfaces. + static void nscbc_report_diagnostics(); + + // Boundary registers (Docs/NSCBC-boundary-registers-design.md): the + // once-per-advance update, and checkpoint/restart of the register array. + + // Warn when THIS refined level's grids touch a domain face that the + // characteristic treatment can reach (Hard/UserBC): the extrapolation + // stencil is level-local, so a fine patch on such a face makes the + // boundary condition level-dependent. Called from post_regrid and + // post_init; no-op on level 0 or when bc_nscbc is off. + void nscbc_check_fine_faces() const; + + // Gate the characteristic fill's periodicity: a ghost cell and its image + // one period away along a periodic tangential direction must agree to the + // last bit, since they are built from the same data by the same arithmetic. + // Runs once, from post_init, and only when a characteristic face has a + // periodic tangential direction. + // Restart from a checkpoint file. void restart( amrex::Amr& papa, std::istream& is, bool bReadSpecial = false) override; diff --git a/Source/PeleC.cpp b/Source/PeleC.cpp index 11f7951ea..fbba65586 100644 --- a/Source/PeleC.cpp +++ b/Source/PeleC.cpp @@ -194,6 +194,21 @@ PeleC::read_params() pp.query("v", verbose); + // Removed parameters. nscbc_adv/nscbc_diff controlled the Fortran + // Ghost-Cell NSCBC deleted in 2d3a6f6; they have had no effect since. + // Abort rather than ignore, so that an inputs file which believes it is + // configuring a characteristic boundary treatment is not silently wrong. + for (const char* removed : {"nscbc_adv", "nscbc_diff"}) { + if (pp.contains(removed)) { + amrex::Abort( + "pelec." + std::string(removed) + + " has been removed: it has had no effect since the Fortran GC-NSCBC " + "was deleted. Remove it from your inputs file. See " + "Docs/sphinx/BoundaryConditions.rst for the currently recommended " + "subsonic inflow/outflow treatment."); + } + } + // Get boundary conditions amrex::Vector lo_bc_char(AMREX_SPACEDIM); amrex::Vector hi_bc_char(AMREX_SPACEDIM); @@ -324,6 +339,88 @@ PeleC::read_params() } } + if (bc_nscbc) { + pc_nscbc::Params probe = nscbc_params(0); + probe.sigma = bc_nscbc_sigma; + std::string why; + if (!probe.validate(why)) { + amrex::Abort("pelec.bc_nscbc_*: " + why); + } + // The treatment only ever acts on ext_dir faces, i.e. Hard/UserBC. + bool any_ext_dir = false; + for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { + for (const auto& b : {lo_bc[dir], hi_bc[dir]}) { + any_ext_dir = any_ext_dir || (b == PCPhysBCType::user_bc) || + (b == PCPhysBCType::inflow); + } + } + if (!any_ext_dir) { + amrex::Abort( + "pelec.bc_nscbc = 1 but no boundary is Hard or UserBC. The " + "characteristic treatment only acts on those faces."); + } + // Chemistry integrated in ghost cells would burn the state the boundary + // condition just wrote there -- a hot inflow target composition would + // ignite in the ghost region. + if (state_nghost > 0) { + amrex::Abort( + "pelec.bc_nscbc = 1 with pelec.state_nghost > 0: reactions are " + "integrated over grown tileboxes, which would advance chemistry in " + "the ghost cells written by the boundary condition."); + } + // The characteristic fill has no EBCellFlag: it detects covered stencil + // cells by their non-positive density. The DEFAULT body state is a + // sampled fluid state (define_body_state) -- positive density, + // undetectable -- so without eb_zero_body_state the fill would silently + // read body values as if they were interior data. Fatal rather than + // silent: that is the standing rule for this boundary condition. + if (eb_in_domain && !eb_zero_body_state) { + amrex::Abort( + "pelec.bc_nscbc = 1 with EB geometry requires " + "pelec.eb_zero_body_state = 1: the characteristic fill detects " + "covered stencil cells by non-positive density, and the default " + "body state is a sampled fluid state it cannot distinguish from " + "interior data."); + } + if (amrex::ParallelDescriptor::IOProcessor()) { + amrex::Print() + << "\n NSCBC: characteristic boundary treatment ENABLED (advective " + "ghost fill).\n" + << " This is the 1-D-normal LODI limit: no transverse terms and no\n" + << " reaction-source correction. Place outflows away from flames,\n" + << " shear layers and composition fronts. See Source/NSCBC.H.\n"; + if (bc_nscbc_pin_farfield) { + amrex::Print() + << " incoming acoustic PINNED to the far field; bc_nscbc_sigma " + "ignored.\n" + << " Anchors to p_target + rho*c*u and does not converge under " + "mesh refinement.\n"; + } else { + amrex::Print() << " sigma = " << bc_nscbc_sigma; + if (bc_nscbc_sigma > 0.0) { + amrex::Print() + << " -> outflow pressure relaxes over " << 1.0 / bc_nscbc_sigma + << " acoustic transit times at low Mach\n" + << " (tau_relax / tau_acoustic = 1 / (sigma (1 - " + "M^2)));\n" + << " it must be much larger than 1 and much smaller " + "than the run time.\n"; + } else { + amrex::Print() << " -> perfectly non-reflecting; the mean pressure " + "is NOT anchored and will drift.\n"; + } + } + amrex::Print() << " relax_u = " << bc_nscbc_relax_u + << ", relax_t = " << bc_nscbc_relax_t + << ", order = " << bc_nscbc_order << "\n"; + amrex::Print() << " L_ref = ("; + for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) { + amrex::Print() << (dir > 0 ? ", " : "") << nscbc_params(dir).L_ref; + } + amrex::Print() << ")\n\n"; + } + } + if (amrex::DefaultGeometry().IsRZ() && (lo_bc[0] != PCPhysBCType::symmetry)) { amrex::Error( "PeleC::read_params: must set r=0 boundary condition to " @@ -1157,6 +1254,19 @@ PeleC::postCoarseTimeStep(amrex::Real cumtime) { BL_PROFILE("PeleC::postCoarseTimeStep()"); AmrLevel::postCoarseTimeStep(cumtime); + + // NSCBC counters must never be silent: without this, a run with the + // default sum_interval = -1 never calls nscbc_report_diagnostics() and a + // boundary counting millions of reversal fills looks identical to one + // counting none (that misread happened; see + // Docs/NSCBC-reversal-branch-defect.md). When the user has not opted + // into periodic reporting, report on a fixed cadence -- the report only + // prints when something actually counted, so healthy runs stay quiet. + if (bc_nscbc && verbose > 0 && sum_interval <= 0 && sum_per <= 0.0) { + if (parent->levelSteps(0) % 100 == 0) { + nscbc_report_diagnostics(); + } + } } void @@ -1176,6 +1286,8 @@ PeleC::post_regrid(int lbase, int /*new_finest*/) if ((do_react) && (use_typical_vals_chem)) { set_typical_values_chem(); } + + nscbc_check_fine_faces(); } void @@ -1183,6 +1295,13 @@ PeleC::post_init(amrex::Real /*stop_time*/) { BL_PROFILE("PeleC::post_init()"); + if (level == 0) { + const int finest = parent->finestLevel(); + for (int lev = 1; lev <= finest; ++lev) { + getLevel(lev).nscbc_check_fine_faces(); + } + } + amrex::Real dtlev = parent->dtLevel(level); amrex::Real cumtime = parent->cumTime(); diff --git a/Source/ProblemSpecificFunctions.H b/Source/ProblemSpecificFunctions.H index ef8262ae3..d6b98a474 100644 --- a/Source/ProblemSpecificFunctions.H +++ b/Source/ProblemSpecificFunctions.H @@ -7,6 +7,7 @@ #include "PeleC.H" #include "IndexDefines.H" +#include "NSCBC.H" struct DefaultProblemSpecificFunctions { @@ -160,6 +161,37 @@ struct DefaultProblemSpecificFunctions { } + // Characteristic (NSCBC) boundary treatment. + // + // Return the target state for one boundary POINT. The default returns + // Type::off, which leaves that point to the ordinary bcnormal() path, so + // pelec.bc_nscbc = 1 is a no-op for a problem that does not override this. + // Because the decision is per point, a single face may mix an inflow, an + // outflow and a wall. + // + // Only the fields relevant to the returned type are read: + // Type::outflow -> p + // Type::inflow -> u[], T, Y[] + // There is deliberately no relaxation coefficient for composition; see + // Source/NSCBC.H. + // + // s_int is the state in the first interior cell along the boundary normal, + // for problems whose target depends on the interior state (a back-pressure + // that tracks a measured plenum, say). Most problems will not need it. + AMREX_GPU_DEVICE + AMREX_FORCE_INLINE + static pc_nscbc::Target bcnormal_nscbc( + const amrex::Real* /*x[AMREX_SPACEDIM]*/, + const amrex::Real* /*s_int[NVAR]*/, + const int /*idir*/, + const int /*sgn*/, + const amrex::Real /*time*/, + amrex::GeometryData const& /*geomdata*/, + ProbParmDevice const& /*prob_parm*/) + { + return pc_nscbc::Target{}; + } + /* TODO static void something_with_ADV_AUX_transport() { diff --git a/Source/SumIQ.cpp b/Source/SumIQ.cpp index 9b3daf376..5c6bfd84d 100644 --- a/Source/SumIQ.cpp +++ b/Source/SumIQ.cpp @@ -11,6 +11,8 @@ PeleC::sum_integrated_quantities() return; } + nscbc_report_diagnostics(); + bool local_flag = true; int finest_level = parent->finestLevel(); From ece2d39d8df235325625dde3986ea335fa38cb40 Mon Sep 17 00:00:00 2001 From: Marc Day Date: Thu, 27 Aug 2026 22:13:57 +0200 Subject: [PATCH 2/3] NSCBC 2/7: transverse terms in the modelled incoming wave The multi-dimensional correction: -(1 - beta) T_in on the incoming amplitude, where T_in carries the tangential convective and dilatational terms in the outward frame (the ghost-cell transplantation of Motheau's T1/T4 -- one expression for lo and hi faces). Tangential derivatives come from the boundary cell's clamped neighbours; at a corner the stencil collapses one-sided or to nothing, which is the measured replacement for the AIAA corner-coupling machinery. The beta blend is exact for a plane wave at incidence theta: beta_opt = 1 - cos(theta)/(1 + cos(theta)) -- 0.5 head-on, toward 1 at grazing -- and both measured optima in the COVO/pulse suite land on it. The response is asymmetric (a wrong beta is worse than an absent one), so the default stays 1: this PR is bit-inert until a problem opts in, and every result of the previous PR is unchanged. The clamp leaves a small measured seam residual where a characteristic face meets a periodic direction, and that trade is guarded rather than assumed: nscbc_check_periodic_wrap() measures the ghost-image mismatch at startup, reports it beside the boundary-row spread (so a vacuous pass is visible), and aborts above 1e-2 -- an order above the worst known-good measurement and an order below what a broken stencil produces. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01XU8M23nucKsu1do2WxXFeq --- Source/BCfill.cpp | 217 +++++++++++++++++- Source/NSCBC.H | 88 ++++++- Source/Params/_cpp_parameters | 22 ++ Source/Params/param_includes/pelec_defaults.H | 1 + Source/Params/param_includes/pelec_params.H | 1 + Source/Params/param_includes/pelec_queries.H | 1 + Source/PeleC.H | 1 + Source/PeleC.cpp | 1 + 8 files changed, 330 insertions(+), 2 deletions(-) diff --git a/Source/BCfill.cpp b/Source/BCfill.cpp index 51b905345..6089534e0 100644 --- a/Source/BCfill.cpp +++ b/Source/BCfill.cpp @@ -188,10 +188,42 @@ struct PCHypFillExtDir // Boundary-register composition (level 0 only; the registers change // once per advance and are read frozen here). The kernel below sees // only the composed Target -- it stays a pure function. + // Tangential neighbours of the boundary cell, for the transverse terms, + // clamped into the domain: at a corner the clamp collapses the stencil + // and inv_dt falls to a one-sided spacing, or to zero if both + // neighbours land on the same cell. + pc_nscbc::Transverse tr; + amrex::Real s_tm[AMREX_SPACEDIM][NVAR], s_tp[AMREX_SPACEDIM][NVAR]; + if (m_nscbc_prm[idir].beta < 1.0) { + for (int d = 0; d < AMREX_SPACEDIM; ++d) { + if (d == idir) { + continue; + } + const int dlo = amrex::max(domlo[d], fab_lo[d]); + const int dhi = amrex::min(domhi[d], fab_hi[d]); + const int jm = amrex::max(ivN[d] - 1, dlo); + const int jp = amrex::min(ivN[d] + 1, dhi); + if (jp == jm) { + continue; + } + amrex::IntVect ivm(ivN), ivp(ivN); + ivm[d] = jm; + ivp[d] = jp; + for (int n = 0; n < NVAR; n++) { + s_tm[d][n] = dest(ivm, n); + s_tp[d][n] = dest(ivp, n); + } + tr.sm[d] = s_tm[d]; + tr.sp[d] = s_tp[d]; + tr.inv_dt[d] = 1.0 / (static_cast(jp - jm) * dx[d]); + tr.valid = true; + } + } + amrex::Real s_ghost[NVAR]; pc_nscbc::apply( s_N, s_Nm1, s_Nm2, n_stencil, dx[idir], idir, sgn, layer, tgt, - m_nscbc_prm[idir], s_ghost, m_nscbc_diag); + m_nscbc_prm[idir], s_ghost, m_nscbc_diag, &tr); for (int n = 0; n < NVAR; n++) { dest(iv, n) = s_ghost[n]; } @@ -483,6 +515,188 @@ PeleC::nscbc_check_fine_faces() const } } +// --------------------------------------------------------------------------- +// Periodic-seam gate. Measures how far the characteristic fill is from +// periodic where the domain is: the worst relative mismatch between ghost +// cells and their images one period away. The clamped tangential stencil +// (see nscbc_fill) leaves a deliberate residual at the seam under amrex's +// corner-strip protocol -- measured 2e-4 on the inert vortex and 1.6e-3 +// with a flame front sitting on the seam corner -- so the gate REPORTS the +// measured value and aborts only above 1e-2, an order of margin above the +// worst known-good state and an order below a broken stencil (the naive +// wrap measured 2.2e-2 here). The report guards its own blind spots: the +// tangential spread of the boundary row is printed beside the mismatch (a pass +// on a boundary-uniform row gates nothing), and a decomposition with no image +// pair in one FAB says NOT CHECKED instead of passing. Exercised by +// NSCBC-COVO/nscbc-wrapgate.inp. +// --------------------------------------------------------------------------- +void +PeleC::nscbc_check_periodic_wrap() +{ + if (!bc_nscbc || (level != 0)) { + return; + } + static bool done = false; + if (done) { + return; + } + + const amrex::Box& dom = geom.Domain(); + auto characteristic = [&](const int dir, const int side) { + const int t = (side == 0) ? phys_bc.lo(dir) : phys_bc.hi(dir); + return (t == PCPhysBCType::inflow) || (t == PCPhysBCType::user_bc); + }; + + bool relevant = false; + for (int idir = 0; idir < AMREX_SPACEDIM; ++idir) { + for (int d = 0; d < AMREX_SPACEDIM; ++d) { + relevant = + relevant || ((characteristic(idir, 0) || characteristic(idir, 1)) && + (d != idir) && geom.isPeriodic(d)); + } + } + if (!relevant) { + return; + } + done = true; + + const int ng = numGrow(); + amrex::MultiFab S(grids, dmap, NVAR, ng, amrex::MFInfo(), Factory()); + FillPatch(*this, S, ng, state[State_Type].curTime(), State_Type, 0, NVAR); + + const amrex::Real big = std::numeric_limits::max(); + + for (int idir = 0; idir < AMREX_SPACEDIM; ++idir) { + for (int side = 0; side < 2; ++side) { + if (!characteristic(idir, side)) { + continue; + } + const int N_pos = (side == 0) ? dom.smallEnd(idir) : dom.bigEnd(idir); + + for (int d = 0; d < AMREX_SPACEDIM; ++d) { + if ((d == idir) || (!geom.isPeriodic(d))) { + continue; + } + const int n_d = dom.length(d); + if (n_d < ng) { + continue; // the image slab would fold over itself + } + + // The ghost cells of this face that sit above the domain in d. The + // image of each is n_d cells below, and is a ghost of this face too. + amrex::Box reg = amrex::grow(dom, ng); + if (side == 0) { + reg.setBig(idir, dom.smallEnd(idir) - 1); + } else { + reg.setSmall(idir, dom.bigEnd(idir) + 1); + } + reg.setSmall(d, dom.bigEnd(d) + 1); + reg.setBig(d, dom.bigEnd(d) + ng); + + const amrex::IntVect img = -n_d * amrex::IntVect::TheDimensionVector(d); + + amrex::ReduceOps< + amrex::ReduceOpMax, amrex::ReduceOpSum, amrex::ReduceOpMax, + amrex::ReduceOpMin> + op; + amrex::ReduceData + rd(op); + using RT = typename decltype(rd)::Type; + + for (amrex::MFIter mfi(S); mfi.isValid(); ++mfi) { + auto const& a = S.const_array(mfi); + const amrex::Box& fbx = mfi.fabbox(); + + // The image pairs this FAB holds both halves of. + amrex::Box sh(fbx); + sh.shift(-img); + const amrex::Box pbx = fbx & reg & sh; + if (!pbx.isEmpty()) { + op.eval(pbx, rd, [=] AMREX_GPU_DEVICE(int i, int j, int k) -> RT { + amrex::ignore_unused(k); // 2-D: AMREX_D_DECL drops it + const amrex::IntVect iv(AMREX_D_DECL(i, j, k)); + const amrex::IntVect iw = iv + img; + amrex::Real e = 0.0; + for (int n = 0; n < NVAR; n++) { + const amrex::Real u = a(iv, n); + const amrex::Real v = a(iw, n); + const amrex::Real s = amrex::max( + amrex::Math::abs(u), + amrex::max( + amrex::Math::abs(v), + std::numeric_limits::min())); + e = amrex::max(e, amrex::Math::abs(u - v) / s); + } + return {e, amrex::Long(1), -big, big}; + }); + } + + // The tangential structure of the boundary row itself: valid data, + // so this measures whether the check above could have failed. + amrex::Box row = dom; + row.setSmall(idir, N_pos); + row.setBig(idir, N_pos); + row &= mfi.validbox(); + if (!row.isEmpty()) { + op.eval(row, rd, [=] AMREX_GPU_DEVICE(int i, int j, int k) -> RT { + const amrex::Real r = a(i, j, k, URHO); + return {0.0, amrex::Long(0), r, r}; + }); + } + } + + auto hv = rd.value(op); + amrex::Real worst = amrex::get<0>(hv); + amrex::Long npairs = amrex::get<1>(hv); + amrex::Real rmax = amrex::get<2>(hv); + amrex::Real rmin = amrex::get<3>(hv); + amrex::ParallelDescriptor::ReduceRealMax(worst); + amrex::ParallelDescriptor::ReduceLongSum(npairs); + amrex::ParallelDescriptor::ReduceRealMax(rmax); + amrex::ParallelDescriptor::ReduceRealMin(rmin); + const amrex::Real spread = + (rmax > -big) + ? (rmax - rmin) / amrex::max( + amrex::Math::abs(rmax), + std::numeric_limits::min()) + : 0.0; + + constexpr amrex::Real tol = 1.0e-2; + if (worst > tol) { + amrex::Abort( + "NSCBC periodic-seam check FAILED on direction " + + std::to_string(idir) + " " + (side == 0 ? "lo" : "hi") + + " with periodic tangential direction " + std::to_string(d) + + ": worst relative mismatch " + std::to_string(worst) + " over " + + std::to_string(npairs) + + " image pairs. The characteristic fill is not periodic where the " + "domain is."); + } + if (amrex::ParallelDescriptor::IOProcessor() && (verbose > 0)) { + amrex::Print() << " NSCBC periodic-seam check: dir " << idir + << (side == 0 ? " lo" : " hi") + << ", periodic tangential dir " << d << " -- "; + if (npairs == 0) { + amrex::Print() + << "NOT CHECKED: no box holds a ghost cell and its image " + "together. Raise amr.max_grid_size in direction " + << d << " to span the domain if you want this gated.\n"; + } else { + amrex::Print() << npairs << " image pairs agree to " << worst + << " (boundary-row density spread " << spread + << ")\n"; + if (spread == 0.0) { + amrex::Print() + << " (that row is uniform along the boundary, so this pass " + "is vacuous -- a clamped stencil would pass it too.)\n"; + } + } + } + } + } + } +} + pc_nscbc::Params PeleC::nscbc_params(const int idir) { @@ -491,6 +705,7 @@ PeleC::nscbc_params(const int idir) p.relax_u = bc_nscbc_relax_u; p.relax_t = bc_nscbc_relax_t; p.order = bc_nscbc_order; + p.beta = bc_nscbc_beta; p.pin_farfield = bc_nscbc_pin_farfield; // Only the ratio sigma/L_ref is physical. L_ref is fixed to the domain // extent along the boundary normal so that sigma keeps the meaning it has diff --git a/Source/NSCBC.H b/Source/NSCBC.H index 57a1bd088..4f187cad0 100644 --- a/Source/NSCBC.H +++ b/Source/NSCBC.H @@ -91,6 +91,20 @@ struct Params // Inflow temperature and tangential-velocity relaxation. amrex::Real relax_t = 0.2; + // Weight on the transverse terms in the modelled incoming acoustic: + // L_in = K (phi - phi_target) - (1 - beta) * T_transverse + // 1 discards them (the 1-D-normal LODI limit, and the default: enabling + // them is an opt-in change), 0 includes them fully, < 0 selects the + // pointwise local Mach number (the legacy Fortran convention). + // + // Exact blend for a plane wave at angle theta: + // beta_opt = 1 - cos(theta)/(1+cos(theta)) -- 0.5 head-on, toward 1 at + // grazing; both measured optima land on it (NSCBC-COVO/README.md). The + // response is asymmetric: beta = 0 is near-unstable and pointwise + // local-Mach is worse than no correction. A wrong beta is more dangerous + // than an absent one: start at 0.5, raise toward 1 for grazing incidence. + amrex::Real beta = 1.0; + // Boundary-normal reference length, probhi[idir] - problo[idir]. Filled per // face by the caller. Only the ratio sigma/L_ref is physical, so this is // not a user knob; it exists so that sigma keeps its literature meaning. @@ -125,6 +139,11 @@ struct Params why = "bc_nscbc_relax_t must be >= 0"; return false; } + if (beta > 1.0) { + why = "bc_nscbc_beta must be <= 1 (1 discards the transverse terms, 0 " + "includes them fully, negative selects the local Mach number)"; + return false; + } if ((order != 1) && (order != 2)) { why = "bc_nscbc_order must be 1 or 2"; return false; @@ -167,6 +186,22 @@ struct Target amrex::Real dudt = 0.0; // d(u_normal target)/dt, LAB frame [cm/s^2] }; +// --------------------------------------------------------------------------- +// Tangential neighbours of the BOUNDARY CELL, for the transverse terms +// T_in = sum_t [ u_t (dp/dt - rho c du_out/dt) ] + gamma p sum_t du_t/dt +// (derivatives at the boundary cell; the outward frame absorbs the lo/hi +// sign flip). Neighbour indexing is the caller's -- see tang_range() in +// BCfill.cpp -- and the spacing arrives as an explicit reciprocal, zero when +// a collapsed corner stencil leaves nothing to difference. +// --------------------------------------------------------------------------- +struct Transverse +{ + const amrex::Real* sm[AMREX_SPACEDIM] = {}; + const amrex::Real* sp[AMREX_SPACEDIM] = {}; + amrex::Real inv_dt[AMREX_SPACEDIM] = {}; + bool valid = false; +}; + // --------------------------------------------------------------------------- // Diagnostic counters. Every fallback is counted; a silent fallback is a // bug that will not be found. Storage is amrex::Long, not int: with the @@ -373,7 +408,8 @@ apply( Target const& tgt, Params const& prm, amrex::Real s_ghost[NVAR], - amrex::Long* diag = nullptr) noexcept + amrex::Long* diag = nullptr, + Transverse const* tr = nullptr) noexcept { auto eos = pele::physics::PhysicsType::eos(); @@ -545,6 +581,55 @@ apply( ? 1.0 : amrex::max(0.0, 1.0 + u_out_N / (1.0e-3 * c_N)); + // Transverse contribution to the modelled incoming wave. Additive on the + // amplitude, which is exactly why the relaxation is written as a gradient + // increment: there is a well-defined slot for it here, and none at all in a + // hard far-field pin of R_-. + amrex::Real beta = prm.beta; + if (beta < 0.0) { + beta = amrex::min(std::abs(mach), 1.0); + } + amrex::Real T_in = 0.0; + if ((tr != nullptr) && tr->valid && (beta < 1.0)) { + amrex::Real gam = 0.0; + eos.RTY2G(qN.rho, qN.T, qN.Y, gam); + amrex::Real div_ut = 0.0; + bool any = false; + for (int d = 0; d < AMREX_SPACEDIM; ++d) { + if ( + (d == idir) || (tr->sm[d] == nullptr) || (tr->sp[d] == nullptr) || + (tr->inv_dt[d] == 0.0)) { + continue; + } + const CellPrim qm = cell_primitives(tr->sm[d]); + const CellPrim qp = cell_primitives(tr->sp[d]); + if (!qm.ok || !qp.ok) { + bump(Diag::transverse_drop); + continue; + } + const amrex::Real dpdt = (qp.p - qm.p) * tr->inv_dt[d]; + // du_out/dt: the normal velocity in the OUTWARD frame, differentiated + // along the tangential direction. The n_sgn here is what turns the + // legacy hi-face T1 and lo-face T4 into one expression. + const amrex::Real duodt = + n_sgn * (qp.u[idir] - qm.u[idir]) * tr->inv_dt[d]; + const amrex::Real dutdt = (qp.u[d] - qm.u[d]) * tr->inv_dt[d]; + T_in += qN.u[d] * (dpdt - rho_c * duodt); + div_ut += dutdt; + any = true; + } + if (any && amrex::Math::isfinite(gam) && (gam > 0.0)) { + T_in += gam * qN.p * div_ut; + } else { + T_in = 0.0; + bump(Diag::transverse_drop); + } + if (!amrex::Math::isfinite(T_in)) { + T_in = 0.0; + bump(Diag::transverse_drop); + } + } + const amrex::Real dRm = 0.0; amrex::Real Rm_g; @@ -571,6 +656,7 @@ apply( // velocity rate through the incoming invariant alone needs twice it. L_in += 2.0 * rho_c * (n_sgn * tgt.dudt); } + L_in -= (1.0 - beta) * T_in; // There is deliberately NO diffusive source term on L_in. The // ghost-cell form does not have the flux-form's viscous-condition gap: // the diffusion operator reads these ghost cells, so a correct ghost diff --git a/Source/Params/_cpp_parameters b/Source/Params/_cpp_parameters index 78558dfba..3c1397804 100644 --- a/Source/Params/_cpp_parameters +++ b/Source/Params/_cpp_parameters @@ -94,6 +94,28 @@ bc_nscbc_relax_u Real 2.0 # PeleC is positive and the internal signs are handled by the kernel. bc_nscbc_relax_t Real 0.2 +# Weight on the transverse terms in the modelled incoming acoustic wave: +# 1 discard them -- the 1-D-normal LODI limit, and the DEFAULT +# 0 include them fully +# <0 use the local Mach number (the legacy Fortran convention; not +# recommended, see below) +# +# The transverse terms are what a boundary needs when a wave arrives obliquely, +# or when a vortex rather than a sound wave crosses it. Measured in +# Exec/RegTests/NSCBC-COVO, against a periodic no-boundary reference: +# +# convected vortex, rms residual: hard 14.0x floor, beta=1 10.0x, +# beta=0.5 3.8x, beta=0.2 16.1x, beta=0 132x +# circular pulse, front amplitude spread: hard 21.4%, beta=1 0.90%, +# beta=0.8 0.066%, beta=0 2.9% +# +# The optimum follows beta = 1 - cos(theta)/(1 + cos(theta)) for a wave meeting +# the boundary at angle theta: 0.5 at normal incidence, rising toward 1 at +# grazing incidence. Too little correction costs a factor of a few; too much +# is catastrophic. RECOMMENDED STARTING VALUE 0.5. The default is 1 only +# because a wrong beta is far more dangerous than an absent one. Must be <= 1. +bc_nscbc_beta Real 1.0 + # Order of the outgoing-invariant extrapolation, 1 or 2. Verification and # debugging knob, not a physics knob; leave at 2. bc_nscbc_order int 2 diff --git a/Source/Params/param_includes/pelec_defaults.H b/Source/Params/param_includes/pelec_defaults.H index 2aeaffe88..efa1a41fc 100644 --- a/Source/Params/param_includes/pelec_defaults.H +++ b/Source/Params/param_includes/pelec_defaults.H @@ -22,6 +22,7 @@ bool PeleC::bc_nscbc = false; amrex::Real PeleC::bc_nscbc_sigma = 0.25; amrex::Real PeleC::bc_nscbc_relax_u = 2.0; amrex::Real PeleC::bc_nscbc_relax_t = 0.2; +amrex::Real PeleC::bc_nscbc_beta = 1.0; int PeleC::bc_nscbc_order = 2; bool PeleC::bc_nscbc_pin_farfield = false; bool PeleC::add_ext_src = false; diff --git a/Source/Params/param_includes/pelec_params.H b/Source/Params/param_includes/pelec_params.H index 970779cd0..c5c6f80e7 100644 --- a/Source/Params/param_includes/pelec_params.H +++ b/Source/Params/param_includes/pelec_params.H @@ -22,6 +22,7 @@ static bool bc_nscbc; static amrex::Real bc_nscbc_sigma; static amrex::Real bc_nscbc_relax_u; static amrex::Real bc_nscbc_relax_t; +static amrex::Real bc_nscbc_beta; static int bc_nscbc_order; static bool bc_nscbc_pin_farfield; static bool add_ext_src; diff --git a/Source/Params/param_includes/pelec_queries.H b/Source/Params/param_includes/pelec_queries.H index 9f850362c..b381f61d0 100644 --- a/Source/Params/param_includes/pelec_queries.H +++ b/Source/Params/param_includes/pelec_queries.H @@ -22,6 +22,7 @@ pp.query("bc_nscbc", bc_nscbc); pp.query("bc_nscbc_sigma", bc_nscbc_sigma); pp.query("bc_nscbc_relax_u", bc_nscbc_relax_u); pp.query("bc_nscbc_relax_t", bc_nscbc_relax_t); +pp.query("bc_nscbc_beta", bc_nscbc_beta); pp.query("bc_nscbc_order", bc_nscbc_order); pp.query("bc_nscbc_pin_farfield", bc_nscbc_pin_farfield); pp.query("add_ext_src", add_ext_src); diff --git a/Source/PeleC.H b/Source/PeleC.H index df7da62b5..ce9c82344 100644 --- a/Source/PeleC.H +++ b/Source/PeleC.H @@ -119,6 +119,7 @@ public: // last bit, since they are built from the same data by the same arithmetic. // Runs once, from post_init, and only when a characteristic face has a // periodic tangential direction. + void nscbc_check_periodic_wrap(); // Restart from a checkpoint file. void restart( diff --git a/Source/PeleC.cpp b/Source/PeleC.cpp index fbba65586..38e99f89e 100644 --- a/Source/PeleC.cpp +++ b/Source/PeleC.cpp @@ -1300,6 +1300,7 @@ PeleC::post_init(amrex::Real /*stop_time*/) for (int lev = 1; lev <= finest; ++lev) { getLevel(lev).nscbc_check_fine_faces(); } + nscbc_check_periodic_wrap(); } amrex::Real dtlev = parent->dtLevel(level); From ecb42d7bd2063fe94839da7e1612afa3c7269780 Mon Sep 17 00:00:00 2001 From: Marc Day Date: Thu, 27 Aug 2026 22:14:51 +0200 Subject: [PATCH 3/3] NSCBC 3/7: reacting boundaries and fronts What a flame does to a characteristic boundary, and the closures that survive it -- all default-off, so the previous two PRs' results are unchanged to the bit until a problem opts in. - The reaction source enters the modelled incoming wave with +(1 - beta_s) S_p (Sutherland-Kennedy): a flame in the boundary cells raises pressure at a rate the relaxation has no model for, and without the term the mean sits off target by S_p/K. The sign is PLUS -- the exact steady flame has dR_-/dn = +S_p/(rho c^2), and the driver PR carries a convention-free check of exactly that. S_p is the compositional dp/dt at fixed (rho, e): closed form for ideal gas, a directional finite difference along the reaction path for SRK. - extrap_temperature closes the lambda_0 family on temperature instead of the entropy invariant. The diffusion operator forms boundary-face fluxes from these ghosts, so the ghost dT/dn IS the heat flux leaving; the entropy closure overstates a flame-like face gradient by 40%. This is the flag that lets a flame survive the outflow. - extrap_material continues the material part of R_- into the ghost, bounded through the acoustically-blind entropy family (extrapolating R_-'s own gradient is positive feedback, measured 16x). For fronts that SIT on the boundary; a transit over-vents with it on. - backflow_material treats FIRM reversal as a local inflow of the Target's reservoir state, ramped over an outward-Mach band continuous with the reversal upwinding underneath -- without it a face in sustained recirculation feeds the domain its own exhaust. Also the structure advisory: material structure sitting in outflow boundary cells is counted, and the counter report prints the measured flame-on-boundary recipe when it fires. Measured guidance and the sitting/transit tables live with the regression cases later in this stack (NSCBC-FlameOutflow, README). Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01XU8M23nucKsu1do2WxXFeq --- Source/BCfill.cpp | 22 +- Source/NSCBC.H | 260 +++++++++++++++++- Source/Params/_cpp_parameters | 49 ++++ Source/Params/param_includes/pelec_defaults.H | 4 + Source/Params/param_includes/pelec_params.H | 4 + Source/Params/param_includes/pelec_queries.H | 4 + 6 files changed, 334 insertions(+), 9 deletions(-) diff --git a/Source/BCfill.cpp b/Source/BCfill.cpp index 6089534e0..5efe03a1b 100644 --- a/Source/BCfill.cpp +++ b/Source/BCfill.cpp @@ -706,6 +706,7 @@ PeleC::nscbc_params(const int idir) p.relax_t = bc_nscbc_relax_t; p.order = bc_nscbc_order; p.beta = bc_nscbc_beta; + p.beta_s = bc_nscbc_beta_s; p.pin_farfield = bc_nscbc_pin_farfield; // Only the ratio sigma/L_ref is physical. L_ref is fixed to the domain // extent along the boundary normal so that sigma keeps the meaning it has @@ -714,6 +715,9 @@ PeleC::nscbc_params(const int idir) // probhi(idir) and was therefore silently wrong for any domain not // anchored at the origin. const auto& geom = amrex::DefaultGeometry(); + p.extrap_temperature = bc_nscbc_extrap_temperature; + p.extrap_material = bc_nscbc_extrap_material; + p.backflow_material = bc_nscbc_backflow_material; p.L_ref = geom.ProbHi(idir) - geom.ProbLo(idir); return p; } @@ -742,7 +746,10 @@ PeleC::nscbc_report_diagnostics() h[pc_nscbc::Diag::eos_failure] + h[pc_nscbc::Diag::floored] + h[pc_nscbc::Diag::transverse_drop] + h[pc_nscbc::Diag::source_drop] + h[pc_nscbc::Diag::target_incomplete]; - if (amrex::ParallelDescriptor::IOProcessor() && (total > 0 || verbose > 1)) { + const amrex::Long structure = h[pc_nscbc::Diag::structure]; + if ( + amrex::ParallelDescriptor::IOProcessor() && + (total > 0 || structure > 0 || verbose > 1)) { amrex::Print() << " NSCBC fallbacks since last report:" << " supersonic " << h[pc_nscbc::Diag::supersonic] << ", flow reversal " << h[pc_nscbc::Diag::reversed] << ", EB body state " @@ -752,6 +759,19 @@ PeleC::nscbc_report_diagnostics() << h[pc_nscbc::Diag::transverse_drop] << ", source dropped " << h[pc_nscbc::Diag::source_drop] << ", target incomplete " << h[pc_nscbc::Diag::target_incomplete] << "\n"; + if (structure > 0) { + // Advisory, not a fallback: a front is in the outflow boundary cells, + // which is the configuration the flame closures exist for. + amrex::Print() + << " NSCBC: material structure (|dS| > 5% of rho per cell) sat in " + << structure << " outflow boundary-cell fills since the last report.\n" + << " A flame or front is on this outflow: set " + "bc_nscbc_extrap_temperature = 1 and bc_nscbc_beta_s = 0 (any " + "sigma then\n" + << " survives a crossing), and keep bc_nscbc_extrap_material " + "off during a transit. See the BCs chapter and " + "NSCBC-FlameOutflow/README.md.\n"; + } } // Settle any counter atomics still in flight on other streams before the // reset; the blocking Gpu::copy above synchronised only its own stream. diff --git a/Source/NSCBC.H b/Source/NSCBC.H index 4f187cad0..3f3f194c6 100644 --- a/Source/NSCBC.H +++ b/Source/NSCBC.H @@ -105,6 +105,66 @@ struct Params // than an absent one: start at 0.5, raise toward 1 for grazing incidence. amrex::Real beta = 1.0; + // Weight on the REACTION source term in the modelled incoming acoustic, + // with the same convention as beta: + // L_in = K (phi - phi_target) - (1 - beta) T + (1 - beta_s) S_reaction + // 1 discards it (the default; correct when no heat release reaches the + // boundary cells), 0 includes it fully. + // + // A flame in the boundary cells raises p at S_p = dp/dt|_react, which the + // relaxation has no model for; without this term the mean sits off target + // by S_p/K. Including it is the Sutherland-Kennedy (JCP 2003) + // cancellation; the sign is PLUS because the exact steady flame has + // dR_-/dn = +S_p/(rho c^2) (convention-free check: + // Verification/NSCBC1D/source_sign_check.py). Guidance, measured in + // NSCBC-FlameOutflow/README.md: sitting flame -- beta_s = 0 halves the + // sigma = 1 error and with extrap_temperature holds hard-outflow level at + // every sigma (at sigma = 16 it overshoots; do not stack both). Crossing + // flame -- with extrap_temperature it cuts the transit disturbance 7-10x + // and lets sigma < 16 survive. Not obtainable from srcq(QPRES); see + // reaction_dpdt(). + amrex::Real beta_s = 1.0; + + // Close the lambda_0 family on TEMPERATURE instead of the entropy + // invariant. The diffusion operator forms boundary-face fluxes from these + // ghost cells, so the ghost dT/dn IS the heat flux leaving; the entropy + // closure sets T as a by-product of the acoustic algebra and overstates a + // flame-like face gradient by 40% (check C8). With this set, T rides the + // same minmod slope as everything else and rho follows from the EOS -- a + // deliberate Neumann condition on T instead of an accidental Dirichlet one. + // Default false (changes every outflow result). Turn it on whenever + // thermal or compositional structure nears the boundary: it is the flag + // that lets a flame survive the outflow (NSCBC-FlameOutflow/README.md). + bool extrap_temperature = false; + + // Continue the MATERIAL part of R_- into the ghost at an outflow, so a + // dilatational gradient crossing the face is not converted into ghost + // pressure (the C9(a) bias, 1/2 rho c l du per layer). The slope cannot be + // read from R_- itself (its gradient contains the boundary's own incoming + // waves; extrapolating them back is positive feedback, 16x on the rate, + // C5), so it is bounded through the acoustically-blind entropy family: + // du_mat = -u_out dS / (rho (1 - M^2)), + // applied slope = minmod(measured dR_-, du_mat (1 + M)). + // For pure acoustics the bound vanishes and the continuation shuts off. + // Quasi-steady bound: conservative for a moving front; keep sigma > 0. + // Default false. For structure that SITS near the outflow; leave it off + // when a front crosses -- a transit over-vents (NSCBC-FlameOutflow/ + // README.md). pin_farfield ignores it. + bool extrap_material = false; + + // Treat FIRM reversal at an outflow as a local inflow of reservoir gas: + // under sustained backflow the lambda_0 ghost content ramps from the + // frozen interior values (the transient-breathing closure) to the ambient + // state the face's Target carries (tgt.T, tgt.Y, tangential tgt.u), over + // an outward-Mach band [1e-3, 1e-2] -- continuous with w_mat's band below + // it, fully reservoir by Mach 0.01. Without this, a face in sustained + // recirculation feeds the domain its own exhaust: a hot domain drawing + // from a cold reservoir never cools (driver C14, the flush test). Only + // active when the Target supplies a physical state (tgt.T > 0); breathing + // reversals (chamber ring-down, radial rarefaction: |M| << 1e-3 at the + // face) never leave the frozen closure. Default false. + bool backflow_material = false; + // Boundary-normal reference length, probhi[idir] - problo[idir]. Filled per // face by the caller. Only the ratio sigma/L_ref is physical, so this is // not a user knob; it exists so that sigma keeps its literature meaning. @@ -139,6 +199,10 @@ struct Params why = "bc_nscbc_relax_t must be >= 0"; return false; } + if (beta_s > 1.0) { + why = "bc_nscbc_beta_s must be <= 1"; + return false; + } if (beta > 1.0) { why = "bc_nscbc_beta must be <= 1 (1 discards the transverse terms, 0 " "includes them fully, negative selects the local Mach number)"; @@ -331,6 +395,84 @@ cell_primitives(const amrex::Real s[NVAR]) noexcept return q; } +// --------------------------------------------------------------------------- +// reaction_dpdt -- dp/dt from chemistry alone at fixed (rho, e). Chemistry +// conserves both, so the effect is purely compositional: +// dp/dt = sum_k (dp/dY_k)|_{rho,e} wdot_k / rho +// = Ru T sum_k wdot_k/W_k - p/(rho T c_v) sum_k e_k wdot_k +// (ideal-gas closed form, gated by C7). SRK/manifold lack the composition +// derivative, so there a single directional finite difference along the +// reaction path is used (sum wdot = 0 keeps sum Y = 1 along it). +// --------------------------------------------------------------------------- +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real +reaction_dpdt(CellPrim const& q, bool& ok) noexcept +{ + auto eos = pele::physics::PhysicsType::eos(); + ok = true; + + amrex::Real wdot[NUM_SPECIES] = {0.0}; // GCC -Wmaybe-uninitialized at -O: + // it cannot see RTY2WDOT fill it + eos.RTY2WDOT(q.rho, q.T, q.Y, wdot); + amrex::Real wmax = 0.0; + for (const amrex::Real wn : wdot) { + if (!amrex::Math::isfinite(wn)) { + ok = false; + return 0.0; + } + wmax = amrex::max(wmax, std::abs(wn)); + } + if (wmax <= 0.0) { + return 0.0; // frozen chemistry here; nothing to correct + } + +#if defined(USE_SRK_EOS) || defined(USE_MANIFOLD_EOS) + // Directional finite difference along the reaction path. + const amrex::Real dY_scale = wmax / q.rho; // [1/s] + const amrex::Real tau = 1.0e-6 / dY_scale; // gives max |dY| ~ 1e-6 + amrex::Real Yp[NUM_SPECIES], sum = 0.0; + for (int n = 0; n < NUM_SPECIES; n++) { + Yp[n] = amrex::max(q.Y[n] + tau * wdot[n] / q.rho, 0.0); + sum += Yp[n]; + } + if (!(sum > 0.0)) { + ok = false; + return 0.0; + } + for (amrex::Real& yp : Yp) { + yp /= sum; + } + amrex::Real Tp = q.T, pp = 0.0; + eos.REY2T(q.rho, q.e, Yp, Tp); + eos.RTY2P(q.rho, Tp, Yp, pp); + if (!amrex::Math::isfinite(pp) || !amrex::Math::isfinite(Tp) || !(Tp > 0.0)) { + ok = false; + return 0.0; + } + return (pp - q.p) / tau; +#else + amrex::Real imw[NUM_SPECIES], ei[NUM_SPECIES], cv = 0.0; + eos.inv_molecular_weight(imw); + eos.RTY2Ei(q.rho, q.T, q.Y, ei); + eos.RTY2Cv(q.rho, q.T, q.Y, cv); + if (!amrex::Math::isfinite(cv) || !(cv > 0.0)) { + ok = false; + return 0.0; + } + amrex::Real mole_term = 0.0, heat_term = 0.0; + for (int n = 0; n < NUM_SPECIES; n++) { + mole_term += wdot[n] * imw[n]; + heat_term += ei[n] * wdot[n]; + } + const amrex::Real dpdt = pele::physics::Constants::RU * q.T * mole_term - + q.p * heat_term / (q.rho * q.T * cv); + if (!amrex::Math::isfinite(dpdt)) { + ok = false; + return 0.0; + } + return dpdt; +#endif +} + // --------------------------------------------------------------------------- // Pack a set of ghost primitives into a conserved state. // @@ -630,7 +772,43 @@ apply( } } - const amrex::Real dRm = 0.0; + // Entropy-family slope at the face: measured for two purposes below. dS + // carries no acoustic content and is an OUTGOING family at an outflow, so + // nothing the incoming model consumes here can feed on the waves it + // launches. Measured on the same stencil, in the same frozen-impedance + // variables, and with the same limiter as R_+. + amrex::Real dRm = 0.0; + if (outflow_face && linear && qNm1.ok && qNm2.ok) { + const amrex::Real inv_c2 = 1.0 / (c_N * c_N); + const amrex::Real S_N = qN.rho - qN.p * inv_c2; + const amrex::Real S_Nm1 = qNm1.rho - qNm1.p * inv_c2; + const amrex::Real S_Nm2 = qNm2.rho - qNm2.p * inv_c2; + const amrex::Real dS = minmod(S_N - S_Nm1, S_Nm1 - S_Nm2); + + // The transit guard: count material structure in the outflow boundary + // cell so nscbc_report_diagnostics() can advise the flame closures + // (extrap_temperature, beta_s). 5% of rho per cell is well above + // acoustic or roundoff content and well below any front. + if ((layer == 1) && (std::abs(dS) > 0.05 * qN.rho)) { + bump(Diag::structure); + } + + // Material-slope continuation of R_- (see Params::extrap_material): the + // measured slope of R_- is used only within the bound the entropy family + // supplies through steady continuity and the momentum-flux pressure + // gradient. + if (prm.extrap_material) { + const amrex::Real Rm_Nm1 = n_sgn * qNm1.u[idir] - qNm1.p / rho_c; + const amrex::Real Rm_Nm2 = n_sgn * qNm2.u[idir] - qNm2.p / rho_c; + const amrex::Real dRm_meas = minmod(Rm_N - Rm_Nm1, Rm_Nm1 - Rm_Nm2); + const amrex::Real du_mat = -u_out_N * dS / (qN.rho * one_m_M2); + dRm = w_mat * minmod(dRm_meas, du_mat * (1.0 + mach)); + if (!amrex::Math::isfinite(dRm)) { + dRm = 0.0; + bump(Diag::floored); + } + } + } amrex::Real Rm_g; if (outflow_face && prm.pin_farfield) { @@ -657,6 +835,15 @@ apply( L_in += 2.0 * rho_c * (n_sgn * tgt.dudt); } L_in -= (1.0 - beta) * T_in; + if (prm.beta_s < 1.0) { + bool src_ok = true; + const amrex::Real S_p = reaction_dpdt(qN, src_ok); + if (src_ok) { + L_in += (1.0 - prm.beta_s) * S_p; + } else { + bump(Diag::source_drop); + } + } // There is deliberately NO diffusive source term on L_in. The // ghost-cell form does not have the flux-form's viscous-condition gap: // the diffusion operator reads these ghost cells, so a correct ghost @@ -723,6 +910,12 @@ apply( const amrex::Real S_g = S_N + fl * dS; rho_g = S_g + p_g * inv_c2; + // Temperature closure: T on the same limited slope; rho from the EOS. + amrex::Real dT = 0.0; + if (prm.extrap_temperature && linear && qNm1.ok && qNm2.ok) { + dT = w_mat * minmod(qN.T - qNm1.T, qNm1.T - qNm2.T); + } + amrex::Real Ysum = 0.0; for (int n = 0; n < NUM_SPECIES; n++) { Y_g[n] = amrex::max(qN.Y[n] + fl * dY[n], 0.0); @@ -746,15 +939,66 @@ apply( } } - const amrex::Real rho_floor = - amrex::max(1.0e-6 * qN.rho, constants::very_small_num()); - if (!amrex::Math::isfinite(rho_g) || (rho_g < rho_floor)) { - rho_g = rho_floor; - bump(Diag::floored); + if (prm.extrap_temperature) { + T_g = qN.T + fl * dT; + // A temperature floor rather than a density one: rho is the derived + // quantity in this closure. + const amrex::Real T_floor = + amrex::max(1.0e-3 * qN.T, constants::very_small_num()); + if (!amrex::Math::isfinite(T_g) || (T_g < T_floor)) { + T_g = T_floor; + bump(Diag::floored); + } + eos.PYT2RE(p_g, Y_g, T_g, rho_g, e_g); + } else { + const amrex::Real rho_floor = + amrex::max(1.0e-6 * qN.rho, constants::very_small_num()); + if (!amrex::Math::isfinite(rho_g) || (rho_g < rho_floor)) { + rho_g = rho_floor; + bump(Diag::floored); + } + eos.RYP2T(rho_g, Y_g, p_g, T_g); + eos.RTY2E(rho_g, T_g, Y_g, e_g); } - eos.RYP2T(rho_g, Y_g, p_g, T_g); - eos.RTY2E(rho_g, T_g, Y_g, e_g); + // Sustained-backflow material (see Params::backflow_material): under + // firm reversal the lambda_0 content ramps from the frozen interior + // values just computed to the reservoir state the Target carries. The + // ramp starts where w_mat's band ends (outward Mach -1e-3) and is + // complete by Mach -1e-2, so breathing reversals never feel it, and the + // acoustic side is untouched -- p_g stays the relaxation's. + if (prm.backflow_material && (u_out_N < 0.0) && (tgt.T > 0.0)) { + const amrex::Real m_lo = 1.0e-3, m_hi = 1.0e-2; + const amrex::Real w_rev = amrex::min( + amrex::max((-mach - m_lo) / (m_hi - m_lo), 0.0), 1.0); + if (w_rev > 0.0) { + amrex::Real Ysum_rev = 0.0; + for (int n = 0; n < NUM_SPECIES; n++) { + Y_g[n] = (1.0 - w_rev) * Y_g[n] + + w_rev * amrex::max(tgt.Y[n], 0.0); + Ysum_rev += Y_g[n]; + } + if (Ysum_rev > constants::very_small_num()) { + const amrex::Real inv = 1.0 / Ysum_rev; + for (amrex::Real& yg : Y_g) { + yg *= inv; + } + } + T_g = (1.0 - w_rev) * T_g + w_rev * tgt.T; + for (int d = 0; d < 3; d++) { + if (d != idir) { + u_g[d] = (1.0 - w_rev) * u_g[d] + w_rev * tgt.u[d]; + } + } + const amrex::Real T_floor = + amrex::max(1.0e-3 * qN.T, constants::very_small_num()); + if (!amrex::Math::isfinite(T_g) || (T_g < T_floor)) { + T_g = T_floor; + bump(Diag::floored); + } + eos.PYT2RE(p_g, Y_g, T_g, rho_g, e_g); + } + } } else { // Inflow. One dimensionless nudge factor shared by temperature and the // tangential velocities -- both ride lambda_0, magnitude |u_out|: diff --git a/Source/Params/_cpp_parameters b/Source/Params/_cpp_parameters index 3c1397804..482dfd5dc 100644 --- a/Source/Params/_cpp_parameters +++ b/Source/Params/_cpp_parameters @@ -116,6 +116,22 @@ bc_nscbc_relax_t Real 0.2 # because a wrong beta is far more dangerous than an absent one. Must be <= 1. bc_nscbc_beta Real 1.0 +# Weight on the REACTION source term in the modelled incoming acoustic, same +# convention as bc_nscbc_beta: 1 discards it (default), 0 includes it fully. +# +# A flame near an outflow raises the pressure in the boundary cell at a rate +# the 1-D relaxation has no model for, so sigma has to absorb it and the mean +# pressure sits off target. This supplies the missing dp/dt. Only worth +# enabling when the reaction zone is close enough to the boundary to matter -- +# the correct first answer to that situation is to move the boundary. +# +# Chemistry conserves mass, and PeleC's constant-volume reactor conserves total +# internal energy because the formation enthalpies live inside e, so the entire +# effect is compositional: dp/dY|_{rho,e}. Exact in closed form for ideal-gas +# EOS; a directional finite difference along the reaction path is used for SRK. +# Must be <= 1. +bc_nscbc_beta_s Real 1.0 + # Order of the outgoing-invariant extrapolation, 1 or 2. Verification and # debugging knob, not a physics knob; leave at 2. bc_nscbc_order int 2 @@ -129,6 +145,39 @@ bc_nscbc_order int 2 # pressure is not the target. bc_nscbc_sigma is ignored when this is set. bc_nscbc_pin_farfield bool false +# Close the lambda_0 family at an outflow on TEMPERATURE rather than on the +# linearised entropy invariant. Matters because PeleC's diffusion operator +# forms the boundary-face conductive and species fluxes from the NSCBC ghost +# cells, and with the entropy closure the ghost temperature is a by-product of +# the acoustic algebra: check C8 in Verification/NSCBC1D measures a 40% +# overstatement of the face temperature gradient on a flame-like ramp. Turn it +# on when a thermal or compositional structure is near the boundary. Default +# false because it changes every outflow result. +bc_nscbc_extrap_temperature bool false + +# Continue the material (non-acoustic) part of the incoming invariant R_- into +# the ghost at an outflow, in addition to the relaxation increment. The +# material slope is bounded through the entropy family (outgoing, carries no +# acoustics), which cancels the ghost-pressure bias of check C9(a) in +# Verification/NSCBC1D without touching reflection, anchoring or the +# relaxation rate. For a front SITTING on the outflow it is worth 5% alone +# and 42% together with bc_nscbc_extrap_temperature; for a front PASSING OUT +# through the face it is a liability -- at small sigma the quasi-steady model +# turns the crossing into a runaway (see NSCBC-FlameOutflow/README.md). Keep +# sigma > 0 with this on. Default false. Ignored by pin_farfield. +bc_nscbc_extrap_material bool false + +# Treat FIRM reversal at an outflow as a local inflow of reservoir gas: under +# sustained backflow the lambda_0 ghost content ramps from the frozen interior +# values to the ambient state the face's bcnormal_nscbc Target carries (T, Y, +# tangential u), over an outward-Mach band [1e-3, 1e-2]. Without it, a face +# held in recirculation feeds the domain its own exhaust (driver check C14, +# the flush test); with it, the entering gas carries the reservoir state. +# Breathing reversals (|M| << 1e-3 at the face) never leave the frozen +# closure either way. Requires the problem's outflow Target to supply T > 0 +# (and Y); inert otherwise. Default false. +bc_nscbc_backflow_material bool false + # NOTE: pelec.nscbc_adv and pelec.nscbc_diff were removed here. They were # declared and queried but never read by any live code path -- the diff --git a/Source/Params/param_includes/pelec_defaults.H b/Source/Params/param_includes/pelec_defaults.H index efa1a41fc..b1b8e7331 100644 --- a/Source/Params/param_includes/pelec_defaults.H +++ b/Source/Params/param_includes/pelec_defaults.H @@ -23,8 +23,12 @@ amrex::Real PeleC::bc_nscbc_sigma = 0.25; amrex::Real PeleC::bc_nscbc_relax_u = 2.0; amrex::Real PeleC::bc_nscbc_relax_t = 0.2; amrex::Real PeleC::bc_nscbc_beta = 1.0; +amrex::Real PeleC::bc_nscbc_beta_s = 1.0; int PeleC::bc_nscbc_order = 2; bool PeleC::bc_nscbc_pin_farfield = false; +bool PeleC::bc_nscbc_extrap_temperature = false; +bool PeleC::bc_nscbc_extrap_material = false; +bool PeleC::bc_nscbc_backflow_material = false; bool PeleC::add_ext_src = false; amrex::GpuArray PeleC::external_forcing = {0.0}; bool PeleC::add_forcing_src = false; diff --git a/Source/Params/param_includes/pelec_params.H b/Source/Params/param_includes/pelec_params.H index c5c6f80e7..338fe2a2d 100644 --- a/Source/Params/param_includes/pelec_params.H +++ b/Source/Params/param_includes/pelec_params.H @@ -23,8 +23,12 @@ static amrex::Real bc_nscbc_sigma; static amrex::Real bc_nscbc_relax_u; static amrex::Real bc_nscbc_relax_t; static amrex::Real bc_nscbc_beta; +static amrex::Real bc_nscbc_beta_s; static int bc_nscbc_order; static bool bc_nscbc_pin_farfield; +static bool bc_nscbc_extrap_temperature; +static bool bc_nscbc_extrap_material; +static bool bc_nscbc_backflow_material; static bool add_ext_src; static amrex::GpuArray external_forcing; static bool add_forcing_src; diff --git a/Source/Params/param_includes/pelec_queries.H b/Source/Params/param_includes/pelec_queries.H index b381f61d0..8884f62a0 100644 --- a/Source/Params/param_includes/pelec_queries.H +++ b/Source/Params/param_includes/pelec_queries.H @@ -23,8 +23,12 @@ pp.query("bc_nscbc_sigma", bc_nscbc_sigma); pp.query("bc_nscbc_relax_u", bc_nscbc_relax_u); pp.query("bc_nscbc_relax_t", bc_nscbc_relax_t); pp.query("bc_nscbc_beta", bc_nscbc_beta); +pp.query("bc_nscbc_beta_s", bc_nscbc_beta_s); pp.query("bc_nscbc_order", bc_nscbc_order); pp.query("bc_nscbc_pin_farfield", bc_nscbc_pin_farfield); +pp.query("bc_nscbc_extrap_temperature", bc_nscbc_extrap_temperature); +pp.query("bc_nscbc_extrap_material", bc_nscbc_extrap_material); +pp.query("bc_nscbc_backflow_material", bc_nscbc_backflow_material); pp.query("add_ext_src", add_ext_src); { amrex::Vector tmp(AMREX_SPACEDIM, 0.0);