Skip to content

GitLab

Commit a2996d87 authored by kraus's avatar kraus
Browse files

Merge branch 'arrayTrack'

parents f345284b 975f8461
Hide whitespace changes
Inline Side-by-side
Showing with 196 additions and 181 deletions
src/Algorithms/ParallelTTracker.cpp
View file @ a2996d87
......@@ -97,12 +97,13 @@ surfaceStatus_m(false),
secondaryFlg_m(false),
mpacflg_m(true),
nEmissionMode_m(false),
zStop_m(0.0),
zStop_m(),
scaleFactor_m(1.0),
vscaleFactor_m(scaleFactor_m),
recpGamma_m(1.0),
rescale_coeff_m(1.0),
dtTrack_m(0.0),
dtCurrentTrack_m(0.0),
dtAllTracks_m(),
surfaceEmissionStop_m(-1),
specifiedNPart_m(N),
minStepforReBin_m(-1),
......@@ -112,7 +113,7 @@ lastVisited_m(-1),
numRefs_m(-1),
gunSubTimeSteps_m(-1),
emissionSteps_m(std::numeric_limits<unsigned int>::max()),
localTrackSteps_m(0),
localTrackSteps_m(),
maxNparts_m(0),
numberOfFieldEmittedParticles_m(std::numeric_limits<size_t>::max()),
bends_m(0),
......@@ -134,9 +135,10 @@ ParallelTTracker::ParallelTTracker(const Beamline &beamline,
const PartData &reference,
bool revBeam,
bool revTrack,
int maxSTEPS,
double zstop,
const std::vector<unsigned long long> &maxSTEPS,
const std::vector<double> &zstop,
int timeIntegrator,
const std::vector<double> &dt,
size_t N):
Tracker(beamline, reference, revBeam, revTrack),
itsBunch(&bunch),
......@@ -159,7 +161,8 @@ scaleFactor_m(itsBunch->getdT() * Physics::c),
vscaleFactor_m(scaleFactor_m),
recpGamma_m(1.0),
rescale_coeff_m(1.0),
dtTrack_m(0.0),
dtCurrentTrack_m(0.0),
dtAllTracks_m(dt),
surfaceEmissionStop_m(-1),
specifiedNPart_m(N),
minStepforReBin_m(-1),
......@@ -195,9 +198,10 @@ ParallelTTracker::ParallelTTracker(const Beamline &beamline,
const PartData &reference,
bool revBeam,
bool revTrack,
int maxSTEPS,
double zstop,
const std::vector<unsigned long long> &maxSTEPS,
const std::vector<double> &zstop,
int timeIntegrator,
const std::vector<double> &dt,
size_t N,
Amr* amrptr_in):
Tracker(beamline, reference, revBeam, revTrack),
......@@ -221,7 +225,8 @@ scaleFactor_m(itsBunch->getdT() * Physics::c),
vscaleFactor_m(scaleFactor_m),
recpGamma_m(1.0),
rescale_coeff_m(1.0),
dtTrack_m(0.0),
dtCurrentTrack_m(0.0),
dtAllTracks_m(dt),
surfaceEmissionStop_m(-1),
specifiedNPart_m(N),
minStepforReBin_m(-1),
......@@ -408,7 +413,7 @@ double ParallelTTracker::schottkyLoop(double rescale_coeff) {
size_t totalParticles_f = 0;
for(; step < localTrackSteps_m; ++step) {
for(; step < localTrackSteps_m.front(); ++step) {
global_EOL = true; // check if any particle hasn't reached the end of the field from the last element
itsOpalBeamline_m.resetStatus();
......@@ -698,8 +703,8 @@ double ParallelTTracker::schottkyLoop(double rescale_coeff) {
/**
Stop simulation if beyond zStop_m
*/
if(sposRef > zStop_m) {
localTrackSteps_m = step;
if(sposRef > zStop_m.front()) {
localTrackSteps_m.front() = step;
}
} else {
INFOMSG("Step " << step << " no emission yet " << " t= " << t << " [s]" << endl);
......@@ -874,7 +879,7 @@ FieldList ParallelTTracker::executeAutoPhaseForSliceTracker() {
for(FieldList::iterator it = monitors.begin(); it != monitors.end(); ++ it) {
double zbegin, zend;
it->getElement()->getDimensions(zbegin, zend);
if(zbegin < zStop_m && zend >= zStop_m) {
if(zbegin < zStop_m.front() && zend >= zStop_m.front()) {
msg << "\033[0;31m"
<< "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n"
<< "% Removing '" << it->getElement()->getName() << "' since it resides in two tracks. %\n"
......@@ -977,7 +982,7 @@ void ParallelTTracker::executeAutoPhase(int numRefs, double zStop) {
Vector_t rmin, rmax;
size_t maxStepsSave = localTrackSteps_m;
size_t maxStepsSave = localTrackSteps_m.front();
size_t step = 0;
double tSave = itsBunch->getT();
......@@ -1031,7 +1036,7 @@ void ParallelTTracker::executeAutoPhase(int numRefs, double zStop) {
double cavity_start = 0.0;
Component *cavity = NULL;
for(; step < localTrackSteps_m * dtfraction; ++step) {
for(; step < localTrackSteps_m.front() * dtfraction; ++step) {
if (itsBunch->getLocalNum() != 0) {
// let's do a drifting step to probe if the particle will reach element in next step
......@@ -1280,7 +1285,7 @@ void ParallelTTracker::executeAutoPhase(int numRefs, double zStop) {
reduce(sposRef,sposRef,OpAddAssign());
if(sposRef > zStop)
localTrackSteps_m = floor(step / dtfraction);
localTrackSteps_m.front() = floor(step / dtfraction);
INFOMSG("step = " << step << ", spos = " << sposRef << " [m], t= " << itsBunch->getT() << " [s], "
......@@ -1302,7 +1307,7 @@ void ParallelTTracker::executeAutoPhase(int numRefs, double zStop) {
bend->resetRecalcRefTrajFlag();
}
localTrackSteps_m = maxStepsSave;
localTrackSteps_m.front() = maxStepsSave;
scaleFactor_m = scaleFactorSave;
itsBunch->setT(tSave);
itsBunch->setdT(dTSave);
......@@ -1347,7 +1352,7 @@ void ParallelTTracker::Tracker_AMR()
const Vector_t vscaleFactor_m = Vector_t(scaleFactor_m);
BorisPusher pusher(itsReference);
secondaryFlg_m = false;
dtTrack_m = itsBunch->getdT();
dtCurrentTrack_m = itsBunch->getdT();
// upper limit of particle number when we do field emission and secondary emission
// simulation. Could be reset to another value in input file with MAXPARTSNUM.
......@@ -1369,7 +1374,7 @@ void ParallelTTracker::Tracker_AMR()
unsigned long long step = itsBunch->getLocalTrackStep();
msg << "Track start at: " << myt1.time() << ", t= " << t << "; zstop at: " << zStop_m << " [m]" << endl;
msg << "Track start at: " << myt1.time() << ", t= " << t << "; zstop at: " << zStop_m.front() << " [m]" << endl;
gunSubTimeSteps_m = 10;
prepareEmission();
......@@ -1377,7 +1382,7 @@ void ParallelTTracker::Tracker_AMR()
doSchottyRenormalization();
msg << "Executing ParallelTTracker, initial DT " << itsBunch->getdT() << " [s];\n"
<< "max integration steps " << localTrackSteps_m << ", next step= " << step << endl;
<< "max integration steps " << localTrackSteps_m.front() << ", next step= " << step << endl;
msg << "Using default (Boris-Buneman) integrator" << endl;
// itsBeamline_m.accept(*this);
......@@ -1405,7 +1410,7 @@ void ParallelTTracker::Tracker_AMR()
itsBunch->Ef = Vector_t(0.0);
itsBunch->Bf = Vector_t(0.0);
for(; step < localTrackSteps_m; ++step)
for(; step < localTrackSteps_m.front(); ++step)
{
bends_m = 0;
numberOfFieldEmittedParticles_m = 0;
......@@ -1545,7 +1550,7 @@ void ParallelTTracker::Tracker_Default() {
const Vector_t vscaleFactor_m = Vector_t(scaleFactor_m);
BorisPusher pusher(itsReference);
secondaryFlg_m = false;
dtTrack_m = itsBunch->getdT();
dtCurrentTrack_m = itsBunch->getdT();
// upper limit of particle number when we do field emission and secondary emission
// simulation. Could be reset to another value in input file with MAXPARTSNUM.
......@@ -1567,7 +1572,7 @@ void ParallelTTracker::Tracker_Default() {
unsigned long long step = itsBunch->getLocalTrackStep();
msg << "Track start at: " << myt1.time() << ", t= " << t << "; zstop at: " << zStop_m << " [m]" << endl;
msg << "Track start at: " << myt1.time() << ", t= " << t << "; zstop at: " << zStop_m.front() << " [m]" << endl;
gunSubTimeSteps_m = 10;
prepareEmission();
......@@ -1575,7 +1580,7 @@ void ParallelTTracker::Tracker_Default() {
doSchottyRenormalization();
msg << "Executing ParallelTTracker, initial DT " << itsBunch->getdT() << " [s];\n"
<< "max integration steps " << localTrackSteps_m << ", next step= " << step << endl;
<< "max integration steps " << localTrackSteps_m.front() << ", next step= " << step << endl;
msg << "Using default (Boris-Buneman) integrator" << endl;
// itsBeamline_m.accept(*this);
......@@ -1599,7 +1604,7 @@ void ParallelTTracker::Tracker_Default() {
wakeStatus_m = false;
surfaceStatus_m = false;
for(; step < localTrackSteps_m; ++step) {
for(; step < localTrackSteps_m.front(); ++step) {
bends_m = 0;
numberOfFieldEmittedParticles_m = 0;
......@@ -1685,7 +1690,7 @@ void ParallelTTracker::handleOverlappingMonitors() {
for(FieldList::iterator it = monitors.begin(); it != monitors.end(); ++ it) {
double zbegin, zend;
it->getElement()->getDimensions(zbegin, zend);
if(zbegin < zStop_m && zend >= zStop_m) {
if(zbegin < zStop_m.front() && zend >= zStop_m.front()) {
msg << "\033[0;31m"
<< "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n"
<< "% Removing '" << it->getElement()->getName() << "' since it resides in two tracks. %\n"
......@@ -2196,7 +2201,7 @@ void ParallelTTracker::timeIntegration2_bgf(BorisPusher & pusher) {
void ParallelTTracker::selectDT() {
double dt = dtTrack_m;
double dt = dtCurrentTrack_m;
itsBunch->setdT(dt);
scaleFactor_m = dt * Physics::c;
vscaleFactor_m = Vector_t(scaleFactor_m);
......@@ -2574,14 +2579,14 @@ void ParallelTTracker::dumpStats(long long step, bool psDump, bool statDump) {
// If we are dealing with field emission and secondary emission, set upper
// limit of particle number in simulation to prevent memory overflow.
if(numParticlesInSimulation_m > maxNparts_m)
localTrackSteps_m = step;
localTrackSteps_m.front() = step;
// ada reset Nimpact_m, does not make sense to integrate this we obtain a rediculus large number !!
Nimpact_m = 0;
}
if(sposRef > zStop_m)
localTrackSteps_m = step;
if(sposRef > zStop_m.front())
localTrackSteps_m.front() = step;
}
......@@ -2912,7 +2917,7 @@ void ParallelTTracker::execute() {
void ParallelTTracker::Tracker_AMTS() {
Inform msg("ParallelTTracker ");
const Vector_t vscaleFactor_m = Vector_t(scaleFactor_m);
dtTrack_m = itsBunch->getdT();
dtCurrentTrack_m = itsBunch->getdT();
// upper limit of particle number when we do field emission and secondary emission
// simulation. Could be reset to another value in input file with MAXPARTSNUM.
......@@ -2926,7 +2931,7 @@ void ParallelTTracker::Tracker_AMTS() {
numParticlesInSimulation_m = itsBunch->getTotalNum();
setTime();
unsigned long long step = itsBunch->getLocalTrackStep();
msg << "Track start at: " << OPALTimer::Timer().time() << ", t = " << itsBunch->getT() << "; zstop at: " << zStop_m << " [m]" << endl;
msg << "Track start at: " << OPALTimer::Timer().time() << ", t = " << itsBunch->getT() << "; zstop at: " << zStop_m.front() << " [m]" << endl;
msg << "Executing ParallelTTracker, next step = " << step << endl;
msg << "Using AMTS (adaptive multiple-time-stepping) integrator" << endl;
itsOpalBeamline_m.print(msg);
......@@ -2972,7 +2977,7 @@ void ParallelTTracker::Tracker_AMTS() {
msg << "AMTS initialization: dt_outer = " << dt_outer << " deltaTau = " << deltaTau << endl;
// AMTS calculation of stopping times
double const tEnd = itsBunch->getT() + double(localTrackSteps_m - step) * dt_inner_target;
double const tEnd = itsBunch->getT() + double(localTrackSteps_m.front() - step) * dt_inner_target;
double const psDumpInterval = double(Options::psDumpFreq) * dt_inner_target;
double const statDumpInterval = double(Options::statDumpFreq) * dt_inner_target;
double const repartInterval = double(repartFreq_m) * dt_inner_target;
......
src/Algorithms/ParallelTTracker.h
View file @ a2996d87
......@@ -124,13 +124,17 @@ public:
// If [b]revBeam[/b] is true, the beam runs from s = C to s = 0.
// If [b]revTrack[/b] is true, we track against the beam.
explicit ParallelTTracker(const Beamline &bl, PartBunch &bunch, DataSink &ds,
const PartData &data, bool revBeam, bool revTrack, int maxSTEPS, double zstop, int timeIntegrator, size_t N);
const PartData &data, bool revBeam, bool revTrack,
const std::vector<unsigned long long> &maxSTEPS, const std::vector<double> &zstop,
int timeIntegrator, const std::vector<double> &dt, size_t N);
/// Constructor
// Amr pointer is taken
#ifdef HAVE_AMR_SOLVER
explicit ParallelTTracker(const Beamline &bl, PartBunch &bunch, DataSink &ds,
const PartData &data, bool revBeam, bool revTrack, int maxSTEPS, double zstop, int timeIntegrator, size_t N, Amr* amrptr_in);
const PartData &data, bool revBeam, bool revTrack,
const std::vector<unsigned long long> &maxSTEPS, const std::vector<double> &zstop,
int timeIntegrator, const std::vector<double> &dt, size_t N, Amr* amrptr_in);
#endif
virtual ~ParallelTTracker();
......@@ -282,7 +286,7 @@ private:
bool nEmissionMode_m;
/// where to stop
double zStop_m;
std::vector<double> zStop_m;
/// The scale factor for dimensionless variables (FIXME: move to PartBunch)
double scaleFactor_m;
......@@ -294,7 +298,8 @@ private:
double rescale_coeff_m;
double dtTrack_m;
double dtCurrentTrack_m;
std::vector<double> dtAllTracks_m;
double surfaceEmissionStop_m;
......@@ -324,7 +329,7 @@ private:
unsigned int emissionSteps_m;
/// The maximal number of steps the system is integrated per TRACK
unsigned long long localTrackSteps_m;
std::vector<unsigned long long> localTrackSteps_m;
size_t maxNparts_m;
size_t numberOfFieldEmittedParticles_m;
......@@ -835,8 +840,8 @@ inline void ParallelTTracker::writePhaseSpace(const long long step, const double
msg << "* Wrote beam statistics." << endl;
itsDataSink_m->writeStatData(*itsBunch, FDext, rmax(2), sposRef, rmin(2), collimatorLosses);
}
// itsBunch->printBinHist();
}
#endif // OPAL_ParallelTTracker_HH
#endif // OPAL_ParallelTTracker_HH
\ No newline at end of file
src/Track/Track.cpp
View file @ a2996d87
......@@ -37,7 +37,7 @@ otherwise a new bunch is allocated in the dictionary.
Track::Track(BeamSequence *u, const PartData &ref, const std::vector<double> & dt,
const std::vector<int> & maxtsteps, int stepsperturn,
const std::vector<unsigned long long> & maxtsteps, int stepsperturn,
const std::vector<double> & zStop, int timeintegrator, int nslices,
double t0, double dtScInit, double deltaTau):
reference(ref),
......
src/Track/Track.h
View file @ a2996d87
......@@ -36,7 +36,7 @@ class Track {
public:
Track(BeamSequence *, const PartData &, const std::vector<double> & dt,
const std::vector<int> & maxtsteps, int stepsperturn,
const std::vector<unsigned long long> & maxtsteps, int stepsperturn,
const std::vector<double> & zStop, int timeintegrator, int nslices,
double t0, double dtScInit, double deltaTau);
~Track();
......@@ -69,7 +69,7 @@ public:
double t0_m;
/// Maximal number of timesteps
std::vector<int> localTimeSteps;
std::vector<unsigned long long> localTimeSteps;
/// The timsteps per revolution period. ONLY available for OPAL-cycl.
int stepsPerTurn;
......
src/Track/TrackCmd.cpp
View file @ a2996d87
......@@ -120,14 +120,19 @@ std::vector<double> TrackCmd::getZSTOP() const {
return zstop;
}
std::vector<int> TrackCmd::getMAXSTEPS() const {
std::vector<unsigned long long> TrackCmd::getMAXSTEPS() const {
std::vector<double> maxsteps_d = Attributes::getRealArray(itsAttr[MAXSTEPS]);
std::vector<int> maxsteps_i;
std::vector<unsigned long long> maxsteps_i;
if (maxsteps_d.size() == 0) {
maxsteps_i.push_back(10);
maxsteps_i.push_back(10ul);
}
for (auto it = maxsteps_d.begin(); it != maxsteps_d.end(); ++ it) {
maxsteps_i.push_back(*it);
if (*it < 0) {
maxsteps_i.push_back(10);
} else {
unsigned long long value = *it;
maxsteps_i.push_back(value);
}
}
return maxsteps_i;
......@@ -166,7 +171,7 @@ void TrackCmd::execute() {
std::vector<double> dt = getDT();
double t0 = getT0();
std::vector<int> maxsteps = getMAXSTEPS();
std::vector<unsigned long long> maxsteps = getMAXSTEPS();
int stepsperturn = getSTEPSPERTURN();
std::vector<double> zstop = getZSTOP();
int timeintegrator = getTIMEINTEGRATOR();
......
src/Track/TrackCmd.h
View file @ a2996d87
......@@ -51,7 +51,7 @@ public:
double getT0() const;
/// Return the maximum timsteps we integrate the system
std::vector<int> getMAXSTEPS() const;
std::vector<unsigned long long> getMAXSTEPS() const;
/// Return the timsteps per revolution period. ONLY available for OPAL-cycl.
/// In OPAL-cycl, timestep is calculated by STEPSPERTURN, rather than given in TRACK command.
......
src/Track/TrackRun.cpp
View file @ a2996d87
......@@ -270,21 +270,21 @@ void TrackRun::execute() {
Beam *beam = Beam::find(Attributes::getString(itsAttr[BEAM]));
if (Attributes::getString(itsAttr[BOUNDARYGEOMETRY]) != "NONE") {
// Ask the dictionary if BoundaryGeometry is allocated.
// If it is allocated use the allocated BoundaryGeometry
if (!OpalData::getInstance()->hasGlobalGeometry()) {
BoundaryGeometry *bg = BoundaryGeometry::find(Attributes::getString(itsAttr[BOUNDARYGEOMETRY]))->
clone(getOpalName() + std::string("_geometry"));
OpalData::getInstance()->setGlobalGeometry(bg);
}
// Ask the dictionary if BoundaryGeometry is allocated.
// If it is allocated use the allocated BoundaryGeometry
if (!OpalData::getInstance()->hasGlobalGeometry()) {
BoundaryGeometry *bg = BoundaryGeometry::find(Attributes::getString(itsAttr[BOUNDARYGEOMETRY]))->
clone(getOpalName() + std::string("_geometry"));
OpalData::getInstance()->setGlobalGeometry(bg);
}
}
fs = FieldSolver::find(Attributes::getString(itsAttr[FIELDSOLVER]));
fs->initCartesianFields();
Track::block->bunch->setSolver(fs);
if (fs->hasPeriodicZ())
Track::block->bunch->setBCForDCBeam();
Track::block->bunch->setBCForDCBeam();
else
Track::block->bunch->setBCAllOpen();
Track::block->bunch->setBCAllOpen();
double charge = SetDistributionParallelT(beam);
......@@ -314,156 +314,156 @@ void TrackRun::execute() {
}
#ifdef HAVE_AMR_SOLVER
Amr* amrptr;
// if (fs->isAMRSolver()) {
*gmsg << "A M R Initialization " << endl;
*gmsg << *Track::block->bunch << endl;
*gmsg << *fs << endl;
// if (fs->isAMRSolver()) {
*gmsg << "A M R Initialization " << endl;
*gmsg << *Track::block->bunch << endl;
*gmsg << *fs << endl;
FieldLayout<3>::iterator_iv locDomBegin = Track::block->bunch->getFieldLayout().begin_iv();
FieldLayout<3>::iterator_iv locDomEnd = Track::block->bunch->getFieldLayout().end_iv();
FieldLayout<3>::iterator_dv globDomBegin = Track::block->bunch->getFieldLayout().begin_rdv();
FieldLayout<3>::iterator_dv globDomEnd = Track::block->bunch->getFieldLayout().end_rdv();
FieldLayout<3>::iterator_iv locDomBegin = Track::block->bunch->getFieldLayout().begin_iv();
FieldLayout<3>::iterator_iv locDomEnd = Track::block->bunch->getFieldLayout().end_iv();
FieldLayout<3>::iterator_dv globDomBegin = Track::block->bunch->getFieldLayout().begin_rdv();
FieldLayout<3>::iterator_dv globDomEnd = Track::block->bunch->getFieldLayout().end_rdv();
NDIndex<3> ipplDom = Track::block->bunch->getFieldLayout().getDomain();
NDIndex<3> ipplDom = Track::block->bunch->getFieldLayout().getDomain();
BoxArray lev0_grids(Ippl::getNodes());
BoxArray lev0_grids(Ippl::getNodes());
Array<int> procMap;
procMap.resize(lev0_grids.size()+1); // +1 is a historical thing, do not ask
Array<int> procMap;
procMap.resize(lev0_grids.size()+1); // +1 is a historical thing, do not ask
// first iterate over the local owned domain(s)
for(FieldLayout<3>::const_iterator_iv v_i = locDomBegin ; v_i != locDomEnd; ++v_i) {
// first iterate over the local owned domain(s)
for(FieldLayout<3>::const_iterator_iv v_i = locDomBegin ; v_i != locDomEnd; ++v_i) {
std::ostringstream stream;
stream << *((*v_i).second);
std::pair<Box,unsigned int> res = getBlGrids(stream.str());
lev0_grids.set(res.second,res.first);
procMap[res.second] = Ippl::myNode();
}
}
// then iterate over the non-local domain(s)
for(FieldLayout<3>::iterator_dv v_i = globDomBegin ; v_i != globDomEnd; ++v_i) {
// then iterate over the non-local domain(s)
for(FieldLayout<3>::iterator_dv v_i = globDomBegin ; v_i != globDomEnd; ++v_i) {
std::ostringstream stream;
stream << *((*v_i).second);
std::pair<Box,unsigned int> res = getBlGrids(stream.str());
lev0_grids.set(res.second,res.first);
procMap[res.second] = res.second;
}
procMap[lev0_grids.size()] = Ippl::myNode();
// This init call will cache the distribution map as determined by procMap
// so that all data will end up on the right processor
RealBox rb;
Array<Real> prob_lo(3);
Array<Real> prob_hi(3);
Vector_t rmin;
Vector_t rmax;
Track::block->bunch->get_bounds(rmin, rmax);
// Set up a problem size with dx = dy = dz and large enough
// to hold the geometric boundary.
prob_lo.set(0,-0.025);
prob_lo.set(1,-0.025);
prob_lo.set(2,-0.050);
prob_hi.set(0, 0.025);
prob_hi.set(1, 0.025);
prob_hi.set(2, 0.050);
rb.setLo(prob_lo);
rb.setHi(prob_hi);
int coord_sys = 0;
Array<int> ncell(3);
ncell[0] = ipplDom[0].length();
ncell[1] = ipplDom[1].length();
ncell[2] = ipplDom[2].length();
std::vector<int > nr(3);
std::vector<double> hr(3);
std::vector<double> prob_lo_in(3);
for (int i = 0; i < 3; i++)
{
nr[i] = ncell[i];
hr[i] = (prob_hi[i] - prob_lo[i]) / ncell[i];
prob_lo_in[i] = prob_lo[i];
}
}
procMap[lev0_grids.size()] = Ippl::myNode();
// This init call will cache the distribution map as determined by procMap
// so that all data will end up on the right processor
RealBox rb;
Array<Real> prob_lo(3);
Array<Real> prob_hi(3);
Vector_t rmin;
Vector_t rmax;
Track::block->bunch->get_bounds(rmin, rmax);
// Set up a problem size with dx = dy = dz and large enough
// to hold the geometric boundary.
prob_lo.set(0,-0.025);
prob_lo.set(1,-0.025);
prob_lo.set(2,-0.050);
prob_hi.set(0, 0.025);
prob_hi.set(1, 0.025);
prob_hi.set(2, 0.050);
rb.setLo(prob_lo);
rb.setHi(prob_hi);
int coord_sys = 0;
Array<int> ncell(3);
ncell[0] = ipplDom[0].length();
ncell[1] = ipplDom[1].length();
ncell[2] = ipplDom[2].length();
std::vector<int > nr(3);
std::vector<double> hr(3);
std::vector<double> prob_lo_in(3);
for (int i = 0; i < 3; i++)
{
nr[i] = ncell[i];
hr[i] = (prob_hi[i] - prob_lo[i]) / ncell[i];
prob_lo_in[i] = prob_lo[i];
}
// We set this to -1 so that we can now control max_lev from the inputs file
int maxLevel = -1;
// We set this to -1 so that we can now control max_lev from the inputs file
int maxLevel = -1;
amrptr = new Amr(&rb,maxLevel,ncell,coord_sys);
amrptr = new Amr(&rb,maxLevel,ncell,coord_sys);
Real strt_time = 0.0;
Real stop_time = 1.0;
Real strt_time = 0.0;
Real stop_time = 1.0;
// This init call will cache the distribution map as determined by procMap // so that all data will end up on the right processor
amrptr->InitializeInit(strt_time, stop_time, &lev0_grids, &procMap);
// This init call will cache the distribution map as determined by procMap // so that all data will end up on the right processor
amrptr->InitializeInit(strt_time, stop_time, &lev0_grids, &procMap);
BoundaryDomain* bd = new BoundaryDomain(nr,hr);
amrptr->setBoundaryGeometry(bd->GetIntersectLoX(), bd->GetIntersectHiX(),
bd->GetIntersectLoY(), bd->GetIntersectHiY());
BoundaryDomain* bd = new BoundaryDomain(nr,hr);
amrptr->setBoundaryGeometry(bd->GetIntersectLoX(), bd->GetIntersectHiX(),
bd->GetIntersectLoY(), bd->GetIntersectHiY());
std::vector<double> x(3);
std::vector<double> attr(11);
std::vector<double> x(3);
std::vector<double> attr(11);
for (size_t i=0; i<Track::block->bunch->getLocalNum(); i++)
{
// X, Y, Z are stored separately from the other attributes
for (unsigned int k=0; k<3; k++)
x[k] = Track::block->bunch->R[i](k);
for (size_t i=0; i<Track::block->bunch->getLocalNum(); i++)
{
// X, Y, Z are stored separately from the other attributes
for (unsigned int k=0; k<3; k++)
x[k] = Track::block->bunch->R[i](k);