Commit 548a5820 authored by frey_m's avatar frey_m
Browse files

simplify

parent c0da8d5d
......@@ -15,17 +15,20 @@ MultiBunchHandler::MultiBunchHandler(PartBunchBase<double, 3> *beam,
const double& eta,
const double& para,
const std::string& mode,
const std::string& binning,
DataSink& ds)
const std::string& binning)
: onebunch_m(OpalData::getInstance()->getInputBasename() + "-onebunch.h5")
, numBunch_m(numBunch)
, eta_m(0.01)
, coeffDBunches_m(para)
, radiusLastTurn_m(0.0)
, radiusThisTurn_m(0.0)
, ds_m(ds)
, bunchCount_m(1)
{
binfo_m.reserve(numBunch);
for (int i = 0; i < numBunch; ++i) {
binfo_m.push_back(beaminfo_t(0.0, 0.0, 0));
}
this->setBinning(binning);
if ( numBunch > 1 ) {
......@@ -140,10 +143,10 @@ void MultiBunchHandler::saveBunch(PartBunchBase<double, 3> *beam,
h5wrapper.writeStep(beam, additionalAttributes);
h5wrapper.close();
// // injection values
// azimuth_m = azimuth;
// time_m = beam->getT();
// spos_m = beam->getLPath();
// injection values
double time = beam->getT();
double lpath = beam->getLPath();
binfo_m.push_back(beaminfo_t(time, azimuth, lpath));
*gmsg << "Done." << endl;
IpplTimings::stopTimer(saveBunchTimer);
......@@ -169,6 +172,10 @@ bool MultiBunchHandler::readBunch(PartBunchBase<double, 3> *beam,
size_t numParticles = h5wrapper.getNumParticles();
const int bunchNum = bunchCount_m - 1;
beam->setTotalNumPerBunch(numParticles, bunchNum);
if ( numParticles == 0 ) {
throw OpalException("MultiBunchHandler::readBunch()",
"No particles in file " + onebunch_m + ".");
......@@ -185,6 +192,8 @@ bool MultiBunchHandler::readBunch(PartBunchBase<double, 3> *beam,
numParticles = lastParticle - firstParticle + 1;
beam->setLocalNumPerBunch(numParticles, bunchNum);
//FIXME
std::unique_ptr<PartBunchBase<double, 3> > tmpBunch = 0;
#ifdef ENABLE_AMR
......@@ -202,8 +211,7 @@ bool MultiBunchHandler::readBunch(PartBunchBase<double, 3> *beam,
h5wrapper.close();
beam->create(numParticles);
const int bunchNum = beam->getNumBunch() - 1;
for(size_t ii = 0; ii < numParticles; ++ ii, ++ localNum) {
beam->R[localNum] = tmpBunch->R[ii];
beam->P[localNum] = tmpBunch->P[ii];
......@@ -231,6 +239,10 @@ short MultiBunchHandler::injectBunch(PartBunchBase<double, 3> *beam,
short result = 0;
if ((bunchCount_m == 1) && (mode_m == MB_MODE::AUTO) && (!flagTransition)) {
// we have still a single bunch
beam->setTotalNumPerBunch(beam->getTotalNum(), 0);
beam->setLocalNumPerBunch(beam->getLocalNum(), 0);
// If all of the following conditions are met, this code will be executed
// to check the distance between two neighboring bunches:
// 1. Only one bunch exists (bunchCount_m == 1)
......@@ -289,6 +301,7 @@ short MultiBunchHandler::injectBunch(PartBunchBase<double, 3> *beam,
case MB_MODE::AUTO:
readBunch(beam, ref);
updateParticleBins(beam);
calcBunchBeamParameters(beam, bunchCount_m - 1);
break;
default:
throw OpalException("MultiBunchHandler::injectBunch()",
......@@ -375,3 +388,145 @@ void MultiBunchHandler::setRadiusTurns(const double& radius) {
// New OPAL 2.0: Init in m -DW
*gmsg << radiusThisTurn_m << " m" << endl;
}
bool MultiBunchHandler::calcBunchBeamParameters(PartBunchBase<double, 3>* beam,
short bunchNr)
{
if ( !OpalData::getInstance()->isInOPALCyclMode() ) {
return false;
}
const unsigned long localNum = beam->getLocalNum();
long int bunchTotalNum = 0;
unsigned long bunchLocalNum = 0;
/* container:
*
* ekin, <x>, <y>, <z>, <p_x>, <p_y>, <p_z>,
* <x^2>, <y^2>, <z^2>, <p_x^2>, <p_y^2>, <p_z^2>,
* <xp_x>, <y_py>, <zp_z>,
* <x^3>, <y^3>, <z^3>, <x^4>, <y^4>, <z^4>
*/
const unsigned int dim = PartBunchBase<double, 3>::Dimension;
std::vector<double> local(7 * dim + 1);
beaminfo_t& binfo = getBunchInfo(bunchNr);
for(unsigned long k = 0; k < localNum; ++ k) {
if ( beam->bunchNum[k] != bunchNr ) { //|| ID[k] == 0 ) {
continue;
}
++bunchTotalNum;
++bunchLocalNum;
// ekin
local[0] += std::sqrt(dot(beam->P[k], beam->P[k]) + 1.0);
for (unsigned int i = 0; i < dim; ++i) {
double r = beam->R[k](i);
double p = beam->P[k](i);
// <x>, <y>, <z>
local[i + 1] += r;
// <p_x>, <p_y, <p_z>
local[i + dim + 1] += p;
// <x^2>, <y^2>, <z^2>
double r2 = r * r;
local[i + 2 * dim + 1] += r2;
// <p_x^2>, <p_y^2>, <p_z^2>
local[i + 3 * dim + 1] += p * p;
// <xp_x>, <y_py>, <zp_z>
local[i + 4 * dim + 1] += r * p;
// <x^3>, <y^3>, <z^3>
local[i + 5 * dim + 1] += r2 * r;
// <x^4>, <y^4>, <z^4>
local[i + 6 * dim + 1] += r2 * r2;
}
}
// inefficient
allreduce(bunchTotalNum, 1, std::plus<long int>());
// here we also update the number of particles of *this* bunch
if (bunchNr >= (short)beam->getNumBunch())
throw OpalException("PartBunchBase::calcBunchBeamParameters()",
"Bunch number " + std::to_string(bunchNr) +
" exceeds bunch index " + std::to_string(beam->getNumBunch() - 1));
beam->setTotalNumPerBunch(bunchTotalNum, bunchNr);
beam->setLocalNumPerBunch(bunchLocalNum, bunchNr);
if ( bunchTotalNum == 0 )
return false;
// ekin
const double m0 = beam->getM() * 1.0e-6;
local[0] -= bunchLocalNum;
local[0] *= m0;
allreduce(local.data(), local.size(), std::plus<double>());
double invN = 1.0 / double(bunchTotalNum);
binfo.ekin = local[0] * invN;
binfo.time = beam->getT() * 1e9; // ns
binfo.nParticles = bunchTotalNum;
for (unsigned int i = 0; i < dim; ++i) {
double w = local[i + 1] * invN;
double pw = local[i + dim + 1] * invN;
double w2 = local[i + 2 * dim + 1] * invN;
double pw2 = local[i + 3 * dim + 1] * invN;
double wpw = local[i + 4 * dim + 1] * invN;
double w3 = local[i + 5 * dim + 1] * invN;
double w4 = local[i + 6 * dim + 1] * invN;
// <x>, <y>, <z>
binfo.mean[i] = w;
// central: <p_w^2> - <p_w>^2 (w = x, y, z)
binfo.prms[i] = pw2 - pw * pw;
if ( binfo.prms[i] < 0 ) {
binfo.prms[i] = 0.0;
}
// central: <wp_w> - <w><p_w>
wpw = wpw - w * pw;
// central: <w^2> - <w>^2 (w = x, y, z)
binfo.rrms[i] = w2 - w * w;
// central: normalized emittance
binfo.emit[i] = (binfo.rrms[i] * binfo.prms[i] - wpw * wpw);
binfo.emit[i] = std::sqrt(std::max(binfo.emit[i], 0.0));
// central: <w^4> - 4 * <w> * <w^3> + 6 * <w>^2 * <w^2> - 3 * <w>^4
double tmp = w4
- 4.0 * w * w3
+ 6.0 * w * w * w2
- 3.0 * w * w * w * w;
binfo.halo[i] = tmp / ( binfo.rrms[i] * binfo.rrms[i] );
// central: sqrt(<w^2> - <w>^2) (w = x, y, z)
binfo.rrms[i] = std::sqrt(binfo.rrms[i]);
// central: sqrt(<p_w^2> - <p_w>^2)
binfo.prms[i] = std::sqrt(binfo.prms[i]);
}
double tmp = 1.0 / std::pow(binfo.ekin / m0 + 1.0, 2.0);
binfo.dEkin = binfo.prms[1] * m0 * std::sqrt(1.0 - tmp);
return true;
}
......@@ -2,7 +2,9 @@
#define OPAL_MULTI_BUNCH_HANDLER_H
#include "Algorithms/PartBunchBase.h"
#include "Structure/DataSink.h"
#include <cassert>
#include <vector>
/* Helper class that stores bunch injection
* information like azimuth, radius etc. of first
......@@ -10,6 +12,31 @@
*/
class MultiBunchHandler {
public:
struct beaminfo_t {
beaminfo_t(double t,
double theta,
double lpath)
: time(t)
, azimuth(theta)
, prevAzimuth(-1.0)
, pathlength(lpath) {
};
double time;
double azimuth;
double radius;
double prevAzimuth;
double pathlength;
long unsigned int nParticles;
double ekin;
double dEkin;
double rrms[3];
double prms[3];
double emit[3];
double mean[3];
double halo[3];
};
// multi-bunch modes
enum class MB_MODE {
FORCE = 0,
......@@ -34,8 +61,7 @@ public:
const double& eta,
const double& para,
const std::string& mode,
const std::string& binning,
DataSink& ds);
const std::string& binning);
void saveBunch(PartBunchBase<double, 3> *beam,
const double& azimuth);
......@@ -75,17 +101,11 @@ public:
bool isForceMode() const;
void setPathLength(const double& pathlength, short bunchNr);
const double& getPathLength(short bunchNr) const;
void setAzimuth(const double& azimuth, short bunchNr);
const double& getAzimuth(short bunchNr) const;
bool calcBunchBeamParameters(PartBunchBase<double, 3>* beam, short bunchNr);
void setTime(const double& time, short bunchNr);
beaminfo_t& getBunchInfo(short bunchNr);
const double& getTime(short bunchNr) const;
const beaminfo_t& getBunchInfo(short bunchNr) const;
private:
// store the data of the beam which are required for injecting a
......@@ -114,10 +134,11 @@ private:
double radiusLastTurn_m;
double radiusThisTurn_m;
DataSink& ds_m;
// record how many bunches has already been injected.
short bunchCount_m;
// each list entry belongs to a bunch
std::vector<beaminfo_t> binfo_m;
};
......@@ -152,44 +173,16 @@ bool MultiBunchHandler::isForceMode() const {
inline
void MultiBunchHandler::setPathLength(const double& pathlength, short bunchNr) {
MultiBunchDump *mbd = ds_m.getMultiBunchWriter(bunchNr);
mbd->setPathLength(pathlength);
}
inline
const double& MultiBunchHandler::getPathLength(short bunchNr) const {
MultiBunchDump *mbd = ds_m.getMultiBunchWriter(bunchNr);
return mbd->getPathLength();
}
inline
void MultiBunchHandler::setAzimuth(const double& azimuth, short bunchNr) {
MultiBunchDump *mbd = ds_m.getMultiBunchWriter(bunchNr);
mbd->setAzimuth(azimuth);
}
inline
const double& MultiBunchHandler::getAzimuth(short bunchNr) const {
MultiBunchDump *mbd = ds_m.getMultiBunchWriter(bunchNr);
return mbd->getAzimuth();
}
inline
void MultiBunchHandler::setTime(const double& time, short bunchNr) {
MultiBunchDump *mbd = ds_m.getMultiBunchWriter(bunchNr);
mbd->setTime(time);
MultiBunchHandler::beaminfo_t& MultiBunchHandler::getBunchInfo(short bunchNr) {
assert(bunchNr < 0 || bunchNr >= (short)binfo_m.size());
return binfo_m[bunchNr];
}
inline
const double& MultiBunchHandler::getTime(short bunchNr) const {
MultiBunchDump *mbd = ds_m.getMultiBunchWriter(bunchNr);
return mbd->getTime();
const MultiBunchHandler::beaminfo_t& MultiBunchHandler::getBunchInfo(short bunchNr) const {
assert(bunchNr < 0 || bunchNr >= (short)binfo_m.size());
return binfo_m[bunchNr];
}
#endif
......@@ -22,6 +22,7 @@
#include <iostream>
#include <limits>
#include <vector>
#include <numeric>
#include "AbstractObjects/Element.h"
#include "AbstractObjects/OpalData.h"
......@@ -130,7 +131,7 @@ ParallelCyclotronTracker::ParallelCyclotronTracker(const Beamline &beamline,
if ( numBunch > 1 ) {
mbHandler_m = std::unique_ptr<MultiBunchHandler>(
new MultiBunchHandler(bunch, numBunch, mbEta,
mbPara, mbMode, mbBinning, ds)
mbPara, mbMode, mbBinning)
);
}
......@@ -212,7 +213,11 @@ void ParallelCyclotronTracker::dumpAngle(const double& theta,
double& azimuth)
{
if ( prevAzimuth < 0.0 ) { // only at first occurrence
azimuth = theta;
double plus = 0.0;
if ( OpalData::getInstance()->inRestartRun() ) {
plus = 360.0 * turnnumber_m;
}
azimuth_m = theta + plus;
} else {
double dtheta = theta - prevAzimuth;
if ( dtheta < 0 ) {
......@@ -227,6 +232,12 @@ void ParallelCyclotronTracker::dumpAngle(const double& theta,
}
double ParallelCyclotronTracker::computeRadius(const Vector_t& meanR) const {
// New OPAL 2.0: m --> mm
return 1000.0 * std::sqrt(meanR(0) * meanR(0) + meanR(1) * meanR(1));
}
double ParallelCyclotronTracker::computePathLengthUpdate(const double& dt,
short bunchNr)
{
......@@ -381,7 +392,9 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) {
*gmsg << "* (This is not an issue for spiral inflectors as they are typically < 100 keV/amu.)" << endl;
*gmsg << endl << "* Note: For now, multi-bunch mode (MBM) needs to be de-activated for spiral inflector" << endl;
*gmsg << "* and space charge needs to be solved every time-step. numBunch_m and scSolveFreq are reset." << endl;
mbHandler_m->setTotalNumBunch(1);
if (isMultiBunch()) {
mbHandler_m = nullptr;
}
}
// Fresh run (no restart):
......@@ -2071,7 +2084,8 @@ bool ParallelCyclotronTracker::deleteParticle(bool flagNeedUpdate){
allreduce(flagNeedUpdate, 1, std::logical_or<bool>());
if(flagNeedUpdate) {
size_t locLostParticleNum = 0;
short bunchCount = itsBunch_m->getNumBunch();
std::vector<size_t> locLostParticleNum(bunchCount, 0);
const int leb = itsBunch_m->getLastemittedBin();
std::unique_ptr<size_t[]> localBinCount;
......@@ -2084,7 +2098,7 @@ bool ParallelCyclotronTracker::deleteParticle(bool flagNeedUpdate){
for(unsigned int i = 0; i < itsBunch_m->getLocalNum(); i++) {
if(itsBunch_m->Bin[i] < 0) {
++locLostParticleNum;
++locLostParticleNum[itsBunch_m->bunchNum[i]];
itsBunch_m->destroy(1, i);
} else if ( isMultiBunch() ) {
/* we need to count the local number of particles
......@@ -2101,26 +2115,54 @@ bool ParallelCyclotronTracker::deleteParticle(bool flagNeedUpdate){
}
}
std::vector<size_t> localnum(bunchCount + 1);
for (size_t i = 0; i < localnum.size() - 1; ++i) {
localnum[i] = itsBunch_m->getLocalNumPerBunch(i) - locLostParticleNum[i];
itsBunch_m->setLocalNumPerBunch(localnum[i], i);
}
/* We need to destroy the particles now
* before we compute the means. We also
* have to update the total number of particles
* otherwise the statistics are wrong.
*/
itsBunch_m->performDestroy(true);
size_t totalnum = 0;
size_t localnum = itsBunch_m->getLocalNum();
allreduce(&localnum, &totalnum, 1, std::plus<size_t>());
itsBunch_m->setTotalNum(totalnum);
/* total number of particles of individual bunches
* last index of vector contains total number over all
* bunches, used as a check
*/
std::vector<size_t> totalnum(bunchCount + 1);
localnum[bunchCount] = itsBunch_m->getLocalNum();
allreduce(localnum.data(), totalnum.data(), localnum.size(), std::plus<size_t>());
itsBunch_m->setTotalNum(totalnum[bunchCount]);
for (short i = 0; i < bunchCount; ++i)
itsBunch_m->setTotalNumPerBunch(totalnum[i], i);
size_t sum = std::accumulate(totalnum.begin(),
totalnum.end() - 1, 0);
if ( sum != totalnum[bunchCount] ) {
throw OpalException("ParallelCyclotronTracker::deleteParticle()",
"Total number of particles " + std::to_string(totalnum[bunchCount]) +
" != " + std::to_string(sum) + " (sum over all bunches)");
}
size_t globLostParticleNum = 0;
reduce(locLostParticleNum, globLostParticleNum, 1, std::plus<size_t>());
size_t locNumLost = std::accumulate(locLostParticleNum.begin(),
locLostParticleNum.end(), 0);
reduce(locNumLost, globLostParticleNum, 1, std::plus<size_t>());
*gmsg << "At step " << step_m
<< ", lost " << globLostParticleNum << " particles "
<< "on stripper, collimator, septum, or out of cyclotron aperture"
<< endl;
if (totalnum == 0) {
if (totalnum[bunchCount] == 0) {
IpplTimings::stopTimer(DelParticleTimer_m);
return flagNeedUpdate;
}
......@@ -2270,6 +2312,9 @@ void ParallelCyclotronTracker::initDistInGlobalFrame() {
}
}
// set the number of particles per bunch
itsBunch_m->countTotalNumPerBunch();
// ------------ Get some Values ---------------------------------------------------------- //
Vector_t const meanR = calcMeanR();
Vector_t const meanP = calcMeanP();
......@@ -2475,6 +2520,8 @@ void ParallelCyclotronTracker::bunchDumpStatData(){
// fix azimuth_m --> make monotonically increasing
dumpAngle(theta, prevAzimuth_m, azimuth_m);
updateAzimuthAndRadius();
if(Options::psDumpFrame != Options::GLOBAL) {
Vector_t meanP = calcMeanP();
......@@ -2491,7 +2538,13 @@ void ParallelCyclotronTracker::bunchDumpStatData(){
globalToLocal(itsBunch_m->P, phi, psi);
}
itsDataSink->writeMultiBunchStatistics(itsBunch_m, azimuth_m);
for (short b = 0; b < mbHandler_m->getNumBunch(); ++b) {
bool isOk = mbHandler_m->calcBunchBeamParameters(itsBunch_m, b);
const MultiBunchHandler::beaminfo_t& binfo = mbHandler_m->getBunchInfo(b);
if (isOk) {
itsDataSink->writeMultiBunchStatistics(itsBunch_m, binfo);
}
}
if(Options::psDumpFrame != Options::GLOBAL) {
localToGlobal(itsBunch_m->R, phi, psi, meanR);
......@@ -2619,7 +2672,7 @@ void ParallelCyclotronTracker::bunchDumpPhaseSpaceData() {
// ---------------- Re-calculate reference values in format of input values ----------------- //
// Position:
// New OPAL 2.0: Init in m (back to mm just for output) -DW
referenceR = 1000.0 * sqrt(meanR(0) * meanR(0) + meanR(1) * meanR(1));
referenceR = computeRadius(meanR); // includes m --> mm conversion
referenceTheta = theta / Physics::deg2rad;
referenceZ = 1000.0 * meanR(2);
......@@ -2775,6 +2828,8 @@ std::tuple<double, double, double> ParallelCyclotronTracker::initializeTracking_
initTrackOrbitFile();
turnnumber_m = 1;
// Get data from h5 file for restart run and reset current step
// to last step of previous simulation
if(OpalData::getInstance()->inRestartRun()) {
......@@ -2782,7 +2837,10 @@ std::tuple<double, double, double> ParallelCyclotronTracker::initializeTracking_
restartStep0_m = itsBunch_m->getLocalTrackStep();
step_m = restartStep0_m;
*gmsg << "* Restart at integration step " << restartStep0_m << endl;
turnnumber_m = step_m / setup_m.stepsPerTurn;
*gmsg << "* Restart at integration step " << restartStep0_m
<< " at turn " << turnnumber_m << endl;
}
setup_m.stepsNextCheck = step_m + setup_m.stepsPerTurn; // Steps to next check for transition
......@@ -3376,9 +3434,8 @@ void ParallelCyclotronTracker::updatePathLength(const double& dt) {
if ( isMultiBunch() ) {
for (short b = 0; b < mbHandler_m->getNumBunch(); ++b) {
double lpath = mbHandler_m->getPathLength(b);
lpath += computePathLengthUpdate(dt, b);
mbHandler_m->setPathLength(lpath, b);
MultiBunchHandler::beaminfo_t& binfo = mbHandler_m->getBunchInfo(b);
binfo.pathlength += computePathLengthUpdate(dt, b);
}
}
}
......@@ -3393,22 +3450,20 @@ void ParallelCyclotronTracker::updateTime(const double& dt) {
if ( isMultiBunch() ) {
for (short b = 0; b < mbHandler_m->getNumBunch(); ++b) {
double time = mbHandler_m->getTime(b);
mbHandler_m->setTime(time + dt, b);
MultiBunchHandler::beaminfo_t& binfo = mbHandler_m->getBunchInfo(b);
binfo.time += dt;
}
}
}
void ParallelCyclotronTracker::updateAzimuth() {
void ParallelCyclotronTracker::updateAzimuthAndRadius() {
for (short b = 0; b < mbHandler_m->getNumBunch(); ++b) {
Vector_t meanR = calcMeanR(b);
double azimuth = calculateAngle(meanR(0), meanR(1)) * Physics::rad2deg;
const double& prevAzimuth = mbHandler_m->getAzimuth();
double newAzimuth = prevAzimuth;
MultiBunchHandler::beaminfo_t& binfo = mbHandler_m->getBunchInfo(b);
dumpAngle(azimuth, prevAzimuth, newAzimuth);
mbHandler_m->setAzimuth(newAzimuth);
binfo.radius = computeRadius(meanR);
double azimuth = calculateAngle(meanR(0), meanR(1)) * Physics::rad2deg;
dumpAngle(azimuth, binfo.prevAzimuth, binfo.azimuth);
}
}