PartBunchBase.hpp 54.5 KB
Newer Older
frey_m's avatar
frey_m committed
1 2 3 4
//
// Class PartBunchBase
//   Base class for representing particle bunches.
//
5
// Copyright (c) 2008 - 2021, Paul Scherrer Institut, Villigen PSI, Switzerland
frey_m's avatar
frey_m committed
6 7 8 9 10 11 12 13 14 15 16 17
// 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/>.
//
18 19 20
#ifndef PART_BUNCH_BASE_HPP
#define PART_BUNCH_BASE_HPP

21
#include <cmath>
frey_m's avatar
frey_m committed
22 23
#include <numeric>

frey_m's avatar
frey_m committed
24 25
#include "Distribution/Distribution.h"

26
#include "AbstractObjects/OpalData.h"
27 28 29
#include "Algorithms/PartBins.h"
#include "Algorithms/PartBinsCyc.h"
#include "Algorithms/PartData.h"
30
#include "Physics/Physics.h"
31
#include "Structure/FieldSolver.h"
32
#include "Utilities/GeneralClassicException.h"
33
#include "Utilities/OpalException.h"
34
#include "Utilities/Options.h"
35
#include "Utilities/SwitcherError.h"
36 37
#include "Utilities/Util.h"

38
extern Inform* gmsg;
39 40

template <class T, unsigned Dim>
41
PartBunchBase<T, Dim>::PartBunchBase(AbstractParticle<T, Dim>* pb, const PartData* ref)
42
    : R(*(pb->R_p)),
43 44 45 46
      ID(*(pb->ID_p)),
      pbin_m(nullptr),
      pmsg_m(nullptr),
      f_stream(nullptr),
47
      fixed_grid(false),
48
      reference(ref),
49 50 51 52 53
      unit_state_(units),
      stateOfLastBoundP_(unitless),
      moments_m(),
      dt_m(0.0),
      t_m(0.0),
54
      spos_m(0.0),
55 56
      globalMeanR_m(Vector_t(0.0, 0.0, 0.0)),
      globalToLocalQuaternion_m(Quaternion_t(1.0, 0.0, 0.0, 0.0)),
57 58 59 60 61 62 63
      rmax_m(0.0),
      rmin_m(0.0),
      hr_m(-1.0),
      nr_m(0),
      fs_m(nullptr),
      couplingConstant_m(0.0),
      qi_m(0.0),
kraus's avatar
kraus committed
64
      massPerParticle_m(0.0),
65 66 67 68 69 70 71 72 73
      distDump_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
74 75
      bunchTotalNum_m(1),
      bunchLocalNum_m(1),
76 77 78 79
      SteptoLastInj_m(0),
      globalPartPerNode_m(nullptr),
      dist_m(nullptr),
      dcBeam_m(false),
80
      periodLength_m(Physics::c / 1e9),
kraus's avatar
kraus committed
81
      pbase_m(pb)
82
{
frey_m's avatar
frey_m committed
83
    setup(pb);
84 85
}

86 87 88 89 90 91
/*
 * Bunch common member functions
 */

template <class T, unsigned Dim>
bool PartBunchBase<T, Dim>::getIfBeamEmitting() {
92 93 94 95 96 97 98 99 100 101 102 103
    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;
}


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


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


121 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
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]);
148
    for (int i = 0; i < numberOfEnergyBins; i++)
149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165
        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) {

166
    if (unit_state_ == unitless)
167 168 169 170
        throw SwitcherError("PartBunch::switchToUnitlessPositions",
                            "Cannot make a unitless PartBunch unitless");

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

173
    for (size_t i = 0; i < getLocalNum(); i++) {
174
        double dt = getdT();
175
        if (use_dt_per_particle)
176 177 178 179 180 181 182
            dt = this->dt[i];

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

    unit_state_ = unitless;

183
    if (hasToReset) R.resetDirtyFlag();
184 185 186 187
}

//FIXME: unify methods, use convention that all particles have own dt
template <class T, unsigned Dim>
188
void PartBunchBase<T, Dim>::switchOffUnitlessPositions(bool use_dt_per_particle) {
189

190
    if (unit_state_ == units)
191 192 193 194
        throw SwitcherError("PartBunch::switchOffUnitlessPositions",
                            "Cannot apply units twice to PartBunch");

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

197
    for (size_t i = 0; i < getLocalNum(); i++) {
198
        double dt = getdT();
199
        if (use_dt_per_particle)
200 201 202 203 204 205 206
            dt = this->dt[i];

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

    unit_state_ = units;

207
    if (hasToReset) R.resetDirtyFlag();
208 209 210 211 212
}


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


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

//    if (Options::cZero)
226
//        dist_m->create(this, addedDistributions, np / 2);
227
//    else
228
//        dist_m->create(this, addedDistributions, np);
229 230 231 232
}


template <class T, unsigned Dim>
233
bool PartBunchBase<T, Dim>::isGridFixed() const {
234 235 236 237
    return fixed_grid;
}


238
template <class T, unsigned Dim>
239
bool PartBunchBase<T, Dim>::hasBinning() const {
240 241 242 243
    return (pbin_m != nullptr);
}


244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266
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 {
267
    if (pbin_m != NULL)
268 269 270 271 272 273 274
        return pbin_m->weHaveBins();
    else
        return false;
}


template <class T, unsigned Dim>
275
void PartBunchBase<T, Dim>::setPBins(PartBins* pbin) {
276 277 278 279 280 281 282
    pbin_m = pbin;
    *gmsg << *pbin_m << endl;
    setEnergyBins(pbin_m->getNBins());
}


template <class T, unsigned Dim>
283
void PartBunchBase<T, Dim>::setPBins(PartBinsCyc* pbin) {
284 285 286 287 288 289 290
    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
291
    return dist_m->emitParticles(this, eZ);
292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313
}


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() {
314
    if (pbin_m != NULL)
315 316 317 318 319 320
        return pbin_m->getLastemittedBin();
    else
        return 0;
}


321 322
template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::setLocalBinCount(size_t num, int bin) {
323
    if (pbin_m != NULL) {
324 325 326 327 328
        pbin_m->setPartNum(bin, num);
    }
}


329 330 331 332 333 334 335
template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::calcGammas() {

    const int emittedBins = dist_m->getNumberOfEnergyBins();

    size_t s = 0;

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

339
    for (unsigned int n = 0; n < getLocalNum(); n++)
340
        bingamma_m[this->Bin[n]] += std::sqrt(1.0 + dot(this->P[n], this->P[n]));
341 342 343 344

    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());
345
    for (int i = 0; i < emittedBins; i++) {
346
        size_t &pInBin = particlesInBin[i];
347
        if (pInBin != 0) {
348 349 350 351 352 353 354 355 356 357 358 359 360 361 362
            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();


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

366 367 368
    if (emittedBins >= 2) {
        for (int i = 1; i < emittedBins; i++) {
            if (binemitted_m[i - 1] != 0 && binemitted_m[i] != 0)
369 370 371 372 373 374 375 376 377 378 379 380
                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();

381
    for (int i = 0; i < emittedBins; i++)
382
        bingamma_m[i] = 0.0;
383
    for (unsigned int n = 0; n < getLocalNum(); n++) {
384
        if ( this->Bin[n] > -1 ) {
385
            bingamma_m[this->Bin[n]] += std::sqrt(1.0 + dot(this->P[n], this->P[n]));
kraus's avatar
kraus committed
386
        }
387
    }
frey_m's avatar
frey_m committed
388 389 390

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

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

}


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) {
412
  this->Q = where(eq(this->Bin, bin), q, 0.0);
413 414 415 416 417
}


template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::setBinCharge(int bin) {
418
  this->Q = where(eq(this->Bin, bin), this->qi_m, 0.0);
419 420 421 422 423 424
}


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

frey_m's avatar
frey_m committed
425
    std::size_t localnum = 0;
426 427
    const Vector_t meanR = get_rmean();

428
    for (unsigned long k = 0; k < getLocalNum(); ++ k)
429 430 431
        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)) {
432

frey_m's avatar
frey_m committed
433
            ++localnum;
434
        }
435

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

frey_m's avatar
frey_m committed
438 439 440
    size_t npOutside = std::accumulate(globalPartPerNode_m.get(),
                                       globalPartPerNode_m.get() + Ippl::getNodes(), 0,
                                       std::plus<size_t>());
441

frey_m's avatar
frey_m committed
442
    return npOutside;
443 444 445
}


