PartBunchBase.hpp 61.8 KB
Newer Older
1 2 3
#ifndef PART_BUNCH_BASE_HPP
#define PART_BUNCH_BASE_HPP

4
#include <cmath>
frey_m's avatar
frey_m committed
5 6
#include <numeric>

frey_m's avatar
frey_m committed
7 8
#include "Distribution/Distribution.h"

9
#include "AbstractObjects/OpalData.h"   // OPAL file
10
#include "Physics/Physics.h"
11
#include "Utilities/GeneralClassicException.h"
12
#include "Utilities/OpalException.h"
13
#include "Utilities/Options.h"
14
#include "Utilities/SwitcherError.h"
15 16 17
#include "Utilities/Util.h"

extern Inform *gmsg;
18 19 20

template <class T, unsigned Dim>
PartBunchBase<T, Dim>::PartBunchBase(AbstractParticle<T, Dim>* pb)
21
    : R(*(pb->R_p)),
22 23
      ID(*(pb->ID_p)),
      pbin_m(nullptr),
24
      lowParticleCount_m(false),
25 26
      pmsg_m(nullptr),
      f_stream(nullptr),
27
      fixed_grid(false),
28 29 30 31 32 33 34 35
//       reference(ref), //FIXME
      unit_state_(units),
      stateOfLastBoundP_(unitless),
      moments_m(),
      dt_m(0.0),
      t_m(0.0),
      eKin_m(0.0),
      dE_m(0.0),
36
      spos_m(0.0),
37 38
      globalMeanR_m(Vector_t(0.0, 0.0, 0.0)),
      globalToLocalQuaternion_m(Quaternion_t(1.0, 0.0, 0.0, 0.0)),
39 40 41 42 43 44 45 46
      rmax_m(0.0),
      rmin_m(0.0),
      rrms_m(0.0),
      prms_m(0.0),
      rmean_m(0.0),
      pmean_m(0.0),
      eps_m(0.0),
      eps_norm_m(0.0),
47
      halo_m(Vector_t(0.0, 0.0, 0.0)),
48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67
      rprms_m(0.0),
      Dx_m(0.0),
      Dy_m(0.0),
      DDx_m(0.0),
      DDy_m(0.0),
      hr_m(-1.0),
      nr_m(0),
      fs_m(nullptr),
      couplingConstant_m(0.0),
      qi_m(0.0),
      distDump_m(0),
      fieldDBGStep_m(0),
      dh_m(1e-12),
      tEmission_m(0.0),
      bingamma_m(nullptr),
      binemitted_m(nullptr),
      stepsPerTurn_m(0),
      localTrackStep_m(0),
      globalTrackStep_m(0),
      numBunch_m(1),
frey_m's avatar
frey_m committed
68 69
      bunchTotalNum_m(1),
      bunchLocalNum_m(1),
70 71 72 73
      SteptoLastInj_m(0),
      globalPartPerNode_m(nullptr),
      dist_m(nullptr),
      dcBeam_m(false),
74
      periodLength_m(Physics::c / 1e9),
75
      pbase(pb)
76
{
frey_m's avatar
frey_m committed
77
    setup(pb);
78 79 80
}

template <class T, unsigned Dim>
frey_m's avatar
frey_m committed
81
PartBunchBase<T, Dim>::PartBunchBase(AbstractParticle<T, Dim>* pb, const PartData *ref)
82
    : PartBunchBase(pb)
83
{
84
    reference = ref;
85
}
86 87 88 89 90 91 92

/*
 * Bunch common member functions
 */

template <class T, unsigned Dim>
bool PartBunchBase<T, Dim>::getIfBeamEmitting() {
93 94 95 96 97 98 99 100 101 102 103 104
    if (dist_m != NULL) {
        size_t isBeamEmitted = dist_m->getIfDistEmitting();
        reduce(isBeamEmitted, isBeamEmitted, OpAddAssign());
        if (isBeamEmitted > 0)
            return true;
        else
            return false;
    } else
        return false;
}


105 106
template <class T, unsigned Dim>
int PartBunchBase<T, Dim>::getLastEmittedEnergyBin() {
107 108 109 110 111 112 113 114 115
    /*
     * Get maximum last emitted energy bin.
     */
    int lastEmittedBin = dist_m->getLastEmittedEnergyBin();
    reduce(lastEmittedBin, lastEmittedBin, OpMaxAssign());
    return lastEmittedBin;
}


116 117
template <class T, unsigned Dim>
size_t PartBunchBase<T, Dim>::getNumberOfEmissionSteps() {
118 119 120 121
    return dist_m->getNumberOfEmissionSteps();
}


122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 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
template <class T, unsigned Dim>
int PartBunchBase<T, Dim>::getNumberOfEnergyBins() {
    return dist_m->getNumberOfEnergyBins();
}


template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::Rebin() {

    size_t isBeamEmitting = dist_m->getIfDistEmitting();
    reduce(isBeamEmitting, isBeamEmitting, OpAddAssign());
    if (isBeamEmitting > 0) {
        *gmsg << "*****************************************************" << endl
              << "Warning: attempted to rebin, but not all distribution" << endl
              << "particles have been emitted. Rebin failed." << endl
              << "*****************************************************" << endl;
    } else {
        if (dist_m->Rebin())
            this->Bin = 0;
    }
}


template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::setEnergyBins(int numberOfEnergyBins) {
    bingamma_m = std::unique_ptr<double[]>(new double[numberOfEnergyBins]);
    binemitted_m = std::unique_ptr<size_t[]>(new size_t[numberOfEnergyBins]);
    for(int i = 0; i < numberOfEnergyBins; i++)
        binemitted_m[i] = 0;
}


template <class T, unsigned Dim>
bool PartBunchBase<T, Dim>::weHaveEnergyBins() {

    if (dist_m != NULL)
        return dist_m->getNumberOfEnergyBins() > 0;
    else
        return false;
}


template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::switchToUnitlessPositions(bool use_dt_per_particle) {

    if(unit_state_ == unitless)
        throw SwitcherError("PartBunch::switchToUnitlessPositions",
                            "Cannot make a unitless PartBunch unitless");

    bool hasToReset = false;
    if(!R.isDirty()) hasToReset = true;

    for(size_t i = 0; i < getLocalNum(); i++) {
        double dt = getdT();
        if(use_dt_per_particle)
            dt = this->dt[i];

        R[i] /= Vector_t(Physics::c * dt);
    }

    unit_state_ = unitless;

    if(hasToReset) R.resetDirtyFlag();
}


//FIXME: unify methods, use convention that all particles have own dt
template <class T, unsigned Dim>
190
void PartBunchBase<T, Dim>::switchOffUnitlessPositions(bool use_dt_per_particle) {
191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214

    if(unit_state_ == units)
        throw SwitcherError("PartBunch::switchOffUnitlessPositions",
                            "Cannot apply units twice to PartBunch");

    bool hasToReset = false;
    if(!R.isDirty()) hasToReset = true;

    for(size_t i = 0; i < getLocalNum(); i++) {
        double dt = getdT();
        if(use_dt_per_particle)
            dt = this->dt[i];

        R[i] *= Vector_t(Physics::c * dt);
    }

    unit_state_ = units;

    if(hasToReset) R.resetDirtyFlag();
}


template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::do_binaryRepart() {
frey_m's avatar
frey_m committed
215
    // do nothing here
216 217 218 219 220 221
}


template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::setDistribution(Distribution *d,
                                            std::vector<Distribution *> addedDistributions,
222
                                            size_t &np)
223 224 225 226
{
    Inform m("setDistribution " );
    dist_m = d;

227
    dist_m->createOpalT(this, addedDistributions, np);
228 229

//    if (Options::cZero)
230
//        dist_m->create(this, addedDistributions, np / 2);
231
//    else
232
//        dist_m->create(this, addedDistributions, np);
233 234 235 236
}


template <class T, unsigned Dim>
237
bool PartBunchBase<T, Dim>::isGridFixed() const {
238 239 240 241
    return fixed_grid;
}


242
template <class T, unsigned Dim>
243
bool PartBunchBase<T, Dim>::hasBinning() const {
244 245 246 247
    return (pbin_m != nullptr);
}


248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295
template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::setTEmission(double t) {
    tEmission_m = t;
}


template <class T, unsigned Dim>
double PartBunchBase<T, Dim>::getTEmission() {
    return tEmission_m;
}


template <class T, unsigned Dim>
bool PartBunchBase<T, Dim>::doEmission() {
    if (dist_m != NULL)
        return dist_m->getIfDistEmitting();
    else
        return false;
}


template <class T, unsigned Dim>
bool PartBunchBase<T, Dim>::weHaveBins() const {

    if(pbin_m != NULL)
        return pbin_m->weHaveBins();
    else
        return false;
}


template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::setPBins(PartBins *pbin) {
    pbin_m = pbin;
    *gmsg << *pbin_m << endl;
    setEnergyBins(pbin_m->getNBins());
}


template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::setPBins(PartBinsCyc *pbin) {
    pbin_m = pbin;
    setEnergyBins(pbin_m->getNBins());
}


template <class T, unsigned Dim>
size_t PartBunchBase<T, Dim>::emitParticles(double eZ) {
frey_m's avatar
frey_m committed
296
    return dist_m->emitParticles(this, eZ);
297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325
}


template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::updateNumTotal() {
    size_t numParticles = getLocalNum();
    reduce(numParticles, numParticles, OpAddAssign());
    setTotalNum(numParticles);
}


template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::rebin() {
    this->Bin = 0;
    pbin_m->resetBins();
    // delete pbin_m; we did not allocate it!
    pbin_m = NULL;
}


template <class T, unsigned Dim>
int PartBunchBase<T, Dim>::getLastemittedBin() {
    if(pbin_m != NULL)
        return pbin_m->getLastemittedBin();
    else
        return 0;
}


326 327 328 329 330 331 332 333
template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::setLocalBinCount(size_t num, int bin) {
    if(pbin_m != NULL) {
        pbin_m->setPartNum(bin, num);
    }
}


