diff --git a/Tests/Particles/SENSEI_Insitu_SOA/CMakeLists.txt b/Tests/Particles/SENSEI_Insitu_SOA/CMakeLists.txt new file mode 100644 index 00000000000..8525a0e6737 --- /dev/null +++ b/Tests/Particles/SENSEI_Insitu_SOA/CMakeLists.txt @@ -0,0 +1,9 @@ +foreach(D IN LISTS AMReX_SPACEDIM) + set(_sources main.cpp) + set(_input_files inputs.rt ) + + setup_test(${D} _sources _input_files NTASKS 2) + + unset(_sources) + unset(_input_files) +endforeach() diff --git a/Tests/Particles/SENSEI_Insitu_SOA/GNUmakefile b/Tests/Particles/SENSEI_Insitu_SOA/GNUmakefile new file mode 100644 index 00000000000..653290723f3 --- /dev/null +++ b/Tests/Particles/SENSEI_Insitu_SOA/GNUmakefile @@ -0,0 +1,24 @@ +AMREX_HOME = ../../../ + +DEBUG = FALSE + +DIM = 3 + +COMP = gcc + +USE_MPI = TRUE +USE_OMP = FALSE +USE_CUDA = FALSE + +TINY_PROFILE = TRUE +USE_PARTICLES = TRUE +USE_SENSEI_INSITU = TRUE + +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)/Src/Extern/SENSEI/Make.package + +include $(AMREX_HOME)/Tools/GNUMake/Make.rules diff --git a/Tests/Particles/SENSEI_Insitu_SOA/inputs.rt b/Tests/Particles/SENSEI_Insitu_SOA/inputs.rt new file mode 100644 index 00000000000..98b6748591e --- /dev/null +++ b/Tests/Particles/SENSEI_Insitu_SOA/inputs.rt @@ -0,0 +1,10 @@ +insitu.size = (32, 64, 64) +insitu.max_grid_size = 32 +insitu.is_periodic = 1 +insitu.num_ppc = 1 +insitu.nlevs = 1 + +insitu.num_runtime_real = 0 +insitu.num_runtime_int = 0 + +particles.do_tiling=1 diff --git a/Tests/Particles/SENSEI_Insitu_SOA/main.cpp b/Tests/Particles/SENSEI_Insitu_SOA/main.cpp new file mode 100644 index 00000000000..8ff8198d3a5 --- /dev/null +++ b/Tests/Particles/SENSEI_Insitu_SOA/main.cpp @@ -0,0 +1,274 @@ +#include +#include +#include +#include + +#if !defined(AMREX_PARTICLES) || !defined(AMREX_NO_SENSEI_AMR_INST) +#error Incompatible AMReX library configuration! This tutorial requires AMREX_PARTICLES and AMREX_NO_SENSEI_AMR_INST +#endif +#include + +using namespace amrex; + +static constexpr int NR = 7; +static constexpr int NI = 4; + +int num_runtime_real = 0; +int num_runtime_int = 0; + +bool remove_negative = true; + +void get_position_unit_cell (Real* r, const IntVect& nppc, int i_part) +{ + int nx = nppc[0]; +#if AMREX_SPACEDIM > 1 + int ny = nppc[1]; +#else + int ny = 1; +#endif +#if AMREX_SPACEDIM > 2 + int nz = nppc[2]; +#else + int nz = 1; +#endif + + int ix_part = i_part/(ny * nz); + int iy_part = (i_part % (ny * nz)) % ny; + int iz_part = (i_part % (ny * nz)) / ny; + + r[0] = (0.5+ix_part)/nx; + r[1] = (0.5+iy_part)/ny; + r[2] = (0.5+iz_part)/nz; +} + +class TestParticleContainer + : public amrex::ParticleContainerPureSoA +{ + +public: + + TestParticleContainer (const Vector & a_geom, + const Vector & a_dmap, + const Vector & a_ba, + const Vector & a_rr) + : amrex::ParticleContainerPureSoA(a_geom, a_dmap, a_ba, a_rr) + { + for (int i = 0; i < num_runtime_real; ++i) + { + AddRealComp(true); + } + for (int i = 0; i < num_runtime_int; ++i) + { + AddIntComp(true); + } + } + + void InitParticles (const amrex::IntVect& a_num_particles_per_cell) + { + BL_PROFILE("InitParticles"); + + const int lev = 0; // only add particles on level 0 + const Real* dx = Geom(lev).CellSize(); + const Real* plo = Geom(lev).ProbLo(); + + const int num_ppc = AMREX_D_TERM( a_num_particles_per_cell[0], + *a_num_particles_per_cell[1], + *a_num_particles_per_cell[2]); + + for(MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi) + { + const Box& tile_box = mfi.tilebox(); + + std::array, NR> host_real; + std::array, NI> host_int; + + std::vector > host_runtime_real(NumRuntimeRealComps()); + std::vector > host_runtime_int(NumRuntimeIntComps()); + + for (IntVect iv = tile_box.smallEnd(); iv <= tile_box.bigEnd(); tile_box.next(iv)) + { + for (int i_part=0; i_part(id)); + host_int[1].push_back(ParallelDescriptor::MyProc()); + host_real[0].push_back(static_cast (plo[0] + (iv[0] + r[0])*dx[0])); +#if AMREX_SPACEDIM > 1 + host_real[1].push_back(static_cast (plo[1] + (iv[1] + r[1])*dx[1])); +#endif +#if AMREX_SPACEDIM > 2 + host_real[2].push_back(static_cast (plo[2] + (iv[2] + r[2])*dx[2])); +#endif + + for (int i = AMREX_SPACEDIM; i < NR; ++i) + host_real[i].push_back(static_cast(id)); + for (int i = 2; i < NI; ++i) + host_int[i].push_back(static_cast(id)); + for (int i = 0; i < NumRuntimeRealComps(); ++i) + host_runtime_real[i].push_back(static_cast(id)); + for (int i = 0; i < NumRuntimeIntComps(); ++i) + host_runtime_int[i].push_back(static_cast(id)); + } + } + + auto& particle_tile = DefineAndReturnParticleTile(lev, mfi.index(), mfi.LocalTileIndex()); + auto old_size = particle_tile.GetArrayOfStructs().size(); + auto new_size = old_size + host_real[0].size(); + particle_tile.resize(new_size); + + auto& soa = particle_tile.GetStructOfArrays(); + for (int i = 0; i < NR; ++i) + { + Gpu::copyAsync(Gpu::hostToDevice, + host_real[i].begin(), + host_real[i].end(), + soa.GetRealData(i).begin() + old_size); + } + + for (int i = 0; i < NI; ++i) + { + Gpu::copyAsync(Gpu::hostToDevice, + host_int[i].begin(), + host_int[i].end(), + soa.GetIntData(i).begin() + old_size); + } + for (int i = 0; i < NumRuntimeRealComps(); ++i) + { + Gpu::copyAsync(Gpu::hostToDevice, + host_runtime_real[i].begin(), + host_runtime_real[i].end(), + soa.GetRealData(NR+i).begin() + old_size); + } + + for (int i = 0; i < NumRuntimeIntComps(); ++i) + { + Gpu::copyAsync(Gpu::hostToDevice, + host_runtime_int[i].begin(), + host_runtime_int[i].end(), + soa.GetIntData(NI+i).begin() + old_size); + } + + Gpu::streamSynchronize(); + } + + Redistribute(); + } +}; + +struct TestParams +{ + IntVect size; + int max_grid_size; + int num_ppc; + int is_periodic; + int nlevs; +}; + +void testRedistribute(); + +int main (int argc, char* argv[]) +{ + amrex::Initialize(argc,argv); + + amrex::Print() << "Running redistribute test \n"; + testRedistribute(); + + amrex::Finalize(); +} + +void get_test_params(TestParams& params, const std::string& prefix) +{ + ParmParse pp(prefix); + pp.get("size", params.size); + pp.get("max_grid_size", params.max_grid_size); + pp.get("num_ppc", params.num_ppc); + pp.get("is_periodic", params.is_periodic); + pp.get("nlevs", params.nlevs); + pp.query("num_runtime_real", num_runtime_real); + pp.query("num_runtime_int", num_runtime_int); +} + +void testRedistribute () +{ + BL_PROFILE("testSENSEI_insitu"); + TestParams params; + get_test_params(params, "insitu"); + + int is_per[BL_SPACEDIM]; + for (int & d : is_per) + d = params.is_periodic; + + Vector rr(params.nlevs-1); + for (int lev = 1; lev < params.nlevs; lev++) + rr[lev-1] = IntVect(AMREX_D_DECL(2,2,2)); + + RealBox real_box; + for (int n = 0; n < BL_SPACEDIM; n++) + { + real_box.setLo(n, 0.0); + real_box.setHi(n, params.size[n]); + } + + IntVect domain_lo(AMREX_D_DECL(0, 0, 0)); + IntVect domain_hi(AMREX_D_DECL(params.size[0]-1,params.size[1]-1,params.size[2]-1)); + const Box base_domain(domain_lo, domain_hi); + + Vector geom(params.nlevs); + geom[0].define(base_domain, &real_box, CoordSys::cartesian, is_per); + for (int lev = 1; lev < params.nlevs; lev++) { + geom[lev].define(amrex::refine(geom[lev-1].Domain(), rr[lev-1]), + &real_box, CoordSys::cartesian, is_per); + } + + Vector ba(params.nlevs); + Vector dm(params.nlevs); + auto lo = IntVect(AMREX_D_DECL(0, 0, 0)); + IntVect size = params.size; + for (int lev = 0; lev < params.nlevs; ++lev) + { + ba[lev].define(Box(lo, lo+params.size-1)); + ba[lev].maxSize(params.max_grid_size); + dm[lev].define(ba[lev]); + lo += size/2; + size *= 2; + } + + TestParticleContainer pc(geom, dm, ba, rr); + + int npc = params.num_ppc; + auto nppc = IntVect(AMREX_D_DECL(npc, npc, npc)); + + amrex::Print() << "About to initialize particles \n"; + + pc.InitParticles(nppc); + + auto *insitu_bridge = new AmrParticleInSituBridge; + + if (insitu_bridge->initialize()) { + amrex::ErrorStream() << "Failed to initialize the in situ bridge." << std::endl; + amrex::Abort(); + } + + // define specifications for fields on particles + std::map> rStructs = {{"u", {0,1,2}}}; + std::map iStructs; + std::map> rArrays; + std::map iArrays; + + if (insitu_bridge->update(&amr, tracers, rStructs)) { + amrex::ErrorStream() << "Failed to update the in situ bridge." << std::endl; + amrex::Abort(); + } + + if (insitu_bridge->finalize()) { + amrex::ErrorStream() << "Failed to finalize the in situ bridge." << std::endl; + } + + delete insitu_bridge; + + // the way this test is set up, if we make it here we pass + amrex::Print() << "pass \n"; +}