Commit 26c45c26 authored by snuverink_j's avatar snuverink_j

cleanup PartBins and PartBinsCyc

parent 823161aa
//
// Class PartBins
// Defines a structure to hold energy bins and their associated data
//
// Copyright (c) 2007-2020, Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
#include "Algorithms/PartBins.h" #include "Algorithms/PartBins.h"
#include "Algorithms/PBunchDefs.h" #include "Algorithms/PBunchDefs.h"
#include "Physics/Physics.h" #include "Physics/Physics.h"
#include "Utility/Inform.h"
#include <cfloat> #include <cfloat>
#include <limits> #include <limits>
#include <vector> #include <vector>
#include <cmath> #include <cmath>
PartBins::PartBins(int bins, int sbins) : PartBins::PartBins(int bins, int sbins) :
gamma_m(1.0), gamma_m(1.0),
bins_m(bins), bins_m(bins),
...@@ -14,17 +32,16 @@ PartBins::PartBins(int bins, int sbins) : ...@@ -14,17 +32,16 @@ PartBins::PartBins(int bins, int sbins) :
xmin_m(0.0), xmin_m(0.0),
xmax_m(0.0), xmax_m(0.0),
hBin_m(0.0), hBin_m(0.0),
nemittedBins_m(0), nemittedBins_m(0) {
nemittedSBins_m(0) {
// number of particles in the bins on the local node // number of particles in the bins on the local node
nBin_m = std::unique_ptr<size_t[]>(new size_t[bins_m]); nBin_m = std::unique_ptr<size_t[]>(new size_t[bins_m]);
xbinmin_m = std::unique_ptr<double[]>(new double[bins_m]); xbinmin_m = std::unique_ptr<double[]>(new double[bins_m]);
xbinmax_m = std::unique_ptr<double[]>(new double[bins_m]); xbinmax_m = std::unique_ptr<double[]>(new double[bins_m]);
// flag whether the bin contain particles or not // flag whether the bin contain particles or not
binsEmitted_m = std::unique_ptr<bool[]>(new bool[bins_m]); binsEmitted_m = std::unique_ptr<bool[]>(new bool[bins_m]);
nDelBin_m = std::unique_ptr<size_t[]>(new size_t[bins_m]); nDelBin_m = std::unique_ptr<size_t[]>(new size_t[bins_m]);
for(int i = 0; i < bins_m; i++) { for(int i = 0; i < bins_m; i++) {
...@@ -57,46 +74,6 @@ size_t PartBins::getTotalNumPerBin(int b) { ...@@ -57,46 +74,6 @@ size_t PartBins::getTotalNumPerBin(int b) {
return s; return s;
} }
void PartBins::updateStatus(const int bunchCount, const size_t partInBin) {
// array index of binsEmitted_m[] starts from 0
// nemittedBins_m and bins_m index starts from 1
binsEmitted_m[bunchCount - 1] = true;
size_t NpartInBin = partInBin;
reduce(NpartInBin, NpartInBin, OpAddAssign());
nBin_m[bunchCount - 1] = NpartInBin;
nemittedBins_m++;
}
void PartBins::updateDeletedPartsInBin(size_t countLost[]) {
Inform msg("updateDeletedPartsInBin ");
const int lastEmittedBin = getLastemittedBin();
reduce(&(countLost[0]), &(countLost[0]) + lastEmittedBin, &(countLost[0]), OpAddAssign());
for(int ii = 0; ii < lastEmittedBin; ii++) {
if(countLost[ii] > 0) {
nDelBin_m[ii] = countLost[ii];
msg << "In Bin: " << ii << ", " << nDelBin_m[ii] << " particle(s) lost" << endl;
}
}
}
void PartBins::updatePartInBin(size_t countLost[]) {
Inform msg0("updatePartInBin ");
for(int ii = 0; ii < nemittedBins_m; ii++) {
msg0 << "In Bin: " << ii << ", " << nBin_m[ii] << " particles " << endl;
}
reduce(&(countLost[0]), &(countLost[0]) + nemittedBins_m, &(countLost[0]), OpAddAssign());
for(int ii = 0; ii < nemittedBins_m; ii++) {
if(countLost[ii] > 0) {
nBin_m[ii] -= countLost[ii];
msg0 << "In Bin: " << ii << ", " << countLost[ii] << " particle(s) lost" << endl;
}
}
}
void PartBins::updatePartInBin_cyc(size_t countLost[]) { void PartBins::updatePartInBin_cyc(size_t countLost[]) {
...@@ -106,16 +83,6 @@ void PartBins::updatePartInBin_cyc(size_t countLost[]) { ...@@ -106,16 +83,6 @@ void PartBins::updatePartInBin_cyc(size_t countLost[]) {
} }
} }
void PartBins::resetPartInBin(size_t newPartNum[]) {
reduce(&(newPartNum[0]), &(newPartNum[0]) + nemittedBins_m, &(newPartNum[0]), OpAddAssign());
for(int ii = 0; ii < nemittedBins_m; ii++) {
nBin_m[ii] = newPartNum[ii];
INFOMSG("After reset Bin: " << ii << ", particle(s): " << newPartNum[ii] << endl);
}
}
void PartBins::resetPartInBin_cyc(size_t newPartNum[], int maxbinIndex) { void PartBins::resetPartInBin_cyc(size_t newPartNum[], int maxbinIndex) {
reduce(maxbinIndex, maxbinIndex, OpMaxAssign()); reduce(maxbinIndex, maxbinIndex, OpMaxAssign());
nemittedBins_m = maxbinIndex + 1; nemittedBins_m = maxbinIndex + 1;
...@@ -127,7 +94,6 @@ void PartBins::resetPartInBin_cyc(size_t newPartNum[], int maxbinIndex) { ...@@ -127,7 +94,6 @@ void PartBins::resetPartInBin_cyc(size_t newPartNum[], int maxbinIndex) {
} }
PartBins::~PartBins() { PartBins::~PartBins() {
tmppart_m.clear(); tmppart_m.clear();
isEmitted_m.clear(); isEmitted_m.clear();
...@@ -171,10 +137,6 @@ void PartBins::sortArray() { ...@@ -171,10 +137,6 @@ void PartBins::sortArray() {
} }
void PartBins::sortArrayT() {
setActualemittedBin(0);
}
void PartBins::calcHBins() { void PartBins::calcHBins() {
for(unsigned int n = 0; n < tmppart_m.size(); n++) for(unsigned int n = 0; n < tmppart_m.size(); n++)
...@@ -182,27 +144,6 @@ void PartBins::calcHBins() { ...@@ -182,27 +144,6 @@ void PartBins::calcHBins() {
calcExtrema(); calcExtrema();
} }
size_t PartBins::getSum() {
size_t s = 0;
for(int n = 0; n < bins_m; n++)
s += nBin_m[n];
return s;
}
void PartBins::calcGlobalExtrema() {
xmin_m = std::numeric_limits<double>::max();
xmax_m = -xmin_m;
for(unsigned int n = 0; n < tmppart_m.size(); n++) {
if(tmppart_m[n][2] <= xmin_m)
xmin_m = tmppart_m[n][2];
if(tmppart_m[n][2] >= xmax_m)
xmax_m = tmppart_m[n][2];
}
double xdiff = 0.01 * (xmax_m - xmin_m);
xmin_m -= xdiff;
xmax_m += xdiff;
}
void PartBins::calcExtrema() { void PartBins::calcExtrema() {
for(unsigned int n = 0; n < tmppart_m.size(); n++) { for(unsigned int n = 0; n < tmppart_m.size(); n++) {
if(xbinmin_m[(int)tmppart_m[n][6]] >= tmppart_m[n][2]) if(xbinmin_m[(int)tmppart_m[n][6]] >= tmppart_m[n][2])
......
/** \file //
* \brief Defines a structure to hold energy bins and their // Class PartBins
* associated data // Defines a structure to hold energy bins and their associated data
* //
* // Copyright (c) 2007-2020, Paul Scherrer Institut, Villigen PSI, Switzerland
* // All rights reserved
* \author Andreas Adelmann //
* \date xx. November 2007 // This file is part of OPAL.
* //
* \warning None. // OPAL is free software: you can redistribute it and/or modify
* \attention // it under the terms of the GNU General Public License as published by
* \bug sure no bug in this code :=) // the Free Software Foundation, either version 3 of the License, or
* \todo // (at your option) any later version.
*/ //
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
#ifndef OPAL_Bins_HH #ifndef OPAL_Bins_HH
#define OPAL_Bins_HH #define OPAL_Bins_HH
#ifndef PartBinTest
#include "Algorithms/PBunchDefs.h"
#else
#include "ranlib.h"
#define Inform ostream
#endif
#include <gsl/gsl_rng.h>
#include <gsl/gsl_histogram.h> #include <gsl/gsl_histogram.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_randist.h>
#include <functional>
#include <memory>
#include <vector>
#include "Algorithms/PBunchDefs.h"
class Inform;
class PartBins { class PartBins {
...@@ -48,14 +48,6 @@ public: ...@@ -48,14 +48,6 @@ public:
virtual ~PartBins(); virtual ~PartBins();
/** \brief How many deleted particles are on one bin */
inline int getDelBinCont(int bin) {
int a = nDelBin_m[bin];
reduce(a, a, OpAddAssign());
return a;
}
/** \brief Add a particle to the temporary container */ /** \brief Add a particle to the temporary container */
void fill(std::vector<double> &p) { void fill(std::vector<double> &p) {
tmppart_m.push_back(p); tmppart_m.push_back(p);
...@@ -73,11 +65,6 @@ public: ...@@ -73,11 +65,6 @@ public:
/** set particles number in given bin */ /** set particles number in given bin */
void setPartNum(int bin, long long num) {nBin_m[bin] = num;} void setPartNum(int bin, long long num) {nBin_m[bin] = num;}
/** assume we emmit in monotinic increasing order */
void setBinEmitted(int bin) {binsEmitted_m[bin] = true;}
bool getBinEmitted(int bin) {return binsEmitted_m[bin];}
/** \brief Is true if we still have particles to emit */ /** \brief Is true if we still have particles to emit */
bool doEmission() {return getNp() > 0;} bool doEmission() {return getNp() > 0;}
...@@ -85,19 +72,6 @@ public: ...@@ -85,19 +72,6 @@ public:
return isEmitted_m[n]; //(isEmitted_m[n][0]==1) && (isEmitted_m[n][1] == bin); return isEmitted_m[n]; //(isEmitted_m[n][0]==1) && (isEmitted_m[n][1] == bin);
} }
void setEmitted(int n, int /*bin*/) {
isEmitted_m[n] = true;
}
void updatePartPos(int n, int /*bin*/, double z) {
tmppart_m[n][2] = z;
}
void updateExtramePos(int bin, double dz) {
xbinmax_m[bin] += dz;
xbinmin_m[bin] += dz;
}
/** assigns the proper position of particle n if it belongs to bin 'bin' */ /** assigns the proper position of particle n if it belongs to bin 'bin' */
bool getPart(size_t n, int bin, std::vector<double> &p); bool getPart(size_t n, int bin, std::vector<double> &p);
...@@ -106,29 +80,19 @@ public: ...@@ -106,29 +80,19 @@ public:
*/ */
void sortArray(); void sortArray();
private:
/** assigns the particles to the bins */ /** assigns the particles to the bins */
void calcHBins(); void calcHBins();
size_t getSum();
void calcGlobalExtrema();
void calcExtrema(); void calcExtrema();
void getExtrema(double &min, double &max) { /** assume we emit in monotonic increasing order */
min = xmin_m; void setBinEmitted(int bin) {binsEmitted_m[bin] = true;}
max = xmax_m;
}
/** update global bin parameters after inject a new bunch */ public:
void updateStatus(int bunchCount, size_t nPartInBin);
/** update particles number in bin after reset Bin ID of PartBunch */
void resetPartInBin(size_t newPartNum[]);
/** update local particles number in bin after reset Bin ID of PartBunch */ /** update local particles number in bin after reset Bin ID of PartBunch */
void resetPartInBin_cyc(size_t newPartNum[], int binID); void resetPartInBin_cyc(size_t newPartNum[], int binID);
/** update particles number in bin after particle deletion */
void updatePartInBin(size_t countLost[]);
/** update local particles number in bin after particle deletion */ /** update local particles number in bin after particle deletion */
void updatePartInBin_cyc(size_t countLost[]); void updatePartInBin_cyc(size_t countLost[]);
void updateDeletedPartsInBin(size_t countLost[]) ;
void setGamma(double gamma) { gamma_m = gamma;} void setGamma(double gamma) { gamma_m = gamma;}
double getGamma() {return gamma_m;} double getGamma() {return gamma_m;}
...@@ -159,13 +123,8 @@ protected: ...@@ -159,13 +123,8 @@ protected:
std::vector< std::vector<double> > tmppart_m; std::vector< std::vector<double> > tmppart_m;
std::vector< bool > isEmitted_m; std::vector< bool > isEmitted_m;
/** holds information whether all particles of a bin are emitted */ /** holds information whether all particles of a bin are emitted */
// std::vector< bool > binsEmitted_m;
std::unique_ptr<bool[]> binsEmitted_m; std::unique_ptr<bool[]> binsEmitted_m;
/**
Here comes the new stuff, t-binning
*/
public: public:
Inform &print(Inform &os); Inform &print(Inform &os);
...@@ -175,35 +134,9 @@ public: ...@@ -175,35 +134,9 @@ public:
/** get the number of used bin */ /** get the number of used bin */
virtual int getNBins() {return gsl_histogram_bins(h_m.get()) / sBins_m; } virtual int getNBins() {return gsl_histogram_bins(h_m.get()) / sBins_m; }
/** Get the total number of sampled bins */
virtual int getNSBins() {return gsl_histogram_bins(h_m.get()); }
int getBinToEmit() {
int save;
if(nemittedBins_m < getNBins()) {
save = nemittedBins_m;
nemittedBins_m++;
return save;
} else
return -1;
}
int getSBinToEmit() {
int save;
if(nemittedSBins_m < getNSBins()) {
save = nemittedSBins_m;
nemittedSBins_m++;
return save;
} else
return -1;
}
/** the last emitted bin is always smaller or equal getNbins */ /** the last emitted bin is always smaller or equal getNbins */
int getLastemittedBin() {return nemittedBins_m; } int getLastemittedBin() {return nemittedBins_m; }
/** set the actual emitted bib */
void setActualemittedBin(int bin) {nemittedBins_m = bin; }
/** \brief If the bunch object rebins we need to call resetBins() */ /** \brief If the bunch object rebins we need to call resetBins() */
void resetBins() { void resetBins() {
h_m.reset(nullptr); h_m.reset(nullptr);
...@@ -213,31 +146,6 @@ public: ...@@ -213,31 +146,6 @@ public:
return h_m != nullptr; return h_m != nullptr;
} }
/** sort the vector of particles according to the bin number */
void sortArrayT();
inline void setHistogram(gsl_histogram *h) { h_m.reset(h);}
/** \brief How many particles are on one bin */
virtual inline size_t getGlobalBinCount(int bin) {
size_t a = gsl_histogram_get(h_m.get(), bin);
reduce(a, a, OpAddAssign());
return a;
}
/** \brief How many particles are on one energy bin */
virtual inline size_t getLocalBinCount(int bin) {
size_t ret = 0;
for(int i = sBins_m * bin; i < sBins_m * (bin + 1); i++) {
ret += gsl_histogram_get(h_m.get(), i);
}
return ret;
}
/** \brief How many particles are in one sampling bin */
inline size_t getLocalSBinCount(int bin) { return gsl_histogram_get(h_m.get(), bin);}
/** \brief How many particles are in all the bins */ /** \brief How many particles are in all the bins */
size_t getTotalNum(); size_t getTotalNum();
...@@ -250,9 +158,6 @@ protected: ...@@ -250,9 +158,6 @@ protected:
/** number of emitted bins */ /** number of emitted bins */
int nemittedBins_m; int nemittedBins_m;
/** Number of total sampled bins emitted */
int nemittedSBins_m;
/** number of particles in the bins, the sum of all the nodes */ /** number of particles in the bins, the sum of all the nodes */
std::unique_ptr<size_t[]> nBin_m; std::unique_ptr<size_t[]> nBin_m;
......
...@@ -26,6 +26,7 @@ ...@@ -26,6 +26,7 @@
// //
#include "Algorithms/PartBinsCyc.h" #include "Algorithms/PartBinsCyc.h"
#include "Physics/Physics.h" #include "Physics/Physics.h"
#include "Utility/Inform.h"
#include <cfloat> #include <cfloat>
#include <vector> #include <vector>
extern Inform *gmsg; extern Inform *gmsg;
...@@ -36,7 +37,7 @@ PartBinsCyc::PartBinsCyc(int specifiedNumBins, int bins, size_t partInBin[]) ...@@ -36,7 +37,7 @@ PartBinsCyc::PartBinsCyc(int specifiedNumBins, int bins, size_t partInBin[])
bins_m = specifiedNumBins; // max bin number bins_m = specifiedNumBins; // max bin number
nemittedBins_m = bins; // the bin number with particles nemittedBins_m = bins; // the bin number with particles
for(int i = 0; i < nemittedBins_m; i++) { for(int i = 0; i < nemittedBins_m; i++) {
nBin_m[i] = partInBin[i]; nBin_m[i] = partInBin[i];
...@@ -51,8 +52,8 @@ PartBinsCyc::PartBinsCyc(int specifiedNumBins, int bins) ...@@ -51,8 +52,8 @@ PartBinsCyc::PartBinsCyc(int specifiedNumBins, int bins)
bins_m = specifiedNumBins; // max bin number bins_m = specifiedNumBins; // max bin number
nemittedBins_m = bins; // the bin number with particles nemittedBins_m = bins; // the bin number with particles
for(int i = 0; i < nemittedBins_m; i++) { for(int i = 0; i < nemittedBins_m; i++) {
binsEmitted_m[i] = true; binsEmitted_m[i] = true;
} }
} }
\ No newline at end of file
...@@ -4,7 +4,7 @@ ...@@ -4,7 +4,7 @@
// associated data for multi-bunch tracking in cyclotrons // associated data for multi-bunch tracking in cyclotrons
// //
// Copyright (c) 2010, Jianjun Yang, Paul Scherrer Institut, Villigen PSI, Switzerland // Copyright (c) 2010, Jianjun Yang, Paul Scherrer Institut, Villigen PSI, Switzerland
// Copyright (c) 2017-2019, Paul Scherrer Institut, Villigen PSI, Switzerland // Copyright (c) 2017-2020, Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved // All rights reserved
// //
// Implemented as part of the PhD thesis // Implemented as part of the PhD thesis
...@@ -27,47 +27,22 @@ ...@@ -27,47 +27,22 @@
#ifndef OPAL_BinsCyc_HH #ifndef OPAL_BinsCyc_HH
#define OPAL_BinsCyc_HH #define OPAL_BinsCyc_HH
#ifndef PartBinTest
#else
#include "ranlib.h"
#define Inform ostream
#endif
#include "Algorithms/PartBins.h" #include "Algorithms/PartBins.h"
#include <gsl/gsl_rng.h>
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_randist.h>
#include <cassert>
class PartBinsCyc: public PartBins { class PartBinsCyc: public PartBins {
public: public:
/** constructer function for cyclotron*/ /** constructor function for cyclotron*/
PartBinsCyc(int bunches, int bins, size_t partInBin[]); PartBinsCyc(int bunches, int bins, size_t partInBin[]);
PartBinsCyc(int specifiedNumBins, int bins); PartBinsCyc(int specifiedNumBins, int bins);
/** get the number of used bin */ /** get the number of used bin */
int getNBins() {return bins_m; } virtual int getNBins() {return bins_m; }
/** \brief How many particles are on one bin */
inline size_t getGlobalBinCount(int bin) {
size_t a = nBin_m[bin];
reduce(a, a, OpAddAssign());
return a;
}
/** \brief How many particles are on one energy bin */
inline size_t getLocalBinCount(int bin) {
assert(bin < bins_m);
return nBin_m[bin];
}
bool weHaveBins() { virtual bool weHaveBins() {
return ( nemittedBins_m > 0 ); return ( nemittedBins_m > 0 );
} }
}; };
......
...@@ -37,6 +37,8 @@ ...@@ -37,6 +37,8 @@
#include "Ippl.h" #include "Ippl.h"
#include <gsl/gsl_randist.h>
#include <cmath> #include <cmath>
#include <iostream> #include <iostream>
#include <fstream> #include <fstream>
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment