Fix issues #233-#239

src/convection/incflo_compute_MAC_projected_velocities.cpp

@@ -208,6 +208,12 @@ incflo::compute_MAC_projected_velocities (
             velBC_MF = make_BC_MF(lev, m_bcrec_velocity_d, "velocity");
         }
 
+        // On a direction_dependent face IncfloVelFill writes the prescribed velocity
+        // into the ghost cell only where it points into the domain; on the outflow
+        // part it copies the interior cell, which is 0 here because
+        // time_dep_inflow_vel was zeroed above.  The inflow tests below must
+        // therefore be strict, otherwise the outflow part of the face is treated as
+        // zero inflow and the extrapolated outflow velocity is lost.
         for (MFIter mfi(time_dep_inflow_vel,false); mfi.isValid(); ++mfi)
         {
             Box const& bx = mfi.validbox();
@@ -222,11 +228,11 @@ incflo::compute_MAC_projected_velocities (
                 int n = 0;
                 const auto bc = HydroBC::getBC(i, j, k, n, domain, bc_vel_d, velbc_arr);
                 if (i == dlo.x && ( bc.lo(0) == BCType::ext_dir ||
-                                   (bc.lo(0) == BCType::direction_dependent && cc_arr(i-1,j,k,0) >= Real(0.0)) ) ) {
+                                   (bc.lo(0) == BCType::direction_dependent && cc_arr(i-1,j,k,0) >  Real(0.0)) ) ) {
                     umac_arr(i,j,k) = cc_arr(i-1,j,k,0);
                 }
                 if (i == dhi.x && ( bc.hi(0) == BCType::ext_dir ||
-                                   (bc.hi(0) == BCType::direction_dependent && cc_arr(i+1,j,k,0) <= Real(0.0)) ) ) {
+                                   (bc.hi(0) == BCType::direction_dependent && cc_arr(i+1,j,k,0) <  Real(0.0)) ) ) {
                     umac_arr(i+1,j,k) = cc_arr(i+1,j,k,0);
                 }
             });
@@ -236,11 +242,11 @@ incflo::compute_MAC_projected_velocities (
                 int n = 1;
                 const auto bc = HydroBC::getBC(i, j, k, n, domain, bc_vel_d, velbc_arr);
                 if (j == dlo.y && ( bc.lo(1) == BCType::ext_dir||
-                                   (bc.lo(1) == BCType::direction_dependent && cc_arr(i,j-1,k,1) >= Real(0.0)) ) ) {
+                                   (bc.lo(1) == BCType::direction_dependent && cc_arr(i,j-1,k,1) >  Real(0.0)) ) ) {
                     vmac_arr(i,j,k) = cc_arr(i,j-1,k,1);
                 }
                 if (j == dhi.y && ( bc.hi(1) == BCType::ext_dir ||
-                                   (bc.hi(1) == BCType::direction_dependent && cc_arr(i,j+1,k,1) <= Real(0.0)) ) ) {
+                                   (bc.hi(1) == BCType::direction_dependent && cc_arr(i,j+1,k,1) <  Real(0.0)) ) ) {
                     vmac_arr(i,j+1,k) = cc_arr(i,j+1,k,1);
                 }
             });
@@ -251,11 +257,11 @@ incflo::compute_MAC_projected_velocities (
                 int n = 2;
                 const auto bc = HydroBC::getBC(i, j, k, n, domain, bc_vel_d, velbc_arr);
                 if (k == dlo.z && ( bc.lo(2) == BCType::ext_dir ||
-                                   (bc.lo(2) == BCType::direction_dependent && cc_arr(i,j,k-1,2) >= Real(0.0)) ) ) {
+                                   (bc.lo(2) == BCType::direction_dependent && cc_arr(i,j,k-1,2) >  Real(0.0)) ) ) {
                     wmac_arr(i,j,k) = cc_arr(i,j,k-1,2);
                 }
                 if (k == dhi.z && ( bc.hi(2) == BCType::ext_dir ||
-                                   (bc.hi(2) == BCType::direction_dependent && cc_arr(i,j,k+1,2) <= Real(0.0)) ) ) {
+                                   (bc.hi(2) == BCType::direction_dependent && cc_arr(i,j,k+1,2) <  Real(0.0)) ) ) {
                     wmac_arr(i,j,k+1) = cc_arr(i,j,k+1,2);
                 }
             });

src/incflo.cpp