446 447 448 449 450 451 452 453 454 455 456
/**
 * \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,
457 458
                                            std::vector<double>& lineDensity,
                                            std::pair<double, double>& meshInfo) {
459 460 461 462
    Vector_t rmin, rmax;
    get_bounds(rmin, rmax);

    if (nBins < 2) {
frey_m's avatar
frey_m committed
463
        Vektor<int, 3>/*NDIndex<3>*/ grid;
464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480
        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;
481
        double tau = (z - zmin) / hz - idx;
482 483 484 485 486 487 488 489 490 491 492

        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
493 494 495 496 497 498 499 500

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

    IpplTimings::startTimer(boundpTimer_m);
501
    //if (!R.isDirty() && stateOfLastBoundP_ == unit_state_) return;
frey_m's avatar
frey_m committed
502 503 504 505
    if ( !(R.isDirty() || ID.isDirty() ) && stateOfLastBoundP_ == unit_state_) return; //-DW

    stateOfLastBoundP_ = unit_state_;

506
    if (!isGridFixed()) {
507
        const int dimIdx = (dcBeam_m? 2: 3);
frey_m's avatar
frey_m committed
508 509

        /**
510 511 512 513
            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
514

515
        this->updateDomainLength(nr_m);
516
        IpplTimings::startTimer(boundpBoundsTimer_m);
frey_m's avatar
frey_m committed
517
        get_bounds(rmin_m, rmax_m);
518
        IpplTimings::stopTimer(boundpBoundsTimer_m);
frey_m's avatar
frey_m committed
519
        Vector_t len = rmax_m - rmin_m;
520

521
        double volume = 1.0;
522
        for (int i = 0; i < dimIdx; i++) {
523 524 525 526 527 528 529
            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
530
            }
531 532 533 534 535 536 537 538 539
            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
540

541 542 543 544 545 546 547
        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;
548

549
                length /= percent;
frey_m's avatar
frey_m committed
550

551 552
                rmax_m[2] += dh_m * length;
                rmin_m[2] -= dh_m * length;
553

554
                hr_m[2] = (rmax_m[2] - rmin_m[2]) / (nr_m[2] - 1);
frey_m's avatar
frey_m committed
555
            }
556
        }
frey_m's avatar
frey_m committed
557

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

563
        if (hr_m[0] * hr_m[1] * hr_m[2] <= 0) {
frey_m's avatar
frey_m committed
564 565
            throw GeneralClassicException("boundp() ", "h<0, can not build a mesh");
        }
566 567 568

        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
569
    }
570
    IpplTimings::startTimer(boundpUpdateTimer_m);
frey_m's avatar
frey_m committed
571
    update();
572
    IpplTimings::stopTimer(boundpUpdateTimer_m);
frey_m's avatar
frey_m committed
573 574 575 576 577 578 579
    R.resetDirtyFlag();

    IpplTimings::stopTimer(boundpTimer_m);
}


template <class T, unsigned Dim>
580
void PartBunchBase<T, Dim>::boundp_destroyCycl() {
frey_m's avatar
frey_m committed
581 582 583 584 585 586

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

    std::unique_ptr<size_t[]> countLost;
587
    if (weHaveBins()) {
frey_m's avatar
frey_m committed
588 589
        const int tempN = pbin_m->getLastemittedBin();
        countLost = std::unique_ptr<size_t[]>(new size_t[tempN]);
590
        for (int ii = 0; ii < tempN; ii++) countLost[ii] = 0;
frey_m's avatar
frey_m committed
591
    }
592

593
    this->updateDomainLength(nr_m);
frey_m's avatar
frey_m committed
594

595
    IpplTimings::startTimer(boundpBoundsTimer_m);
frey_m's avatar
frey_m committed
596
    get_bounds(rmin_m, rmax_m);
frey_m's avatar
frey_m committed
597
    IpplTimings::stopTimer(boundpBoundsTimer_m);
598

frey_m's avatar
frey_m committed
599 600 601 602 603 604 605 606
    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
kraus's avatar
kraus committed
607 608 609
        Vector_t rmean = momentsComputer_m.getMeanPosition();
        Vector_t rrms = momentsComputer_m.getStandardDeviationPosition();
        if(checkfactor < 0) {
frey_m's avatar
frey_m committed
610
            checkfactor *= -1;
kraus's avatar
kraus committed
611 612 613 614 615
            if (len[0] > checkfactor * rrms[0] ||
                len[1] > checkfactor * rrms[1] ||
                len[2] > checkfactor * rrms[2])
            {
                for(unsigned int ii = 0; ii < this->getLocalNum(); ii++) {
616
                    /* delete the particle if the distance to the beam center
frey_m's avatar
frey_m committed
617 618
                     * is larger than 8 times of beam's rms size
                     */
kraus's avatar
kraus committed
619 620 621 622
                    if (std::abs(R[ii](0) - rmean(0)) > checkfactor * rrms[0] ||
                        std::abs(R[ii](1) - rmean(1)) > checkfactor * rrms[1] ||
                        std::abs(R[ii](2) - rmean(2)) > checkfactor * rrms[2])
                    {
frey_m's avatar
frey_m committed
623 624 625 626 627 628 629 630 631 632 633
                        // 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;);
                         */
                    }
                }
            }
