fix issues 171 176 181 188 189 191

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);

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);
+    }
 }

src/incflo.H

@@ -462,6 +462,13 @@ private:
 
     std::unique_ptr<Hydro::MacProjector> 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

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<int>(std::round((m_cur_time-m_dt) / a_plot_per_approx));
-        int num_per_new = static_cast<int>(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<int>(std::floor((m_cur_time-m_dt) / a_plot_per_approx));
+        int num_per_new = static_cast<int>(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

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<Hydro::MacProjector>(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<Hydro::MacProjector>(Geom(0,finest_level));
+#endif
+    }
+    return macproj.get();
+}
+
 // Delete level data
 // overrides the pure virtual function in AmrCore
 void incflo::ClearLevel (int lev)

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;

src/projection/incflo_apply_cc_projection.cpp

@@ -240,15 +240,6 @@ void incflo::ApplyCCProjection (Vector<MultiFab const*> 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<MultiFab const*> 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
     // ***************************************************************************************

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