334 335 336 337 338 339 340 341 342 343 344
template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::calcGammas() {

    const int emittedBins = dist_m->getNumberOfEnergyBins();

    size_t s = 0;

    for(int i = 0; i < emittedBins; i++)
        bingamma_m[i] = 0.0;

    for(unsigned int n = 0; n < getLocalNum(); n++)
345
        bingamma_m[this->Bin[n]] += std::sqrt(1.0 + dot(this->P[n], this->P[n]));
346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387

    std::unique_ptr<size_t[]> particlesInBin(new size_t[emittedBins]);
    reduce(bingamma_m.get(), bingamma_m.get() + emittedBins, bingamma_m.get(), OpAddAssign());
    reduce(binemitted_m.get(), binemitted_m.get() + emittedBins, particlesInBin.get(), OpAddAssign());
    for(int i = 0; i < emittedBins; i++) {
        size_t &pInBin = particlesInBin[i];
        if(pInBin != 0) {
            bingamma_m[i] /= pInBin;
            INFOMSG(level2 << "Bin " << std::setw(3) << i
                           << " gamma = " << std::setw(8) << std::scientific
                           << std::setprecision(5) << bingamma_m[i]
                           << "; NpInBin= " << std::setw(8)
                           << std::setfill(' ') << pInBin << endl);
        } else {
            bingamma_m[i] = 1.0;
            INFOMSG(level2 << "Bin " << std::setw(3) << i << " has no particles " << endl);
        }
        s += pInBin;
    }
    particlesInBin.reset();


    if(s != getTotalNum() && !OpalData::getInstance()->hasGlobalGeometry())
        ERRORMSG("sum(Bins)= " << s << " != sum(R)= " << getTotalNum() << endl;);

    if(emittedBins >= 2) {
        for(int i = 1; i < emittedBins; i++) {
            if(binemitted_m[i - 1] != 0 && binemitted_m[i] != 0)
                INFOMSG(level2 << "d(gamma)= " << 100.0 * std::abs(bingamma_m[i - 1] - bingamma_m[i]) / bingamma_m[i] << " [%] "
                        << "between bin " << i - 1 << " and " << i << endl);
        }
    }
}


template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::calcGammas_cycl() {

    const int emittedBins = pbin_m->getLastemittedBin();

    for(int i = 0; i < emittedBins; i++)
        bingamma_m[i] = 0.0;
388 389
    for(unsigned int n = 0; n < getLocalNum(); n++) {
        if ( this->Bin[n] > -1 ) {
390
            bingamma_m[this->Bin[n]] += std::sqrt(1.0 + dot(this->P[n], this->P[n]));
kraus's avatar
kraus committed
391
        }
392
    }
frey_m's avatar
frey_m committed
393 394 395

    allreduce(*bingamma_m.get(), emittedBins, std::plus<double>());

396 397 398 399 400
    for(int i = 0; i < emittedBins; i++) {
        if(pbin_m->getTotalNumPerBin(i) > 0)
            bingamma_m[i] /= pbin_m->getTotalNumPerBin(i);
        else
            bingamma_m[i] = 0.0;
frey_m's avatar
frey_m committed
401 402
        INFOMSG("Bin " << i << " : particle number = " << pbin_m->getTotalNumPerBin(i)
                       << " gamma = " << bingamma_m[i] << endl);
403 404 405 406 407 408 409 410 411 412 413 414 415
    }

}


template <class T, unsigned Dim>
double PartBunchBase<T, Dim>::getBinGamma(int bin) {
    return bingamma_m[bin];
}


template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::setBinCharge(int bin, double q) {
416
  this->Q = where(eq(this->Bin, bin), q, 0.0);
417 418 419 420 421
}


template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::setBinCharge(int bin) {
422
  this->Q = where(eq(this->Bin, bin), this->qi_m, 0.0);
423 424 425 426 427 428
}


template <class T, unsigned Dim>
size_t PartBunchBase<T, Dim>::calcNumPartsOutside(Vector_t x) {

frey_m's avatar
frey_m committed
429
    std::size_t localnum = 0;
430 431 432
    const Vector_t meanR = get_rmean();

    for(unsigned long k = 0; k < getLocalNum(); ++ k)
433 434 435
        if (std::abs(R[k](0) - meanR(0)) > x(0) ||
            std::abs(R[k](1) - meanR(1)) > x(1) ||
            std::abs(R[k](2) - meanR(2)) > x(2)) {
436

frey_m's avatar
frey_m committed
437
            ++localnum;
438
        }
439

frey_m's avatar
frey_m committed
440
    gather(&localnum, &globalPartPerNode_m[0], 1);
441

frey_m's avatar
frey_m committed
442 443 444
    size_t npOutside = std::accumulate(globalPartPerNode_m.get(),
                                       globalPartPerNode_m.get() + Ippl::getNodes(), 0,
                                       std::plus<size_t>());
445

frey_m's avatar
frey_m committed
446
    return npOutside;
447 448 449
}


450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467
/**
 * \method calcLineDensity()
 * \brief calculates the 1d line density (not normalized) and append it to a file.
 * \see ParallelTTracker
 * \warning none yet
 *
 * DETAILED TODO
 *
 */