@@ -170,6 +170,13 @@ void incflo::Evolve()
         {
             if (m_verbose > 0) amrex::Print() << "Regridding...\n";
             regrid(0, m_cur_time);
+#ifdef INCFLO_USE_PARTICLES
+            // This must be done after regrid() returns: inside RemakeLevel /
+            // MakeNewLevelFromCoarse the ParGDB still sees the old BoxArray,
+            // DistributionMapping and finest_level, so a Redistribute there would
+            // map the particles onto the pre-regrid layout.
+            particleData.Redistribute();
+#endif
             if (m_verbose > 0 && ParallelDescriptor::IOProcessor()) {
                 printGridSummary(amrex::OutStream(), 0, finest_level);
             }

src/incflo_apply_predictor.cpp

@@ -229,7 +229,10 @@ void incflo::ApplyPredictor (bool incremental_projection)
     // **************************************************************************************
     // Update the particle positions
     // **************************************************************************************
-    if (m_advection_type != "MOL") {
+    // Not during the initial pressure iterations: the fields are reset after each
+    // of those, but the particle positions would not be, so the particles would
+    // start the run m_initial_iterations*dt ahead of the fluid.
+    if (m_advection_type != "MOL" && !incremental_projection) {
         evolveTracerParticles(AMREX_D_DECL(GetVecOfConstPtrs(u_mac), GetVecOfConstPtrs(v_mac),
                                            GetVecOfConstPtrs(w_mac)));
     }

src/incflo_compute_dt.cpp

@@ -52,6 +52,51 @@ void incflo::ComputeDt (int initialization, bool explicit_diffusion, double cur_
 
        compute_vel_forces_on_level (lev, vel_forces, vel, rho, tra_o, tra);
 
+       // Explicit-diffusion bound: the largest diffusivity that is applied
+       // explicitly, divided by rho.  We must use the strain-rate dependent
+       // viscosity here (not just the constant m_mu) and we must cover the tracer
+       // and temperature diffusivities as well, since those are updated explicitly
+       // with the same m_diff_type switch.  Covered cells are set to zero.
+       MultiFab nu;
+       if (explicit_diffusion) {
+           nu.define(grids[lev], dmap[lev], 1, 0, MFInfo(), Factory(lev));
+           compute_viscosity_at_level(lev, &nu, &m_leveldata[lev]->density,
+                                      &m_leveldata[lev]->velocity, geom[lev],
+                                      m_cur_time, 0);
+           Real mu_s_max = Real(0.0);
+           if (m_advect_tracer) {
+               for (int n = 0; n < m_ntrac; ++n) {
+                   mu_s_max = amrex::max(mu_s_max, m_mu_s[n]);
+               }
+           }
+           Real mu_T_eff = m_use_temperature ? m_mu_T / m_cp : Real(0.0);
+#ifdef AMREX_USE_EB
+           auto const& dt_flags = EBFactory(lev).getMultiEBCellFlagFab();
+#endif
+#ifdef _OPENMP
+#pragma omp parallel if (Gpu::notInLaunchRegion())
+#endif
+           for (MFIter mfi(nu,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+               Box const& bx = mfi.tilebox();
+               Array4<Real> const& nu_a = nu.array(mfi);
+               Array4<Real const> const& r = rho.const_array(mfi);
+#ifdef AMREX_USE_EB
+               Array4<EBCellFlag const> const& f = dt_flags.const_array(mfi);
+#endif
+               ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+               {
+#ifdef AMREX_USE_EB
+                   if (f(i,j,k).isCovered()) { nu_a(i,j,k) = Real(0.0); return; }
+#endif
+                   Real rinv = Real(1.0)/r(i,j,k);
+                   // Velocity and temperature diffuse with eta/rho and mu_T/(rho cp);
+                   // a non-conservative tracer diffuses with mu_s itself.
+                   nu_a(i,j,k) = amrex::max(amrex::max(nu_a(i,j,k), mu_T_eff)*rinv,
+                                            mu_s_max*amrex::max(Real(1.0), rinv));
+               });
+           }
+       }
+
 #ifdef AMREX_USE_EB
         if (!vel.isAllRegular()) {
             auto const& flag = EBFactory(lev).getMultiEBCellFlagFab();
@@ -72,22 +117,7 @@ void incflo::ComputeDt (int initialization, bool explicit_diffusion, double cur_
                            return mx;
                        });
             if (explicit_diffusion) {
-                diff_lev = amrex::ReduceMax(rho, flag, 0,
-                           [=] AMREX_GPU_HOST_DEVICE (Box const& b,
-                                                      Array4<Real const> const& r,
-                                                      Array4<EBCellFlag const> const& f) -> Real
-                          {
-                              Real mx = Real(-1.0);
-                              amrex::Loop(b, [=,&mx] (int i, int j, int k) noexcept
-                              {
-                                  if (!f(i,j,k).isCovered()) {
-                                      Real rho_inv = Real(1.0)/r(i,j,k);
-                                      mx = amrex::max(rho_inv, mx);
-                                  }
-                              });
-                              return mx;
-                          });
-                diff_lev *= m_mu;
+                diff_lev = nu.max(0, 0, true);
             }
 
             // Forcing term -- old way of computing
@@ -131,19 +161,7 @@ void incflo::ComputeDt (int initialization, bool explicit_diffusion, double cur_
                        });
 
             if (explicit_diffusion) {
-                diff_lev = amrex::ReduceMax(rho, 0,
-                           [=] AMREX_GPU_HOST_DEVICE (Box const& b,
-                                                      Array4<Real const> const& r) -> Real
-                           {
-                               Real mx = Real(-1.0);
-                               amrex::Loop(b, [=,&mx] (int i, int j, int k) noexcept
-                               {
-                                   Real rho_inv = Real(1.0)/r(i,j,k);
-                                   mx = amrex::max(rho_inv, mx);
-                               });
-                               return mx;
-                           });
-                diff_lev *= m_mu;
+                diff_lev = nu.max(0, 0, true);
             }
 
             // Forcing term -- old way of computing

src/incflo_regrid.cpp

@@ -59,10 +59,6 @@ void incflo::MakeNewLevelFromCoarse (int lev,
 #else
     macproj = std::make_unique<Hydro::MacProjector>(Geom(0,lev));
 #endif
-
-#ifdef INCFLO_USE_PARTICLES
-    particleData.Redistribute();
-#endif
 }
 
 // Remake an existing level using provided BoxArray and DistributionMapping and
@@ -120,10 +116,6 @@ void incflo::RemakeLevel (int lev, Real time, const BoxArray& ba,
 #else
     macproj = std::make_unique<Hydro::MacProjector>(Geom(0,finest_level));
 #endif
-
-#ifdef INCFLO_USE_PARTICLES
-    particleData.Redistribute();
-#endif
 }
 
 // Rebuild macproj on demand.  Note that this must not be done inside ClearLevel: