Skip to content

GitLab

Commit 8c4e57c7 authored by frey_m's avatar frey_m
Browse files

save injection values differently

parent 354de3c5
Hide whitespace changes
Inline Side-by-side
Showing with 109 additions and 61 deletions
src/Algorithms/MultiBunchHandler.cpp
View file @ 8c4e57c7
......@@ -23,13 +23,10 @@ MultiBunchHandler::MultiBunchHandler(PartBunchBase<double, 3> *beam,
, radiusLastTurn_m(0.0)
, radiusThisTurn_m(0.0)
, bunchCount_m(1)
, injTime_m(0.0)
, injPathlength_m(0.0)
, injAzimuth_m(0.0)
{
binfo_m.reserve(numBunch);
for (int i = 0; i < beam->getNumBunch(); ++i) {
binfo_m.push_back(beaminfo_t(beam->getT(), beam->getLPath(), 0.0));
binfo_m.push_back(beaminfo_t());
}
this->setBinning(binning);
......@@ -80,8 +77,7 @@ MultiBunchHandler::MultiBunchHandler(PartBunchBase<double, 3> *beam,
}
void MultiBunchHandler::saveBunch(PartBunchBase<double, 3> *beam,
const double& azimuth)
void MultiBunchHandler::saveBunch(PartBunchBase<double, 3> *beam)
{
if ( numBunch_m < 2 )
return;
......@@ -146,11 +142,6 @@ void MultiBunchHandler::saveBunch(PartBunchBase<double, 3> *beam,
h5wrapper.writeStep(beam, additionalAttributes);
h5wrapper.close();
// injection values
injTime_m = beam->getT();
injPathlength_m = beam->getLPath();
injAzimuth_m = azimuth;
*gmsg << "Done." << endl;
IpplTimings::stopTimer(saveBunchTimer);
}
......@@ -227,7 +218,7 @@ bool MultiBunchHandler::readBunch(PartBunchBase<double, 3> *beam,
beam->boundp();
binfo_m.push_back(beaminfo_t(injTime_m, injPathlength_m, injAzimuth_m));
binfo_m.push_back(beaminfo_t(injection_m));
*gmsg << "Done." << endl;
......@@ -238,8 +229,7 @@ bool MultiBunchHandler::readBunch(PartBunchBase<double, 3> *beam,
short MultiBunchHandler::injectBunch(PartBunchBase<double, 3> *beam,
const PartData& ref,
bool& flagTransition,
const double& azimuth)
bool& flagTransition)
{
short result = 0;
if ((bunchCount_m == 1) && (mode_m == MB_MODE::AUTO) && (!flagTransition)) {
......@@ -272,7 +262,7 @@ short MultiBunchHandler::injectBunch(PartBunchBase<double, 3> *beam,
// start multi-bunch simulation, fill current phase space to initialR and initialP arrays
if ((radiusThisTurn_m - radiusLastTurn_m) < coeffDBunches_m * XYrms) {
// since next turn, start multi-bunches
saveBunch(beam, azimuth);
saveBunch(beam);
flagTransition = true;
}
......
src/Algorithms/MultiBunchHandler.h
View file @ 8c4e57c7
......@@ -12,14 +12,27 @@
*/
class MultiBunchHandler {
public:
struct injection_t {
injection_t()
: time(0.0)
, pathlength(0.0)
, azimuth(0.0)
, radius(0.0)
{ };
double time; // ns
double pathlength; // m
double azimuth; // deg
double radius; // mm
};
struct beaminfo_t {
beaminfo_t(const double& t,
const double& lpath,
const double& theta)
: time(t)
, azimuth(theta)
beaminfo_t(const injection_t& injection = injection_t())
: time(injection.time)
, azimuth(injection.azimuth)
, radius(injection.radius)
, prevAzimuth(-1.0)
, pathlength(lpath)
, pathlength(injection.pathlength)
{ };
double time;
......@@ -64,8 +77,7 @@ public:
const std::string& mode,
const std::string& binning);
void saveBunch(PartBunchBase<double, 3> *beam,
const double& azimuth);
void saveBunch(PartBunchBase<double, 3> *beam);
bool readBunch(PartBunchBase<double, 3> *beam,
const PartData& ref);
......@@ -77,8 +89,7 @@ public:
*/
short injectBunch(PartBunchBase<double, 3> *beam,
const PartData& ref,
bool& flagTransition,
const double& azimuth);
bool& flagTransition);
void updateParticleBins(PartBunchBase<double, 3> *beam);
......@@ -108,6 +119,8 @@ public:
const beaminfo_t& getBunchInfo(short bunchNr) const;
injection_t& getInjectionValues();
private:
// store the data of the beam which are required for injecting a
// new bunch for multibunch filename
......@@ -141,10 +154,8 @@ private:
// each list entry belongs to a bunch
std::vector<beaminfo_t> binfo_m;
// injection values
double injTime_m;
double injPathlength_m;
double injAzimuth_m;
// global attributes of injection
injection_t injection_m;
};
......@@ -191,4 +202,10 @@ const MultiBunchHandler::beaminfo_t& MultiBunchHandler::getBunchInfo(short bunch
return binfo_m[bunchNr];
}
inline
MultiBunchHandler::injection_t& MultiBunchHandler::getInjectionValues() {
return injection_m;
}
#endif
src/Algorithms/ParallelCyclotronTracker.cpp
View file @ 8c4e57c7
......@@ -238,38 +238,41 @@ double ParallelCyclotronTracker::computeRadius(const Vector_t& meanR) const {
}
double ParallelCyclotronTracker::computePathLengthUpdate(const double& dt,
short bunchNr)
void ParallelCyclotronTracker::computePathLengthUpdate(std::vector<double>& dl,
const double& dt)
{
double dotP = 0.0;
// the last element in dotP is the dot-product over all particles
std::vector<double> dotP(dl.size());
if ( Options::psDumpFrame == Options::BUNCH_MEAN || isMultiBunch()) {
for(unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) {
// take all particles if bunchNr <= -1
if ( bunchNr > -1 && itsBunch_m->bunchNum[i] != bunchNr)
continue;
dotP += dot(itsBunch_m->P[i], itsBunch_m->P[i]);
for(unsigned int i = 0; i < itsBunch_m->getLocalNum(); ++i) {
dotP[itsBunch_m->bunchNum[i]] += dot(itsBunch_m->P[i], itsBunch_m->P[i]);
}
allreduce(dotP, 1, std::plus<double>());
allreduce(dotP.data(), dotP.size(), std::plus<double>());
size_t n = itsBunch_m->getTotalNum();
if ( bunchNr > -1 )
n = itsBunch_m->getTotalNumPerBunch(bunchNr);
// dot-product over all particles
double sum = std::accumulate(dotP.begin(), dotP.end() - 1, 0);
dotP.back() = sum / double(itsBunch_m->getTotalNum());
dotP /= double(n);
// bunch specific --> multi-bunches only
for (short b = 0; b < (short)dotP.size() - 1; ++b) {
dotP[b] /= double(itsBunch_m->getTotalNumPerBunch(b));
}
} else if ( itsBunch_m->getLocalNum() == 0 ) {
// here we are in Options::GLOBAL mode
return 0.0;
dotP[0] = 0.0;
} else {
// here we are in Options::GLOBAL mode
dotP = dot(itsBunch_m->P[0], itsBunch_m->P[0]);
dotP[0] = dot(itsBunch_m->P[0], itsBunch_m->P[0]);
}
double const gamma = std::sqrt(1.0 + dotP);
double const beta = std::sqrt(dotP) / gamma;
return c_mmtns * dt * 1.0e-3 * beta; // unit: m
for (size_t i = 0; i < dotP.size(); ++i) {
double const gamma = std::sqrt(1.0 + dotP[i]);
double const beta = std::sqrt(dotP[i]) / gamma;
dl[i] = c_mmtns * dt * 1.0e-3 * beta; // unit: m
}
}
......@@ -1154,8 +1157,11 @@ void ParallelCyclotronTracker::execute() {
for the integrators
*/
step_m = 0;
restartStep0_m = 0;
step_m = 0;
restartStep0_m = 0;
turnnumber_m = 1;
azimuth_m = -1.0;
prevAzimuth_m = -1.0;
// Record how many bunches have already been injected. ONLY FOR MPM
if (isMultiBunch())
......@@ -2281,7 +2287,7 @@ void ParallelCyclotronTracker::initDistInGlobalFrame() {
// Backup initial distribution if multi bunch mode
if ((initialTotalNum_m > 2) && isMultiBunch() && mbHandler_m->isForceMode()) {
mbHandler_m->saveBunch(itsBunch_m, azimuth_m);
mbHandler_m->saveBunch(itsBunch_m);
}
// Else: Restart from the distribution in the h5 file
......@@ -2377,6 +2383,9 @@ void ParallelCyclotronTracker::initDistInGlobalFrame() {
itsBunch_m->calcBeamParameters();
// multi-bunch simulation only
saveInjectValues();
*gmsg << endl << "* *********************** Bunch information in global frame: ***********************";
*gmsg << *itsBunch_m << endl;
......@@ -2833,8 +2842,6 @@ 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()) {
......@@ -2855,10 +2862,6 @@ std::tuple<double, double, double> ParallelCyclotronTracker::initializeTracking_
if (isMultiBunch())
mbHandler_m->updateParticleBins(itsBunch_m);
turnnumber_m = 1;
azimuth_m = -1.0;
prevAzimuth_m = -1.0;
// --- Output to user --- //
*gmsg << "* Beginning of this run is at t = " << t << " [ns]" << endl;
*gmsg << "* The time step is set to dt = " << dt << " [ns]" << endl;
......@@ -3405,7 +3408,7 @@ void ParallelCyclotronTracker::injectBunch(bool& flagTransition) {
}
const short result = mbHandler_m->injectBunch(itsBunch_m, itsReference,
flagTransition, azimuth_m);
flagTransition);
switch ( result ) {
case 0: {
......@@ -3414,6 +3417,7 @@ void ParallelCyclotronTracker::injectBunch(bool& flagTransition) {
}
case 1: {
// bunch got saved
saveInjectValues();
setup_m.stepsNextCheck += setup_m.stepsPerTurn;
if ( flagTransition ) {
*gmsg << "* MBM: Saving beam distribution at turn " << turnnumber_m << endl;
......@@ -3433,15 +3437,49 @@ void ParallelCyclotronTracker::injectBunch(bool& flagTransition) {
}
void ParallelCyclotronTracker::saveInjectValues() {
if ( !isMultiBunch() ) {
return;
}
Vector_t meanR = calcMeanR();
// Bunch (global) angle w.r.t. x-axis (cylinder coordinates)
double theta = calculateAngle(meanR(0), meanR(1)) * Physics::rad2deg;
// fix azimuth_m --> make monotonically increasing
dumpAngle(theta, prevAzimuth_m, azimuth_m);
const double radius = computeRadius(meanR);
MultiBunchHandler::injection_t& inj = mbHandler_m->getInjectionValues();
inj.time = itsBunch_m->getT();
inj.pathlength = itsBunch_m->getLPath();
inj.azimuth = azimuth_m;
inj.radius = radius;
}
void ParallelCyclotronTracker::updatePathLength(const double& dt) {
pathLength_m += computePathLengthUpdate(dt);
/* the last element includes all particles,
* all other elements are bunch number specific
*/
std::vector<double> lpaths(1);
if ( isMultiBunch() ) {
lpaths.resize(mbHandler_m->getNumBunch() + 1);
}
computePathLengthUpdate(lpaths, dt);
pathLength_m += lpaths[0];
itsBunch_m->setLPath(pathLength_m);
if ( isMultiBunch() ) {
for (short b = 0; b < mbHandler_m->getNumBunch(); ++b) {
MultiBunchHandler::beaminfo_t& binfo = mbHandler_m->getBunchInfo(b);
binfo.pathlength += computePathLengthUpdate(dt, b);
binfo.pathlength += lpaths[b];
}
}
}
......
src/Algorithms/ParallelCyclotronTracker.h
View file @ 8c4e57c7
......@@ -286,8 +286,7 @@ private:
double computeRadius(const Vector_t& meanR) const;
// take all particles if bunchNr <= -1
double computePathLengthUpdate(const double& dt, short bunchNr = -1);
void computePathLengthUpdate(std::vector<double>& dl, const double& dt);
// external field arrays for dumping
Vector_t FDext_m[2], extE_m, extB_m;
......@@ -501,6 +500,8 @@ private:
void injectBunch(bool& flagTransition);
void saveInjectValues();
bool isMultiBunch() const;
bool hasMultiBunch() const;
......
src/Structure/MultiBunchDump.cpp
View file @ 8c4e57c7
......@@ -25,6 +25,7 @@ void MultiBunchDump::fillHeader() {
columns_m.addColumn("t", "double", "ns", "Time");
columns_m.addColumn("s", "double", "m", "Path length");
columns_m.addColumn("azimuth", "double", "deg", "Azimuth in global coordinates");
columns_m.addColumn("radius", "double", "mm", "Radius in global coordinates");
columns_m.addColumn("numParticles", "long", "1", "Number of Macro Particles");
columns_m.addColumn("energy", "double", "MeV", "Mean Bunch Energy");
columns_m.addColumn("dE", "double", "MeV", "energy spread of the beam");
......@@ -88,6 +89,7 @@ void MultiBunchDump::write(PartBunchBase<double, 3>* beam,
columns_m.addColumnValue("t", binfo.time);
columns_m.addColumnValue("s", binfo.pathlength);
columns_m.addColumnValue("azimuth", binfo.azimuth);
columns_m.addColumnValue("radius", binfo.radius);
columns_m.addColumnValue("numParticles", binfo.nParticles);
columns_m.addColumnValue("energy", binfo.ekin);
columns_m.addColumnValue("dE", binfo.dEkin);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment