Merge branch 'main' of github.com:AMReX-FHD/FHDeX into main

AMReX-FHD · May 22, 2024 · 0382c06 · 0382c06
2 parents 5c8cf2a + dd4d543
commit 0382c06
Show file tree

Hide file tree

Showing 3 changed files with 38 additions and 19 deletions.
diff --git 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
@@ -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
@@ -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<Real const>& 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<Real> phi_arr = phi_for_reflux.array(gid);
 
             amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE (int i)