kraus's avatar
kraus committed
634 635 636 637 638 639
        }
        else {
            if (len[0] > checkfactor * rrms[0] ||
                len[2] > checkfactor * rrms[2])
            {
                for(unsigned int ii = 0; ii < this->getLocalNum(); ii++) {
640
                    /* delete the particle if the distance to the beam center
frey_m's avatar
frey_m committed
641 642
                     * is larger than 8 times of beam's rms size
                     */
kraus's avatar
kraus committed
643 644 645
                    if (std::abs(R[ii](0) - rmean(0)) > checkfactor * rrms[0] ||
                        std::abs(R[ii](2) - rmean(2)) > checkfactor * rrms[2])
                    {
frey_m's avatar
frey_m committed
646 647 648 649 650 651 652 653 654 655 656 657 658 659
                        // 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;);
                         */
                    }
                }
            }
        }
    }

660
    for (int i = 0; i < dimIdx; i++) {
frey_m's avatar
frey_m committed
661 662 663 664 665 666 667
        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
668
    this->updateFields(hr_m, rmin_m);
669

670
    if (weHaveBins()) {
frey_m's avatar
frey_m committed
671 672
        pbin_m->updatePartInBin_cyc(countLost.get());
    }
673

frey_m's avatar
frey_m committed
674 675 676 677 678
    /* we also need to update the number of particles per bunch
     * expensive since does an allreduce!
     */
    countTotalNumPerBunch();

679
    IpplTimings::startTimer(boundpUpdateTimer_m);
frey_m's avatar
frey_m committed
680
    update();
681
    IpplTimings::stopTimer(boundpUpdateTimer_m);
682

frey_m's avatar
frey_m committed
683 684 685 686 687 688 689
    IpplTimings::stopTimer(boundpTimer_m);
}


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

690
    this->updateDomainLength(nr_m);
frey_m's avatar
frey_m committed
691

692
    std::vector<size_t> tmpbinemitted;
frey_m's avatar
frey_m committed
693 694 695 696 697 698

    boundp();

    size_t ne = 0;
    const size_t localNum = getLocalNum();

699 700 701 702 703 704 705 706 707 708 709 710
    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]]++;
frey_m's avatar
frey_m committed
711 712 713
        }
    }

714 715
    boundp();

frey_m's avatar
frey_m committed
716 717 718
    calcBeamParameters();
    gatherLoadBalanceStatistics();

719
    if (haveEnergyBins) {
frey_m's avatar
frey_m committed
720
        const int lastBin = dist_m->getLastEmittedEnergyBin();
721
        for (int i = 0; i < lastBin; i++) {
frey_m's avatar
frey_m committed
722 723 724 725 726 727 728 729 730 731 732 733 734 735 736
            binemitted_m[i] = tmpbinemitted[i];
        }
    }
    reduce(ne, ne, OpAddAssign());
    return ne;
}


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

    std::unique_ptr<size_t[]> tmpbinemitted;

    const size_t localNum = getLocalNum();
    const size_t totalNum = getTotalNum();