template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::calcLineDensity(unsigned int nBins,
                                            std::vector<double> &lineDensity,
                                            std::pair<double, double> &meshInfo)
{
    Vector_t rmin, rmax;
    get_bounds(rmin, rmax);

    if (nBins < 2) {
frey_m's avatar
frey_m committed
468
        Vektor<int, 3>/*NDIndex<3>*/ grid;
469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485
        this->updateDomainLength(grid);
        nBins = grid[2];
    }

    double length = rmax(2) - rmin(2);
    double zmin = rmin(2) - dh_m * length, zmax = rmax(2) + dh_m * length;
    double hz = (zmax - zmin) / (nBins - 2);
    double perMeter = 1.0 / hz;//(zmax - zmin);
    zmin -= hz;

    lineDensity.resize(nBins, 0.0);
    std::fill(lineDensity.begin(), lineDensity.end(), 0.0);

    const unsigned int lN = getLocalNum();
    for (unsigned int i = 0; i < lN; ++ i) {
        const double z = R[i](2) - 0.5 * hz;
        unsigned int idx = (z - zmin) / hz;
486
        double tau = (z - zmin) / hz - idx;
487 488 489 490 491 492 493 494 495 496 497

        lineDensity[idx] += Q[i] * (1.0 - tau) * perMeter;
        lineDensity[idx + 1] += Q[i] * tau * perMeter;
    }

    reduce(&(lineDensity[0]), &(lineDensity[0]) + nBins, &(lineDensity[0]), OpAddAssign());

    meshInfo.first = zmin;
    meshInfo.second = hz;
}

frey_m's avatar
frey_m committed
498 499 500 501 502 503 504 505 506 507 508 509 510 511 512


template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::boundp() {
    /*
      Assume rmin_m < 0.0
     */

    IpplTimings::startTimer(boundpTimer_m);
    //if(!R.isDirty() && stateOfLastBoundP_ == unit_state_) return;
    if ( !(R.isDirty() || ID.isDirty() ) && stateOfLastBoundP_ == unit_state_) return; //-DW

    stateOfLastBoundP_ = unit_state_;

    if(!isGridFixed()) {
513
        const int dimIdx = (dcBeam_m? 2: 3);
frey_m's avatar
frey_m committed
514 515

        /**
516 517 518 519
            In case of dcBeam_m && hr_m < 0
            this is the first call to boundp and we
            have to set hr completely i.e. x,y and z.
         */
frey_m's avatar
frey_m committed
520

521
        this->updateDomainLength(nr_m);
522
        IpplTimings::startTimer(boundpBoundsTimer_m);
frey_m's avatar
frey_m committed
523
        get_bounds(rmin_m, rmax_m);
524
        IpplTimings::stopTimer(boundpBoundsTimer_m);
frey_m's avatar
frey_m committed
525
        Vector_t len = rmax_m - rmin_m;
526

527 528 529 530 531 532 533 534 535
        double volume = 1.0;
        for(int i = 0; i < dimIdx; i++) {
            double length = std::abs(rmax_m[i] - rmin_m[i]);
            if (length < 1e-10) {
                rmax_m[i] += 1e-10;
                rmin_m[i] -= 1e-10;
            } else {
                rmax_m[i] += dh_m * length;
                rmin_m[i] -= dh_m * length;
frey_m's avatar
frey_m committed
536
            }
537 538 539 540 541 542 543 544 545
            hr_m[i]    = (rmax_m[i] - rmin_m[i]) / (nr_m[i] - 1);
        }
        if (dcBeam_m) {
            rmax_m[2] = rmin_m[2] + periodLength_m;
            hr_m[2] = periodLength_m / (nr_m[2] - 1);
        }
        for (int i = 0; i < dimIdx; ++ i) {
            volume *= std::abs(rmax_m[i] - rmin_m[i]);
        }
frey_m's avatar
frey_m committed
546

547 548 549 550 551 552 553
        if (getIfBeamEmitting() && dist_m != NULL) {
            // keep particles per cell ratio high, don't spread a hand full particles across the whole grid
            double percent = std::max(1.0 / (nr_m[2] - 1), dist_m->getPercentageEmitted());
            double length  = std::abs(rmax_m[2] - rmin_m[2]) / (1.0 + 2 * dh_m);
            if (percent < 1.0 && percent > 0.0) {
                rmax_m[2] -= dh_m * length;
                rmin_m[2] = rmax_m[2] - length / percent;
554

555
                length /= percent;
frey_m's avatar
frey_m committed
556

557 558
                rmax_m[2] += dh_m * length;
                rmin_m[2] -= dh_m * length;
559

560
                hr_m[2] = (rmax_m[2] - rmin_m[2]) / (nr_m[2] - 1);
frey_m's avatar
frey_m committed
561
            }
562
        }
frey_m's avatar
frey_m committed
563

564
        if (volume < 1e-21 && getTotalNum() > 1 && std::abs(sum(Q)) > 0.0) {
565
            WARNMSG(level1 << "!!! Extremely high particle density detected !!!" << endl);
frey_m's avatar
frey_m committed
566
        }
567
        //INFOMSG("It is a full boundp hz= " << hr_m << " rmax= " << rmax_m << " rmin= " << rmin_m << endl);
568

569
        if(hr_m[0] * hr_m[1] * hr_m[2] <= 0) {
frey_m's avatar
frey_m committed
570 571
            throw GeneralClassicException("boundp() ", "h<0, can not build a mesh");
        }
572 573 574

        Vector_t origin = rmin_m - Vector_t(hr_m[0] / 2.0, hr_m[1] / 2.0, hr_m[2] / 2.0);
        this->updateFields(hr_m, origin);
frey_m's avatar
frey_m committed
575
    }
576
    IpplTimings::startTimer(boundpUpdateTimer_m);
frey_m's avatar
frey_m committed
577
    update();
578
    IpplTimings::stopTimer(boundpUpdateTimer_m);
frey_m's avatar
frey_m committed
579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599
    R.resetDirtyFlag();

    IpplTimings::stopTimer(boundpTimer_m);
}


template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::boundp_destroy() {

    Inform gmsgAll("boundp_destroy ", INFORM_ALL_NODES);

    Vector_t len;
    const int dimIdx = 3;
    IpplTimings::startTimer(boundpTimer_m);

    std::unique_ptr<size_t[]> countLost;
    if(weHaveBins()) {
        const int tempN = pbin_m->getLastemittedBin();
        countLost = std::unique_ptr<size_t[]>(new size_t[tempN]);
        for(int ii = 0; ii < tempN; ii++) countLost[ii] = 0;
    }
600

601
    this->updateDomainLength(nr_m);
frey_m's avatar
frey_m committed
602

603

604
    IpplTimings::startTimer(boundpBoundsTimer_m);
frey_m's avatar
frey_m committed
605
    get_bounds(rmin_m, rmax_m);
frey_m's avatar
frey_m committed
606
    IpplTimings::stopTimer(boundpBoundsTimer_m);
607

frey_m's avatar
frey_m committed
608 609 610 611 612 613 614 615 616 617 618 619 620 621 622
    len = rmax_m - rmin_m;

    calcBeamParameters();

    int checkfactor = Options::remotePartDel;
    if (checkfactor != 0) {
        //INFOMSG("checkfactor= " << checkfactor << endl);
        // check the bunch if its full size is larger than checkfactor times of its rms size
        if(checkfactor < 0) {
            checkfactor *= -1;
            if (len[0] > checkfactor * rrms_m[0] ||
                len[1] > checkfactor * rrms_m[1] ||
                len[2] > checkfactor * rrms_m[2])
            {
                for(unsigned int ii = 0; ii < this->getLocalNum(); ii++) {
623
                    /* delete the particle if the distance to the beam center
frey_m's avatar
frey_m committed
624 625
                     * is larger than 8 times of beam's rms size
                     */
626 627 628
                    if (std::abs(R[ii](0) - rmean_m(0)) > checkfactor * rrms_m[0] ||
                        std::abs(R[ii](1) - rmean_m(1)) > checkfactor * rrms_m[1] ||
                        std::abs(R[ii](2) - rmean_m(2)) > checkfactor * rrms_m[2])
frey_m's avatar
frey_m committed
629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646
                    {
                        // put particle onto deletion list
                        destroy(1, ii);
                        //update bin parameter
                        if (weHaveBins())
                            countLost[Bin[ii]] += 1 ;
                        /* INFOMSG("REMOTE PARTICLE DELETION: ID = " << ID[ii] << ", R = " << R[ii]
                         * << ", beam rms = " << rrms_m << endl;);
                         */
                    }
                }
            }
        }
        else {
            if (len[0] > checkfactor * rrms_m[0] ||
                len[2] > checkfactor * rrms_m[2])
            {
                for(unsigned int ii = 0; ii < this->getLocalNum(); ii++) {
647
                    /* delete the particle if the distance to the beam center
frey_m's avatar
frey_m committed
648 649
                     * is larger than 8 times of beam's rms size
                     */
650 651
                    if (std::abs(R[ii](0) - rmean_m(0)) > checkfactor * rrms_m[0] ||
                        std::abs(R[ii](2) - rmean_m(2)) > checkfactor * rrms_m[2])
frey_m's avatar
frey_m committed
652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674
                    {
                        // put particle onto deletion list
                        destroy(1, ii);
                        //update bin parameter
                        if (weHaveBins())
                            countLost[Bin[ii]] += 1 ;
                        /* INFOMSG("REMOTE PARTICLE DELETION: ID = " << ID[ii] << ", R = " << R[ii]
                         * << ", beam rms = " << rrms_m << endl;);
                         */
                    }
                }
            }
        }
    }

    for(int i = 0; i < dimIdx; i++) {
        double length = std::abs(rmax_m[i] - rmin_m[i]);
        rmax_m[i] += dh_m * length;
        rmin_m[i] -= dh_m * length;
        hr_m[i]    = (rmax_m[i] - rmin_m[i]) / (nr_m[i] - 1);
    }

    // rescale mesh
675
    this->updateFields(hr_m, rmin_m);
676

frey_m's avatar
frey_m committed
677 678 679
    if(weHaveBins()) {
        pbin_m->updatePartInBin_cyc(countLost.get());
    }
680

frey_m's avatar
frey_m committed
681 682 683 684 685
    /* we also need to update the number of particles per bunch
     * expensive since does an allreduce!
     */
    countTotalNumPerBunch();

686
    IpplTimings::startTimer(boundpUpdateTimer_m);
frey_m's avatar
frey_m committed
687
    update();
688
    IpplTimings::stopTimer(boundpUpdateTimer_m);
689

frey_m's avatar
frey_m committed
690 691 692 693 694 695 696 697 698
    IpplTimings::stopTimer(boundpTimer_m);
}


template <class T, unsigned Dim>
size_t PartBunchBase<T, Dim>::boundp_destroyT() {

    const unsigned int minNumParticlesPerCore = getMinimumNumberOfParticlesPerCore();

699
    this->updateDomainLength(nr_m);
frey_m's avatar
frey_m committed
700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759

    std::unique_ptr<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);
            }
        }
        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();
    }
    calcBeamParameters();
    gatherLoadBalanceStatistics();

    if(weHaveEnergyBins()) {
        const int lastBin = dist_m->getLastEmittedEnergyBin();
        for(int i = 0; i < lastBin; i++) {
            binemitted_m[i] = tmpbinemitted[i];
        }
    }
    reduce(ne, ne, OpAddAssign());
    return ne;
}


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();
    const size_t totalNum = getTotalNum();
