Skip to content

Commit

Permalink
Port InitRandom for pure SoA (#3325)
Browse files Browse the repository at this point in the history 
This is useful for porting and testing other features, like Checkpoint /
Restart. I've also added a test.

The proposed changes:
- [x] fix a bug or incorrect behavior in AMReX
- [x] add new capabilities to AMReX
- [ ] changes answers in the test suite to more than roundoff level
- [ ] are likely to significantly affect the results of downstream AMReX
users
- [ ] include documentation in the code and/or rst files, if appropriate
  • Loading branch information
@atmyers
atmyers authored May 23, 2023
1 parent b625d63 commit 4167cbe
Show file tree
Hide file tree
Showing 6 changed files with 373 additions and 53 deletions.
180 changes: 127 additions & 53 deletions Src/Particle/AMReX_ParticleInit.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1046,17 +1046,13 @@ InitRandom (Long icount,

for (Long j = 0; j < icount; j++)
{
ParticleType p;
Particle<0, 0> ptest;

for (int i = 0; i < AMREX_SPACEDIM; i++) {
p.pos(i) = pos[j*AMREX_SPACEDIM + i];
ptest.pos(i) = pos[j*AMREX_SPACEDIM + i];
}

for (int i = 0; i < NStructReal; i++) {
p.rdata(i) = static_cast<ParticleReal>(pdata.real_struct_data[i]);
}

if (!Where(p, pld)) {
if (!Where(ptest, pld)) {
amrex::Abort("ParticleContainer::InitRandom(): invalid particle");
}

Expand All @@ -1067,27 +1063,58 @@ InitRandom (Long icount,

if (who == MyProc) {

// We own it. Add it at the appropriate level.
p.id() = ParticleType::NextID();
p.cpu() = MyProc;
if constexpr(!ParticleType::is_soa_particle)
{
ParticleType p;
for (int i = 0; i < AMREX_SPACEDIM; i++) {
p.pos(i) = pos[j*AMREX_SPACEDIM + i];
}

for (int i = 0; i < NStructInt; i++) {
p.idata(i) = pdata.int_struct_data[i];
}
// We own it. Add it at the appropriate level.
p.id() = ParticleType::NextID();
p.cpu() = MyProc;

// add the struct
host_particles[pld.m_lev][ind].push_back(p);
for (int i = 0; i < NStructInt; i++) {
p.idata(i) = pdata.int_struct_data[i];
}

// add the real...
for (int i = 0; i < NArrayReal; i++) {
host_real_attribs[pld.m_lev][ind][i].push_back(static_cast<ParticleReal>(pdata.real_array_data[i]));
}
for (int i = 0; i < NStructReal; i++) {
p.rdata(i) = static_cast<ParticleReal>(pdata.real_struct_data[i]);
}

// ... and int array data
for (int i = 0; i < NArrayInt; i++) {
host_int_attribs[pld.m_lev][ind][i].push_back(pdata.int_array_data[i]);
// add the struct
host_particles[pld.m_lev][ind].push_back(p);

// add the real...
for (int i = 0; i < NArrayReal; i++) {
host_real_attribs[pld.m_lev][ind][i].push_back(static_cast<ParticleReal>(pdata.real_array_data[i]));
}

// ... and int array data
for (int i = 0; i < NArrayInt; i++) {
host_int_attribs[pld.m_lev][ind][i].push_back(pdata.int_array_data[i]);
}
} else {
for (int i = 0; i < AMREX_SPACEDIM; i++) {
host_real_attribs[pld.m_lev][ind][i].push_back(pos[j*AMREX_SPACEDIM+i]);
}

host_int_attribs[pld.m_lev][ind][0].push_back(ParticleType::NextID());
host_int_attribs[pld.m_lev][ind][1].push_back(MyProc);

host_particles[pld.m_lev][ind];

// add the real...
for (int i = AMREX_SPACEDIM; i < NArrayReal; i++) {
host_real_attribs[pld.m_lev][ind][i].push_back(static_cast<ParticleReal>(pdata.real_array_data[i]));
}

// ... and int array data
for (int i = 2; i < NArrayInt; i++) {
host_int_attribs[pld.m_lev][ind][i].push_back(pdata.int_array_data[i]);
}
}
}
}
}

for (int host_lev = 0; host_lev < static_cast<int>(host_particles.size()); ++host_lev)
Expand All @@ -1099,11 +1126,20 @@ InitRandom (Long icount,

auto& dst_tile = GetParticles(host_lev)[std::make_pair(grid,tile)];
auto old_size = dst_tile.GetArrayOfStructs().size();
auto new_size = old_size + src_tile.size();
auto new_size = old_size;
if constexpr(!ParticleType::is_soa_particle)
{
new_size += src_tile.size();
} else {
new_size += host_real_attribs[host_lev][std::make_pair(grid,tile)][0].size();
}
dst_tile.resize(new_size);

Gpu::copyAsync(Gpu::hostToDevice, src_tile.begin(), src_tile.end(),
dst_tile.GetArrayOfStructs().begin() + old_size);
if constexpr(!ParticleType::is_soa_particle)
{
Gpu::copyAsync(Gpu::hostToDevice, src_tile.begin(), src_tile.end(),
dst_tile.GetArrayOfStructs().begin() + old_size);
}

for (int i = 0; i < NArrayReal; ++i) {
Gpu::copyAsync(Gpu::hostToDevice,
Expand Down Expand Up @@ -1148,50 +1184,79 @@ InitRandom (Long icount,
host_int_attribs.resize(finestLevel()+1);

for (Long icnt = 0; icnt < M; icnt++) {
ParticleType p;
Particle<0, 0> ptest;
for (int i = 0; i < AMREX_SPACEDIM; i++) {
do {
r = amrex::Random();
x = geom.ProbLo(i) + (r * len[i]);
}
while (static_cast<ParticleReal>(x) < static_cast<ParticleReal>(xlo[i]) || static_cast<ParticleReal>(x) >= static_cast<ParticleReal>(xhi[i]));

p.pos(i) = static_cast<ParticleReal>(x);
ptest.pos(i) = static_cast<ParticleReal>(x);

AMREX_ASSERT(p.pos(i) < geom.ProbHi(i));
AMREX_ASSERT(ptest.pos(i) < geom.ProbHi(i));
}

for (int i = 0; i < NStructReal; i++) {
p.rdata(i) = static_cast<ParticleReal>(pdata.real_struct_data[i]);
}

// the int struct data
p.id() = ParticleType::NextID();
p.cpu() = ParallelDescriptor::MyProc();

for (int i = 0; i < NStructInt; i++) {
p.idata(i) = pdata.int_struct_data[i];
}
ptest.id() = ParticleType::NextID();
ptest.cpu() = ParallelDescriptor::MyProc();

// locate the particle
if (!Where(p, pld))
if (!Where(ptest, pld))
{
amrex::Abort("ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt>::InitRandom(): invalid particle");
}
AMREX_ASSERT(pld.m_lev >= 0 && pld.m_lev <= finestLevel());
std::pair<int, int> ind(pld.m_grid, pld.m_tile);

// add the struct
host_particles[pld.m_lev][ind].push_back(p);
if constexpr(!ParticleType::is_soa_particle)
{
ParticleType p;
for (int i = 0; i < AMREX_SPACEDIM; i++) {
p.pos(i) = ptest.pos(i);;
}

// add the real...
for (int i = 0; i < NArrayReal; i++) {
host_real_attribs[pld.m_lev][ind][i].push_back(static_cast<ParticleReal>(pdata.real_array_data[i]));
}
p.id() = ptest.id();
p.cpu() = ptest.cpu();

// ... and int array data
for (int i = 0; i < NArrayInt; i++) {
host_int_attribs[pld.m_lev][ind][i].push_back(pdata.int_array_data[i]);
for (int i = 0; i < NStructReal; i++) {
p.rdata(i) = static_cast<ParticleReal>(pdata.real_struct_data[i]);
}

for (int i = 0; i < NStructInt; i++) {
p.idata(i) = pdata.int_struct_data[i];
}

// add the struct
host_particles[pld.m_lev][ind].push_back(p);

// add the real...
for (int i = 0; i < NArrayReal; i++) {
host_real_attribs[pld.m_lev][ind][i].push_back(static_cast<ParticleReal>(pdata.real_array_data[i]));
}

// ... and int array data
for (int i = 0; i < NArrayInt; i++) {
host_int_attribs[pld.m_lev][ind][i].push_back(pdata.int_array_data[i]);
}
} else {
for (int i = 0; i < AMREX_SPACEDIM; i++) {
host_real_attribs[pld.m_lev][ind][i].push_back(ptest.pos(i));
}

host_int_attribs[pld.m_lev][ind][0].push_back(ptest.id());
host_int_attribs[pld.m_lev][ind][1].push_back(ptest.cpu());

host_particles[pld.m_lev][ind];

// add the real...
for (int i = AMREX_SPACEDIM; i < NArrayReal; i++) {
host_real_attribs[pld.m_lev][ind][i].push_back(static_cast<ParticleReal>(pdata.real_array_data[i]));
}

// ... and int array data
for (int i = 2; i < NArrayInt; i++) {
host_int_attribs[pld.m_lev][ind][i].push_back(pdata.int_array_data[i]);
}
}
}

Expand All @@ -1204,11 +1269,20 @@ InitRandom (Long icount,

auto& dst_tile = GetParticles(host_lev)[std::make_pair(grid,tile)];
auto old_size = dst_tile.GetArrayOfStructs().size();
auto new_size = old_size + src_tile.size();
auto new_size = old_size;
if constexpr(!ParticleType::is_soa_particle)
{
new_size += src_tile.size();
} else {
new_size += host_real_attribs[host_lev][std::make_pair(grid,tile)][0].size();
}
dst_tile.resize(new_size);

Gpu::copyAsync(Gpu::hostToDevice, src_tile.begin(), src_tile.end(),
dst_tile.GetArrayOfStructs().begin() + old_size);
if constexpr(!ParticleType::is_soa_particle)
{
Gpu::copyAsync(Gpu::hostToDevice, src_tile.begin(), src_tile.end(),
dst_tile.GetArrayOfStructs().begin() + old_size);
}

for (int i = 0; i < NArrayReal; ++i) {
Gpu::copyAsync(Gpu::hostToDevice,
Expand Down
14 changes: 14 additions & 0 deletions Tests/Particles/InitRandom/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,14 @@
# This tests requires particle support
if (NOT AMReX_PARTICLES)
return()
endif ()

foreach(D IN LISTS AMReX_SPACEDIM)
set(_sources main.cpp)
set(_input_files inputs )

setup_test(${D} _sources _input_files)

unset(_sources)
unset(_input_files)
endforeach()
31 changes: 31 additions & 0 deletions Tests/Particles/InitRandom/GNUmakefile
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,31 @@
AMREX_HOME = ../../../

# DEBUG = TRUE
DEBUG = FALSE

DIM = 3

COMP = gnu

PRECISION = DOUBLE

USE_MPI = TRUE
MPI_THREAD_MULTIPLE = FALSE

USE_OMP = FALSE

TINY_PROFILE = TRUE

USE_PARTICLES = TRUE

###################################################

EBASE = main

include $(AMREX_HOME)/Tools/GNUMake/Make.defs

include ./Make.package
include $(AMREX_HOME)/Src/Base/Make.package
include $(AMREX_HOME)/Src/Particle/Make.package

include $(AMREX_HOME)/Tools/GNUMake/Make.rules
1 change: 1 addition & 0 deletions Tests/Particles/InitRandom/Make.package
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
CEXE_sources += main.cpp
13 changes: 13 additions & 0 deletions Tests/Particles/InitRandom/inputs
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,13 @@
# Domain size
ncells = 64

# Maximum allowable size of each subdomain in the problem domain;
# this is used to decompose the domain for parallel calculations.
max_grid_size = 8

# Number of levels
nlevs = 1

# Number of particles per cell
nppc = 2

Loading

0 comments on commit 4167cbe

Please sign in to comment.