Skip to content

Commit

Permalink
fix particle initialization so all particles are actually in the doma…
Browse files Browse the repository at this point in the history 
…in; also add the ability to specify where particles are created by defining a box in physical space (#1514)
  • Loading branch information
@asalmgren
asalmgren authored Mar 20, 2024
1 parent c8677ef commit 24e4a9a
Show file tree
Hide file tree
Showing 5 changed files with 195 additions and 41 deletions.
6 changes: 5 additions & 1 deletion Exec/RegTests/ParticlesOverWoA/inputs
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
# ------------------ INPUTS TO MAIN PROGRAM -------------------
max_step = 40
max_step = 80

amrex.fpe_trap_invalid = 1

Expand All @@ -25,6 +25,10 @@ zhi.type = "SlipWall"

# PARTICLES
erf.use_tracer_particles = 1
tracer_particles.initial_distribution_type = box
tracer_particles.particle_box_lo = 3.95 -1.0 -1.0
tracer_particles.particle_box_hi = 4.00 2.0 1.0
tracer_particles.place_randomly_in_cells = true

# TIME STEP CONTROL
erf.fixed_dt = 1E-3
Expand Down
7 changes: 6 additions & 1 deletion Source/Particles/ERFPC.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,6 +43,7 @@ namespace ERFParticleInitializations
{
/* list of particle initializations */
const std::string init_default = "default";
const std::string init_box = "box";
const std::string init_uniform = "uniform";

}
Expand Down Expand Up @@ -215,12 +216,17 @@ class ERFPC : public amrex::ParticleContainer< ERFParticlesRealIdxAoS::ncomps,
void initializeParticlesDefaultHydro (const std::unique_ptr<amrex::MultiFab>& a_ptr = nullptr);
/*! Default particle initialization */
void initializeParticlesUniformDistribution (const std::unique_ptr<amrex::MultiFab>& a_ptr = nullptr);
void initializeParticlesUniformDistributionInBox (const std::unique_ptr<amrex::MultiFab>& a_ptr,
const amrex::RealBox& particle_box);

protected:

bool m_advect_w_flow; /*!< advect with flow velocity */
bool m_advect_w_gravity; /*!< advect under gravitational force */

amrex::RealBox particle_box;
bool place_randomly_in_cells;

std::string m_name; /*!< name of this particle species */

std::string m_initialization_type; /*!< initial particle distribution type */
Expand All @@ -233,7 +239,6 @@ class ERFPC : public amrex::ParticleContainer< ERFParticlesRealIdxAoS::ncomps,

/*! Default particle initialization */
void initializeParticlesDefault (const std::unique_ptr<amrex::MultiFab>& a_ptr = nullptr);

};