760
    size_t ne = 0;
frey_m's avatar
frey_m committed
761 762 763 764 765 766 767 768 769

    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) {
                destroy(1, i);
770
                ++ ne;
frey_m's avatar
frey_m committed
771 772 773 774 775 776
            } else
                tmpbinemitted[Bin[i]]++;
        }
    } else {
        Inform dmsg("destroy: ", INFORM_ALL_NODES);
        for(size_t i = 0; i < localNum; i++) {
777 778 779 780 781
            if((Bin[i] < 0)) {
                if ((localNum - ne) > minNumParticlesPerCore) {   // need in minimum x particles per node
                    ne++;
                    destroy(1, i);
                }
frey_m's avatar
frey_m committed
782 783 784 785
            }
        }
        lowParticleCount_m = ((localNum - ne) <= minNumParticlesPerCore);
        reduce(lowParticleCount_m, lowParticleCount_m, OpOr());
786 787 788 789
    }

    if (ne > 0) {
        performDestroy(true);
frey_m's avatar
frey_m committed
790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808
    }

    calcBeamParameters();
    gatherLoadBalanceStatistics();

    if (weHaveEnergyBins()) {
        const int lastBin = dist_m->getLastEmittedEnergyBin();
        for(int i = 0; i < lastBin; i++) {
            binemitted_m[i] = tmpbinemitted[i];
        }
    }
    size_t newTotalNum = getLocalNum();
    reduce(newTotalNum, newTotalNum, OpAddAssign());

    setTotalNum(newTotalNum);

    return totalNum - newTotalNum;
}

809
template <class T, unsigned Dim>
810
double PartBunchBase<T, Dim>::getPx(int /*i*/) {
811 812 813 814 815
    return 0.0;
}


template <class T, unsigned Dim>
816
double PartBunchBase<T, Dim>::getPy(int) {
817 818 819 820 821
    return 0.0;
}


template <class T, unsigned Dim>
822
double PartBunchBase<T, Dim>::getPz(int) {
823 824 825 826 827
    return 0.0;
}


template <class T, unsigned Dim>
828
double PartBunchBase<T, Dim>::getPx0(int) {
829 830 831 832 833
    return 0.0;
}


template <class T, unsigned Dim>
834
double PartBunchBase<T, Dim>::getPy0(int) {
835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861
    return 0;
}


//ff
template <class T, unsigned Dim>
double PartBunchBase<T, Dim>::getX(int i) {
    return this->R[i](0);
}


//ff
template <class T, unsigned Dim>
double PartBunchBase<T, Dim>::getY(int i) {
    return this->R[i](1);
}


//ff
template <class T, unsigned Dim>
double PartBunchBase<T, Dim>::getZ(int i) {
    return this->R[i](2);
}


//ff
template <class T, unsigned Dim>
862
double PartBunchBase<T, Dim>::getX0(int /*i*/) {
863 864 865 866 867 868
    return 0.0;
}


//ff
template <class T, unsigned Dim>
869
double PartBunchBase<T, Dim>::getY0(int /*i*/) {
870 871 872 873 874
    return 0.0;
}


template <class T, unsigned Dim>
875
void PartBunchBase<T, Dim>::setZ(int /*i*/, double /*zcoo*/) {
876 877 878 879 880
};


template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::get_bounds(Vector_t &rmin, Vector_t &rmax) {
881

frey_m's avatar
frey_m committed
882
    this->getLocalBounds(rmin, rmax);
883

kraus's avatar
kraus committed
884
    double min[2*Dim];
885

886
    for (unsigned int i = 0; i < Dim; ++i) {
kraus's avatar
kraus committed
887 888
        min[2*i] = rmin[i];
        min[2*i + 1] = -rmax[i];
889
    }
890

kraus's avatar
kraus committed
891
    allreduce(min, 2*Dim, std::less<double>());
892

893
    for (unsigned int i = 0; i < Dim; ++i) {
kraus's avatar
kraus committed
894 895
        rmin[i] = min[2*i];
        rmax[i] = -min[2*i + 1];
896
    }
897 898 899 900 901 902
}


template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::getLocalBounds(Vector_t &rmin, Vector_t &rmax) {
    const size_t localNum = getLocalNum();
frey_m's avatar
frey_m committed
903
    if (localNum == 0) {
904 905 906
        double maxValue = 1e8;
        rmin = Vector_t(maxValue, maxValue, maxValue);
        rmax = Vector_t(-maxValue, -maxValue, -maxValue);
kraus's avatar
kraus committed
907
        return;
frey_m's avatar
frey_m committed
908
    }
909

910 911 912 913 914
    rmin = R[0];
    rmax = R[0];
    for (size_t i = 1; i < localNum; ++ i) {
        for (unsigned short d = 0; d < 3u; ++ d) {
            if (rmin(d) > R[i](d)) rmin(d) = R[i](d);
kraus's avatar
kraus committed
915
            else if (rmax(d) < R[i](d)) rmax(d) = R[i](d);
916 917 918 919 920 921 922 923 924 925 926 927
        }
    }
}


