Commit e574dc13 authored by frey_m's avatar frey_m
Browse files

Merge branch '486-replace-mpi_allreduce-with-allreduce' into 'master'

Resolve "Replace MPI_Allreduce with allreduce"

Closes #486

See merge request !322
parents d31834f0 a42c2c76
...@@ -1141,7 +1141,7 @@ protected: ...@@ -1141,7 +1141,7 @@ protected:
} }
//reduce message count so every node knows how many messages to receive //reduce message count so every node knows how many messages to receive
MPI_Allreduce(&(msgsend[0]), &(msgrecv[0]), N, MPI_INT, MPI_SUM, Ippl::getComm()); allreduce(msgsend.data(), msgrecv.data(), N, std::plus<int>());
int tag = Ippl::Comm->next_tag(P_SPATIAL_TRANSFER_TAG,P_LAYOUT_CYCLE); int tag = Ippl::Comm->next_tag(P_SPATIAL_TRANSFER_TAG,P_LAYOUT_CYCLE);
...@@ -1273,7 +1273,7 @@ protected: ...@@ -1273,7 +1273,7 @@ protected:
} }
//reduce message count so every node knows how many messages to receive //reduce message count so every node knows how many messages to receive
MPI_Allreduce(&(msgsend[0]), &(msgrecv[0]), N, MPI_INT, MPI_SUM, Ippl::getComm()); allreduce(msgsend.data(), msgrecv.data(), N, std::plus<int>());
int tag = Ippl::Comm->next_tag(P_SPATIAL_TRANSFER_TAG,P_LAYOUT_CYCLE); int tag = Ippl::Comm->next_tag(P_SPATIAL_TRANSFER_TAG,P_LAYOUT_CYCLE);
......
...@@ -288,7 +288,7 @@ void ParallelSliceTracker::execute() { ...@@ -288,7 +288,7 @@ void ParallelSliceTracker::execute() {
//reduce(&globalEOL_m, &globalEOL_m, OpBitwiseOrAssign()); //reduce(&globalEOL_m, &globalEOL_m, OpBitwiseOrAssign());
//reduce(&globalEOL_m, &globalEOL_m + 1, &globalEOL_m, OpBitwiseAndAssign()); //reduce(&globalEOL_m, &globalEOL_m + 1, &globalEOL_m, OpBitwiseAndAssign());
MPI_Allreduce(MPI_IN_PLACE, &globalEOL_m, 1, MPI_CXX_BOOL, MPI_LAND, Ippl::getComm()); allreduce(&globalEOL_m, 1, std::logical_and<bool>());
computeSpaceChargeFields(); computeSpaceChargeFields();
timeIntegration(); timeIntegration();
......
...@@ -48,7 +48,6 @@ extern Inform *gmsg; ...@@ -48,7 +48,6 @@ extern Inform *gmsg;
// Hack allows odeint in rk.C to be called with a class member function // Hack allows odeint in rk.C to be called with a class member function
static EnvelopeBunch *thisBunch = NULL; // pointer to access calling bunch static EnvelopeBunch *thisBunch = NULL; // pointer to access calling bunch
static void Gderivs(double t, double Y[], double dYdt[]) { thisBunch->derivs(t, Y, dYdt); } static void Gderivs(double t, double Y[], double dYdt[]) { thisBunch->derivs(t, Y, dYdt); }
//static double Gcur(double z) { return thisBunch->currentProfile_m->get(z, itype_lin); }
static double rootValue = 0.0; static double rootValue = 0.0;
// used in setLShape for Gaussian // used in setLShape for Gaussian
...@@ -123,10 +122,6 @@ void EnvelopeBunch::calcBeamParameters() { ...@@ -123,10 +122,6 @@ void EnvelopeBunch::calcBeamParameters() {
Irms_m = Q_m * nValid_m * Physics::c / (zRms * sqrt(Physics::two_pi) * numSlices_m); Irms_m = Q_m * nValid_m * Physics::c / (zRms * sqrt(Physics::two_pi) * numSlices_m);
Px_m = Px / Physics::c; Px_m = Px / Physics::c;
Py_m = Py / Physics::c; Py_m = Py / Physics::c;
// dx02_m = dx0_m;
// dy02_m = dy0_m;
// dfi_x_m = 180.0 * dfi_x / Physics::pi;
// dfi_y_m = 180.0 * dfi_y / Physics::pi;
Ez_m = 1.0e-6 * Efz; Ez_m = 1.0e-6 * Efz;
Bz_m = Bfz; Bz_m = Bfz;
...@@ -140,7 +135,6 @@ void EnvelopeBunch::calcBeamParameters() { ...@@ -140,7 +135,6 @@ void EnvelopeBunch::calcBeamParameters() {
//minX in the T-T sense is -RxMax_m //minX in the T-T sense is -RxMax_m
minX_m = Vector_t(-RxMax_m, -RyMax_m, zMin); minX_m = Vector_t(-RxMax_m, -RyMax_m, zMin);
minP_m = Vector_t(-PxMax * E_m * sqrt((g0 + 1.0) / (g0 - 1.0)) / Physics::c * Physics::pi, -PyMax * E_m * sqrt((g0 + 1.0) / (g0 - 1.0)) / Physics::c * Physics::pi, PzMin); minP_m = Vector_t(-PxMax * E_m * sqrt((g0 + 1.0) / (g0 - 1.0)) / Physics::c * Physics::pi, -PyMax * E_m * sqrt((g0 + 1.0) / (g0 - 1.0)) / Physics::c * Physics::pi, PzMin);
//minP_m = Vector_t(-maxP_m[0], -maxP_m[1], PzMin);
/* PxRms is the rms of the divergence in x direction in [rad]: x' /* PxRms is the rms of the divergence in x direction in [rad]: x'
* x' = Px/P0 -> Px = x' * P0 = x' * \beta/c*E (E is the particle total * x' = Px/P0 -> Px = x' * P0 = x' * \beta/c*E (E is the particle total
...@@ -234,7 +228,7 @@ void EnvelopeBunch::runStats(EnvelopeBunchParameter sp, double *xAvg, double *xM ...@@ -234,7 +228,7 @@ void EnvelopeBunch::runStats(EnvelopeBunchParameter sp, double *xAvg, double *xM
} }
int nVTot = nV; int nVTot = nV;
MPI_Allreduce(MPI_IN_PLACE, &nVTot, 1, MPI_INT, MPI_SUM, Ippl::getComm()); allreduce(&nVTot, 1, std::plus<int>());
if(nVTot <= 0) { if(nVTot <= 0) {
*xAvg = 0.0; *xAvg = 0.0;
*xMax = 0.0; *xMax = 0.0;
...@@ -259,10 +253,10 @@ void EnvelopeBunch::runStats(EnvelopeBunchParameter sp, double *xAvg, double *xM ...@@ -259,10 +253,10 @@ void EnvelopeBunch::runStats(EnvelopeBunchParameter sp, double *xAvg, double *xM
} }
} }
MPI_Allreduce(MPI_IN_PLACE, &M1, 1, MPI_DOUBLE, MPI_SUM, Ippl::getComm()); allreduce(&M1, 1, std::plus<double>());
MPI_Allreduce(MPI_IN_PLACE, &M2, 1, MPI_DOUBLE, MPI_SUM, Ippl::getComm()); allreduce(&M2, 1, std::plus<double>());
MPI_Allreduce(MPI_IN_PLACE, &maxv, 1, MPI_DOUBLE, MPI_MAX, Ippl::getComm()); allreduce(&maxv, 1, std::greater<double>());
MPI_Allreduce(MPI_IN_PLACE, &minv, 1, MPI_DOUBLE, MPI_MIN, Ippl::getComm()); allreduce(&minv, 1, std::less<double>());
*xAvg = M1 / nVTot; *xAvg = M1 / nVTot;
*xMax = maxv; *xMax = maxv;
...@@ -576,7 +570,6 @@ void EnvelopeBunch::setBinnedLShape(EnvelopeBunchShape shape, double z0, double ...@@ -576,7 +570,6 @@ void EnvelopeBunch::setBinnedLShape(EnvelopeBunchShape shape, double z0, double
case bsRect: case bsRect:
bunch_width = Physics::c * emission_time_s * slices_m[0]->p[SLI_beta]; bunch_width = Physics::c * emission_time_s * slices_m[0]->p[SLI_beta];
// std::cout << "bunch_width = " << bunch_width << " SLI_beta= " << slices_m[0]->p[SLI_beta] << std::endl;
for(int i = 0; i < numMySlices_m; i++) { for(int i = 0; i < numMySlices_m; i++) {
slices_m[i]->p[SLI_z] = -(((numSlices_m - 1) - (mySliceStartOffset_m + i)) * bunch_width) / numSlices_m; slices_m[i]->p[SLI_z] = -(((numSlices_m - 1) - (mySliceStartOffset_m + i)) * bunch_width) / numSlices_m;
} }
...@@ -629,17 +622,6 @@ void EnvelopeBunch::setBinnedLShape(EnvelopeBunchShape shape, double z0, double ...@@ -629,17 +622,6 @@ void EnvelopeBunch::setBinnedLShape(EnvelopeBunchShape shape, double z0, double
} }
} }
////XXX: debug output
/*
for(unsigned int i = 0; i < bins_m.size(); i++) {
if(bins_m[i].size() > 0) {
std::cout << Ippl::Comm->myNode() << ": Bin " << i << ": ";
for(unsigned int j = 0; j < bins_m[i].size(); j++)
std::cout << " " << mySliceStartOffset_m + bins_m[i][j] << "(" << slices_m[j]->p[SLI_z] << ")";
std::cout << std::endl;
}
}
*/
//find first bin containing slices //find first bin containing slices
firstBinWithValue_m = 0; firstBinWithValue_m = 0;
for(; firstBinWithValue_m < nebin_m; firstBinWithValue_m++) for(; firstBinWithValue_m < nebin_m; firstBinWithValue_m++)
...@@ -649,19 +631,12 @@ void EnvelopeBunch::setBinnedLShape(EnvelopeBunchShape shape, double z0, double ...@@ -649,19 +631,12 @@ void EnvelopeBunch::setBinnedLShape(EnvelopeBunchShape shape, double z0, double
} }
void EnvelopeBunch::setTShape(double enx, double eny, double rx, double ry, double /*b0*/) { void EnvelopeBunch::setTShape(double enx, double eny, double rx, double ry, double /*b0*/) {
/*
Inform msg("setTshape");
msg << "set SLI_x to " << rx / 2.0 << endl;
msg << "set SLI_y to " << ry / 2.0 << endl;
*/
// set emittances // set emittances
emtnx0_m = enx; emtnx0_m = enx;
emtny0_m = eny; emtny0_m = eny;
//FIXME: is rx here correct? //FIXME: is rx here correct?
emtbx0_m = Physics::q_e * rx * rx * Bz0_m / (8.0 * Physics::EMASS * Physics::c); emtbx0_m = Physics::q_e * rx * rx * Bz0_m / (8.0 * Physics::EMASS * Physics::c);
emtby0_m = Physics::q_e * ry * ry * Bz0_m / (8.0 * Physics::EMASS * Physics::c); emtby0_m = Physics::q_e * ry * ry * Bz0_m / (8.0 * Physics::EMASS * Physics::c);
//msg << "set emtbx0 to " << emtbx0 << endl;
//msg << "set emtby0 to " << emtby0 << endl;
//XXX: rx = 2*rms_x //XXX: rx = 2*rms_x
for(int i = 0; i < numMySlices_m; i++) { for(int i = 0; i < numMySlices_m; i++) {
...@@ -707,8 +682,8 @@ void EnvelopeBunch::synchronizeSlices() { ...@@ -707,8 +682,8 @@ void EnvelopeBunch::synchronizeSlices() {
z_m[mySliceStartOffset_m+i] = slices_m[i]->p[SLI_z]; z_m[mySliceStartOffset_m+i] = slices_m[i]->p[SLI_z];
} }
MPI_Allreduce(MPI_IN_PLACE, &(z_m[0]), numSlices_m, MPI_DOUBLE, MPI_SUM, Ippl::getComm()); allreduce(z_m.data(), numSlices_m, std::plus<double>());
MPI_Allreduce(MPI_IN_PLACE, &(b_m[0]), numSlices_m, MPI_DOUBLE, MPI_SUM, Ippl::getComm()); allreduce(b_m.data(), numSlices_m, std::plus<double>());
} }
void EnvelopeBunch::calcI() { void EnvelopeBunch::calcI() {
...@@ -746,7 +721,6 @@ void EnvelopeBunch::calcI() { ...@@ -746,7 +721,6 @@ void EnvelopeBunch::calcI() {
if(n1 < 2) { if(n1 < 2) {
msg << "n1 (= " << n1 << ") < 2" << endl; msg << "n1 (= " << n1 << ") < 2" << endl;
//throw OpalException("EnvelopeBunch", "Insufficient points to calculate the current (n1)");
currentProfile_m = std::unique_ptr<Profile>(new Profile(0.0)); currentProfile_m = std::unique_ptr<Profile>(new Profile(0.0));
return; return;
} }
...@@ -805,7 +779,7 @@ void EnvelopeBunch::calcI() { ...@@ -805,7 +779,7 @@ void EnvelopeBunch::calcI() {
} }
} }
MPI_Allreduce(MPI_IN_PLACE, &(I1[0]), n1, MPI_DOUBLE, MPI_SUM, Ippl::getComm()); allreduce(I1.data(), n1, std::plus<double>());
for(int i = 1; i < n1 - 1; i++) { for(int i = 1; i < n1 - 1; i++) {
if(I1[i] == 0.0) if(I1[i] == 0.0)
I1[i] = I1[i-1]; I1[i] = I1[i-1];
...@@ -839,7 +813,6 @@ void EnvelopeBunch::calcI() { ...@@ -839,7 +813,6 @@ void EnvelopeBunch::calcI() {
// following values // following values
int k = 0; int k = 0;
for(int i = 1; i < n1; i++) { for(int i = 1; i < n1; i++) {
//for(int i = my_start; i <= my_end; i++) {
// add new points // add new points
int j = 0; int j = 0;
while(((i + j) < n1) && ((z1[i+j] - z1[i]) <= dz)) { while(((i + j) < n1) && ((z1[i+j] - z1[i]) <= dz)) {
...@@ -861,8 +834,6 @@ void EnvelopeBunch::calcI() { ...@@ -861,8 +834,6 @@ void EnvelopeBunch::calcI() {
} }
++j; ++j;
} }
//z2_temp[i] = Mz1 / np;
//I2_temp[i] = MI1 / np;
z2[i-k] = Mz1 / np; z2[i-k] = Mz1 / np;
I2[i-k] = MI1 / np; I2[i-k] = MI1 / np;
...@@ -873,24 +844,6 @@ void EnvelopeBunch::calcI() { ...@@ -873,24 +844,6 @@ void EnvelopeBunch::calcI() {
} }
} }
//MPI_Allreduce(MPI_IN_PLACE, &(z2_temp[0]), n1, MPI_DOUBLE, MPI_SUM, Ippl::getComm());
//MPI_Allreduce(MPI_IN_PLACE, &(I2_temp[0]), n1, MPI_DOUBLE, MPI_SUM, Ippl::getComm());
////FIXME: we dont need copy of z2 and I2: z2[i-k] = z2[i];
//int k = 0;
//for(int i = 1; i < n1; i++) {
//z2[i-k] = z2_temp[i];
//I2[i-k] = I2_temp[i];
//// make sure there are no duplicate z coordinates
//if(z2[i-k] <= z2[i-k-1]) {
//I2[i-k-1] = 0.5 * (I2[i-k] + I2[i-k-1]);
//k++;
//}
//}
//free(z2_temp);
//free(I2_temp);
int n2 = n1 - k; int n2 = n1 - k;
if(n2 < 1) { if(n2 < 1) {
msg << "Insufficient points to calculate the current (m = " << n2 << ")" << endl; msg << "Insufficient points to calculate the current (m = " << n2 << ")" << endl;
...@@ -945,8 +898,6 @@ void EnvelopeBunch::cSpaceCharge() { ...@@ -945,8 +898,6 @@ void EnvelopeBunch::cSpaceCharge() {
double H1 = sqrt((1 - zeta / L) * (1 - zeta / L) + A * A) - sqrt((zeta / L) * (zeta / L) + A * A) - std::abs(1 - zeta / L) + std::abs(zeta / L); double H1 = sqrt((1 - zeta / L) * (1 - zeta / L) + A * A) - sqrt((zeta / L) * (zeta / L) + A * A) - std::abs(1 - zeta / L) + std::abs(zeta / L);
double H2 = sqrt((1 - xi / L) * (1 - xi / L) + A * A) - sqrt((xi / L) * (xi / L) + A * A) - std::abs(1 - xi / L) + std::abs(xi / L); double H2 = sqrt((1 - xi / L) * (1 - xi / L) + A * A) - sqrt((xi / L) * (xi / L) + A * A) - std::abs(1 - xi / L) + std::abs(xi / L);
//FIXME: Q_act or Q?
//double Q = activeSlices_m * Q_m / numSlices_m;
Esct_m[i] = (Q_m / 2 / Physics::pi / Physics::epsilon_0 / R / R) * (H1 - icON * H2); Esct_m[i] = (Q_m / 2 / Physics::pi / Physics::epsilon_0 / R / R) * (H1 - icON * H2);
double G1 = (1 - zeta / L) / sqrt((1 - zeta / L) * (1 - zeta / L) + A * A) + (zeta / L) / sqrt((zeta / L) * (zeta / L) + A * A); double G1 = (1 - zeta / L) / sqrt((1 - zeta / L) * (1 - zeta / L) + A * A) + (zeta / L) / sqrt((zeta / L) * (zeta / L) + A * A);
double G2 = (1 - xi / L) / sqrt((1 - xi / L) * (1 - xi / L) + A * A) + (xi / L) / sqrt((xi / L) * (xi / L) + A * A); double G2 = (1 - xi / L) / sqrt((1 - xi / L) * (1 - xi / L) + A * A) + (xi / L) / sqrt((xi / L) * (xi / L) + A * A);
...@@ -986,14 +937,14 @@ void EnvelopeBunch::cSpaceCharge() { ...@@ -986,14 +937,14 @@ void EnvelopeBunch::cSpaceCharge() {
} }
int nVTot = nV; int nVTot = nV;
MPI_Allreduce(MPI_IN_PLACE, &nVTot, 1, MPI_INT, MPI_SUM, Ippl::getComm()); allreduce(&nVTot, 1, std::plus<int>());
if(nVTot < 2) { if(nVTot < 2) {
msg << "Exiting, to few nV slices" << endl; msg << "Exiting, to few nV slices" << endl;
return; return;
} }
MPI_Allreduce(MPI_IN_PLACE, &xi[0], numSlices_m, MPI_DOUBLE, MPI_SUM, Ippl::getComm()); allreduce(xi.data(), numSlices_m, std::plus<double>());
MPI_Allreduce(MPI_IN_PLACE, &sm, 1, MPI_DOUBLE, MPI_SUM, Ippl::getComm()); allreduce(&sm, 1, std::plus<double>());
A0 = sm / nVTot; A0 = sm / nVTot;
double dzMin = 5.0 * Physics::c * Q_m / (Imax * numSlices_m); double dzMin = 5.0 * Physics::c * Q_m / (Imax * numSlices_m);
...@@ -1115,8 +1066,6 @@ void EnvelopeBunch::derivs(double /*tc*/, double Y[], double dYdt[]) { ...@@ -1115,8 +1066,6 @@ void EnvelopeBunch::derivs(double /*tc*/, double Y[], double dYdt[]) {
#ifdef USE_HOMDYN_SC_MODEL #ifdef USE_HOMDYN_SC_MODEL
//FIXME: Q_act or Q?
//double Q = activeSlices_m * Q_m / numSlices_m;
double kpc = 0.5 * Physics::c * Physics::c * (Y[SLI_beta] * Physics::c) * activeSlices_m * Q_m / numSlices_m / (curZHead_m - curZTail_m) / Physics::Ia; double kpc = 0.5 * Physics::c * Physics::c * (Y[SLI_beta] * Physics::c) * activeSlices_m * Q_m / numSlices_m / (curZHead_m - curZTail_m) / Physics::Ia;
/// \f[ \dot{\sigma} = p \f] /// \f[ \dot{\sigma} = p \f]
...@@ -1182,11 +1131,6 @@ void EnvelopeBunch::computeSpaceCharge() { ...@@ -1182,11 +1131,6 @@ void EnvelopeBunch::computeSpaceCharge() {
cSpaceCharge(); cSpaceCharge();
IpplTimings::stopTimer(spaceChargeTimer_m); IpplTimings::stopTimer(spaceChargeTimer_m);
//if(numSlices_m > 40) {
// smooth Esct to get rid of numerical noise
//double Esc = 0;
//sgSmooth(Esct,n,n/20,n/20,0,4);
//}
} else { } else {
currentProfile_m = std::unique_ptr<Profile>(new Profile(0.0)); currentProfile_m = std::unique_ptr<Profile>(new Profile(0.0));
} }
...@@ -1212,35 +1156,6 @@ void EnvelopeBunch::timeStep(double tStep, double _zCat) { ...@@ -1212,35 +1156,6 @@ void EnvelopeBunch::timeStep(double tStep, double _zCat) {
//FIXME: HAVE A PROBLEM WITH EMISSION (EMISSION OFF GIVES GOOD RESULTS) //FIXME: HAVE A PROBLEM WITH EMISSION (EMISSION OFF GIVES GOOD RESULTS)
activeSlices_m = numSlices_m; activeSlices_m = numSlices_m;
//if(lastEmittedBin_m < nebin_m) {
//time_step_s = emission_time_step_;
//size_t nextBin = nebin_m - 1 - lastEmittedBin_m;
//if(Ippl::Comm->myNode() == 0) cout << "Emitting bin " << nextBin << endl;
//double gz0 = 0.0;
//if(Ippl::Comm->myNode() == 0)
//gz0 = slices_m[0]->p[SLI_z];
//MPI_Allreduce(MPI_IN_PLACE, &gz0, 1, MPI_DOUBLE, MPI_SUM, Ippl::getComm());
////XXX: bin holds local slice number
//for(int j = 0; j < bins_m[nextBin].size(); j++) {
//slices_m[bins_m[nextBin][j]]->p[SLI_z] -= gz0;
//slices_m[bins_m[nextBin][j]]->p[SLI_z] -= hbin_m * nextBin;
//slices_m[bins_m[nextBin][j]]->emitted = true;
//cout << "\tEmitting slice " << mySliceStartOffset_m + bins_m[nextBin][j] << " at " << slices_m[bins_m[nextBin][j]]->p[SLI_z] << endl;
//}
//lastEmittedBin_m++;
//activeSlices_m += bins_m[nextBin].size();
//}
//activeSlices_m = min(activeSlices_m+1, numSlices_m);
//if(activeSlices_m != numSlices_m)
//time_step_s = emission_time_step_;
//else
//time_step_s = tStep;
curZHead_m = zHead(); curZHead_m = zHead();
curZTail_m = zTail(); curZTail_m = zTail();
...@@ -1249,12 +1164,6 @@ void EnvelopeBunch::timeStep(double tStep, double _zCat) { ...@@ -1249,12 +1164,6 @@ void EnvelopeBunch::timeStep(double tStep, double _zCat) {
currentSlice_m = i; currentSlice_m = i;
std::shared_ptr<EnvelopeSlice> sp = slices_m[i]; std::shared_ptr<EnvelopeSlice> sp = slices_m[i];
//if(i + mySliceStartOffset_m == numSlices_m - activeSlices_m) {
//sp->emitted = true;
//sp->p[SLI_z] = 0.0;
//std::cout << "emitting " << i + mySliceStartOffset_m << std::endl;
//}
// Assign values of K for certain slices // Assign values of K for certain slices
KRsl_m = KR[i]; KRsl_m = KR[i];
KTsl_m = KT[i]; KTsl_m = KT[i];
...@@ -1413,7 +1322,7 @@ double EnvelopeBunch::AvBField() { ...@@ -1413,7 +1322,7 @@ double EnvelopeBunch::AvBField() {
} }
} }
MPI_Allreduce(MPI_IN_PLACE, &bf, 1, MPI_DOUBLE, MPI_SUM, Ippl::getComm()); allreduce(&bf, 1, std::plus<double>());
return bf / numSlices_m; return bf / numSlices_m;
} }
...@@ -1425,7 +1334,7 @@ double EnvelopeBunch::AvEField() { ...@@ -1425,7 +1334,7 @@ double EnvelopeBunch::AvEField() {
} }
} }
MPI_Allreduce(MPI_IN_PLACE, &ef, 1, MPI_DOUBLE, MPI_SUM, Ippl::getComm()); allreduce(&ef, 1, std::plus<double>());
return ef / numSlices_m; return ef / numSlices_m;
} }
...@@ -1438,8 +1347,8 @@ double EnvelopeBunch::Eavg() { ...@@ -1438,8 +1347,8 @@ double EnvelopeBunch::Eavg() {
nValid++; nValid++;
} }
} }
MPI_Allreduce(MPI_IN_PLACE, &nValid, 1, MPI_INT, MPI_SUM, Ippl::getComm()); allreduce(&nValid, 1, std::plus<int>());
MPI_Allreduce(MPI_IN_PLACE, &sum, 1, MPI_DOUBLE, MPI_SUM, Ippl::getComm()); allreduce(&sum, 1, std::plus<double>());
sum /= nValid; sum /= nValid;
return (nValid > 0 ? ((Physics::EMASS * Physics::c * Physics::c / Physics::q_e) * (sum - 1.0)) : 0.0); return (nValid > 0 ? ((Physics::EMASS * Physics::c * Physics::c / Physics::q_e) * (sum - 1.0)) : 0.0);
} }
...@@ -1455,8 +1364,8 @@ double EnvelopeBunch::get_sPos() { ...@@ -1455,8 +1364,8 @@ double EnvelopeBunch::get_sPos() {
} }
} }
MPI_Allreduce(MPI_IN_PLACE, &count, 1, MPI_UNSIGNED_LONG, MPI_SUM, Ippl::getComm()); allreduce(&count, 1, std::plus<size_t>());
MPI_Allreduce(MPI_IN_PLACE, &refpos, 1, MPI_DOUBLE, MPI_SUM, Ippl::getComm()); allreduce(&refpos, 1, std::plus<double>());
return refpos / count; return refpos / count;
} }
...@@ -1470,14 +1379,13 @@ double EnvelopeBunch::zAvg() { ...@@ -1470,14 +1379,13 @@ double EnvelopeBunch::zAvg() {
} }
} }
MPI_Allreduce(MPI_IN_PLACE, &nV, 1, MPI_INT, MPI_SUM, Ippl::getComm()); allreduce(&nV, 1, std::plus<int>());
if(nV < 1) { if(nV < 1) {
isValid_m = false; isValid_m = false;
return -1; return -1;
//throw OpalException("EnvelopeBunch", "EnvelopeBunch::zAvg() no valid slices left");
} }
MPI_Allreduce(MPI_IN_PLACE, &sum, 1, MPI_DOUBLE, MPI_SUM, Ippl::getComm()); allreduce(&sum, 1, std::plus<double>());
return (sum / nV); return (sum / nV);
} }
...@@ -1487,7 +1395,6 @@ double EnvelopeBunch::zTail() { ...@@ -1487,7 +1395,6 @@ double EnvelopeBunch::zTail() {
int i = 0; int i = 0;
while((i < numMySlices_m) && (!slices_m[i]->is_valid())) while((i < numMySlices_m) && (!slices_m[i]->is_valid()))
i++; i++;
//throw OpalException("EnvelopeBunch", "EnvelopeBunch::zTail() no valid slices left");
if(i == numMySlices_m) if(i == numMySlices_m)
isValid_m = false; isValid_m = false;
else else
...@@ -1497,8 +1404,7 @@ double EnvelopeBunch::zTail() { ...@@ -1497,8 +1404,7 @@ double EnvelopeBunch::zTail() {
if((slices_m[i]->p[SLI_z] < min) && (slices_m[i]->is_valid())) if((slices_m[i]->p[SLI_z] < min) && (slices_m[i]->is_valid()))
min = slices_m[i]->p[SLI_z]; min = slices_m[i]->p[SLI_z];
//reduce(min, min, OpMinAssign()); allreduce(&min, 1, std::less<double>());
MPI_Allreduce(MPI_IN_PLACE, &min, 1, MPI_DOUBLE, MPI_MIN, Ippl::getComm());
return min; return min;
} }
...@@ -1508,7 +1414,6 @@ double EnvelopeBunch::zHead() { ...@@ -1508,7 +1414,6 @@ double EnvelopeBunch::zHead() {
int i = 0; int i = 0;
while((i < numMySlices_m) && (slices_m[i]->is_valid() == 0)) while((i < numMySlices_m) && (slices_m[i]->is_valid() == 0))
i++; i++;
//throw OpalException("EnvelopeBunch", "EnvelopeBunch::zHead() no valid slices left");
if(i == numMySlices_m) if(i == numMySlices_m)
isValid_m = false; isValid_m = false;
else else
...@@ -1517,8 +1422,7 @@ double EnvelopeBunch::zHead() { ...@@ -1517,8 +1422,7 @@ double EnvelopeBunch::zHead() {
for(i = 1; i < numMySlices_m; i++) for(i = 1; i < numMySlices_m; i++)
if(slices_m[i]->p[SLI_z] > max) max = slices_m[i]->p[SLI_z]; if(slices_m[i]->p[SLI_z] > max) max = slices_m[i]->p[SLI_z];
//reduce(max, max, OpMaxAssign()); allreduce(&max, 1, std::greater<double>());
MPI_Allreduce(MPI_IN_PLACE, &max, 1, MPI_DOUBLE, MPI_MAX, Ippl::getComm());
return max; return max;
} }
...@@ -1526,41 +1430,13 @@ double EnvelopeBunch::zHead() { ...@@ -1526,41 +1430,13 @@ double EnvelopeBunch::zHead() {
Inform &EnvelopeBunch::slprint(Inform &os) { Inform &EnvelopeBunch::slprint(Inform &os) {
if(this->getTotalNum() != 0) { // to suppress Nan's if(this->getTotalNum() != 0) { // to suppress Nan's
os << "* ************** S L B U N C H ***************************************************** " << endl; os << "* ************** S L B U N C H ***************************************************** " << endl;
os << "* NSlices= " << this->getTotalNum() << " Qtot= " << Q_m << endl; //" [nC] Qi= " << std::abs(qi_m) << " [C]" << endl; os << "* NSlices= " << this->getTotalNum() << " Qtot= " << Q_m << endl;
os << "* Emean= " << get_meanKineticEnergy() * 1e-6 << " [MeV]" << endl; os << "* Emean= " << get_meanKineticEnergy() * 1e-6 << " [MeV]" << endl;
os << "* dT= " << this->getdT() << " [s]" << endl; os << "* dT= " << this->getdT() << " [s]" << endl;
os << "* spos= " << this->zAvg() << " [m]" << endl; os << "* spos= " << this->zAvg() << " [m]" << endl;
//os << "* rmax= " << rmax_m << " [m]" << endl;
//os << "* rmin= " << rmin_m << " [m]" << endl;
//os << "* rms beam size= " << rrms_m << " [m]" << endl;
//os << "* rms momenta= " << prms_m << " [beta gamma]" << endl;
//os << "* mean position= " << rmean_m << " [m]" << endl;
//os << "* mean momenta= " << pmean_m << " [beta gamma]" << endl;
//os << "* rms emmitance= " << eps_m << " (not normalized)" << endl;
//os << "* rms correlation= " << rprms_m << endl;
//os << "* tEmission= " << getTEmission() << " [s]" << endl;
os << "* ********************************************************************************** " << endl; os << "* ********************************************************************************** " << endl;
//Inform::FmtFlags_t ff = os.flags();
//os << scientific;
//os << "* ************** B U N C H ********************************************************* " << endl;
//os << "* NSlices= " << this->getTotalNum() << " Qtot= " << abs(sum(Q) * 1.0E9) << " [nC]" << endl;
//os << "* Ekin = " << setw(12) << setprecision(5) << eKin_m << " [MeV] dEkin= " << dE_m << " [MeV]" << endl;
//os << "* rmax = " << setw(12) << setprecision(5) << rmax_m << " [m]" << endl;
//os << "* rmin = " << setw(12) << setprecision(5) << rmin_m << " [m]" << endl;
//os << "* rms beam size = " << setw(12) << setprecision(5) << rrms_m << " [m]" << endl;
//os << "* rms momenta = " << setw(12) << setprecision(5) << prms_m << " [beta gamma]" << endl;
//os << "* mean position = " << setw(12) << setprecision(5) << rmean_m << " [m]" << endl;
//os << "* mean momenta = " << setw(12) << setprecision(5) << pmean_m << " [beta gamma]" << endl;
//os << "* rms emittance = " << setw(12) << setprecision(5) << eps_m << " (not normalized)" << endl;
//os << "* rms correlation = " << setw(12) << setprecision(5) << rprms_m << endl;
//os << "* hr = " << setw(12) << setprecision(5) << hr_m << " [m]" << endl;
//os << "* dh = " << setw(12) << setprecision(5) << dh_m << " [m]" << endl;
//os << "* t = " << setw(12) << setprecision(5) << getT() << " [s] dT= " << getdT() << " [s]" << endl;
//os << "* spos = " << setw(12) << setprecision(5) << get_sPos() << " [m]" << endl;
//os << "* ********************************************************************************** " << endl;
//os.flags(ff);