Resolve 'Delete particles that move in the opposite direction than the bunch.'

ippl/src/Particle/AbstractParticle.h

@@ -59,9 +59,6 @@ public:
     virtual void setTotalNum(size_t n) = 0;
     virtual void setLocalNum(size_t n) = 0;
 
-    virtual unsigned int getMinimumNumberOfParticlesPerCore() const = 0;
-    virtual void setMinimumNumberOfParticlesPerCore(unsigned int n) = 0;
-
     virtual Layout_t& getLayout() = 0;
     virtual const Layout_t& getLayout() const = 0;
 

ippl/src/Particle/IpplParticleBase.h

@@ -143,9 +143,6 @@ public:
     typedef attrib_container_t::iterator      attrib_iterator;
     typedef ParticleAttribBase::SortList_t    SortList_t;
 
-    // useful constants
-    unsigned int MIN_NUM_PART_PER_CORE;
-
     // our position, and our global ID's
     ParticlePos_t   R;
     ParticleIndex_t ID;
@@ -155,7 +152,6 @@ public:
     // If this constructor is used, the user must call 'initialize' with
     // a layout object in order to use this.
     IpplParticleBase() :
-        MIN_NUM_PART_PER_CORE(0),
         Layout(NULL),
         TotalNum(0),
         LocalNum(0),
@@ -165,7 +161,6 @@ public:
 
     // constructor 2: arguments = layout to use.
     IpplParticleBase(PLayout *layout) :
-        MIN_NUM_PART_PER_CORE(0),
         Layout(layout),
         TotalNum(0),
         LocalNum(0),
@@ -204,10 +199,6 @@ public:
     void setTotalNum(size_t n) { TotalNum = n; }
     void setLocalNum(size_t n) { LocalNum = n; }
 
-    unsigned int getMinimumNumberOfParticlesPerCore() const { return MIN_NUM_PART_PER_CORE; };
-    void setMinimumNumberOfParticlesPerCore(unsigned int n) { MIN_NUM_PART_PER_CORE=n; };
-
-
     // get the layout manager
     PLayout& getLayout() { return *Layout; }
     const PLayout& getLayout() const { return *Layout; }
@@ -362,7 +353,7 @@ protected:
     // Each destroy is stored as a pair of unsigned ints, the particle
     // index I to start at and the number of particles M to destroy.
     std::vector< std::pair<size_t,size_t> > DestroyList;
-    
+
 private:
     // our layout object, which we delete in our destructor
     PLayout *Layout;
@@ -380,7 +371,7 @@ private:
 
     // unique particle ID number generation value
     unsigned NextID;
-    
+
     //
     // private methods
     //
@@ -401,4 +392,4 @@ private:
  * $RCSfile: IpplParticleBase.h,v $   $Author: adelmann $
  * $Revision: 1.1.1.1 $   $Date: 2003/01/23 07:40:28 $
  * IPPL_VERSION_ID: $Id: IpplParticleBase.h,v 1.1.1.1 2003/01/23 07:40:28 adelmann Exp $
- ***************************************************************************/
+ ***************************************************************************/
\ No newline at end of file

ippl/src/Particle/ParticleSpatialLayout.h

@@ -1094,7 +1094,6 @@ protected:
 
         std::multimap<unsigned, unsigned> p2n; //<node ID, particle ID>
 
-        int minParticlesPerNode = PData.getMinimumNumberOfParticlesPerCore();
         int particlesLeft = LocalNum;
         bool responsibleNodeNotFound = false;
         for (unsigned int ip=0; ip<LocalNum; ++ip)
@@ -1113,9 +1112,6 @@ protected:
             if (found)
                 continue;
 
-            if (particlesLeft <= minParticlesPerNode)
-	        break; //leave atleast minimum number of particles per core
-
             typename RegionLayout<T,Dim,Mesh>::touch_range_dv touchingVN = RLayout.touch_range_rdv(pLoc);
 
             //external location
@@ -1224,7 +1220,6 @@ protected:
 
         std::multimap<unsigned, unsigned> p2n; //<node ID, particle ID>
 
-        int minParticlesPerNode = PData.getMinimumNumberOfParticlesPerCore();
         int particlesLeft = LocalNum;
         bool responsibleNodeNotFound = false;
         for (unsigned int ip=0; ip<LocalNum; ++ip)
@@ -1246,9 +1241,6 @@ protected:
             if (found)
                 continue;
 
-            if (particlesLeft <= minParticlesPerNode)
-                continue; //leave atleast minimum number of particles per core
-
             typename RegionLayout<T,Dim,Mesh>::touch_range_dv touchingVN = RLayout.touch_range_rdv(pLoc);
 
             //external location