template <class T, unsigned Dim>
std::pair<Vector_t, double> PartBunchBase<T, Dim>::getBoundingSphere() {
    Vector_t rmin, rmax;
    get_bounds(rmin, rmax);

    std::pair<Vector_t, double> sphere;
    sphere.first = 0.5 * (rmin + rmax);
928
    sphere.second = std::sqrt(dot(rmax - sphere.first, rmax - sphere.first));
929 930 931 932 933 934 935 936 937 938 939 940

    return sphere;
}


template <class T, unsigned Dim>
std::pair<Vector_t, double> PartBunchBase<T, Dim>::getLocalBoundingSphere() {
    Vector_t rmin, rmax;
    getLocalBounds(rmin, rmax);

    std::pair<Vector_t, double> sphere;
    sphere.first = 0.5 * (rmin + rmax);
941
    sphere.second = std::sqrt(dot(rmax - sphere.first, rmax - sphere.first));
942 943 944 945 946 947

    return sphere;
}


template <class T, unsigned Dim>
948
void PartBunchBase<T, Dim>::push_back(OpalParticle p) {
949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978
    Inform msg("PartBunch ");

    create(1);
    size_t i = getTotalNum();

    R[i](0) = p[0];
    R[i](1) = p[1];
    R[i](2) = p[2];

    P[i](0) = p[3];
    P[i](1) = p[4];
    P[i](2) = p[5];

    update();
    msg << "Created one particle i= " << i << endl;
}


template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::set_part(FVector<double, 6> z, int ii) {
    R[ii](0) = z[0];
    P[ii](0) = z[1];
    R[ii](1) = z[2];
    P[ii](1) = z[3];
    R[ii](2) = z[4];
    P[ii](2) = z[5];
}


template <class T, unsigned Dim>
979
void PartBunchBase<T, Dim>::set_part(OpalParticle p, int ii) {
980 981 982 983 984 985 986 987 988 989
    R[ii](0) = p[0];
    P[ii](0) = p[1];
    R[ii](1) = p[2];
    P[ii](1) = p[3];
    R[ii](2) = p[4];
    P[ii](2) = p[5];
}


template <class T, unsigned Dim>
990 991
OpalParticle PartBunchBase<T, Dim>::get_part(int ii) {
    OpalParticle part;
992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005
    part[0] = R[ii](0);
    part[1] = P[ii](0);
    part[2] = R[ii](1);
    part[3] = P[ii](1);
    part[4] = R[ii](2);
    part[5] = P[ii](2);
    return part;
}


template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::maximumAmplitudes(const FMatrix<double, 6, 6> &D,
                                              double &axmax, double &aymax) {
    axmax = aymax = 0.0;
1006
    OpalParticle part;
1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048

    for(unsigned int ii = 0; ii < getLocalNum(); ii++) {

        part = get_part(ii);

        double xnor =
            D(0, 0) * part.x()  + D(0, 1) * part.px() + D(0, 2) * part.y() +
            D(0, 3) * part.py() + D(0, 4) * part.t()  + D(0, 5) * part.pt();
        double pxnor =
            D(1, 0) * part.x()  + D(1, 1) * part.px() + D(1, 2) * part.y() +
            D(1, 3) * part.py() + D(1, 4) * part.t()  + D(1, 5) * part.pt();
        double ynor =
            D(2, 0) * part.x()  + D(2, 1) * part.px() + D(2, 2) * part.y() +
            D(2, 3) * part.py() + D(2, 4) * part.t()  + D(2, 5) * part.pt();
        double pynor =
            D(3, 0) * part.x()  + D(3, 1) * part.px() + D(3, 2) * part.y() +
            D(3, 3) * part.py() + D(3, 4) * part.t()  + D(3, 5) * part.pt();

        axmax = std::max(axmax, (xnor * xnor + pxnor * pxnor));
        aymax = std::max(aymax, (ynor * ynor + pynor * pynor));
    }
}


template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::setdT(double dt) {
    dt_m = dt;
}


template <class T, unsigned Dim>
double PartBunchBase<T, Dim>::getdT() const {
    return dt_m;
}


template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::setT(double t) {
    t_m = t;
}


1049 1050 1051 1052 1053 1054
template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::incrementT() {
    t_m += dt_m;
}


1055 1056 1057 1058 1059 1060 1061
template <class T, unsigned Dim>
double PartBunchBase<T, Dim>::getT() const {
    return t_m;
}


template <class T, unsigned Dim>
1062
double PartBunchBase<T, Dim>::get_sPos() const {
1063
    return spos_m;
1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134
}


template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::set_sPos(double s) {
    spos_m = s;
}


template <class T, unsigned Dim>
double PartBunchBase<T, Dim>::get_gamma() const {
    return eKin_m / (getM()*1e-6) + 1.0;
}


template <class T, unsigned Dim>
double PartBunchBase<T, Dim>::get_meanKineticEnergy() const {
    return eKin_m;
}


template <class T, unsigned Dim>
Vector_t PartBunchBase<T, Dim>::get_origin() const {
    return rmin_m;
}


template <class T, unsigned Dim>
Vector_t PartBunchBase<T, Dim>::get_maxExtent() const {
    return rmax_m;
}


template <class T, unsigned Dim>
Vector_t PartBunchBase<T, Dim>::get_centroid() const {
    return rmean_m;
}


template <class T, unsigned Dim>
Vector_t PartBunchBase<T, Dim>::get_rrms() const {
    return rrms_m;
}