737
    size_t ne = 0;
frey_m's avatar
frey_m committed
738

739
    if (weHaveEnergyBins()) {
frey_m's avatar
frey_m committed
740
        tmpbinemitted = std::unique_ptr<size_t[]>(new size_t[getNumberOfEnergyBins()]);
741
        for (int i = 0; i < getNumberOfEnergyBins(); i++) {
frey_m's avatar
frey_m committed
742 743
            tmpbinemitted[i] = 0;
        }
744
        for (unsigned int i = 0; i < localNum; i++) {
frey_m's avatar
frey_m committed
745 746
            if (Bin[i] < 0) {
                destroy(1, i);
747
                ++ ne;
frey_m's avatar
frey_m committed
748 749 750 751 752
            } else
                tmpbinemitted[Bin[i]]++;
        }
    } else {
        Inform dmsg("destroy: ", INFORM_ALL_NODES);
753 754
        for (size_t i = 0; i < localNum; i++) {
            if ((Bin[i] < 0)) {
755 756
                ne++;
                destroy(1, i);
frey_m's avatar
frey_m committed
757 758
            }
        }
759 760 761 762
    }

    if (ne > 0) {
        performDestroy(true);
frey_m's avatar
frey_m committed
763 764 765 766 767 768 769
    }

    calcBeamParameters();
    gatherLoadBalanceStatistics();

    if (weHaveEnergyBins()) {
        const int lastBin = dist_m->getLastEmittedEnergyBin();
770
        for (int i = 0; i < lastBin; i++) {
frey_m's avatar
frey_m committed
771 772 773 774 775 776 777 778 779 780 781
            binemitted_m[i] = tmpbinemitted[i];
        }
    }
    size_t newTotalNum = getLocalNum();
    reduce(newTotalNum, newTotalNum, OpAddAssign());

    setTotalNum(newTotalNum);

    return totalNum - newTotalNum;
}

782

783
template <class T, unsigned Dim>
784
double PartBunchBase<T, Dim>::getPx(int /*i*/) {
785 786 787 788
    return 0.0;
}

template <class T, unsigned Dim>
789
double PartBunchBase<T, Dim>::getPy(int) {
790 791 792 793
    return 0.0;
}

template <class T, unsigned Dim>
794
double PartBunchBase<T, Dim>::getPz(int) {
795 796 797 798
    return 0.0;
}

template <class T, unsigned Dim>
799
double PartBunchBase<T, Dim>::getPx0(int) {
800 801 802 803
    return 0.0;
}

template <class T, unsigned Dim>
804
double PartBunchBase<T, Dim>::getPy0(int) {
805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827
    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>
828
double PartBunchBase<T, Dim>::getX0(int /*i*/) {
829 830 831 832 833
    return 0.0;
}

//ff
template <class T, unsigned Dim>
834
double PartBunchBase<T, Dim>::getY0(int /*i*/) {
835 836 837 838 839
    return 0.0;
}


template <class T, unsigned Dim>
840
void PartBunchBase<T, Dim>::setZ(int /*i*/, double /*zcoo*/) {
841 842 843 844
};


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

frey_m's avatar
frey_m committed
847
    this->getLocalBounds(rmin, rmax);
848

kraus's avatar
kraus committed
849
    double min[2*Dim];
850

851
    for (unsigned int i = 0; i < Dim; ++i) {
kraus's avatar
kraus committed
852 853
        min[2*i] = rmin[i];
        min[2*i + 1] = -rmax[i];
854
    }
855

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

858
    for (unsigned int i = 0; i < Dim; ++i) {
kraus's avatar
kraus committed
859 860
        rmin[i] = min[2*i];
        rmax[i] = -min[2*i + 1];
861
    }
862 863 864 865
}


template <class T, unsigned Dim>
866
void PartBunchBase<T, Dim>::getLocalBounds(Vector_t& rmin, Vector_t& rmax) {
867
    const size_t localNum = getLocalNum();
frey_m's avatar
frey_m committed
868
    if (localNum == 0) {
869 870 871
        double maxValue = 1e8;
        rmin = Vector_t(maxValue, maxValue, maxValue);
        rmax = Vector_t(-maxValue, -maxValue, -maxValue);
kraus's avatar
kraus committed
872
        return;
frey_m's avatar
frey_m committed
873
    }
874

875 876 877 878 879
    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
880
            else if (rmax(d) < R[i](d)) rmax(d) = R[i](d);
881 882 883 884 885 886 887 888 889 890 891 892
        }
    }
}


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);
893
    sphere.second = std::sqrt(dot(rmax - sphere.first, rmax - sphere.first));
