fix issue 187

src/particles/incflo_PC.H

@@ -98,7 +98,8 @@ class incflo_PC : public amrex::ParticleContainer<  incflo_ParticlesRealIdxAoS::
 #endif
                                          );
 
-        /*! Evolve particles for one time step */
+        /*! Evolve particles on one level for one time step; the caller is
+         *  responsible for calling Redistribute() once all levels are done */
         virtual void EvolveParticles (int,
                                       amrex::Real,
                                       AMREX_D_DECL(const amrex::MultiFab* a_u_mac,

src/particles/incflo_PCEvolve.cpp

@@ -18,7 +18,12 @@ void incflo_PC::EvolveParticles ( int                                        a_l
         AdvectWithFlow(a_lev, a_dt_lev, AMREX_D_DECL(a_umac, a_vmac, a_wmac));
     }
 
-    Redistribute();
+    //
+    // NOTE: we do *not* Redistribute here -- the caller must do that once all
+    //       levels have been advected.  Redistributing after each level would
+    //       hand a particle that crossed from level lev into a level lev+1 grid
+    //       to level lev+1, which then advects it a second time in the same step.
+    //
     return;
 }
 
@@ -96,15 +101,21 @@ void incflo_PC::AdvectWithFlow (int                                 a_lev,
                         if (!is_periodic[dim] && p.pos(dim) < plo[dim])
                         {
                             p.pos(dim) = 2.0*plo[dim] - p.pos(dim);
-                            v_ptr[dim][i] *= -1.0;
+                            //
+                            // Here v_ptr holds the position at the start of the step,
+                            // which must be mirrored about the boundary we crossed --
+                            // not negated as in the corrector pass below, where v_ptr
+                            // holds the velocity.
+                            //
+                            v_ptr[dim][i] = static_cast<ParticleReal>(2.0*plo[dim]) - v_ptr[dim][i];
                         }
                         //
                         // Reflect off high domain boundaries if not periodic
                         //
                         if (!is_periodic[dim] && p.pos(dim) > phi[dim])
                         {
                             p.pos(dim) = 2.0*phi[dim] - p.pos(dim);
-                            v_ptr[dim][i] *= -1.0;
+                            v_ptr[dim][i] = static_cast<ParticleReal>(2.0*phi[dim]) - v_ptr[dim][i];
                         }
                     }
                 } else {
@@ -157,7 +168,7 @@ void incflo_PC::AdvectWithFlow (int                                 a_lev,
 
         int cyl_direction;
         pp.get("direction",cyl_direction);
-        AMREX_ALWAYS_ASSERT(cyl_direction >= 0 && cyl_direction <= AMREX_SPACEDIM);
+        AMREX_ALWAYS_ASSERT(cyl_direction >= 0 && cyl_direction <= 2);
 
         // Remove particles that are outside of the cylindner
         for (ParIterType pti(*this, a_lev); pti.isValid(); ++pti)
@@ -177,16 +188,26 @@ void incflo_PC::AdvectWithFlow (int                                 a_lev,
                 Real z =  p.pos(2) - z_ctr;
 #endif
 
-                Real r;
+                //
+                // Distance from the axis of the cylinder; in 2D directions 0 and 1
+                // define a slab, matching amrex::EB2::CylinderIF in 2D.
+                //
+                Real r = 0.;
                 if (cyl_direction == 2) {
                     r =  std::sqrt(x*x + y*y);
 #if (AMREX_SPACEDIM == 3)
                 } else if (cyl_direction == 1) {
                     r =  std::sqrt(x*x + z*z);
-                } else if (cyl_direction == 0) {
+                } else {
                     r =  std::sqrt(y*y + z*z);
-#endif
                 }
+#else
+                } else if (cyl_direction == 1) {
+                    r =  std::abs(x);
+                } else {
+                    r =  std::abs(y);
+                }
+#endif
 
                 if (r > cyl_radius) {
                     p.id() = -1;

src/particles/incflo_Tracers.cpp

@@ -76,6 +76,13 @@ void incflo::evolveTracerParticles (AMREX_D_DECL(Vector<MultiFab const*> const&
             particleData[incfloParticleNames::tracers]->EvolveParticles(lev, m_dt,
                                                                         AMREX_D_DECL(u_mac[lev],v_mac[lev],w_mac[lev]));
         }
+
+        //
+        // Redistribute only once all levels have been advected -- Redistribute
+        // assigns each particle to the finest level containing it, so doing this
+        // inside the loop above would advect a coarse-to-fine crosser twice.
+        //
+        particleData.Redistribute();
     }
 }
 #endif