template <class T, unsigned Dim>
Vector_t PartBunchBase<T, Dim>::get_rprms() const {
    return rprms_m;
}


template <class T, unsigned Dim>
Vector_t PartBunchBase<T, Dim>::get_rmean() const {
    return rmean_m;
}


template <class T, unsigned Dim>
Vector_t PartBunchBase<T, Dim>::get_prms() const {
    return prms_m;
}


template <class T, unsigned Dim>
Vector_t PartBunchBase<T, Dim>::get_pmean() const {
    return pmean_m;
}


template <class T, unsigned Dim>
Vector_t PartBunchBase<T, Dim>::get_pmean_Distribution() const {
1135
    if (dist_m)// && dist_m->getType() != DistrTypeT::FROMFILE)
1136 1137 1138
        return dist_m->get_pmean();

    double gamma = 0.1 / getM() + 1; // set default 0.1 eV
1139
    return Vector_t(0, 0, std::sqrt(std::pow(gamma, 2) - 1));
1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154
}


template <class T, unsigned Dim>
Vector_t PartBunchBase<T, Dim>::get_emit() const {
    return eps_m;
}


template <class T, unsigned Dim>
Vector_t PartBunchBase<T, Dim>::get_norm_emit() const {
    return eps_norm_m;
}


frey_m's avatar
frey_m committed
1155
template <class T, unsigned Dim>
1156 1157
Vector_t PartBunchBase<T, Dim>::get_halo() const {
    return halo_m;
frey_m's avatar
frey_m committed
1158 1159 1160
}


1161 1162 1163 1164 1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186 1187 1188 1189 1190 1191 1192 1193 1194 1195 1196
template <class T, unsigned Dim>
Vector_t PartBunchBase<T, Dim>::get_hr() const {
    return hr_m;
}


template <class T, unsigned Dim>
double PartBunchBase<T, Dim>::get_Dx() const {
    return Dx_m;
}


template <class T, unsigned Dim>
double PartBunchBase<T, Dim>::get_Dy() const {
    return Dy_m;
}


template <class T, unsigned Dim>
double PartBunchBase<T, Dim>::get_DDx() const {
    return DDx_m;
}


template <class T, unsigned Dim>
double PartBunchBase<T, Dim>::get_DDy() const {
    return DDy_m;
}


template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::set_meshEnlargement(double dh) {
    dh_m = dh;
}


1197 1198 1199 1200 1201 1202
template <class T, unsigned Dim>
double PartBunchBase<T, Dim>::getMeshEnlargement() const {
    return dh_m;
}


1203 1204 1205 1206
template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::gatherLoadBalanceStatistics() {

    for(int i = 0; i < Ippl::getNodes(); i++)
frey_m's avatar
frey_m committed
1207
        globalPartPerNode_m[i] = 0;
1208

frey_m's avatar
frey_m committed
1209 1210
    std::size_t localnum = getLocalNum();
    gather(&localnum, &globalPartPerNode_m[0], 1);
1211 1212 1213 1214 1215 1216 1217 1218 1219 1220 1221 1222 1223 1224 1225 1226 1227 1228 1229 1230 1231 1232 1233
}


template <class T, unsigned Dim>
size_t PartBunchBase<T, Dim>::getLoadBalance(int p) const {
    return globalPartPerNode_m[p];
}


template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::get_PBounds(Vector_t &min, Vector_t &max) const {
    bounds(this->P, min, max);
}


template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::calcBeamParameters() {

    Vector_t eps2, fac, rsqsum, psqsum, rpsum;

    IpplTimings::startTimer(statParamTimer_m);

    const size_t locNp = getLocalNum();
1234
    const size_t totalNum = getTotalNum();
1235 1236 1237 1238
    const double zero = 0.0;

    get_bounds(rmin_m, rmax_m);

1239
    if(totalNum == 0) {
1240 1241 1242 1243 1244 1245 1246 1247 1248 1249 1250 1251 1252 1253
        for(unsigned int i = 0 ; i < Dim; i++) {
            rmean_m(i) = 0.0;
            pmean_m(i) = 0.0;
            rrms_m(i) = 0.0;
            prms_m(i) = 0.0;
            eps_norm_m(i)  = 0.0;
        }
        rprms_m = 0.0;
        eKin_m = 0.0;
        eps_m = 0.0;
        IpplTimings::stopTimer(statParamTimer_m);
        return;
    }

1254 1255
    const size_t intN = calcMoments();
    const double N = static_cast<double>(intN);
1256 1257 1258 1259 1260 1261 1262 1263 1264 1265 1266 1267 1268 1269

    for(unsigned int i = 0 ; i < Dim; i++) {
        rmean_m(i) = centroid_m[2 * i] / N;
        pmean_m(i) = centroid_m[(2 * i) + 1] / N;
        rsqsum(i) = moments_m(2 * i, 2 * i) - N * rmean_m(i) * rmean_m(i);
        psqsum(i) = moments_m((2 * i) + 1, (2 * i) + 1) - N * pmean_m(i) * pmean_m(i);
        if(psqsum(i) < 0)
            psqsum(i) = 0;
        rpsum(i) = moments_m((2 * i), (2 * i) + 1) - N * rmean_m(i) * pmean_m(i);
    }
    eps2 = (rsqsum * psqsum - rpsum * rpsum) / (N * N);
    rpsum /= N;

    for