#endif
Expand Down
2 changes: 1 addition & 1 deletion Source/Particles/ERFPCEvolve.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,7 +64,7 @@ void ERFPC::AdvectWithFlow ( MultiFab* a_umac,
{
for (int i = 0; i < AMREX_SPACEDIM; i++)
{
int ng = a_umac[i].nGrow();
IntVect ng = a_umac[i].nGrowVect();
raii_umac[i] = std::make_unique<MultiFab>
(convert(m_gdb->ParticleBoxArray(a_lev), IntVect::TheDimensionVector(i)),
m_gdb->ParticleDistributionMap(a_lev), a_umac[i].nComp(), ng);
Expand Down
213 changes: 179 additions & 34 deletions Source/Particles/ERFPCInitializations.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,30 @@ void ERFPC::readInputs ()
m_initialization_type = ERFParticleInitializations::init_default;
pp.query("initial_distribution_type", m_initialization_type);

if (m_initialization_type == ERFParticleInitializations::init_box)
{
Vector<Real> particle_box_lo(AMREX_SPACEDIM);
Vector<Real> particle_box_hi(AMREX_SPACEDIM);

// Defaults
for (int i = 0; i < AMREX_SPACEDIM; i++) { particle_box_lo[i] = Geom(0).ProbLo(i); }
for (int i = 0; i < AMREX_SPACEDIM; i++) { particle_box_hi[i] = Geom(0).ProbHi(i); }

pp.queryAdd("particle_box_lo", particle_box_lo, AMREX_SPACEDIM);
AMREX_ASSERT(particle_box_lo.size() == AMREX_SPACEDIM);

pp.queryAdd("particle_box_hi", particle_box_hi, AMREX_SPACEDIM);
AMREX_ASSERT(particle_box_hi.size() == AMREX_SPACEDIM);

particle_box.setLo(particle_box_lo);
particle_box.setHi(particle_box_hi);

// We default to placing the particles randomly within each cell,
// but can override this for regression testing
place_randomly_in_cells = true;
pp.query("place_randomly_in_cells", place_randomly_in_cells);
}

m_ppc_init = 1;
pp.query("initial_particles_per_cell", m_ppc_init);

Expand All @@ -35,6 +59,8 @@ void ERFPC::InitializeParticles (const std::unique_ptr<MultiFab>& a_height_ptr)

if (m_initialization_type == ERFParticleInitializations::init_default) {
initializeParticlesDefault( a_height_ptr );
} else if (m_initialization_type == ERFParticleInitializations::init_box) {
initializeParticlesUniformDistributionInBox( a_height_ptr , particle_box );
} else if (m_initialization_type == ERFParticleInitializations::init_uniform) {
initializeParticlesUniformDistribution( a_height_ptr );
} else {
Expand Down Expand Up @@ -76,7 +102,7 @@ void ERFPC::initializeParticlesDefaultTracersWoA (const std::unique_ptr<MultiFab
auto num_particles_arr = num_particles[mfi].array();
ParallelFor(tile_box, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
if (i == 3) {
if (i == 100) {
num_particles_arr(i,j,k) = 1;
}
});
Expand Down Expand Up @@ -126,15 +152,17 @@ void ERFPC::initializeParticlesDefaultTracersWoA (const std::unique_ptr<MultiFab
const auto height_arr = (*a_height_ptr)[mfi].array();
ParallelFor(tile_box, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
if (i == 3) {
Real r[3] = {0.5, 0.5, 0.5}; // this means place at cell center
if (i == 100) {
Real r[3] = {0.5, 0.5, 0.25}; // this means place at cell center in x/y, 1/4 way up cell in z
Real v[3] = {0.0, 0.0, 0.0}; // with 0 initial velocity

Real x = plo[0] + (i + r[0])*dx[0];
Real y = plo[1] + (j + r[1])*dx[1];
Real z = height_arr(i,j,k)
+ r[2] * ( height_arr(i,j,k+1)
- height_arr(i,j,k) );
Real zlo = 0.25 * (height_arr(i,j ,k ) + height_arr(i+1,j ,k ) +
height_arr(i,j+1,k ) + height_arr(i+1,j+1,k ) );
Real zhi = 0.25 * (height_arr(i,j ,k+1) + height_arr(i+1,j ,k+1) +
height_arr(i,j+1,k+1) + height_arr(i+1,j+1,k ) );
Real z = zlo + r[2] * (zhi - zlo);

auto& p = aos[offset_arr(i,j,k)];
p.id() = pid + offset_arr(i,j,k);
Expand Down Expand Up @@ -264,9 +292,11 @@ void ERFPC::initializeParticlesDefaultHydro (const std::unique_ptr<MultiFab>& a_

Real x = plo[0] + (i + r[0])*dx[0];
Real y = plo[1] + (j + r[1])*dx[1];
Real z = height_arr(i,j,k)
+ r[2] * ( height_arr(i,j,k+1)
- height_arr(i,j,k) );
Real zlo = 0.25 * (height_arr(i,j ,k ) + height_arr(i+1,j ,k ) +
height_arr(i,j+1,k ) + height_arr(i+1,j+1,k ) );
Real zhi = 0.25 * (height_arr(i,j ,k+1) + height_arr(i+1,j ,k+1) +
height_arr(i,j+1,k+1) + height_arr(i+1,j+1,k ) );
Real z = zlo + r[2] * (zhi - zlo);

auto& p = aos[offset_arr(i,j,k)];
p.id() = pid + offset_arr(i,j,k);
Expand Down Expand Up @@ -325,7 +355,15 @@ void ERFPC::initializeParticlesDefaultHydro (const std::unique_ptr<MultiFab>& a_
* by "initial_particles_per_cell", and they are randomly distributed. */
void ERFPC::initializeParticlesUniformDistribution (const std::unique_ptr<MultiFab>& a_height_ptr)
{
BL_PROFILE("ERFPC::initializeParticlesUniformDistribution");
initializeParticlesUniformDistributionInBox ( a_height_ptr, Geom(0).ProbDomain() );
}

/*! Uniform distribution: the number of particles per grid cell is specified
* by "initial_particles_per_cell", and they are randomly distributed. */
void ERFPC::initializeParticlesUniformDistributionInBox (const std::unique_ptr<MultiFab>& a_height_ptr,
const RealBox& particle_init_domain)
{
BL_PROFILE("ERFPC::initializeParticlesUniformDistributionInBox");

const int lev = 0;
const auto dx = Geom(lev).CellSizeArray();
Expand All @@ -340,10 +378,32 @@ void ERFPC::initializeParticlesUniformDistribution (const std::unique_ptr<MultiF
for(MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) {
const Box& tile_box = mfi.tilebox();
auto num_particles_arr = num_particles[mfi].array();
ParallelFor(tile_box, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
num_particles_arr(i,j,k) = particles_per_cell;
});
if (a_height_ptr) {
const auto height_arr = (*a_height_ptr)[mfi].array();
ParallelFor(tile_box, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
Real x = plo[0] + (i + 0.5)*dx[0];
Real y = plo[1] + (j + 0.5)*dx[1];
Real z = 0.125 * (height_arr(i,j ,k ) + height_arr(i+1,j ,k ) +
height_arr(i,j+1,k ) + height_arr(i+1,j+1,k ) +
height_arr(i,j ,k+1) + height_arr(i+1,j ,k+1) +
height_arr(i,j+1,k+1) + height_arr(i+1,j+1,k ) );
if (particle_init_domain.contains(RealVect(x,y,z))) {
num_particles_arr(i,j,k) = particles_per_cell;
}
});

} else {
ParallelFor(tile_box, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
Real x = plo[0] + (i + 0.5)*dx[0];
Real y = plo[1] + (j + 0.5)*dx[1];
Real z = plo[2] + (k + 0.5)*dx[2];
if (particle_init_domain.contains(RealVect(x,y,z))) {
num_particles_arr(i,j,k) = particles_per_cell;
}
});
}
}

iMultiFab offsets( ParticleBoxArray(lev),
Expand Down Expand Up @@ -376,6 +436,8 @@ void ERFPC::initializeParticlesUniformDistribution (const std::unique_ptr<MultiF
auto* vz_ptr = soa.GetRealData(ERFParticlesRealIdxSoA::vz).data();
auto* mass_ptr = soa.GetRealData(ERFParticlesRealIdxSoA::mass).data();

const auto num_particles_arr = num_particles[mfi].array();

auto my_proc = ParallelDescriptor::MyProc();
Long pid;
{
Expand All @@ -385,47 +447,106 @@ void ERFPC::initializeParticlesUniformDistribution (const std::unique_ptr<MultiF
AMREX_ALWAYS_ASSERT_WITH_MESSAGE( static_cast<Long>(pid + np) < LastParticleID,
"Error: overflow on particle id numbers!" );

if (a_height_ptr) {
if (a_height_ptr && place_randomly_in_cells) {

const auto height_arr = (*a_height_ptr)[mfi].array();

const auto height_arr = (*a_height_ptr)[mfi].array();
ParallelForRNG(tile_box, [=] AMREX_GPU_DEVICE (int i, int j, int k,
const RandomEngine& rnd_engine) noexcept
{
int start = offset_arr(i,j,k);
for (int n = start; n < start+particles_per_cell; n++) {
for (int n = start; n < start+num_particles_arr(i,j,k); n++) {
Real r[3] = {Random(rnd_engine), Random(rnd_engine), Random(rnd_engine)};
Real v[3] = {0.0, 0.0, 0.0};

Real x = plo[0] + (i + r[0])*dx[0];
Real y = plo[1] + (j + r[1])*dx[1];
Real z = height_arr(i,j,k)
+ r[2] * ( height_arr(i,j,k+1)
- height_arr(i,j,k) );

Real sx[] = { amrex::Real(1.) - r[0], r[0]};
Real sy[] = { amrex::Real(1.) - r[1], r[1]};

Real height_at_pxy_lo = 0.;
for (int ii = 0; ii < 2; ++ii) {
for (int jj = 0; jj < 2; ++jj) {
height_at_pxy_lo += sx[ii] * sy[jj] * height_arr(i+ii,j+jj,k);
}
}
Real height_at_pxy_hi = 0.;
for (int ii = 0; ii < 2; ++ii) {
for (int jj = 0; jj < 2; ++jj) {
height_at_pxy_hi += sx[ii] * sy[jj] * height_arr(i+ii,j+jj,k+1);
}
}

Real z = height_at_pxy_lo + r[2] * (height_at_pxy_hi - height_at_pxy_lo);

auto& p = aos[n];
p.id() = pid + n;
p.cpu() = my_proc;
p.pos(0) = x;
p.pos(1) = y;
p.pos(2) = z;

p.pos(0) = x; p.pos(1) = y; p.pos(2) = z;

p.idata(ERFParticlesIntIdxAoS::k) = k;

vx_ptr[n] = v[0];
vy_ptr[n] = v[1];
vz_ptr[n] = v[2];
vx_ptr[n] = v[0]; vy_ptr[n] = v[1]; vz_ptr[n] = v[2];

mass_ptr[n] = 1.0e-6;
}
});

} else {
} else if (a_height_ptr && !place_randomly_in_cells) {

const auto height_arr = (*a_height_ptr)[mfi].array();

ParallelFor(tile_box, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
int start = offset_arr(i,j,k);
for (int n = start; n < start+num_particles_arr(i,j,k); n++) {
Real r[3] = {0.3, 0.7, 0.25};
Real v[3] = {0.0, 0.0, 0.0};

Real x = plo[0] + (i + r[0])*dx[0];
Real y = plo[1] + (j + r[1])*dx[1];

Real sx[] = { amrex::Real(1.) - r[0], r[0]};
Real sy[] = { amrex::Real(1.) - r[1], r[1]};

Real height_at_pxy_lo = 0.;
for (int ii = 0; ii < 2; ++ii) {
for (int jj = 0; jj < 2; ++jj) {
height_at_pxy_lo += sx[ii] * sy[jj] * height_arr(i+ii,j+jj,k);
}
}
Real height_at_pxy_hi = 0.;
for (int ii = 0; ii < 2; ++ii) {
for (int jj = 0; jj < 2; ++jj) {
height_at_pxy_hi += sx[ii] * sy[jj] * height_arr(i+ii,j+jj,k+1);
}
}

Real z = height_at_pxy_lo + r[2] * (height_at_pxy_hi - height_at_pxy_lo);

auto& p = aos[n];
p.id() = pid + n;
p.cpu() = my_proc;

p.pos(0) = x; p.pos(1) = y; p.pos(2) = z;

p.idata(ERFParticlesIntIdxAoS::k) = k;

vx_ptr[n] = v[0]; vy_ptr[n] = v[1]; vz_ptr[n] = v[2];

mass_ptr[n] = 1.0e-6;
}
});

} else if (!a_height_ptr && place_randomly_in_cells) {

ParallelForRNG(tile_box, [=] AMREX_GPU_DEVICE (int i, int j, int k,
const RandomEngine& rnd_engine) noexcept
{
int start = offset_arr(i,j,k);
for (int n = start; n < start+particles_per_cell; n++) {
for (int n = start; n < start+num_particles_arr(i,j,k); n++) {
Real r[3] = {Random(rnd_engine), Random(rnd_engine), Random(rnd_engine)};
Real v[3] = {0.0, 0.0, 0.0};

Expand All @@ -436,15 +557,39 @@ void ERFPC::initializeParticlesUniformDistribution (const std::unique_ptr<MultiF
auto& p = aos[n];
p.id() = pid + n;
p.cpu() = my_proc;
p.pos(0) = x;
p.pos(1) = y;
p.pos(2) = z;

p.pos(0) = x; p.pos(1) = y; p.pos(2) = z;

p.idata(ERFParticlesIntIdxAoS::k) = k;

vx_ptr[n] = v[0]; vy_ptr[n] = v[1]; vz_ptr[n] = v[2];

mass_ptr[n] = 1.0e-6;
}
});

} else { // if (!a_height_ptr && !place_randomly_in_cells) {

ParallelFor(tile_box, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
int start = offset_arr(i,j,k);
for (int n = start; n < start+num_particles_arr(i,j,k); n++) {
Real r[3] = {0.3, 0.7, 0.25};
Real v[3] = {0.0, 0.0, 0.0};

Real x = plo[0] + (i + r[0])*dx[0];
Real y = plo[1] + (j + r[1])*dx[1];
Real z = plo[2] + (k + r[2])*dx[2];

auto& p = aos[n];
p.id() = pid + n;
p.cpu() = my_proc;

p.pos(0) = x; p.pos(1) = y; p.pos(2) = z;

p.idata(ERFParticlesIntIdxAoS::k) = k;

vx_ptr[n] = v[0];
vy_ptr[n] = v[1];
vz_ptr[n] = v[2];
vx_ptr[n] = v[0]; vy_ptr[n] = v[1]; vz_ptr[n] = v[2];

mass_ptr[n] = 1.0e-6;
}
Expand Down
Loading

0 comments on commit 24e4a9a

Please sign in to comment.