@@ -1338,4 +1330,4 @@ protected:
 
 #include "Particle/ParticleSpatialLayout.hpp"
 
-#endif // PARTICLE_SPATIAL_LAYOUT_H
+#endif // PARTICLE_SPATIAL_LAYOUT_H
\ No newline at end of file

src/Algorithms/ParallelCyclotronTracker.cpp

@@ -3345,7 +3345,7 @@ void ParallelCyclotronTracker::computeSpaceChargeFields_m() {
         //itsBunch_m->R *= Vector_t(0.001); // mm --> m
 
         if ((step_m + 1) % Options::boundpDestroyFreq == 0)
-            itsBunch_m->boundp_destroy();
+            itsBunch_m->boundp_destroyCycl();
         else
             itsBunch_m->boundp();
 

src/Classic/Algorithms/PartBunchBase.h

@@ -180,7 +180,7 @@ public:
     virtual void boundp();
 
     /** delete particles which are too far away from the center of beam*/
-    void boundp_destroy();
+    void boundp_destroyCycl();
 
     /** This is only temporary in order to get the collimator and pepperpot working */
     size_t boundp_destroyT();
@@ -514,9 +514,6 @@ public:
     size_t getDestroyNum() const;
     size_t getGhostNum() const;
 
-    unsigned int getMinimumNumberOfParticlesPerCore() const;
-    void setMinimumNumberOfParticlesPerCore(unsigned int n);
-
     ParticleLayout<T, Dim>& getLayout();
     const ParticleLayout<T, Dim>& getLayout() const;
 
@@ -607,9 +604,6 @@ public:
     // For AMTS integrator in OPAL-T
     double dtScInit_m, deltaTau_m;
 
-    /// if a local node has less than 2 particles  lowParticleCount_m == true
-    bool lowParticleCount_m;
-
     // get 2nd order momentum matrix
     FMatrix<double, 2 * Dim, 2 * Dim> getSigmaMatrix();
 
@@ -744,4 +738,4 @@ typename PartBunchBase<T, Dim>::ConstIterator end(PartBunchBase<T, Dim> const& b
 
 #include "PartBunchBase.hpp"
 
-#endif
\ No newline at end of file
+#endif

src/Classic/Algorithms/PartBunchBase.hpp

@@ -42,7 +42,6 @@ PartBunchBase<T, Dim>::PartBunchBase(AbstractParticle<T, Dim>* pb, const PartDat
     : R(*(pb->R_p)),
       ID(*(pb->ID_p)),
       pbin_m(nullptr),
-      lowParticleCount_m(false),
       pmsg_m(nullptr),
       f_stream(nullptr),
       fixed_grid(false),
@@ -578,9 +577,7 @@ void PartBunchBase<T, Dim>::boundp() {
 
 
 template <class T, unsigned Dim>
-void PartBunchBase<T, Dim>::boundp_destroy() {
-
-    Inform gmsgAll("boundp_destroy ", INFORM_ALL_NODES);
+void PartBunchBase<T, Dim>::boundp_destroyCycl() {
 
     Vector_t len;
     const int dimIdx = 3;
@@ -690,52 +687,36 @@ void PartBunchBase<T, Dim>::boundp_destroy() {
 template <class T, unsigned Dim>
 size_t PartBunchBase<T, Dim>::boundp_destroyT() {
 
-    const unsigned int minNumParticlesPerCore = getMinimumNumberOfParticlesPerCore();
-
     this->updateDomainLength(nr_m);
 
-    std::unique_ptr<size_t[]> tmpbinemitted;
+    std::vector<size_t> tmpbinemitted;
 
     boundp();
 
     size_t ne = 0;
     const size_t localNum = getLocalNum();
 
-    if (weHaveEnergyBins()) {
-        tmpbinemitted = std::unique_ptr<size_t[]>(new size_t[getNumberOfEnergyBins()]);
-        for (int i = 0; i < getNumberOfEnergyBins(); i++) {
-            tmpbinemitted[i] = 0;
-        }
-        for (unsigned int i = 0; i < localNum; i++) {
-            if (Bin[i] < 0) {
-                ne++;
-                destroy(1, i);
-            } else {
-                tmpbinemitted[Bin[i]]++;
-            }
-        }
-    } else {
-        for (unsigned int i = 0; i < localNum; i++) {
-            if ((Bin[i] < 0) && ((localNum - ne) > minNumParticlesPerCore)) {   // need in minimum x particles per node
-                ne++;
-                destroy(1, i);
-            }
+    double rzmean = momentsComputer_m.getMeanPosition()(2);
+    double rzrms = momentsComputer_m.getStandardDeviationPosition()(2);
+    const bool haveEnergyBins = weHaveEnergyBins();
+    if (haveEnergyBins) {
+        tmpbinemitted.resize(getNumberOfEnergyBins(), 0.0);
+    }
+    for (unsigned int i = 0; i < localNum; i++) {
+        if (Bin[i] < 0 || (Options::remotePartDel > 0 && std::abs(R[i](2) - rzmean) < Options::remotePartDel * rzrms)) {
+            ne++;
+            destroy(1, i);
+        } else if (haveEnergyBins) {
+            tmpbinemitted[Bin[i]]++;
         }
-        lowParticleCount_m = ((localNum - ne) <= minNumParticlesPerCore);
-        reduce(lowParticleCount_m, lowParticleCount_m, OpOr());
     }
 
-    if (lowParticleCount_m) {
-        Inform m ("boundp_destroyT a) ", INFORM_ALL_NODES);
-        m << level3 << "Warning low number of particles on some cores localNum= "
-          << localNum << " ne= " << ne << " NLocal= " << getLocalNum() << endl;
-    } else {
-        boundp();
-    }
+    boundp();
+
     calcBeamParameters();
     gatherLoadBalanceStatistics();
 
-    if (weHaveEnergyBins()) {
+    if (haveEnergyBins) {
         const int lastBin = dist_m->getLastEmittedEnergyBin();
         for (int i = 0; i < lastBin; i++) {
             binemitted_m[i] = tmpbinemitted[i];
@@ -749,7 +730,6 @@ size_t PartBunchBase<T, Dim>::boundp_destroyT() {
 template <class T, unsigned Dim>
 size_t PartBunchBase<T, Dim>::destroyT() {
 
-    const unsigned int minNumParticlesPerCore = getMinimumNumberOfParticlesPerCore();
     std::unique_ptr<size_t[]> tmpbinemitted;
 
     const size_t localNum = getLocalNum();
@@ -772,14 +752,10 @@ size_t PartBunchBase<T, Dim>::destroyT() {
         Inform dmsg("destroy: ", INFORM_ALL_NODES);
         for (size_t i = 0; i < localNum; i++) {
             if ((Bin[i] < 0)) {
-                if ((localNum - ne) > minNumParticlesPerCore) {   // need in minimum x particles per node
-                    ne++;
-                    destroy(1, i);
-                }
+                ne++;
+                destroy(1, i);
             }
         }
-        lowParticleCount_m = ((localNum - ne) <= minNumParticlesPerCore);
-        reduce(lowParticleCount_m, lowParticleCount_m, OpOr());
     }
 
     if (ne > 0) {
@@ -1889,9 +1865,6 @@ void PartBunchBase<T, Dim>::setup(AbstractParticle<T, Dim>* pb) {
     globalPartPerNode_m = std::unique_ptr<size_t[]>(new size_t[Ippl::getNodes()]);
 
     pmsg_m.release();
-
-    // set the default IPPL behaviour
-    setMinimumNumberOfParticlesPerCore(0);
 }
 
 
@@ -1926,16 +1899,6 @@ void PartBunchBase<T, Dim>::setLocalNum(size_t n) {
     pbase_m->setLocalNum(n);
 }
 
-template <class T, unsigned Dim>
-unsigned int PartBunchBase<T, Dim>::getMinimumNumberOfParticlesPerCore() const {
-    return pbase_m->getMinimumNumberOfParticlesPerCore();
-}
-
-template <class T, unsigned Dim>
-void PartBunchBase<T, Dim>::setMinimumNumberOfParticlesPerCore(unsigned int n) {
-    pbase_m->setMinimumNumberOfParticlesPerCore(n);
-}
-
 template <class T, unsigned Dim>
 ParticleLayout<T, Dim> & PartBunchBase<T, Dim>::getLayout() {
     return pbase_m->getLayout();

src/Classic/Solvers/ScatteringPhysics.cpp

@@ -547,10 +547,8 @@ void ScatteringPhysics::copyFromBunch(PartBunchBase<double, 3>* bunch,
 
     size_t ne = 0;
     std::set<size_t> partsToDel;
-    const unsigned int minNumOfParticlesPerCore = bunch->getMinimumNumberOfParticlesPerCore();
     for (size_t i = 0; i < nL; ++ i) {
         if ((bunch->Bin[i] == -1 || bunch->Bin[i] == 1) &&
-            ((nL - ne) > minNumOfParticlesPerCore) &&
             hitTester_m->checkHit(bunch->R[i]))
         {
             // adjust the time step for those particles that enter the material