PartBins.cpp 7.21 KB
Newer Older
gsell's avatar
gsell committed
1
#include "Algorithms/PartBins.h"
2
#include "Algorithms/PBunchDefs.h"
gsell's avatar
gsell committed
3
#include "Physics/Physics.h"
kraus's avatar
kraus committed
4 5 6
#include <cfloat>
#include <vector>

gsell's avatar
gsell committed
7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134
extern Inform *gmsg;

PartBins::PartBins(int bins, int sbins) :
    gamma_m(1.0),
    bins_m(bins),
    sBins_m(sbins),
    xmin_m(0.0),
    xmax_m(0.0),
    hBin_m(0.0),
    nemittedBins_m(0),
    nemittedSBins_m(0) {

    nBin_m = new size_t[bins_m];
    xbinmin_m = new double[bins_m];
    xbinmax_m = new double[bins_m];
    binsEmitted_m = new bool[bins_m];
    nDelBin_m = new size_t[bins_m];

    for(int i = 0; i < bins_m; i++) {
        nDelBin_m[i] = nBin_m[i] = 0;
        xbinmax_m[i] = -DBL_MAX;
        xbinmin_m[i] = DBL_MAX;
        binsEmitted_m[i] = false;
    }
}


size_t PartBins::getTotalNum() {
    size_t s = 0;
    size_t sd = 0;
    size_t gs = 0;

    for(int i = 0; i < getLastemittedBin(); i++) {
        s  += nBin_m[i];
        sd += nDelBin_m[i];
    }
    gs = s - sd;
    reduce(gs, gs, OpAddAssign());
    return gs;
}

size_t PartBins::getTotalNumPerBin(int b) {
    size_t s = 0;
    s  = nBin_m[b];
    reduce(s, s, OpAddAssign());
    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 m2all("updateDeletedPartsInBin ", INFORM_ALL_NODES);

    for(int ii = 0; ii < getLastemittedBin(); ii++) {
        bool flagNeedUpdate = false;
        flagNeedUpdate = (countLost[ii] > 0);

        reduce(flagNeedUpdate, flagNeedUpdate, OpOr());

        if(flagNeedUpdate) {
            reduce(countLost[ii], countLost[ii], OpAddAssign());
            nDelBin_m[ii] = countLost[ii];
        }
        m2all << "In Bin: " << ii << ", " << nDelBin_m[ii] << " particle(s) lost" << endl;
    }
}


void PartBins::updatePartInBin(size_t countLost[]) {

    Inform msg0("updatePartInBin ");
    //  for (int ii=0; ii < bins_m; ii++){
    for(int ii = 0; ii < nemittedBins_m; ii++) {
        msg0 << "In Bin: " << ii << ", " << nBin_m[ii] << " particles " << endl;
    }
    for(int ii = 0; ii < nemittedBins_m; ii++) {
        bool flagNeedUpdate = false;
        if(countLost[ii] > 0)
            flagNeedUpdate = true;
        reduce(&flagNeedUpdate, &flagNeedUpdate + 1, &flagNeedUpdate, OpBitwiseOrAssign());
        //          reduce(flagNeedUpdate, flagNeedUpdate, OpOr());
        if(flagNeedUpdate) {
            reduce(countLost[ii], countLost[ii], OpAddAssign());
            nBin_m[ii] -= countLost[ii];
            msg0 << "In Bin: " << ii << ", " << countLost[ii] << " particle(s) lost" << endl;
        }
    }
}

void PartBins::updatePartInBin_cyc(size_t countLost[]) {

  for(int ii = 0; ii < nemittedBins_m; ii++) {
    if(countLost[ii] > 0)
      nBin_m[ii] -= countLost[ii];
  }
}


void PartBins::resetPartInBin(size_t newPartNum[]) {
    for(int ii = 0; ii < nemittedBins_m; ii++) {
        reduce(newPartNum[ii], newPartNum[ii], OpAddAssign());
        nBin_m[ii] = newPartNum[ii];
        INFOMSG("After reset Bin: " << ii << ", particle(s): " << newPartNum[ii] << endl);
    }
}


