From 19021dfa05bb38f0e0690e7958721f9afa71d340 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Fri, 21 Aug 2026 11:15:07 -0700 Subject: [PATCH 1/2] fix issues 171 176 181 188 189 191 --- ...ncflo_compute_MAC_projected_velocities.cpp | 2 +- src/embedded_boundaries/eb_twocylinders.cpp | 19 ++++++---- src/incflo.H | 7 ++++ src/incflo.cpp | 18 +++++---- src/incflo_regrid.cpp | 21 ++++++++++ src/particles/incflo_PCInit.cpp | 20 ++++++++-- src/projection/incflo_apply_cc_projection.cpp | 38 ++++++++++++++----- src/setup/init.cpp | 6 +++ 8 files changed, 102 insertions(+), 29 deletions(-) diff --git a/src/convection/incflo_compute_MAC_projected_velocities.cpp b/src/convection/incflo_compute_MAC_projected_velocities.cpp index 1919d979..a40bb292 100644 --- a/src/convection/incflo_compute_MAC_projected_velocities.cpp +++ b/src/convection/incflo_compute_MAC_projected_velocities.cpp @@ -93,7 +93,7 @@ incflo::compute_MAC_projected_velocities ( // // Initialize (or redefine the beta in) the MacProjector // - if (macproj->needInitialization()) + if (get_mac_projector()->needInitialization()) { LPInfo lp_info; lp_info.setMaxCoarseningLevel(m_mac_mg_max_coarsening_level); diff --git a/src/embedded_boundaries/eb_twocylinders.cpp b/src/embedded_boundaries/eb_twocylinders.cpp index 07bb3413..0665d4c9 100644 --- a/src/embedded_boundaries/eb_twocylinders.cpp +++ b/src/embedded_boundaries/eb_twocylinders.cpp @@ -61,14 +61,19 @@ void incflo::make_eb_twocylinders() // Build the implicit function as a union of two cylinders EB2::CylinderIF cyl1(radius1, direction1, center1, false); EB2::CylinderIF cyl2(radius2, direction2, center2, false); - auto twocylinders = inside ? EB2::makeComplement(EB2::makeUnion(cyl1, cyl2)) - : EB2::makeUnion(cyl1, cyl2); - - // Generate GeometryShop - auto gshop = EB2::makeShop(twocylinders); - // Build index space int max_level_here = 0; int max_coarsening_level = 100; - EB2::Build(gshop, geom.back(), max_level_here, max_level_here + max_coarsening_level); + + // NOTE: this must not be written as a ternary -- a conditional expression has + // a single type, and ComplementIF's converting constructor is not explicit, so + // the "else" branch would be silently wrapped in a complement as well, giving + // internal-flow geometry for both settings. + if (inside) { + auto gshop = EB2::makeShop(EB2::makeComplement(EB2::makeUnion(cyl1, cyl2))); + EB2::Build(gshop, geom.back(), max_level_here, max_level_here + max_coarsening_level); + } else { + auto gshop = EB2::makeShop(EB2::makeUnion(cyl1, cyl2)); + EB2::Build(gshop, geom.back(), max_level_here, max_level_here + max_coarsening_level); + } } diff --git a/src/incflo.H b/src/incflo.H index a4ff18f4..a2a96850 100644 --- a/src/incflo.H +++ b/src/incflo.H @@ -462,6 +462,13 @@ private: std::unique_ptr macproj; + // Lazily (re)builds macproj if it is null. ClearLevel resets it, and + // AmrCore::regrid runs ClearLevel *after* the RemakeLevel / + // MakeNewLevelFromCoarse calls, so a regrid that lowers finest_level would + // otherwise leave macproj null until the next projection dereferenced it. + // Same lazy pattern as get_diffusion_tensor_op / get_diffusion_scalar_op. + Hydro::MacProjector* get_mac_projector (); + int m_mac_mg_max_coarsening_level = 100; #ifdef AMREX_USE_FLOAT diff --git a/src/incflo.cpp b/src/incflo.cpp index 78e4749e..8122dff1 100644 --- a/src/incflo.cpp +++ b/src/incflo.cpp @@ -89,11 +89,6 @@ void incflo::InitData () WriteSmallPlotFile(); m_last_smallplt = 0; } - if (m_KE_int > 0) - { - amrex::Abort("xxxxx m_KE_int todo"); -// amrex::Print() << "Time, Kinetic Energy: " << m_cur_time << ", " << ComputeKineticEnergy() << "\n"; - } } else { @@ -134,7 +129,11 @@ void incflo::Evolve() { BL_PROFILE("incflo::Evolve()"); - bool do_not_evolve = ((m_max_step == 0) || ((m_stop_time >= 0.) && (m_cur_time > m_stop_time)) || + // The stop_time test here mirrors the loop-exit test below, tolerance included, so + // that restarting from the final checkpoint of a completed run does not take an + // extra step past stop_time. + bool do_not_evolve = ((m_max_step == 0) || + ((m_stop_time > 0.) && (m_cur_time >= m_stop_time - (1.e-12 * m_dt))) || ((m_stop_time <= 0.) && (m_max_step <= 0)) || (m_max_step >= 0 && m_nstep >= m_max_step) ) && !m_steady_state; @@ -302,8 +301,11 @@ incflo::writeNow(int a_plot_int, Real a_plot_per_approx, Real a_plot_per_exact) // the number of intervals that have elapsed for both the current // time and the time at the beginning of this timestep. - int num_per_old = static_cast(std::round((m_cur_time-m_dt) / a_plot_per_approx)); - int num_per_new = static_cast(std::round((m_cur_time ) / a_plot_per_approx)); + // Note: this must truncate, not round -- rounding would detect a crossing at + // half-multiples of the period (plotting up to half a period early) and would + // also make the two epsilon corrections below unreachable. + int num_per_old = static_cast(std::floor((m_cur_time-m_dt) / a_plot_per_approx)); + int num_per_new = static_cast(std::floor((m_cur_time ) / a_plot_per_approx)); // Before using these, however, we must test for the case where we're // within machine epsilon of the next interval. In that case, increment diff --git a/src/incflo_regrid.cpp b/src/incflo_regrid.cpp index c8a8a867..11ef44c7 100644 --- a/src/incflo_regrid.cpp +++ b/src/incflo_regrid.cpp @@ -126,6 +126,27 @@ void incflo::RemakeLevel (int lev, Real time, const BoxArray& ba, #endif } +// Rebuild macproj on demand. Note that this must not be done inside ClearLevel: +// AmrCore::regrid assigns finest_level = new_finest only *after* its ClearLevel +// loop, so finest_level is still stale there and a multi-level shrink would build +// a projector with too many levels. Deferring to the first use gets the final +// finest_level, and also covers the regrid_on_restart path. +Hydro::MacProjector* +incflo::get_mac_projector () +{ + if (!macproj) { +#ifdef AMREX_USE_EB + macproj = std::make_unique(Geom(0,finest_level), + MLMG::Location::FaceCentroid, // Location of mac_vec + MLMG::Location::FaceCentroid, // Location of beta + MLMG::Location::CellCenter ); // Location of solution variable phi +#else + macproj = std::make_unique(Geom(0,finest_level)); +#endif + } + return macproj.get(); +} + // Delete level data // overrides the pure virtual function in AmrCore void incflo::ClearLevel (int lev) diff --git a/src/particles/incflo_PCInit.cpp b/src/particles/incflo_PCInit.cpp index f1561b4d..1a28141e 100644 --- a/src/particles/incflo_PCInit.cpp +++ b/src/particles/incflo_PCInit.cpp @@ -113,11 +113,9 @@ void incflo_PC::initializeParticlesUniformDistributionInBox ( const RealBox& par #ifdef AMREX_USE_EB if (vf_arr(i,j,k) > 0.) { #endif - if (i%2 == 0 && j%2 == 1) { if (particle_init_domain.contains(RealVect(AMREX_D_DECL(x,y,z)))) { num_particles_arr(i,j,k) = particles_per_cell; } - } #ifdef AMREX_USE_EB } #endif @@ -236,6 +234,15 @@ void incflo_PC::initializeParticlesUniformDistributionInBox ( const RealBox& par pp.get("center",cyl_center); Real x_ctr = cyl_center[0]; Real y_ctr = cyl_center[1]; +#if (AMREX_SPACEDIM == 3) + Real z_ctr = cyl_center[2]; +#endif + + // The cylinder axis matters here exactly as it does in AdvectWithFlow; + // assuming a z-parallel axis culls against the wrong axis for direction 0/1. + int cyl_direction; + pp.get("direction",cyl_direction); + AMREX_ALWAYS_ASSERT(cyl_direction >= 0 && cyl_direction < AMREX_SPACEDIM); // Remove particles that are outside of the cylinder for (ParIterType pti(*this, lev); pti.isValid(); ++pti) @@ -251,8 +258,15 @@ void incflo_PC::initializeParticlesUniformDistributionInBox ( const RealBox& par Real x = p.pos(0) - x_ctr; Real y = p.pos(1) - y_ctr; - +#if (AMREX_SPACEDIM == 3) + Real z = p.pos(2) - z_ctr; + Real r = (cyl_direction == 0) ? std::sqrt(y*y + z*z) + : (cyl_direction == 1) ? std::sqrt(x*x + z*z) + : std::sqrt(x*x + y*y); +#else + // In 2D the only cylinder axis is the out-of-plane one Real r = std::sqrt(x*x + y*y); +#endif if (r > cyl_radius) { p.id() = -1; diff --git a/src/projection/incflo_apply_cc_projection.cpp b/src/projection/incflo_apply_cc_projection.cpp index 7b5aa8cb..5503822c 100644 --- a/src/projection/incflo_apply_cc_projection.cpp +++ b/src/projection/incflo_apply_cc_projection.cpp @@ -240,15 +240,6 @@ void incflo::ApplyCCProjection (Vector density, } } - // Define "vel" to be U^* - U^n rather than U^* - if (proj_for_small_dt || incremental) - { - for (int lev = 0; lev <= finest_level; ++lev) { - MultiFab::Subtract(m_leveldata[lev]->velocity, - m_leveldata[lev]->velocity_o, 0, 0, AMREX_SPACEDIM, 0); - } - } - auto bclo = get_projection_bc(Orientation::low); auto bchi = get_projection_bc(Orientation::high); @@ -258,6 +249,33 @@ void incflo::ApplyCCProjection (Vector density, fillpatch_velocity(lev, m_t_new[lev], *vel[lev], 1); } + // Define "vel" to be U^* - U^n rather than U^*. + // + // This must happen *after* the fillpatch above, and must include the ghost + // cells: the field being projected is a difference, so its inflow ghosts must + // not carry the full physical inflow value (U^n already carries that). The + // nodal twin achieves this by suppressing the inflow fill altogether + // (set_inflow_bc = !proj_for_small_dt && !incremental). Here we instead fill + // U^n's ghosts the same way as U^*, at t_old, and subtract including one ghost, + // so the inflow ghosts of the difference are the inflow *increment* -- zero for + // steady inflow. Otherwise MOL::ExtrapVelToFacesBox would copy the full + // u_inflow onto the inflow faces and the MAC solve would see a spurious + // O(u_inflow) divergence source. + if (proj_for_small_dt || incremental) + { + for (int lev = 0; lev <= finest_level; ++lev) { + // A temporary, so that fillpatch_velocity is not writing into one of + // its own source MultiFabs (it interpolates between velocity_o and + // velocity), and so velocity_o itself is left untouched. + MultiFab vel_old(m_leveldata[lev]->velocity_o.boxArray(), + m_leveldata[lev]->velocity_o.DistributionMap(), + AMREX_SPACEDIM, 1, MFInfo(), Factory(lev)); + fillpatch_velocity(lev, m_t_old[lev], vel_old, 1); + MultiFab::Subtract(m_leveldata[lev]->velocity, vel_old, + 0, 0, AMREX_SPACEDIM, 1); + } + } + // *************************************************************************************** // START OF MAC STUFF // *************************************************************************************** @@ -293,7 +311,7 @@ void incflo::ApplyCCProjection (Vector density, // // Initialize (or redefine the beta in) the MacProjector - if (macproj->needInitialization()) + if (get_mac_projector()->needInitialization()) { LPInfo lp_info; lp_info.setMaxCoarseningLevel(m_mac_mg_max_coarsening_level); diff --git a/src/setup/init.cpp b/src/setup/init.cpp index dcc18a7d..f2536970 100644 --- a/src/setup/init.cpp +++ b/src/setup/init.cpp @@ -27,6 +27,12 @@ void incflo::ReadParameters () pp.query("refine_particles", m_refine_particles); #endif pp.query("KE_int", m_KE_int); + if (m_KE_int > 0) { + // ComputeKineticEnergy() is not implemented (its body is #if 0'd and it + // returns 0), so refuse the option here rather than printing a zero as if + // it were a real diagnostic. Checked here so that restarts refuse it too. + amrex::Abort("incflo.KE_int > 0: ComputeKineticEnergy() is not implemented yet"); + } } // end prefix amr From ea9c98e47f631c5c62ad71cfdb24c9accc0b1589 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Fri, 21 Aug 2026 11:20:07 -0700 Subject: [PATCH 2/2] update to previous commit --- src/projection/incflo_apply_cc_projection.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/projection/incflo_apply_cc_projection.cpp b/src/projection/incflo_apply_cc_projection.cpp index 5503822c..8b29c463 100644 --- a/src/projection/incflo_apply_cc_projection.cpp +++ b/src/projection/incflo_apply_cc_projection.cpp @@ -311,7 +311,7 @@ void incflo::ApplyCCProjection (Vector density, // // Initialize (or redefine the beta in) the MacProjector - if (get_mac_projector()->needInitialization()) + if (macproj->needInitialization()) { LPInfo lp_info; lp_info.setMaxCoarseningLevel(m_mac_mg_max_coarsening_level);