Commit 2b6e0805 authored by frey_m's avatar frey_m
Browse files

Merge branch 'master' into 'master'

Master

See merge request !59
parents 5acc4c14 2f148336
......@@ -264,6 +264,7 @@ public:
void storeArguments(int argc, char *argv[]);
std::vector<std::string> getArguments();
private:
static bool isInstantiated;
......@@ -282,4 +283,4 @@ private:
//extern OpalData OPAL;
#endif // OPAL_OpalData_HH
\ No newline at end of file
#endif // OPAL_OpalData_HH
......@@ -1444,6 +1444,10 @@ void ParallelCyclotronTracker::GenericTracker() {
// Apply the plugin elements: probe, collimator, stripper, septum once before first step
applyPluginElements(dt);
// Destroy particles if they are marked as Bin = -1 in the plugin elements
// or out of global aperture
deleteParticle();
/********************************
* Main integration loop *
********************************/
......@@ -2207,10 +2211,7 @@ void ParallelCyclotronTracker::borisExternalFields(double h) {
// apply the plugin elements: probe, collimator, stripper, septum
applyPluginElements(h);
// destroy particles if they are marked as Bin=-1 in the plugin elements or out of global apeture
bool const flagNeedUpdate = deleteParticle();
if(itsBunch_m->weHaveBins() && flagNeedUpdate)
updateParticleBins_m();
deleteParticle();
}
......@@ -2249,13 +2250,34 @@ bool ParallelCyclotronTracker::deleteParticle(){
if(flagNeedUpdate) {
size_t locLostParticleNum = 0;
const int leb = itsBunch_m->getLastemittedBin();
std::unique_ptr<size_t[]> localBinCount;
if ( itsBunch_m->weHaveBins() ) {
localBinCount = std::unique_ptr<size_t[]>(new size_t[leb]);
for (int i = 0; i < leb; ++i)
localBinCount[i] = 0;
}
for(unsigned int i = 0; i < itsBunch_m->getLocalNum(); i++) {
if(itsBunch_m->Bin[i] < 0) {
++locLostParticleNum;
itsBunch_m->destroy(1, i);
} else if ( itsBunch_m->weHaveBins() ) {
/* we need to count the local number of particles
* per energy bin
*/
++localBinCount[itsBunch_m->Bin[i]];
}
}
if ( itsBunch_m->weHaveBins() ) {
// set the local bin count
for (int i = 0; i < leb; ++i) {
itsBunch_m->setLocalBinCount(localBinCount[i], i);
}
}
/* We need to destroy the particles now
* before we compute the means. We also
* have to update the total number of particles
......@@ -2302,6 +2324,10 @@ bool ParallelCyclotronTracker::deleteParticle(){
localToGlobal(itsBunch_m->R, phi, psi, meanR);
localToGlobal(itsBunch_m->P, phi, psi, Vector_t(0.0));
// If particles were deleted, recalculate bingamma and reset BinID for remaining particles
if ( itsBunch_m->weHaveBins() )
updateParticleBins_m();
}
return flagNeedUpdate;
}
......@@ -2428,7 +2454,7 @@ void ParallelCyclotronTracker::initDistInGlobalFrame() {
if(multiBunchMode_m == MB_MODE::AUTO) {
RLastTurn_m = sqrt(meanR[0] * meanR[0] + meanR[1] * meanR[1]);
RLastTurn_m = std::sqrt(meanR[0] * meanR[0] + meanR[1] * meanR[1]); // [m]
RThisTurn_m = RLastTurn_m;
if(OpalData::getInstance()->inRestartRun()) {
......@@ -2440,7 +2466,7 @@ void ParallelCyclotronTracker::initDistInGlobalFrame() {
*gmsg << "Initial radial position = ";
}
// New OPAL 2.0: Init in m -DW
*gmsg << 1000.0 * RThisTurn_m << " mm" << endl;
*gmsg << RThisTurn_m << " m" << endl;
}
// Do boundp and repartition in the local frame at beginning of this run
......@@ -2624,7 +2650,35 @@ void ParallelCyclotronTracker::bunchDumpStatData(){
// dump stat file per bin in case of multi-bunch mode
if ( multiBunchMode_m != MB_MODE::NONE ) {
bunchDumpStatDataPerBin();
double phi = 0.0, psi = 0.0;
Vector_t meanR = calcMeanR();
// Bunch (global) angle w.r.t. x-axis (cylinder coordinates)
double theta = calculateAngle(meanR(0), meanR(1)) * Physics::rad2deg;
if(Options::psDumpFrame != Options::GLOBAL) {
Vector_t meanP = calcMeanP();
// Bunch (local) azimuth at meanR w.r.t. y-axis
phi = calculateAngle(meanP(0), meanP(1)) - 0.5 * pi;
// Bunch (local) elevation at meanR w.r.t. xy plane
psi = 0.5 * pi - acos(meanP(2) / sqrt(dot(meanP, meanP)));
// Rotate so Pmean is in positive y direction. No shift, so that normalized emittance and
// unnormalized emittance as well as centroids are calculated correctly in
// PartBunch::calcBeamParameters()
globalToLocal(itsBunch_m->R, phi, psi, meanR);
globalToLocal(itsBunch_m->P, phi, psi);
}
bunchDumpStatDataPerBin(theta);
if(Options::psDumpFrame != Options::GLOBAL) {
localToGlobal(itsBunch_m->R, phi, psi, meanR);
localToGlobal(itsBunch_m->P, phi, psi);
}
return;
}
......@@ -2655,6 +2709,10 @@ void ParallelCyclotronTracker::bunchDumpStatData(){
}
double phi = 0;
double psi = 0;
// Bunch (global) angle w.r.t. x-axis (cylinder coordinates)
double azimuth = calculateAngle(meanR(0), meanR(1)) * Physics::rad2deg;
// -------------- Calculate the external fields at the center of the bunch ----------------- //
beamline_list::iterator DumpSindex = FieldDimensions.begin();
......@@ -2688,7 +2746,7 @@ void ParallelCyclotronTracker::bunchDumpStatData(){
FDext_m[1] = extE_m; // kV/mm? -DW
// Save the stat file
itsDataSink->writeStatData(itsBunch_m, FDext_m, E);
itsDataSink->writeStatData(itsBunch_m, FDext_m, E, azimuth);
//itsBunch_m->R *= Vector_t(1000.0); // m -> mm
......@@ -2702,7 +2760,7 @@ void ParallelCyclotronTracker::bunchDumpStatData(){
}
void ParallelCyclotronTracker::bunchDumpStatDataPerBin() {
void ParallelCyclotronTracker::bunchDumpStatDataPerBin(const double& azimuth) {
IpplTimings::startTimer(DumpTimer_m);
for (short b = 0; b < BunchCount_m; ++b) {
......@@ -2710,6 +2768,7 @@ void ParallelCyclotronTracker::bunchDumpStatDataPerBin() {
MultiBunchDump::beaminfo_t binfo;
if ( itsBunch_m->calcBunchBeamParameters(binfo, b) ) {
binfo.azimuth = azimuth;
itsMBDump_m->writeData(binfo, b);
}
}
......@@ -2911,7 +2970,6 @@ std::tuple<double, double, double> ParallelCyclotronTracker::initializeTracking_
double oldReferenceTheta = referenceTheta * Physics::deg2rad; // init here, reset each step
setup_m.deltaTheta = pi / (setup_m.stepsPerTurn); // half of the average angle per step
setup_m.stepsNextCheck = step_m + setup_m.stepsPerTurn; // Steps to next check for transition
//int stepToLastInj = itsBunch_m->getSteptoLastInj(); // TODO: Do we need this? -DW
// Record how many bunches have already been injected. ONLY FOR MBM
......@@ -2929,6 +2987,8 @@ std::tuple<double, double, double> ParallelCyclotronTracker::initializeTracking_
*gmsg << "* Restart at integration step " << restartStep0_m << endl;
}
setup_m.stepsNextCheck = step_m + setup_m.stepsPerTurn; // Steps to next check for transition
initDistInGlobalFrame();
if ( multiBunchMode_m != MB_MODE::NONE ) {
......@@ -3183,7 +3243,7 @@ void ParallelCyclotronTracker::singleMode_m(double& t, const double dt,
void ParallelCyclotronTracker::bunchMode_m(double& t, const double dt, bool& dumpEachTurn) {
// Flag for transition from single-bunch to multi-bunches mode
bool flagTransition = false;
static bool flagTransition = false;
// single particle dumping
if(step_m % Options::sptDumpFreq == 0)
......@@ -3258,11 +3318,7 @@ void ParallelCyclotronTracker::bunchMode_m(double& t, const double dt, bool& dum
// Destroy particles if they are marked as Bin = -1 in the plugin elements
// or out of global aperture
bool flagNeedUpdate = deleteParticle();
// If particles were deleted, recalculate bingamma and reset BinID for remaining particles
if(itsBunch_m->weHaveBins() && flagNeedUpdate)
updateParticleBins_m();
deleteParticle();
// Recalculate bingamma and reset the BinID for each particles according to its current gamma
if (itsBunch_m->weHaveBins() && BunchCount_m > 1 && step_m % Options::rebinFreq == 0)
......@@ -3439,8 +3495,36 @@ void ParallelCyclotronTracker::computeSpaceChargeFields_m() {
getQuaternionTwoVectors(PreviousMeanP, yaxis, quaternionToYAxis);
#if COORD_PRINT
if ( BunchCount_m > 1 ) {
std::stringstream ss;
ss << "global-coordinates-step-" << step_m;
std::ofstream out(ss.str());
for (uint i = 0; i < itsBunch_m->getLocalNum(); ++i) {
out << itsBunch_m->R[i](0) << " "
<< itsBunch_m->R[i](1) << " "
<< itsBunch_m->R[i](2) << std::endl;
}
out.close();
}
#endif
globalToLocal(itsBunch_m->R, quaternionToYAxis, meanR);
#if COORD_PRINT
if ( BunchCount_m > 1 ) {
std::stringstream ss;
ss << "local-coordinates-step-" << step_m;
std::ofstream out(ss.str());
for (uint i = 0; i < itsBunch_m->getLocalNum(); ++i) {
out << itsBunch_m->R[i](0) << " "
<< itsBunch_m->R[i](1) << " "
<< itsBunch_m->R[i](2) << std::endl;
}
out.close();
}
#endif
//itsBunch_m->R *= Vector_t(0.001); // mm --> m
if ((step_m + 1) % Options::boundpDestroyFreq == 0)
......@@ -3532,11 +3616,11 @@ void ParallelCyclotronTracker::injectBunch_m(bool& flagTransition) {
itsBunch_m->calcBeamParameters();
//itsBunch_m->R *= Vector_t(1000.0); // m --> mm
Vector_t Rmean = itsBunch_m->get_centroid() * 1000.0; // m --> mm
Vector_t Rmean = itsBunch_m->get_centroid(); // m
RThisTurn_m = sqrt(pow(Rmean[0], 2.0) + pow(Rmean[1], 2.0));
Vector_t Rrms = itsBunch_m->get_rrms() * 1000.0; // m --> mm
Vector_t Rrms = itsBunch_m->get_rrms(); // m
double XYrms = sqrt(pow(Rrms[0], 2.0) + pow(Rrms[1], 2.0));
......@@ -3552,9 +3636,9 @@ void ParallelCyclotronTracker::injectBunch_m(bool& flagTransition) {
setup_m.stepsNextCheck += setup_m.stepsPerTurn;
*gmsg << "* MBM: RLastTurn = " << RLastTurn_m << " [mm]" << endl;
*gmsg << "* MBM: RThisTurn = " << RThisTurn_m << " [mm]" << endl;
*gmsg << "* MBM: XYrms = " << XYrms << " [mm]" << endl;
*gmsg << "* MBM: RLastTurn = " << RLastTurn_m << " [m]" << endl;
*gmsg << "* MBM: RThisTurn = " << RThisTurn_m << " [m]" << endl;
*gmsg << "* MBM: XYrms = " << XYrms << " [m]" << endl;
RLastTurn_m = RThisTurn_m;
}
......
......@@ -464,7 +464,8 @@ private:
void bunchDumpStatData();
void bunchDumpStatDataPerBin();
// @param azimuth (global) [deg] of dump
void bunchDumpStatDataPerBin(const double& azimuth);
void bunchDumpPhaseSpaceData();
......
......@@ -20,14 +20,16 @@ AmrBoxLib::AmrBoxLib(const AmrDomain_t& domain,
const AmrIntArray_t& nGridPts,
int maxLevel,
AmrPartBunch* bunch_p)
: AmrObject(),
amrex::AmrMesh(&domain, maxLevel, nGridPts, 0 /* cartesian */),
bunch_mp(bunch_p),
layout_mp(static_cast<AmrLayout_t*>(&bunch_p->getLayout())),
rho_m(maxLevel + 1),
phi_m(maxLevel + 1),
efield_m(maxLevel + 1),
meshScaling_m(Vector_t(1.0, 1.0, 1.0))
: AmrObject()
, amrex::AmrMesh(&domain, maxLevel, nGridPts, 0 /* cartesian */)
, bunch_mp(bunch_p)
, layout_mp(static_cast<AmrLayout_t*>(&bunch_p->getLayout()))
, rho_m(maxLevel + 1)
, phi_m(maxLevel + 1)
, efield_m(maxLevel + 1)
, meshScaling_m(Vector_t(1.0, 1.0, 1.0))
, isFirstTagging_m(maxLevel + 1, true)
, isPoissonSolved_m(false)
{
/*
* The layout needs to know how many levels we can make.
......@@ -533,6 +535,8 @@ void AmrBoxLib::computeSelfFields_cycl(int bin) {
ytWriter.writeBunch(bunch_mp, time, scalefactor);
}
isPoissonSolved_m = true;
}
......@@ -825,15 +829,15 @@ void AmrBoxLib::tagForChargeDensity_m(int lev, TagBoxArray_t& tags,
void AmrBoxLib::tagForPotentialStrength_m(int lev, TagBoxArray_t& tags,
AmrReal_t time, int ngrow)
{
/* Perform a single level solve a level lev and tag all cells for refinement
/* Tag all cells for refinement
* where the value of the potential is higher than 75 percent of the maximum potential
* value of this level.
*/
if ( !phi_m[lev]->ok() || (time == 0 && !(phi_m[lev]->norm0(0) > 0)) ) {
if ( !isPoissonSolved_m || !phi_m[lev]->ok() || this->isFirstTagging_m[lev] ) {
*gmsg << level2 << "* Level " << lev << ": We need to perform "
<< "charge tagging in the first time step" << endl;
this->tagForChargeDensity_m(lev, tags, time, ngrow);
this->isFirstTagging_m[lev] = false;
return;
}
......@@ -878,16 +882,15 @@ void AmrBoxLib::tagForPotentialStrength_m(int lev, TagBoxArray_t& tags,
void AmrBoxLib::tagForEfield_m(int lev, TagBoxArray_t& tags,
AmrReal_t time, int ngrow)
{
/* Perform a single level solve a level lev and tag all cells for refinement
/* Tag all cells for refinement
* where the value of the efield is higher than 75 percent of the maximum efield
* value of this level.
*/
if ( !efield_m[lev][0]->ok() || (time == 0 && !(efield_m[lev][0]->norm0(0) > 0)) ) {
if ( !isPoissonSolved_m || !efield_m[lev][0]->ok() || this->isFirstTagging_m[lev] ) {
*gmsg << level2 << "* Level " << lev << ": We need to perform "
<< "charge tagging in the first time step" << endl;
this->tagForChargeDensity_m(lev, tags, time, ngrow);
this->isFirstTagging_m[lev] = false;
return;
}
......
......@@ -365,6 +365,11 @@ private:
/// in particle rest frame, the longitudinal length enlarged
Vector_t meshScaling_m;
// only used in case of potential and efield tagging
std::vector<bool> isFirstTagging_m;
bool isPoissonSolved_m;
};
#endif
......@@ -226,7 +226,7 @@ public:
}
/** \brief How many particles are on one energy bin */
inline size_t getLocalBinCount(int bin) {
virtual inline size_t getLocalBinCount(int bin) {
size_t ret = 0;
for(int i = sBins_m * bin; i < sBins_m * (bin + 1); i++) {
ret += gsl_histogram_get(h_m.get(), i);
......
......@@ -28,6 +28,8 @@
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_randist.h>
#include <cassert>
class PartBinsCyc: public PartBins {
public:
......@@ -47,6 +49,12 @@ public:
return a;
}
/** \brief How many particles are on one energy bin */
inline size_t getLocalBinCount(int bin) {
assert(bin < bins_m);
return nBin_m[bin];
}
bool weHaveBins() {
return ( nemittedBins_m > 0 );
}
......
......@@ -120,6 +120,8 @@ public:
int getLastemittedBin();
void setLocalBinCount(size_t num, int bin);
/** \brief Compute the gammas of all bins */
void calcGammas();
......
......@@ -356,6 +356,14 @@ int PartBunchBase<T, Dim>::getLastemittedBin() {
}
template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::setLocalBinCount(size_t num, int bin) {
if(pbin_m != NULL) {
pbin_m->setPartNum(bin, num);
}
}
template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::calcGammas() {
......@@ -410,8 +418,11 @@ void PartBunchBase<T, Dim>::calcGammas_cycl() {
for(int i = 0; i < emittedBins; i++)
bingamma_m[i] = 0.0;
for(unsigned int n = 0; n < getLocalNum(); n++)
bingamma_m[this->Bin[n]] += sqrt(1.0 + dot(this->P[n], this->P[n]));
for(unsigned int n = 0; n < getLocalNum(); n++) {
if ( this->Bin[n] > -1 ) {
bingamma_m[this->Bin[n]] += sqrt(1.0 + dot(this->P[n], this->P[n]));
}
}
allreduce(*bingamma_m.get(), emittedBins, std::plus<double>());
......@@ -1794,7 +1805,13 @@ bool PartBunchBase<T, Dim>::resetPartBinID2(const double eta) {
for(unsigned long int n = 0; n < getLocalNum(); n++) {
double temp_betagamma = sqrt(pow(P[n](0), 2) + pow(P[n](1), 2));
#ifdef ENABLE_AMR
// FIXME: Restrict to 1 GeV proton
if ( temp_betagamma > 1.808 ) {
Bin[n] = -1;
continue;
}
#endif
int itsBinID = floor((asinh(temp_betagamma) - asinh0) / eta + 1.0E-6);
Bin[n] = itsBinID;
if(maxbinIndex < itsBinID) {
......@@ -2435,4 +2452,4 @@ FMatrix<double, 2 * Dim, 2 * Dim> PartBunchBase<T, Dim>::getSigmaMatrix() {
return sigmaMatrix;
}
#endif
\ No newline at end of file
#endif
......@@ -88,6 +88,21 @@ namespace {
}
}
bool checkInitAmrFlag(int argc, char* argv[]) {
std::string noamr = "noInitAMR";
bool initAMR = true;
for (int i = 0; i < argc; ++i) {
std::string sargv = std::string(argv[i]);
if ( sargv.find(noamr) != std::string::npos ) {
initAMR = false;
break;
}
}
return initAMR;
}
int main(int argc, char *argv[]) {
Ippl *ippl = new Ippl(argc, argv);
gmsg = new Inform("OPAL");
......@@ -95,8 +110,11 @@ int main(int argc, char *argv[]) {
namespace fs = boost::filesystem;
#ifdef ENABLE_AMR
// false: build no parmparse, we use the OPAL parser instead.
amrex::Initialize(argc, argv, false, Ippl::getComm());
bool initAMR = checkInitAmrFlag(argc, argv);
if ( initAMR ) {
// false: build no parmparse, we use the OPAL parser instead.
amrex::Initialize(argc, argv, false, Ippl::getComm());
}
#endif
OPALTimer::Timer simtimer;
......@@ -269,6 +287,8 @@ int main(int argc, char *argv[]) {
}
INFOMSG(header << options << endl);
exit(0);
} else if ( argStr.find("noInitAMR") != std::string::npos) {
// do nothing here
} else if (argStr == std::string("-help") ||
argStr == std::string("--help")) {
IpplInfo::printHelp(argv);
......@@ -379,10 +399,13 @@ int main(int argc, char *argv[]) {
Fieldmap::clearDictionary();
OpalData::deleteInstance();
delete gmsg;
#ifdef ENABLE_AMR
amrex::Finalize(true);
if ( initAMR ) {
amrex::Finalize(true);
}
#endif
delete ippl;
delete Ippl::Info;
delete Ippl::Warn;
......
......@@ -22,6 +22,7 @@ set (HDRS
SampleIndividual.h
SamplePilot.h
Sampler.h
SampleRandomizedSequence.h
SampleSequence.h
SampleWorker.h
SamplingMethod.h
......
......@@ -11,6 +11,7 @@
#include "Sample/SampleGaussianSequence.h"
#include "Sample/FromFile.h"
#include "Sample/LatinHyperCube.h"
#include "Sample/SampleRandomizedSequence.h"
// Class OpalSample
......@@ -25,6 +26,7 @@ namespace {
FNAME, // file to read from sampling points
N,
RANDOM,
STEP,
SIZE
};
}
......@@ -52,6 +54,9 @@ OpalSample::OpalSample():
itsAttr[RANDOM] = Attributes::makeBool
("RANDOM", "Whether sequence should be sampled randomly (default: false)", false);
itsAttr[STEP] = Attributes::makeReal
("STEP", "Increment for randomized sequences (default: 1)", 1.0);
registerOwnership(AttributeHandler::STATEMENT);
}
......@@ -96,6 +101,7 @@ void OpalSample::initialize(const std::string &dvarName,
int seed = Attributes::getReal(itsAttr[SEED]);
size_m = Attributes::getReal(itsAttr[N]);
double step = Attributes::getReal(itsAttr[STEP]);
bool random = Attributes::getBool(itsAttr[RANDOM]);
......@@ -143,6 +149,26 @@ void OpalSample::initialize(const std::string &dvarName,
} else {
sampleMethod_m.reset( new LatinHyperCube(lower, upper) );
}
} else if (type == "RANDOM_SEQUENCE_UNIFORM_INT") {
if (Attributes::getReal(itsAttr[SEED])) {
sampleMethod_m.reset(
new SampleRandomizedSequence<int>(lower, upper, step, seed)
);
} else {
sampleMethod_m.reset(
new SampleRandomizedSequence<int>(lower, upper, step)
);
}
} else if (type == "RANDOM_SEQUENCE_UNIFORM") {
if (Attributes::getReal(itsAttr[SEED])) {
sampleMethod_m.reset(
new SampleRandomizedSequence<double>(lower, upper, step, seed)
);
} else {
sampleMethod_m.reset(
new SampleRandomizedSequence<double>(lower, upper, step)
);
}
} else {
throw OpalException("OpalSample::initialize()",
"Unknown sampling method: '" + type + "'.");
......@@ -153,4 +179,4 @@ void OpalSample::initialize(const std::string &dvarName,
std::string OpalSample::getVariable() const {
return Attributes::getString(itsAttr[VARIABLE]);
}
\ No newline at end of file
}
#ifndef OPAL_SAMPLE_WEIGHTED_SEQUENCE_H
#define OPAL_SAMPLE_WEIGHTED_SEQUENCE_H
#include "Sample/Uniform.h"
template <typename T>
class SampleRandomizedSequence : public SamplingMethod
{
public:
SampleRandomizedSequence(T lower, T upper, double step)
: unif_m(0, size_t((upper - lower) / step))