aboutsummaryrefslogtreecommitdiff
path: root/Source
diff options
context:
space:
mode:
Diffstat (limited to 'Source')
-rw-r--r--Source/Particles/Collision/CollisionType.cpp40
-rw-r--r--Source/Particles/PhysicalParticleContainer.cpp3
-rw-r--r--Source/Particles/Resampling/CMakeLists.txt1
-rw-r--r--Source/Particles/Resampling/LevelingThinning.H27
-rw-r--r--Source/Particles/Resampling/LevelingThinning.cpp96
-rw-r--r--Source/Particles/Resampling/Make.package1
-rw-r--r--Source/Particles/Resampling/Resampling.H7
-rw-r--r--Source/Particles/Resampling/Resampling.cpp4
-rw-r--r--Source/Utils/CMakeLists.txt1
-rw-r--r--Source/Utils/Make.package1
-rw-r--r--Source/Utils/ParticleUtils.H31
-rw-r--r--Source/Utils/ParticleUtils.cpp52
12 files changed, 217 insertions, 47 deletions
diff --git a/Source/Particles/Collision/CollisionType.cpp b/Source/Particles/Collision/CollisionType.cpp
index 8c5fac875..fc30f151a 100644
--- a/Source/Particles/Collision/CollisionType.cpp
+++ b/Source/Particles/Collision/CollisionType.cpp
@@ -7,7 +7,8 @@
#include "CollisionType.H"
#include "ShuffleFisherYates.H"
#include "ElasticCollisionPerez.H"
-#include <WarpX.H>
+#include "Utils/ParticleUtils.H"
+#include "WarpX.H"
CollisionType::CollisionType(
const std::vector<std::string>& species_names,
@@ -47,42 +48,7 @@ using ParticleTileType = WarpXParticleContainer::ParticleTileType;
using ParticleBins = DenseBins<ParticleType>;
using index_type = ParticleBins::index_type;
-namespace {
-
- /* Find the particles and count the particles that are in each cell.
- Note that this does *not* rearrange particle arrays */
- ParticleBins
- findParticlesInEachCell( int const lev, MFIter const& mfi,
- ParticleTileType const& ptile) {
-
- // Extract particle structures for this tile
- int const np = ptile.numParticles();
- ParticleType const* particle_ptr = ptile.GetArrayOfStructs()().data();
-
- // Extract box properties
- Geometry const& geom = WarpX::GetInstance().Geom(lev);
- Box const& cbx = mfi.tilebox(IntVect::TheZeroVector()); //Cell-centered box
- const auto lo = lbound(cbx);
- const auto dxi = geom.InvCellSizeArray();
- const auto plo = geom.ProbLoArray();
-
- // Find particles that are in each cell;
- // results are stored in the object `bins`.
- ParticleBins bins;
- bins.build(np, particle_ptr, cbx,
- // Pass lambda function that returns the cell index
- [=] AMREX_GPU_HOST_DEVICE (const ParticleType& p) noexcept -> IntVect
- {
- return IntVect(AMREX_D_DECL(
- static_cast<int>((p.pos(0)-plo[0])*dxi[0] - lo.x),
- static_cast<int>((p.pos(1)-plo[1])*dxi[1] - lo.y),
- static_cast<int>((p.pos(2)-plo[2])*dxi[2] - lo.z)));
- });
-
- return bins;
- }
-
-}
+using namespace ParticleUtils;
/** Perform all binary collisions within a tile
*
diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp
index 68900463b..1ac048cec 100644
--- a/Source/Particles/PhysicalParticleContainer.cpp
+++ b/Source/Particles/PhysicalParticleContainer.cpp
@@ -2003,11 +2003,12 @@ void PhysicalParticleContainer::resample (const Resampling& resampler, const int
if (resampler.triggered(timestep, global_numparts))
{
+ amrex::Print() << "Resampling " << species_name << ".\n";
for (int lev = 0; lev <= maxLevel(); lev++)
{
for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti)
{
- resampler(pti);
+ resampler(pti, lev, this);
}
}
}
diff --git a/Source/Particles/Resampling/CMakeLists.txt b/Source/Particles/Resampling/CMakeLists.txt
index 7bba14166..0b44f7a5d 100644
--- a/Source/Particles/Resampling/CMakeLists.txt
+++ b/Source/Particles/Resampling/CMakeLists.txt
@@ -2,4 +2,5 @@ target_sources(WarpX
PRIVATE
Resampling.cpp
ResamplingTrigger.cpp
+ LevelingThinning.cpp
)
diff --git a/Source/Particles/Resampling/LevelingThinning.H b/Source/Particles/Resampling/LevelingThinning.H
index 77d1b967a..920fa07c9 100644
--- a/Source/Particles/Resampling/LevelingThinning.H
+++ b/Source/Particles/Resampling/LevelingThinning.H
@@ -9,15 +9,32 @@
#include "Resampling.H"
-#include <AMReX_Print.H>
-
/**
- * \brief This class will soon implement the leveling thinning algorithm, but for now it is just
- * a dummy class that prints "Resampling." used to test the resampling structure.
+ * \brief This class implements the leveling thinning algorithm as defined in Muraviev, A., et al.
+ * arXiv:2006.08593 (2020).
+ * The main steps of the algorithm are the following: for every cell we calculate a level weight,
+ * defined by the average weight of the species particles in that cell multiplied by the target
+ * ratio. Then, particles with a weight lower than the level weight are either removed, with a
+ * probability 1 - particle_weight/level_weight, or have their weight set to the level weight.
*/
class LevelingThinning: public ResamplingAlgorithm {
public:
- void operator() (WarpXParIter& /*pti*/) override final {amrex::Print() << "Resampling. \n";}
+ /**
+ * \brief Constructor of the LevelingThinning class
+ */
+ LevelingThinning ();
+
+ /**
+ * \brief A method that performs leveling thinning for the considered species.
+ *
+ * @param[in] pti WarpX particle iterator of the particles to resample.
+ * @param[in] lev the index of the refinement level.
+ * @param[in] pc a pointer to the particle container.
+ */
+ void operator() (WarpXParIter& pti, const int lev, WarpXParticleContainer * const pc) const override final;
+
+private:
+ amrex::Real m_target_ratio = amrex::Real(1.5);
};
diff --git a/Source/Particles/Resampling/LevelingThinning.cpp b/Source/Particles/Resampling/LevelingThinning.cpp
new file mode 100644
index 000000000..d3ffc9530
--- /dev/null
+++ b/Source/Particles/Resampling/LevelingThinning.cpp
@@ -0,0 +1,96 @@
+/* Copyright 2019-2020 Neil Zaim
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+#include "LevelingThinning.H"
+#include "Utils/ParticleUtils.H"
+
+#include <AMReX_Particles.H>
+
+LevelingThinning::LevelingThinning ()
+{
+ using namespace amrex::literals;
+
+ amrex::ParmParse pp("resampling_algorithm");
+ pp.query("target_ratio", m_target_ratio);
+ AMREX_ALWAYS_ASSERT_WITH_MESSAGE( m_target_ratio > 0._rt,
+ "Resampling target ratio should be strictly greater than 0");
+ if (m_target_ratio <= 1._rt)
+ {
+ amrex::Warning("WARNING: target ratio for leveling thinning is smaller or equal to one."
+ " It is possible that no particle will be removed during resampling");
+ }
+}
+
+void LevelingThinning::operator() (WarpXParIter& pti, const int lev,
+ WarpXParticleContainer * const pc) const
+{
+ using namespace amrex::literals;
+
+ auto& ptile = pc->ParticlesAt(lev, pti);
+ auto& soa = ptile.GetStructOfArrays();
+ amrex::ParticleReal * const AMREX_RESTRICT w = soa.GetRealData(PIdx::w).data();
+ WarpXParticleContainer::ParticleType * const AMREX_RESTRICT
+ particle_ptr = ptile.GetArrayOfStructs()().data();
+
+ // Using this function means that we must loop over the cells in the ParallelFor. In the case
+ // of the leveling thinning algorithm, it would have possibly been more natural and more
+ // efficient to directly loop over the particles. Nevertheless, this structure with a loop over
+ // the cells is more general and can be readily used to implement almost any other resampling
+ // algorithm.
+ auto bins = ParticleUtils::findParticlesInEachCell(lev, pti, ptile);
+
+ const int n_cells = bins.numBins();
+ const auto indices = bins.permutationPtr();
+ const auto cell_offsets = bins.offsetsPtr();
+
+ const amrex::Real target_ratio = m_target_ratio;
+
+ // Loop over cells
+ amrex::ParallelFor( n_cells,
+ [=] AMREX_GPU_DEVICE (int i_cell) noexcept
+ {
+ // The particles that are in the cell `i_cell` are
+ // given by the `indices[cell_start:cell_stop]`
+ const auto cell_start = cell_offsets[i_cell];
+ const auto cell_stop = cell_offsets[i_cell+1];
+ const int cell_numparts = cell_stop - cell_start;
+
+ // do nothing for cells without particles
+ if (cell_numparts == 0)
+ return;
+ amrex::Real average_weight = 0._rt;
+
+ // First loop over cell particles to compute average particle weight in the cell
+ for (int i = cell_start; i < cell_stop; ++i)
+ {
+ average_weight += w[indices[i]];
+ }
+ average_weight /= cell_numparts;
+
+ const amrex::Real level_weight = average_weight*target_ratio;
+
+
+ // Second loop over cell particles to perform the thinning
+ for (int i = cell_start; i < cell_stop; ++i)
+ {
+ // Particles with weight greater than level_weight are left unchanged
+ if (w[indices[i]] > level_weight) {continue;}
+
+ amrex::Real const random_number = amrex::Random();
+ // Remove particle with probability 1 - particle_weight/level_weight
+ if (random_number > w[indices[i]]/level_weight)
+ {
+ particle_ptr[indices[i]].id() = -1;
+ }
+ // Set particle weight to level weight otherwise
+ else
+ {
+ w[indices[i]] = level_weight;
+ }
+ }
+ }
+ );
+}
diff --git a/Source/Particles/Resampling/Make.package b/Source/Particles/Resampling/Make.package
index 89470af55..e70e3446c 100644
--- a/Source/Particles/Resampling/Make.package
+++ b/Source/Particles/Resampling/Make.package
@@ -1,4 +1,5 @@
CEXE_sources += Resampling.cpp
CEXE_sources += ResamplingTrigger.cpp
+CEXE_sources += LevelingThinning.cpp
VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles/Resampling/
diff --git a/Source/Particles/Resampling/Resampling.H b/Source/Particles/Resampling/Resampling.H
index 83b3ccfd5..de919bfcb 100644
--- a/Source/Particles/Resampling/Resampling.H
+++ b/Source/Particles/Resampling/Resampling.H
@@ -10,6 +10,7 @@
#include "ResamplingTrigger.H"
class WarpXParIter; // forward declaration
+class WarpXParticleContainer; // forward declaration
/**
* \brief An empty base class from which specific resampling algorithms are derived.
@@ -19,7 +20,7 @@ struct ResamplingAlgorithm
/**
* \brief Virtual operator() of the abstract ResamplingAlgorithm class
*/
- virtual void operator() (WarpXParIter& /*pti*/) = 0;
+ virtual void operator() (WarpXParIter& /*pti*/, const int /*lev*/, WarpXParticleContainer */*pc*/) const = 0;
/**
* \brief Virtual destructor of the abstract ResamplingAlgorithm class
@@ -56,8 +57,10 @@ public:
* \brief A method that uses the ResamplingAlgorithm object to perform resampling.
*
* @param[in] pti WarpX particle iterator of the particles to resample.
+ * @param[in] lev the index of the refinement level.
+ * @param[in] pc a pointer to the particle container.
*/
- void operator() (WarpXParIter& pti) const;
+ void operator() (WarpXParIter& pti, const int lev, WarpXParticleContainer * const pc) const;
private:
ResamplingTrigger m_resampling_trigger;
diff --git a/Source/Particles/Resampling/Resampling.cpp b/Source/Particles/Resampling/Resampling.cpp
index a08145915..65332d191 100644
--- a/Source/Particles/Resampling/Resampling.cpp
+++ b/Source/Particles/Resampling/Resampling.cpp
@@ -26,7 +26,7 @@ bool Resampling::triggered (const int timestep, const amrex::Real global_numpart
return m_resampling_trigger.triggered(timestep, global_numparts);
}
-void Resampling::operator() (WarpXParIter& pti) const
+void Resampling::operator() (WarpXParIter& pti, const int lev, WarpXParticleContainer * const pc) const
{
- (*m_resampling_algorithm)(pti);
+ (*m_resampling_algorithm)(pti, lev, pc);
}
diff --git a/Source/Utils/CMakeLists.txt b/Source/Utils/CMakeLists.txt
index ce785cc33..c9da2ae47 100644
--- a/Source/Utils/CMakeLists.txt
+++ b/Source/Utils/CMakeLists.txt
@@ -4,6 +4,7 @@ target_sources(WarpX
CoarsenMR.cpp
Interpolate.cpp
IntervalsParser.cpp
+ ParticleUtils.cpp
RelativeCellPosition.cpp
WarpXAlgorithmSelection.cpp
WarpXMovingWindow.cpp
diff --git a/Source/Utils/Make.package b/Source/Utils/Make.package
index 08b9aa3a9..61ced5d7f 100644
--- a/Source/Utils/Make.package
+++ b/Source/Utils/Make.package
@@ -7,5 +7,6 @@ CEXE_sources += CoarsenMR.cpp
CEXE_sources += Interpolate.cpp
CEXE_sources += IntervalsParser.cpp
CEXE_sources += RelativeCellPosition.cpp
+CEXE_sources += ParticleUtils.cpp
VPATH_LOCATIONS += $(WARPX_HOME)/Source/Utils
diff --git a/Source/Utils/ParticleUtils.H b/Source/Utils/ParticleUtils.H
new file mode 100644
index 000000000..cb25607a3
--- /dev/null
+++ b/Source/Utils/ParticleUtils.H
@@ -0,0 +1,31 @@
+/* Copyright 2019-2020 Neil Zaim, Yinjian Zhao
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+#ifndef WARPX_PARTICLE_UTILS_H_
+#define WARPX_PARTICLE_UTILS_H_
+
+#include "Particles/WarpXParticleContainer.H"
+#include <AMReX_DenseBins.H>
+
+namespace ParticleUtils {
+
+ /**
+ * \brief Find the particles and count the particles that are in each cell. More specifically
+ * this function returns an amrex::DenseBins object containing an offset array and a permutation
+ * array which can be used to loop over all the cells in a tile and apply an algorithm to
+ * particles of a given species present in each cell.
+ * Note that this does *not* rearrange particle arrays.
+ *
+ * @param[in] lev the index of the refinement level.
+ * @param[in] mfi the MultiFAB iterator.
+ * @param[in] ptile the particle tile.
+ */
+ amrex::DenseBins<WarpXParticleContainer::ParticleType>
+ findParticlesInEachCell( int const lev, amrex::MFIter const& mfi,
+ WarpXParticleContainer::ParticleTileType const& ptile);
+}
+
+#endif // WARPX_PARTICLE_UTILS_H_
diff --git a/Source/Utils/ParticleUtils.cpp b/Source/Utils/ParticleUtils.cpp
new file mode 100644
index 000000000..92c03bc30
--- /dev/null
+++ b/Source/Utils/ParticleUtils.cpp
@@ -0,0 +1,52 @@
+/* Copyright 2019-2020 Neil Zaim, Yinjian Zhao
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+#include "ParticleUtils.H"
+#include "WarpX.H"
+
+namespace ParticleUtils {
+
+ using namespace amrex;
+ // Define shortcuts for frequently-used type names
+ using ParticleType = WarpXParticleContainer::ParticleType;
+ using ParticleTileType = WarpXParticleContainer::ParticleTileType;
+ using ParticleBins = DenseBins<ParticleType>;
+ using index_type = ParticleBins::index_type;
+
+ /* Find the particles and count the particles that are in each cell.
+ Note that this does *not* rearrange particle arrays */
+ ParticleBins
+ findParticlesInEachCell( int const lev, MFIter const& mfi,
+ ParticleTileType const& ptile) {
+
+ // Extract particle structures for this tile
+ int const np = ptile.numParticles();
+ ParticleType const* particle_ptr = ptile.GetArrayOfStructs()().data();
+
+ // Extract box properties
+ Geometry const& geom = WarpX::GetInstance().Geom(lev);
+ Box const& cbx = mfi.tilebox(IntVect::TheZeroVector()); //Cell-centered box
+ const auto lo = lbound(cbx);
+ const auto dxi = geom.InvCellSizeArray();
+ const auto plo = geom.ProbLoArray();
+
+ // Find particles that are in each cell;
+ // results are stored in the object `bins`.
+ ParticleBins bins;
+ bins.build(np, particle_ptr, cbx,
+ // Pass lambda function that returns the cell index
+ [=] AMREX_GPU_HOST_DEVICE (const ParticleType& p) noexcept -> IntVect
+ {
+ return IntVect(AMREX_D_DECL(
+ static_cast<int>((p.pos(0)-plo[0])*dxi[0] - lo.x),
+ static_cast<int>((p.pos(1)-plo[1])*dxi[1] - lo.y),
+ static_cast<int>((p.pos(2)-plo[2])*dxi[2] - lo.z)));
+ });
+
+ return bins;
+ }
+
+}
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-11 21:02:38 +0000