894 895 896 897 898 899 900 901 902 903 904 905

    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);
906
    sphere.second = std::sqrt(dot(rmax - sphere.first, rmax - sphere.first));
907 908 909 910 911 912

    return sphere;
}


template <class T, unsigned Dim>
kraus's avatar
kraus committed
913
void PartBunchBase<T, Dim>::push_back(OpalParticle const& particle) {
914 915
    Inform msg("PartBunch ");

kraus's avatar
kraus committed
916
    size_t i = getLocalNum();
917 918
    create(1);

kraus's avatar
kraus committed
919 920
    R[i] = particle.getR();
    P[i] = particle.getP();
921 922 923 924 925 926 927

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


template <class T, unsigned Dim>
kraus's avatar
kraus committed
928
void PartBunchBase<T, Dim>::setParticle(FVector<double, 6> z, int ii) {
929 930 931 932 933 934 935 936 937 938
    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>
kraus's avatar
kraus committed
939 940 941
void PartBunchBase<T, Dim>::setParticle(OpalParticle const& particle, int ii) {
    R[ii] = particle.getR();
    P[ii] = particle.getP();
942 943 944 945
}


template <class T, unsigned Dim>
kraus's avatar
kraus committed
946 947 948 949 950
OpalParticle PartBunchBase<T, Dim>::getParticle(int ii) {
    OpalParticle particle(ID[ii],
                          Vector_t(R[ii](0), R[ii](1), 0), P[ii],
                          R[ii](2), Q[ii], M[ii]);
    return particle;
951 952 953 954
}


template <class T, unsigned Dim>
955 956
void PartBunchBase<T, Dim>::maximumAmplitudes(const FMatrix<double, 6, 6>& D,
                                              double& axmax, double& aymax) {
957
    axmax = aymax = 0.0;
kraus's avatar
kraus committed
958
    OpalParticle particle;
959

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

kraus's avatar
kraus committed
962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982
        particle = getParticle(ii);
        FVector<double, 6> vec({particle.getX(), particle.getPx(),
                                particle.getY(), particle.getPy(),
                                particle.getZ(), particle.getPz()});
        FVector<double, 6> result;
        result = D * vec;
        // double xnor =
        //     D(0, 0) * part.getX()  + D(0, 1) * part.getPx() + D(0, 2) * part.getY() +
        //     D(0, 3) * part.getPy() + D(0, 4) * part.getL()  + D(0, 5) * part.getPLon();
        // double pxnor =
        //     D(1, 0) * part.getX()  + D(1, 1) * part.getPx() + D(1, 2) * part.getY() +
        //     D(1, 3) * part.getPy() + D(1, 4) * part.getL()  + D(1, 5) * part.getPLon();
        // double ynor =
        //     D(2, 0) * part.getX()  + D(2, 1) * part.getPx() + D(2, 2) * part.getY() +
        //     D(2, 3) * part.getPy() + D(2, 4) * part.getL()  + D(2, 5) * part.getPLon();
        // double pynor =
        //     D(3, 0) * part.getX()  + D(3, 1) * part.getPx() + D(3, 2) * part.getY() +
        //     D(3, 3) * part.getPy() + D(3, 4) * part.getL()  + D(3, 5) * part.getPLon();

        axmax = std::max(axmax, (std::pow(result[0], 2) + std::pow(result[1], 2)));
        aymax = std::max(aymax, (std::pow(result[2], 2) + std::pow(result[3], 2)));
983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004
    }
}


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;
}


1005 1006 1007 1008 1009 1010
template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::incrementT() {
    t_m += dt_m;
}


1011 1012 1013 1014 1015 1016 1017
template <class T, unsigned Dim>
double PartBunchBase<T, Dim>::getT() const {
    return t_m;
}


template <class T, unsigned Dim>