diff --git a/src_deankow/multilevel/AmrCoreAdv.cpp b/src_deankow/multilevel/AmrCoreAdv.cpp index 395ff57a..45ddea8d 100644 --- a/src_deankow/multilevel/AmrCoreAdv.cpp +++ b/src_deankow/multilevel/AmrCoreAdv.cpp @@ -267,6 +267,7 @@ AmrCoreAdv::RemakeLevel (int lev, Real time, const BoxArray& ba, const int ncomp = phi_new[lev].nComp(); const int ng = phi_new[lev].nGrow(); + BoxArray old_fine_ba = phi_old[1].boxArray(); amrex::Print() << " REGRIDDING: NEW GRIDS AT LEVEL " << lev << " " << ba << std::endl; if (lev == 1) { @@ -291,7 +292,7 @@ AmrCoreAdv::RemakeLevel (int lev, Real time, const BoxArray& ba, #ifdef AMREX_PARTICLES if (lev == 1) { - particleData.regrid_particles(grown_fba); + particleData.regrid_particles(grown_fba, old_fine_ba, phi_new[1]); } #endif } diff --git a/src_deankow/multilevel/ParticleData.H b/src_deankow/multilevel/ParticleData.H index c3e5af28..acc0524e 100644 --- a/src_deankow/multilevel/ParticleData.H +++ b/src_deankow/multilevel/ParticleData.H @@ -76,7 +76,7 @@ struct ParticleData { stochastic_particles->RemoveParticlesNotInBA(phi_fine.boxArray()); // ********************************************************************** - // Redistribute so every particle is now in the "correct" place + // Redistribute so every particle is now in the "correct" place // ********************************************************************** stochastic_particles->Redistribute(); amrex::Print() << "Initialized " << stochastic_particles->TotalNumberOfParticles() << @@ -93,12 +93,14 @@ struct ParticleData { } } - void regrid_particles(const amrex::BoxArray& grown_fba) + void regrid_particles(const amrex::BoxArray& new_fine_ba, const amrex::BoxArray& old_fine_ba, amrex::MultiFab& phi_fine) { - amrex::DistributionMapping grown_dm { grown_fba, amrex::ParallelDescriptor::NProcs() }; - stochastic_particles->SetParticleDistributionMap(1,grown_dm); - stochastic_particles->SetParticleBoxArray(1,grown_fba); + amrex::DistributionMapping new_dm { new_fine_ba, amrex::ParallelDescriptor::NProcs() }; + stochastic_particles->SetParticleDistributionMap(1, new_dm); + stochastic_particles->SetParticleBoxArray(1, new_fine_ba); stochastic_particles->Redistribute(); + stochastic_particles->RemoveParticlesAtLevel(0); + stochastic_particles->AddParticles(phi_fine, old_fine_ba); } void Checkpoint (const std::string& filename) const @@ -117,18 +119,18 @@ struct ParticleData { } } - void advance_particles(int lev, amrex::Real dt_lev, amrex::Real cell_vol, + void advance_particles(int lev, amrex::Real dt_lev, amrex::Real cell_vol, amrex::MultiFab& phi_crse_old, amrex::MultiFab& phi_crse_new, amrex::MultiFab& phi_fine) { // Update stochastic particles if (use_particles) { - AMREX_ALWAYS_ASSERT(lev == 1); + AMREX_ALWAYS_ASSERT(lev == 1); const amrex::Geometry& geom_fine = stochastic_particles->GetParGDB()->Geom(lev); // ********************************************************************** - // Fill a temporary phi_grown with phi from "coarse" level + // Fill a temporary phi_grown with phi from "coarse" level // ********************************************************************** amrex::MultiFab phi_grown(stochastic_particles->ParticleBoxArray(1), stochastic_particles->ParticleDistributionMap(1),1,0); @@ -163,7 +165,7 @@ struct ParticleData { phi_grown.setVal(0.); stochastic_particles->RefluxFineToCrse(phi_fine.boxArray(), phi_grown); phi_grown.mult(1./cell_vol); - phi_crse_new.ParallelAdd(phi_grown,0,0,1); + phi_crse_new.ParallelAdd(phi_grown,0,0,1); // Remove particles that are not in the valid boxArray used by phi_fine stochastic_particles->RemoveParticlesNotInBA(phi_fine.boxArray()); @@ -176,7 +178,7 @@ struct ParticleData { stochastic_particles->RefluxCrseToFine(phi_fine.boxArray(), phi_grown2); phi_grown2.SumBoundary(geom_fine.periodicity()); phi_grown2.mult(1./cell_vol); - phi_crse_new.ParallelAdd(phi_grown2,0,0,1); + phi_crse_new.ParallelAdd(phi_grown2,0,0,1); // ********************************************************************** // Define phi at level 1 from the particle count diff --git a/src_deankow/multilevel/StochasticPC.cpp b/src_deankow/multilevel/StochasticPC.cpp index 9677778b..3626d589 100644 --- a/src_deankow/multilevel/StochasticPC.cpp +++ b/src_deankow/multilevel/StochasticPC.cpp @@ -36,6 +36,13 @@ StochasticPC:: AddParticles (MultiFab& phi_fine, const BoxArray& ba_to_exclude) if (ba_to_exclude.contains(tile_box)) {continue;} + if (!m_reflux_particle_locator.isValid(ba_to_exclude)) { + m_reflux_particle_locator.build(ba_to_exclude, Geom(lev)); + } + m_reflux_particle_locator.setGeometry(Geom(lev)); + + auto assign_grid = m_reflux_particle_locator.getGridAssignor(); + const Array4& phi_arr = phi_fine.const_array(mfi); // count the number of particles to create in each cell @@ -45,6 +52,7 @@ StochasticPC:: AddParticles (MultiFab& phi_fine, const BoxArray& ba_to_exclude) amrex::ParallelForRNG(tile_box, [=] AMREX_GPU_DEVICE (int i, int j, int k, amrex::RandomEngine const& engine) noexcept { + if (assign_grid(IntVect(AMREX_D_DECL(i, j, k))) >= 0) {return;} Real rannum = amrex::Random(engine); int npart_in_cell = int(phi_arr(i,j,k,0)*cell_vol+rannum); pcount[flat_index(i, j, k)] += npart_in_cell; @@ -129,13 +137,22 @@ StochasticPC::RemoveParticlesNotInBA (const BoxArray& ba_to_keep) const int np = aos.numParticles(); auto *pstruct = aos().data(); - if (!ba_to_keep.contains(pti.tilebox())) { - amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE (int i) - { - ParticleType& p = pstruct[i]; + if (ba_to_keep.contains(pti.tilebox())) {continue;} + + if (!m_reflux_particle_locator.isValid(ba_to_keep)) { + m_reflux_particle_locator.build(ba_to_keep, Geom(lev)); + } + m_reflux_particle_locator.setGeometry(Geom(lev)); + + auto assign_grid = m_reflux_particle_locator.getGridAssignor(); + + amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE (int i) + { + ParticleType& p = pstruct[i]; + if (assign_grid(p) < 0) { p.id() = -1; - }); -} + } + }); } Redistribute(); } @@ -249,8 +266,7 @@ StochasticPC::RefluxCrseToFine (const BoxArray& ba_to_keep, MultiFab& phi_for_re const int np = aos.numParticles(); auto *pstruct = aos().data(); - if (ba_to_keep.contains(pti.tilebox())) - { + if (ba_to_keep.contains(pti.tilebox())) { Array4 phi_arr = phi_for_reflux.array(gid); amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE (int i)