From 12869a959edf30d6c64c833fe08b1ffcdbc8ca16 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Mon, 29 Apr 2024 10:44:12 -0700 Subject: [PATCH] make sure this compiles in 2D and add comments re cylinder --- src/particles/incflo_PC.H | 12 ++-- src/particles/incflo_PCEvolve.cpp | 112 +++++++++++++++++++----------- src/particles/incflo_PCInit.cpp | 2 +- 3 files changed, 80 insertions(+), 46 deletions(-) diff --git a/src/particles/incflo_PC.H b/src/particles/incflo_PC.H index 4cdf0eb1..3def0aef 100644 --- a/src/particles/incflo_PC.H +++ b/src/particles/incflo_PC.H @@ -101,9 +101,9 @@ class incflo_PC : public amrex::ParticleContainer< incflo_ParticlesRealIdxAoS:: /*! Evolve particles for one time step */ virtual void EvolveParticles (int, amrex::Real, - const amrex::MultiFab* a_u_mac, - const amrex::MultiFab* a_v_mac, - const amrex::MultiFab* a_w_mac); + AMREX_D_DECL(const amrex::MultiFab* a_u_mac, + const amrex::MultiFab* a_v_mac, + const amrex::MultiFab* a_w_mac)); /*! Get real-type particle attribute names */ virtual amrex::Vector varNames () const @@ -122,9 +122,9 @@ class incflo_PC : public amrex::ParticleContainer< incflo_ParticlesRealIdxAoS:: /*! Uses midpoint method to advance particles using flow velocity. */ virtual void AdvectWithFlow ( int, amrex::Real, - const amrex::MultiFab* a_u_mac, - const amrex::MultiFab* a_v_mac, - const amrex::MultiFab* a_w_mac); + AMREX_D_DECL(const amrex::MultiFab* a_u_mac, + const amrex::MultiFab* a_v_mac, + const amrex::MultiFab* a_w_mac)); /*! Compute mass density */ virtual void massDensity ( amrex::MultiFab&, const int&, const int& a_comp = 0) const; diff --git a/src/particles/incflo_PCEvolve.cpp b/src/particles/incflo_PCEvolve.cpp index f237604a..0689097b 100644 --- a/src/particles/incflo_PCEvolve.cpp +++ b/src/particles/incflo_PCEvolve.cpp @@ -27,9 +27,9 @@ void incflo_PC::EvolveParticles ( int a_l // void incflo_PC::AdvectWithFlow (int a_lev, Real a_dt, - const MultiFab* a_umac, const MultiFab* a_vmac, const MultiFab* a_wmac) + AMREX_D_DECL(const MultiFab* a_umac, const MultiFab* a_vmac, const MultiFab* a_wmac)) { - BL_PROFILE("incflo_PC::AdvectWithUmac()"); + BL_PROFILE("incflo_PC::AdvectWithFlow()"); AMREX_ASSERT(OK(a_lev, a_lev, a_umac[0].nGrow()-1)); AMREX_ASSERT(a_lev >= 0 && a_lev < GetParticles().size()); @@ -43,9 +43,9 @@ void incflo_PC::AdvectWithFlow (int a_lev, const auto phi = geom.ProbHiArray(); const auto dxi = geom.InvCellSizeArray(); - bool is_periodic_in_x = geom.isPeriodic(0); - bool is_periodic_in_y = geom.isPeriodic(1); - bool is_periodic_in_z = geom.isPeriodic(2); + GpuArray is_periodic = {AMREX_D_DECL(geom.isPeriodic(0), + geom.isPeriodic(1), + geom.isPeriodic(2))}; for (int ipass = 0; ipass < 2; ipass++) { @@ -90,14 +90,18 @@ void incflo_PC::AdvectWithFlow (int a_lev, v_ptr[dim][i] = p.pos(dim); p.pos(dim) += static_cast(ParticleReal(0.5)*a_dt*v[dim]); - // Reflect off bottom boundary - if (dim == 2 && !is_periodic_in_z && p.pos(dim) < plo[2]) + // + // Reflect off low domain boundaries if not periodic + // + 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; } - // Reflect off top boundary - if (dim == 2 && !is_periodic_in_z && p.pos(dim) > phi[2]) + // + // 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; @@ -109,14 +113,18 @@ void incflo_PC::AdvectWithFlow (int a_lev, p.pos(dim) = v_ptr[dim][i] + static_cast(a_dt*v[dim]); v_ptr[dim][i] = v[dim]; - // Reflect off bottom boundary - if (dim == 2 && !is_periodic_in_z && p.pos(dim) < plo[2]) + // + // Reflect off low domain boundaries if not periodic + // + 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; } - // Reflect off top boundary - if (dim == 2 && !is_periodic_in_z && p.pos(dim) > phi[2]) + // + // 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; @@ -129,38 +137,64 @@ void incflo_PC::AdvectWithFlow (int a_lev, ParmParse pp("cylinder"); - Real cyl_radius; - pp.get("radius",cyl_radius); - - Array cyl_center; - pp.get("center",cyl_center); - Real x_ctr = cyl_center[0]; - Real y_ctr = cyl_center[1]; - Real z_ctr = cyl_center[2]; + Real cyl_radius = -1.; + pp.query("radius",cyl_radius); + + // + // The code below is for the specific problem of flow inside a cylinder + // + // This is just one possible problem set-up but demonstrates how to remove + // particles if they leave a specified region. + // + if (cyl_radius > 0.) { + Array cyl_center; + 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 - // Remove particles that are outside of the cylindner - for (ParIterType pti(*this, a_lev); pti.isValid(); ++pti) - { - auto& ptile = ParticlesAt(a_lev, pti); - auto& aos = ptile.GetArrayOfStructs(); - auto& soa = ptile.GetStructOfArrays(); - const int n = aos.numParticles(); - auto *p_pbox = aos().data(); + int cyl_direction; + pp.get("direction",cyl_direction); + AMREX_ALWAYS_ASSERT(cyl_direction >= 0 && cyl_direction <= AMREX_SPACEDIM); - ParallelFor(n, [=] AMREX_GPU_DEVICE (int i) + // Remove particles that are outside of the cylindner + for (ParIterType pti(*this, a_lev); pti.isValid(); ++pti) { - ParticleType& p = p_pbox[i]; + auto& ptile = ParticlesAt(a_lev, pti); + auto& aos = ptile.GetArrayOfStructs(); + auto& soa = ptile.GetStructOfArrays(); + const int n = aos.numParticles(); + auto *p_pbox = aos().data(); - Real x = p.pos(0) - x_ctr; - Real y = p.pos(1) - y_ctr; + ParallelFor(n, [=] AMREX_GPU_DEVICE (int i) + { + ParticleType& p = p_pbox[i]; - Real r = std::sqrt(x*x + y*y); + 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; +#endif - if (r > cyl_radius) { - p.id() = -1; - } - }); - } + Real r; + 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) { + r = std::sqrt(y*y + z*z); +#endif + } + + if (r > cyl_radius) { + p.id() = -1; + } + }); + } + } // cyl_radius > 0 if (m_verbose > 1) { diff --git a/src/particles/incflo_PCInit.cpp b/src/particles/incflo_PCInit.cpp index 12d2e1b8..ada57ac1 100644 --- a/src/particles/incflo_PCInit.cpp +++ b/src/particles/incflo_PCInit.cpp @@ -114,7 +114,7 @@ void incflo_PC::initializeParticlesUniformDistributionInBox ( const RealBox& par if (vf_arr(i,j,k) > 0.) { #endif if (i%2 == 0 && j%2 == 1) { - if (particle_init_domain.contains(RealVect(x,y,z))) { + if (particle_init_domain.contains(RealVect(AMREX_D_DECL(x,y,z)))) { num_particles_arr(i,j,k) = particles_per_cell; } }