Skip to content

Commit

Permalink
make sure this compiles in 2D and add comments re cylinder
Browse files Browse the repository at this point in the history
  • Loading branch information
@asalmgren
asalmgren committed Apr 29, 2024
1 parent 2830d38 commit 12869a9
Show file tree
Hide file tree
Showing 3 changed files with 80 additions and 46 deletions.
12 changes: 6 additions & 6 deletions src/particles/incflo_PC.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<std::string> varNames () const
Expand All @@ -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;
Expand Down
112 changes: 73 additions & 39 deletions src/particles/incflo_PCEvolve.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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());

Expand All @@ -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<const int, AMREX_SPACEDIM> is_periodic = {AMREX_D_DECL(geom.isPeriodic(0),
geom.isPeriodic(1),
geom.isPeriodic(2))};

for (int ipass = 0; ipass < 2; ipass++)
{
Expand Down Expand Up @@ -90,14 +90,18 @@ void incflo_PC::AdvectWithFlow (int a_lev,
v_ptr[dim][i] = p.pos(dim);
p.pos(dim) += static_cast<ParticleReal>(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;
Expand All @@ -109,14 +113,18 @@ void incflo_PC::AdvectWithFlow (int a_lev,
p.pos(dim) = v_ptr[dim][i] + static_cast<ParticleReal>(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;
Expand All @@ -129,38 +137,64 @@ void incflo_PC::AdvectWithFlow (int a_lev,

ParmParse pp("cylinder");

Real cyl_radius;
pp.get("radius",cyl_radius);

Array<Real,3> 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<Real,AMREX_SPACEDIM> 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)
{
Expand Down
2 changes: 1 addition & 1 deletion src/particles/incflo_PCInit.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
}
}
Expand Down

0 comments on commit 12869a9

Please sign in to comment.