void PartBins::resetPartInBin_cyc(size_t newPartNum[], int maxbinIndex) {
    reduce(maxbinIndex, maxbinIndex, OpMaxAssign());
    nemittedBins_m =  maxbinIndex + 1;

    for(int ii = 0; ii < nemittedBins_m; ii++) {
        nBin_m[ii] = newPartNum[ii]; // only count particles on the local node
        setBinEmitted(ii);  // set true for this bin
    }
}



PartBins::~PartBins() {
135

gsell's avatar
gsell committed
136
    if(nBin_m) {
137 138 139 140
        delete [] nBin_m;
        delete [] xbinmax_m;
        delete [] xbinmin_m;
        delete [] binsEmitted_m;
gsell's avatar
gsell committed
141 142 143
    }
    tmppart_m.clear();
    isEmitted_m.clear();
144

gsell's avatar
gsell committed
145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258
}


bool PartBins::getPart(size_t n, int bin, std::vector<double> &p) {

    if(tmppart_m[n][6] == bin) {
        p = tmppart_m[n];
        return true;
    } else
        return false;
}

/** /brief There is only a local sort, no global yet */
void PartBins::sortArray() {
    /** sort the vector of particles such that position of the particles decrease with increasing index.
        Then push the particles back by 1e-13 s * beta * c (approximately one step).
        In order that the method getBin(double x) works xmin_m has to be lowered a bit more.
    */

    double sshift = sqrt(1. - (1. / (gamma_m * gamma_m))) * Physics::c * 1e-13;
    std::sort(tmppart_m.begin(), tmppart_m.end(), DescendingLocationSort(2));
    xmax_m = tmppart_m[0][2];
    xmin_m = tmppart_m.back()[2];

    for(unsigned int n = 0; n < tmppart_m.size(); n++)
        tmppart_m[n][2] -= xmax_m + sshift; /* push particles back */

    xmin_m -= xmax_m + 0.0001 * (xmax_m - xmin_m) + sshift; /* lower the limits */
    xmax_m = -sshift;

    reduce(xmin_m, xmin_m, OpMinAssign());
    reduce(xmax_m, xmax_m, OpMaxAssign());

    hBin_m = (fabs(xmax_m - xmin_m)) / (bins_m);
    calcHBins();
    for(int n = 0; n < bins_m; n++)
        if(nBin_m[n] == 0) setBinEmitted(n);
}


void PartBins::sortArrayT() {
    setActualemittedBin(0);
}

void PartBins::calcHBins() {

    for(unsigned int n = 0; n < tmppart_m.size(); n++)
        tmppart_m[n][6] = getBin(tmppart_m[n][2]);
    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 = DBL_MAX;
    xmax_m = -DBL_MAX;
    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() {
    for(unsigned int n = 0; n < tmppart_m.size(); n++) {
        if(xbinmin_m[(int)tmppart_m[n][6]] >= tmppart_m[n][2])
            xbinmin_m[(int)tmppart_m[n][6]] = tmppart_m[n][2];

        if(xbinmax_m[(int)tmppart_m[n][6]] <= tmppart_m[n][2])
            xbinmax_m[(int)tmppart_m[n][6]] = tmppart_m[n][2];
    }
}

Inform &PartBins::print(Inform &os) {

    os << "-----------------------------------------" << endl;
    os << "     CREATE BINNED GAUSS DISTRIBUTION DONE        " << endl;

    os << "Bins= " << bins_m << " hBin= " << hBin_m << " Particle vector length " << tmppart_m.size() << endl;

    //for(int i = 0; i < gsl_histogram_bins(h_m); i++)
    //os << "Bin # " << i << " val " << gsl_histogram_get(h_m, i) << endl;
    for(int i = 0; i < bins_m; i++) {
        size_t msum = 0;
        for(int j = 0; j < sBins_m; j++)
            msum += gsl_histogram_get(h_m, i * sBins_m + j);
        os << "Bin # " << i << " val " << msum << endl;
    }

    if(getLastemittedBin() >= 0)
        os << "Last emitted bin is " << getLastemittedBin() << endl;
    else
        os << "No bin is emitted !" << endl;
    return os;
}

int PartBins::getBin(double x) {
    /**
       returns the index of the bin to which the particle with z = 'x' belongs.
       If getBin returns b < 0 || b >= bins_m, then is x out of range!
    */
    int b = (int) floor(fabs(xmax_m - x) / hBin_m);
    nBin_m[b]++;
    return b;
}