From 24436f448424d0cef83cfaca0fca300432ef7e7b Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Fri, 21 Aug 2026 12:57:08 -0700 Subject: [PATCH] fix issue 187 --- src/particles/incflo_PC.H | 3 ++- src/particles/incflo_PCEvolve.cpp | 35 ++++++++++++++++++++++++------- src/particles/incflo_Tracers.cpp | 7 +++++++ 3 files changed, 37 insertions(+), 8 deletions(-) diff --git a/src/particles/incflo_PC.H b/src/particles/incflo_PC.H index 011e1d21..2ab62341 100644 --- a/src/particles/incflo_PC.H +++ b/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, diff --git a/src/particles/incflo_PCEvolve.cpp b/src/particles/incflo_PCEvolve.cpp index 979f71ec..8283feab 100644 --- a/src/particles/incflo_PCEvolve.cpp +++ b/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,7 +101,13 @@ 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(2.0*plo[dim]) - v_ptr[dim][i]; } // // Reflect off high domain boundaries if not periodic @@ -104,7 +115,7 @@ void incflo_PC::AdvectWithFlow (int a_lev, 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(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; diff --git a/src/particles/incflo_Tracers.cpp b/src/particles/incflo_Tracers.cpp index 9a3eb7b3..81cd031a 100644 --- a/src/particles/incflo_Tracers.cpp +++ b/src/particles/incflo_Tracers.cpp @@ -76,6 +76,13 @@ void incflo::evolveTracerParticles (AMREX_D_DECL(Vector 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