TrackRun.cpp 45.1 KB
Newer Older
gsell's avatar
gsell committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
// -----------------------------------------------------------------------
// $RCSfile: TrackRun.cpp,v $
// ------------------------------------------------------------------------
// $Revision: 1.1.1.1.4.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: TrackRun
//   The class for the OPAL RUN command.
//
// ------------------------------------------------------------------------
//
// $Date: 2004/11/12 20:10:11 $
// $Author: adelmann $
//
// ------------------------------------------------------------------------

#include "Track/TrackRun.h"
#include "AbstractObjects/BeamSequence.h"
#include "AbstractObjects/OpalData.h"
#include "AbstractObjects/ObjectFunction.h"
23
#include "Algorithms/Tracker.h"
gsell's avatar
gsell committed
24 25 26
#include "Algorithms/ThinTracker.h"
#include "Algorithms/ThickTracker.h"

27 28
#include "Algorithms/bet/EnvelopeBunch.h"

gsell's avatar
gsell committed
29 30 31
#include "Algorithms/ParallelTTracker.h"
#include "Algorithms/ParallelSliceTracker.h"
#include "Algorithms/ParallelCyclotronTracker.h"
32
#include "Algorithms/AutophaseTracker.h"
gsell's avatar
gsell committed
33 34 35 36 37 38 39 40 41 42 43 44 45

#include "Attributes/Attributes.h"
#include "Beamlines/TBeamline.h"

#include "BasicActions/Option.h"

#include "Elements/OpalBeamBeam3D.h"
#include "Track/Track.h"
#include "Utilities/OpalException.h"
#include "Utilities/Round.h"
#include "Structure/Beam.h"
#include "Structure/FieldSolver.h"
#include "Structure/DataSink.h"
46 47 48 49
#include "Structure/H5PartWrapper.h"
#include "Structure/H5PartWrapperForPT.h"
#include "Structure/H5PartWrapperForPC.h"
#include "Structure/H5PartWrapperForPS.h"
gsell's avatar
gsell committed
50
#include "Distribution/Distribution.h"
51
#include "Structure/BoundaryGeometry.h"
gsell's avatar
gsell committed
52

53 54 55 56 57 58 59
#ifdef HAVE_AMR_SOLVER
#define DIM 3
#include <Amr.H>
#include <ParallelDescriptor.H>
#include "Solvers/amr/BoundaryDomain.h"
#endif

gsell's avatar
gsell committed
60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79
#include <fstream>
#include <iomanip>

extern Inform *gmsg;

using namespace Physics;

// ------------------------------------------------------------------------

namespace {

    // The attributes of class TrackRun.
    enum {
        METHOD,       // Tracking method to use.
        FNAME,        // The name of file to be written.
        TURNS,        // The number of turns to be tracked.
        MBMODE,       // The working way for multi-bunch mode for OPAL-cycl: "FORCE" or "AUTO"
        PARAMB,       // The control parameter for "AUTO" mode of multi-bunch
        BEAM,         // The beam to track
        FIELDSOLVER,  // The field solver attached
80
        BOUNDARYGEOMETRY, // The boundary geometry
gsell's avatar
gsell committed
81 82 83 84 85 86 87
        DISTRIBUTION, // The particle distribution
        MULTIPACTING, // MULTIPACTING flag
        // THE INTEGRATION TIMESTEP IN SEC
        SIZE
    };
}

88 89
const std::string TrackRun::defaultDistribution("DISTRIBUTION");

gsell's avatar
gsell committed
90 91 92
TrackRun::TrackRun():
    Action(SIZE, "RUN",
           "The \"RUN\" sub-command tracks the defined particles through "
93 94 95 96 97
           "the given lattice."),
    itsTracker(NULL),
    dist(NULL),
    distrs_m(),
    fs(NULL),
98 99
    ds(NULL),
    phaseSpaceSink_m(NULL) {
gsell's avatar
gsell committed
100 101
    itsAttr[METHOD] = Attributes::makeString
                      ("METHOD", "Name of tracking algorithm to use:\n"
102
                       "\t\t\t\"THIN\" (default) or \"THICK,PARALLEL-T,CYCLOTRON-T,PARALLEL-SLICE,AUTOPHASE\".", "THIN");
gsell's avatar
gsell committed
103 104 105 106 107 108 109 110 111 112
    itsAttr[TURNS] = Attributes::makeReal
                     ("TURNS", "Number of turns to be tracked; Number of neighboring bunches to be tracked in cyclotron", 1.0);

    itsAttr[MBMODE] = Attributes::makeString
                      ("MBMODE", "The working way for multi-bunch mode for OPAL-cycl: FORCE or AUTO ", "FORCE");

    itsAttr[PARAMB] = Attributes::makeReal
                      ("PARAMB", " Control parameter to define when to start multi-bunch mode, only available in \"AUTO\" mode ", 5.0);

    itsAttr[FNAME] = Attributes::makeString
kraus's avatar
kraus committed
113
                     ("FILE", "Name of file to be written", "TRACK");
gsell's avatar
gsell committed
114 115 116 117 118

    itsAttr[BEAM] = Attributes::makeString
                    ("BEAM", "Name of beam ", "BEAM");
    itsAttr[FIELDSOLVER] = Attributes::makeString
                           ("FIELDSOLVER", "Field solver to be used ", "FIELDSOLVER");
119 120
    itsAttr[BOUNDARYGEOMETRY] = Attributes::makeString
                           ("BOUNDARYGEOMETRY", "Boundary geometry to be used NONE (default)", "NONE");
121 122
    itsAttr[DISTRIBUTION] = Attributes::makeStringArray
                             ("DISTRIBUTION", "List of particle distributions to be used ");
gsell's avatar
gsell committed
123 124 125 126 127 128
    itsAttr[MULTIPACTING] = Attributes::makeBool
                            ("MULTIPACTING", "Multipacting flag, default: false. Set true to initialize primary particles according to BoundaryGeometry", false);
    OPAL = OpalData::getInstance();
}


129
TrackRun::TrackRun(const std::string &name, TrackRun *parent):
130 131 132 133 134
    Action(name, parent),
    itsTracker(NULL),
    dist(NULL),
    distrs_m(),
    fs(NULL),
135 136
    ds(NULL),
    phaseSpaceSink_m(NULL) {
gsell's avatar
gsell committed
137 138 139 140 141
    OPAL = OpalData::getInstance();
}


TrackRun::~TrackRun()
142 143 144 145
{
    delete phaseSpaceSink_m;
    phaseSpaceSink_m = NULL;
}
gsell's avatar
gsell committed
146 147


148
TrackRun *TrackRun::clone(const std::string &name) {
gsell's avatar
gsell committed
149 150 151 152 153 154 155
    return new TrackRun(name, this);
}


void TrackRun::execute() {

    // Get algorithm to use.
156
    std::string method = Attributes::getString(itsAttr[METHOD]);
gsell's avatar
gsell committed
157 158 159 160 161 162 163 164 165 166 167
    if(method == "THIN") {
        //std::cerr << "  method == \"THIN\"" << std::endl;
        itsTracker = new ThinTracker(*Track::block->use->fetchLine(),
                                     *Track::block->bunch, Track::block->reference,
                                     false, false);
    } else if(method == "THICK") {
        //std::cerr << "  method == \"THICK\"" << std::endl;
        itsTracker = new ThickTracker(*Track::block->use->fetchLine(),
                                      *Track::block->bunch, Track::block->reference,
                                      false, false);
    } else if(method == "PARALLEL-SLICE") {
168 169 170 171 172
        setupSliceTracker();
    } else if(method == "PARALLEL-T") {
        setupTTracker();
    } else if(method == "PARALLEL-Z") {
        *gmsg << "  method == \"PARALLEL-Z\"" << endl;
173

174 175
    } else if(method == "CYCLOTRON-T") {
        setupCyclotronTracker();
176 177 178 179
    } else if(method == "AUTOPHASE") {
        executeAutophaseTracker();

        return;
180 181 182 183
    } else {
        throw OpalException("TrackRun::execute()",
                            "Method name \"" + method + "\" unknown.");
    }
184

185 186 187 188 189 190 191 192 193 194
    if(method == "THIN" || method == "THICK") {
        /*
          OLD SERIAL STUFF
        */
        // Open output file.
        std::string file = Attributes::getString(itsAttr[FNAME]);
        std::ofstream os(file.c_str());
        if(os.bad()) {
            throw OpalException("TrackRun::execute()",
                                "Unable to open output file \"" + file + "\".");
gsell's avatar
gsell committed
195 196
        }

197 198 199
        // Print initial conditions.
        os << "\nInitial particle positions:\n"
           << itsTracker->getBunch() << std::endl;
200

201 202 203 204 205 206 207 208 209
        int turns = int(Round(Attributes::getReal(itsAttr[TURNS])));
        // Track for the all but last turn.
        for(int turn = 1; turn < turns; ++turn) {
            itsTracker->execute();
            os << "\nParticle positions after turn " << turn << ":\n"
               << itsTracker->getBunch() << std::endl;
        }
        // Track last turn, with statistics.
        itsTracker->execute();
210

211 212 213 214 215 216 217 218
        // Print final conditions.
        os << "Particle positions after turn " << turns << ":\n"
           << itsTracker->getBunch() << std::endl;
        //    Track::block->bunch = itsTracker->getBunch();
    } else {
        itsTracker->execute();
        OPAL->setRestartRun(false);
    }
gsell's avatar
gsell committed
219

220 221 222
    OPAL->bunchIsAllocated();
    if(method == "PARALLEL-SLICE")
        OPAL->slbunchIsAllocated();
gsell's avatar
gsell committed
223

224 225
    delete itsTracker;
}
gsell's avatar
gsell committed
226

227 228
void TrackRun::setupSliceTracker() {
    OpalData::getInstance()->setInOPALEnvMode();
229
    bool isFollowupTrack = OPAL->hasSLBunchAllocated() && !Options::scan;
230 231
    if(!OPAL->hasSLBunchAllocated()) {
        *gmsg << "* ********************************************************************************** " << endl;
232
        *gmsg << "* Selected Tracking Method == PARALLEL-SLICE, NEW TRACK" << endl;
233
        *gmsg << "* ********************************************************************************** " << endl;
234
    } else if(isFollowupTrack) {
235
        *gmsg << "* ********************************************************************************** " << endl;
236
        *gmsg << "* Selected Tracking Method == PARALLEL-SLICE, FOLLOWUP TRACK" << endl;
237 238 239
        *gmsg << "* ********************************************************************************** " << endl;
    } else if(OPAL->hasSLBunchAllocated() && Options::scan) {
        *gmsg << "* ********************************************************************************** " << endl;
240
        *gmsg << "* Selected Tracking Method == PARALLEL-SLICE, SCAN TRACK" << endl;
241 242
        *gmsg << "* ********************************************************************************** " << endl;
    }
gsell's avatar
gsell committed
243

244
    Beam   *beam = Beam::find(Attributes::getString(itsAttr[BEAM]));
gsell's avatar
gsell committed
245

246 247 248 249 250
    if(OPAL->inRestartRun()) {
        phaseSpaceSink_m = new H5PartWrapperForPS(OPAL->getInputBasename() + std::string(".h5"),
                                                  OPAL->getRestartStep(),
                                                  OpalData::getInstance()->getRestartFileName(),
                                                  H5_O_WRONLY);
251 252 253 254 255
    } else if (isFollowupTrack) {
        phaseSpaceSink_m = new H5PartWrapperForPS(OPAL->getInputBasename() + std::string(".h5"),
                                                  -1,
                                                  OPAL->getInputBasename() + std::string(".h5"),
                                                  H5_O_WRONLY);
256 257 258 259 260
    } else {
        phaseSpaceSink_m = new H5PartWrapperForPS(OPAL->getInputBasename() + std::string(".h5"),
                                                  H5_O_WRONLY);
    }

261 262 263 264 265 266 267
    std::vector<std::string> distr_str = Attributes::getStringArray(itsAttr[DISTRIBUTION]);
    const size_t numberOfDistributions = distr_str.size();
    if (numberOfDistributions == 0) {
        dist = Distribution::find(defaultDistribution);
    } else {
        dist = Distribution::find(distr_str.at(0));
        dist->setNumberOfDistributions(numberOfDistributions);
gsell's avatar
gsell committed
268

269 270 271 272
        if(numberOfDistributions > 1) {
            *gmsg << "Found more than one distribution: ";
            for(size_t i = 1; i < numberOfDistributions; ++ i) {
                Distribution *d = Distribution::find(distr_str.at(i));
gsell's avatar
gsell committed
273

274 275
                d->setNumberOfDistributions(numberOfDistributions);
                distrs_m.push_back(d);
gsell's avatar
gsell committed
276

277
                *gmsg << " " << distr_str.at(i);
278
            }
279
            *gmsg << endl;
280
        }
281
    }
adelmann's avatar
adelmann committed
282

283 284
    fs = FieldSolver::find(Attributes::getString(itsAttr[FIELDSOLVER]));
    fs->initCartesianFields();
gsell's avatar
gsell committed
285

286
    double charge = 0.0;
gsell's avatar
gsell committed
287

288 289
    if(!OPAL->hasSLBunchAllocated()) {
        if(!OPAL->inRestartRun()) {
adelmann's avatar
adelmann committed
290

291 292
            dist->CreateOpalE(beam, distrs_m, Track::block->slbunch, 0.0, 0.0);
            OPAL->setGlobalPhaseShift(0.5 * dist->GetTEmission());
adelmann's avatar
adelmann committed
293

gsell's avatar
gsell committed
294
        } else {
295 296 297
            /***
                reload slice distribution
            */
gsell's avatar
gsell committed
298

299
            dist->DoRestartOpalE(*Track::block->slbunch, beam->getNumberOfParticles(), OPAL->getRestartStep(), phaseSpaceSink_m);
gsell's avatar
gsell committed
300
        }
301 302 303
    } else {
        charge = 1.0;
    }
304

305 306 307 308 309 310 311 312
    Track::block->slbunch->setdT(Track::block->dT.front());
    // set the total charge
    charge = beam->getCharge() * beam->getCurrent() / beam->getFrequency();
    Track::block->slbunch->setCharge(charge);
    // set coupling constant
    double coefE = 1.0 / (4 * pi * epsilon_0);
    Track::block->slbunch->setCouplingConstant(coefE);
    //Track::block->slbunch->calcBeamParameters();
313 314


315
    if(!OPAL->inRestartRun()) {
316 317 318 319 320 321 322 323 324
        if(!OPAL->hasDataSinkAllocated()) {
            OPAL->setDataSink(new DataSink(phaseSpaceSink_m));
        } else {
            ds = OPAL->getDataSink();
            ds->changeH5Wrapper(phaseSpaceSink_m);
        }
    } else {
        OPAL->setDataSink(new DataSink(phaseSpaceSink_m, -1));
    }
325

326
    ds = OPAL->getDataSink();
327

328 329 330 331 332 333 334 335
    if(!OPAL->hasBunchAllocated())
        *gmsg << *dist << endl;
    *gmsg << *beam << endl;
    *gmsg << *Track::block->slbunch  << endl;
    *gmsg << "Phase space dump frequency is set to " << Options::psDumpFreq
          << " Inputfile is " << OPAL->getInputFn() << endl;


336 337
    findPhasesForMaxEnergy();

338 339 340 341 342 343
    itsTracker = new ParallelSliceTracker(*Track::block->use->fetchLine(),
                                          dynamic_cast<EnvelopeBunch &>(*Track::block->slbunch),
                                          *ds,
                                          Track::block->reference,
                                          false, false,
                                          Track::block->localTimeSteps.front(),
344
                                          Track::block->zstop.front());
345
}
346

347 348 349
void TrackRun::setupTTracker(){
    OpalData::getInstance()->setInOPALTMode();
    bool isFollowupTrack = (OPAL->hasBunchAllocated() && !Options::scan);
350

351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374
    if(OPAL->hasBunchAllocated()) {
        if (Options::scan) {
            *gmsg << "* ********************************************************************************** " << endl;
            *gmsg << "* Selected Tracking Method == PARALLEL-T, FOLLOWUP TRACK in SCAN MODE" << endl;
            *gmsg << "* ********************************************************************************** " << endl;
            Track::block->bunch->setLocalTrackStep(0);
            Track::block->bunch->setGlobalTrackStep(0);
        } else {
            *gmsg << "* ********************************************************************************** " << endl;
            *gmsg << "* Selected Tracking Method == PARALLEL-T, FOLLOWUP TRACK" << endl;
            *gmsg << "* ********************************************************************************** " << endl;
            Track::block->bunch->setLocalTrackStep(0);
        }
    } else {
        if (Options::scan) {
            *gmsg << "* ********************************************************************************** " << endl;
            *gmsg << "* Selected Tracking Method == PARALLEL-T, NEW TRACK in SCAN MODE" << endl;
            *gmsg << "* ********************************************************************************** " << endl;
        } else {
            *gmsg << "* ********************************************************************************** " << endl;
            *gmsg << "* Selected Tracking Method == PARALLEL-T, NEW TRACK" << endl;
            *gmsg << "* ********************************************************************************** " << endl;
        }
    }
375

376 377 378 379 380 381 382 383
    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);
384
        }
385
    }
386

387
    setupFieldsolver();
388

389 390 391 392 393
    if(OPAL->inRestartRun()) {
        phaseSpaceSink_m = new H5PartWrapperForPT(OPAL->getInputBasename() + std::string(".h5"),
                                                  OPAL->getRestartStep(),
                                                  OpalData::getInstance()->getRestartFileName(),
                                                  H5_O_WRONLY);
394 395 396 397 398
    } else if (isFollowupTrack) {
        phaseSpaceSink_m = new H5PartWrapperForPT(OPAL->getInputBasename() + std::string(".h5"),
                                                  -1,
                                                  OPAL->getInputBasename() + std::string(".h5"),
                                                  H5_O_WRONLY);
399 400 401 402 403
    } else {
        phaseSpaceSink_m = new H5PartWrapperForPT(OPAL->getInputBasename() + std::string(".h5"),
                                                  H5_O_WRONLY);
    }

404
    double charge = setDistributionParallelT(beam);
405

406 407 408
    Track::block->bunch->setdT(Track::block->dT.front());
    Track::block->bunch->dtScInit_m = Track::block->dtScInit;
    Track::block->bunch->deltaTau_m = Track::block->deltaTau;
409

410 411
    if (!isFollowupTrack && !OPAL->inRestartRun())
        Track::block->bunch->setT(Track::block->t0_m);
412

413 414 415 416 417 418
    bool mpacflg = Attributes::getBool(itsAttr[MULTIPACTING]);
    if(!mpacflg) {
        Track::block->bunch->setCharge(charge);
        // set coupling constant
        double coefE = 1.0 / (4 * pi * epsilon_0);
        Track::block->bunch->setCouplingConstant(coefE);
419

420

421 422 423 424 425 426
        // statistical data are calculated (rms, eps etc.)
        Track::block->bunch->calcBeamParameters();
    } else {
        Track::block->bunch->setChargeZeroPart(charge);// just set bunch->qi_m=charge, don't set bunch->Q[] as we have not initialized any particle yet.
        Track::block->bunch->calcBeamParametersInitial();// we have not initialized any particle yet.
    }
427

428 429
#ifdef HAVE_AMR_SOLVER
    Amr *amrptr = setupAMRSolver();
430
#endif
431 432 433

    if(!OPAL->inRestartRun()) {
        if(!OPAL->hasDataSinkAllocated() && !Options::scan) {
434
            OPAL->setDataSink(new DataSink(phaseSpaceSink_m));
435 436 437 438
        } else if(Options::scan) {
            ds = OPAL->getDataSink();
            if(ds)
                delete ds;
439 440 441 442
            OPAL->setDataSink(new DataSink(phaseSpaceSink_m));
        } else {
            ds = OPAL->getDataSink();
            ds->changeH5Wrapper(phaseSpaceSink_m);
gsell's avatar
gsell committed
443
        }
444
    } else {
445
        OPAL->setDataSink(new DataSink(phaseSpaceSink_m, -1));//OPAL->getRestartStep()));
446
    }
gsell's avatar
gsell committed
447

448
    ds = OPAL->getDataSink();
gsell's avatar
gsell committed
449

450 451
    if(OPAL->hasBunchAllocated() && Options::scan)
        ds->reset();
gsell's avatar
gsell committed
452

453 454 455 456 457
    if(!OPAL->hasBunchAllocated() || Options::scan) {
        if(!mpacflg) {
            *gmsg << *dist << endl;
        } else {
            *gmsg << "* Multipacting flag is true. The particle distribution in the run command will be ignored " << endl;
gsell's avatar
gsell committed
458
        }
459
    }
gsell's avatar
gsell committed
460

461 462 463
    *gmsg << *beam << endl;
    *gmsg << *fs   << endl;
    *gmsg << *Track::block->bunch  << endl;
gsell's avatar
gsell committed
464

465 466
    findPhasesForMaxEnergy();

kraus's avatar
kraus committed
467 468
    *gmsg << level2
          << "Phase space dump frequency " << Options::psDumpFreq << " and "
469
          << "statistics dump frequency " << Options::statDumpFreq << " w.r.t. the time step." << endl;
470
#ifdef HAVE_AMR_SOLVER
471 472 473 474 475
    itsTracker = new ParallelTTracker(*Track::block->use->fetchLine(),
                                      dynamic_cast<PartBunch &>(*Track::block->bunch), *ds,
                                      Track::block->reference, false, false, Track::block->localTimeSteps,
                                      Track::block->zstop, Track::block->timeIntegrator, Track::block->dT,
                                      beam->getNumberOfParticles(),amrptr);
476
#else
477 478 479 480 481
    itsTracker = new ParallelTTracker(*Track::block->use->fetchLine(),
                                      dynamic_cast<PartBunch &>(*Track::block->bunch), *ds,
                                      Track::block->reference, false, false, Track::block->localTimeSteps,
                                      Track::block->zstop, Track::block->timeIntegrator, Track::block->dT,
                                      beam->getNumberOfParticles());
482
#endif
483 484
    itsTracker->setMpacflg(mpacflg); // set multipacting flag in ParallelTTracker
}
gsell's avatar
gsell committed
485

486 487 488
void TrackRun::setupCyclotronTracker(){
    OpalData::getInstance()->setInOPALCyclMode();
    Beam *beam = Beam::find(Attributes::getString(itsAttr[BEAM]));
gsell's avatar
gsell committed
489

490
    if (Attributes::getString(itsAttr[BOUNDARYGEOMETRY]) != "NONE") {
491
        // Ask the dictionary if BoundaryGeometry is allocated.
492
        // If it is allocated use the allocated BoundaryGeometry
493
        if (!OpalData::getInstance()->hasGlobalGeometry()) {
494
            BoundaryGeometry *bg = BoundaryGeometry::find(Attributes::getString(itsAttr[BOUNDARYGEOMETRY]))->
495
                clone(getOpalName() + std::string("_geometry"));
496 497
            OpalData::getInstance()->setGlobalGeometry(bg);
        }
498
    }
adelmann's avatar
adelmann committed
499

500
    setupFieldsolver();
gsell's avatar
gsell committed
501

kraus's avatar
kraus committed
502
    Track::block->bunch->PType = ParticleType::REGULAR;
gsell's avatar
gsell committed
503

504 505 506 507 508 509
    std::vector<std::string> distr_str = Attributes::getStringArray(itsAttr[DISTRIBUTION]);
    if (distr_str.size() == 0) {
        dist = Distribution::find(defaultDistribution);
    } else {
        dist = Distribution::find(distr_str.at(0));
    }
gsell's avatar
gsell committed
510

511 512 513
    // set macromass and charge for simulation particles
    double macromass = 0.0;
    double macrocharge = 0.0;
gsell's avatar
gsell committed
514

515
    const int specifiedNumBunch = int(std::abs(Round(Attributes::getReal(itsAttr[TURNS]))));
gsell's avatar
gsell committed
516

517 518 519 520 521
    if(OPAL->inRestartRun()) {
        phaseSpaceSink_m = new H5PartWrapperForPC(OPAL->getInputBasename() + std::string(".h5"),
                                                  OPAL->getRestartStep(),
                                                  OpalData::getInstance()->getRestartFileName(),
                                                  H5_O_WRONLY);
522 523 524 525 526
    } else if (OPAL->hasBunchAllocated() && !Options::scan) {
        phaseSpaceSink_m = new H5PartWrapperForPC(OPAL->getInputBasename() + std::string(".h5"),
                                                  -1,
                                                  OPAL->getInputBasename() + std::string(".h5"),
                                                  H5_O_WRONLY);
527 528 529 530 531
    } else {
        phaseSpaceSink_m = new H5PartWrapperForPC(OPAL->getInputBasename() + std::string(".h5"),
                                                  H5_O_WRONLY);
    }

532 533 534 535 536 537 538
    if(beam->getNumberOfParticles() < 3 || beam->getCurrent() == 0.0) {
        macrocharge = beam->getCharge() * q_e;
        macromass = beam->getMass();
        dist->CreateOpalCycl(*Track::block->bunch,
                             beam->getNumberOfParticles(),
                             beam->getCurrent(),*Track::block->use->fetchLine(),
                             Options::scan);
gsell's avatar
gsell committed
539

540
    } else {
gsell's avatar
gsell committed
541

542 543 544
        /**
           getFrequency() gets RF frequency [MHz], NOT isochronous  revolution frequency of particle!
           getCurrent() gets beamcurrent [A]
gsell's avatar
gsell committed
545

546 547
        */
        macrocharge = beam->getCurrent() / (beam->getFrequency() * 1.0e6); // [MHz]-->[Hz]
gsell's avatar
gsell committed
548

549 550
        if(!OPAL->hasBunchAllocated()) {
            if(!OPAL->inRestartRun()) {
gsell's avatar
gsell committed
551 552
                macrocharge /= beam->getNumberOfParticles();
                macromass = beam->getMass() * macrocharge / (beam->getCharge() * q_e);
553 554 555 556 557
                dist->CreateOpalCycl(*Track::block->bunch,
                                     beam->getNumberOfParticles(),
                                     beam->getCurrent(),
                                     *Track::block->use->fetchLine(),
                                     Options::scan);
558 559

            } else {
560 561 562 563 564
                dist->DoRestartOpalCycl(*Track::block->bunch,
                                        beam->getNumberOfParticles(),
                                        OPAL->getRestartStep(),
                                        specifiedNumBunch,
                                        phaseSpaceSink_m);
565 566
                macrocharge /= beam->getNumberOfParticles();
                macromass = beam->getMass() * macrocharge / (beam->getCharge() * q_e);
gsell's avatar
gsell committed
567
            }
568 569 570
        } else if(OPAL->hasBunchAllocated() && Options::scan) {
            macrocharge /= beam->getNumberOfParticles();
            macromass = beam->getMass() * macrocharge / (beam->getCharge() * q_e);
571 572 573 574 575
            dist->CreateOpalCycl(*Track::block->bunch,
                                 beam->getNumberOfParticles(),
                                 beam->getCurrent(),
                                 *Track::block->use->fetchLine(),
                                 Options::scan);
gsell's avatar
gsell committed
576
        }
577 578 579
    }
    Track::block->bunch->setMass(macromass); // set the Mass per macro-particle, [GeV/c^2]
    Track::block->bunch->setCharge(macrocharge);  // set the charge per macro-particle, [C]
gsell's avatar
gsell committed
580

581 582
    *gmsg << "* Mass of simulation particle= " << macromass << " GeV/c^2" << endl;
    *gmsg << "* Charge of simulation particle= " << macrocharge << " [C]" << endl;
gsell's avatar
gsell committed
583

584 585 586 587 588 589 590 591 592
    try {
        throw (beam->getNumberOfParticles() != Track::block->bunch->getTotalNum());
    }
    catch (bool notEqual) {
        if (notEqual && !OPAL->inRestartRun()) {
            throw OpalException("TrackRun::execute CYCLOTRON_T",
                                "Number of macro particles and NPART on BEAM are not equal");
        }
    }
gsell's avatar
gsell committed
593

594 595
    Track::block->bunch->setdT(1.0 / (Track::block->stepsPerTurn * beam->getFrequency() * 1.0e6));
    Track::block->bunch->setStepsPerTurn(Track::block->stepsPerTurn);
gsell's avatar
gsell committed
596

597 598 599
    // set coupling constant
    double coefE = 1.0 / (4 * pi * epsilon_0);
    Track::block->bunch->setCouplingConstant(coefE);
gsell's avatar
gsell committed
600

601 602 603 604 605
    // statistical data are calculated (rms, eps etc.)
    Track::block->bunch->calcBeamParameters_cycl();

    if(!OPAL->inRestartRun())
        if(!OPAL->hasDataSinkAllocated()) {
606
            ds = new DataSink(phaseSpaceSink_m);
gsell's avatar
gsell committed
607
            OPAL->setDataSink(ds);
608
        } else {
609
            ds = OPAL->getDataSink();
610 611
            ds->changeH5Wrapper(phaseSpaceSink_m);
        }
612
    else {
613
        ds = new DataSink(phaseSpaceSink_m, -1);
614 615
        OPAL->setDataSink(ds);
    }
616

617 618 619 620 621
    if(OPAL->hasBunchAllocated() && Options::scan)
        ds->reset();

    if(!OPAL->hasBunchAllocated() && !Options::scan) {
        *gmsg << "* ********************************************************************************** " << endl;
622
        *gmsg << "* Selected Tracking Method == CYCLOTRON-T, NEW TRACK" << endl;
623 624
        *gmsg << "* ********************************************************************************** " << endl;
    } else if(OPAL->hasBunchAllocated() && !Options::scan) {
gsell's avatar
gsell committed
625
        *gmsg << "* ********************************************************************************** " << endl;
626
        *gmsg << "* Selected Tracking Method == CYCLOTRON-T, FOLLOWUP TRACK" << endl;
627 628 629
        *gmsg << "* ********************************************************************************** " << endl;
    } else if(OPAL->hasBunchAllocated() && Options::scan) {
        *gmsg << "* ********************************************************************************** " << endl;
630
        *gmsg << "* Selected Tracking Method == CYCLOTRON-T, SCAN TRACK" << endl;
631 632 633 634 635 636 637 638
        *gmsg << "* ********************************************************************************** " << endl;
    }
    *gmsg << "* Number of neighbour bunches= " << specifiedNumBunch << endl;
    *gmsg << "* DT                         = " << Track::block->dT.front() << endl;
    *gmsg << "* MAXSTEPS                   = " << Track::block->localTimeSteps.front() << endl;
    *gmsg << "* Phase space dump frequency = " << Options::psDumpFreq << endl;
    *gmsg << "* Statistics dump frequency  = " << Options::statDumpFreq << " w.r.t. the time step." << endl;
    *gmsg << "* ********************************************************************************** " << endl;
gsell's avatar
gsell committed
639

640 641 642
    itsTracker = new ParallelCyclotronTracker(*Track::block->use->fetchLine(),
                                              dynamic_cast<PartBunch &>(*Track::block->bunch), *ds, Track::block->reference,
                                              false, false, Track::block->localTimeSteps.front(), Track::block->timeIntegrator);
gsell's avatar
gsell committed
643

644
    itsTracker->setNumBunch(specifiedNumBunch);
gsell's avatar
gsell committed
645

646
    if(OPAL->inRestartRun()) {
647

648 649
        H5PartWrapperForPC *h5pw = static_cast<H5PartWrapperForPC*>(phaseSpaceSink_m);
        itsTracker->setBeGa(h5pw->getMeanMomentum());
650

651 652 653
        itsTracker->setPr(h5pw->getReferencePr());
        itsTracker->setPt(h5pw->getReferencePt());
        itsTracker->setPz(h5pw->getReferencePz());
654

655 656 657
        itsTracker->setR(h5pw->getReferenceR());
        itsTracker->setTheta(h5pw->getReferenceT());
        itsTracker->setZ(h5pw->getReferenceZ());
658

659
        // The following is for restarts in local frame
660 661 662
        itsTracker->setPhi(h5pw->getAzimuth());
        itsTracker->setPsi(h5pw->getElevation());
        itsTracker->setPreviousH5Local(h5pw->getPreviousH5Local());
663
    }
664

665 666
    // statistical data are calculated (rms, eps etc.)
    Track::block->bunch->calcBeamParameters_cycl();
667

668 669 670 671
    *gmsg << *dist << endl;
    *gmsg << *beam << endl;
    *gmsg << *fs   << endl;
    // *gmsg << *Track::block->bunch  << endl;
672

673
    if(specifiedNumBunch > 1) {
gsell's avatar
gsell committed
674

675 676 677
        // only for regular  run of multi bunches, instantiate the  PartBins class
        // note that for restart run of multi bunches, PartBins class is instantiated in function doRestart_cycl()
        if(!OPAL->inRestartRun()) {
gsell's avatar
gsell committed
678

679 680 681 682
            // already exist bins number initially
            const int BinCount = 1;
            //allowed maximal bin
            const int MaxBinNum = 1000;
gsell's avatar
gsell committed
683

684 685 686 687
            // initialize particles number for each bin (both existed and not yet emmitted)
            size_t partInBin[MaxBinNum];
            for(int ii = 0; ii < MaxBinNum; ii++) partInBin[ii] = 0;
            partInBin[0] =  beam->getNumberOfParticles();
gsell's avatar
gsell committed
688

689 690 691
            Track::block->bunch->setPBins(new PartBinsCyc(MaxBinNum, BinCount, partInBin));
            // the allowed maximal bin number is set to 100
            //Track::block->bunch->setPBins(new PartBins(100));
gsell's avatar
gsell committed
692

693
        }
gsell's avatar
gsell committed
694

695 696 697 698 699
        // mode of generating new bunches:
        // "FORCE" means generating one bunch after each revolution, until get "TURNS" bunches.
        // "AUTO" means only when the distance between two neighbor bunches is bellow the limitation,
        //        then starts to generate new bunches after each revolution,until get "TURNS" bunches;
        //        otherwise, run single bunch track
gsell's avatar
gsell committed
700

701
        *gmsg << "***---------------------------- MULTI-BUNCHES MULTI-ENERGY-BINS MODE------ ----------------------------*** " << endl;
gsell's avatar
gsell committed
702

703 704
        double paraMb = Attributes::getReal(itsAttr[PARAMB]);
        itsTracker->setParaAutoMode(paraMb);
gsell's avatar
gsell committed
705

706
        if(OPAL->inRestartRun()) {
gsell's avatar
gsell committed
707

708
            itsTracker->setLastDumpedStep(OPAL->getRestartStep());
gsell's avatar
gsell committed
709

710 711 712
            if(Track::block->bunch->pbin_m->getLastemittedBin() < 2) {
                itsTracker->setMultiBunchMode(2);
                *gmsg << "In this restart job, the multi-bunches mode is forcely set to AUTO mode." << endl;
gsell's avatar
gsell committed
713
            } else {
714 715 716 717 718 719 720 721 722
                itsTracker->setMultiBunchMode(1);
                *gmsg << "In this restart job, the multi-bunches mode is forcely set to FORCE mode." << endl
                      << "If the existing bunch number is less than the specified number of TURN, readin the phase space of STEP#0 from h5 file consecutively" << endl;
            }
        } else {
            //////
            if((Attributes::getString(itsAttr[MBMODE])) == std::string("FORCE")) {
                itsTracker->setMultiBunchMode(1);
                *gmsg << "FORCE mode: The multi bunches will be injected consecutively after each revolution, until get \"TURNS\" bunches." << endl;
gsell's avatar
gsell committed
723 724 725


            }
726 727
            //////
            else if((Attributes::getString(itsAttr[MBMODE])) == std::string("AUTO")) {
gsell's avatar
gsell committed
728 729


730
                itsTracker->setMultiBunchMode(2);
gsell's avatar
gsell committed
731

732 733 734 735 736 737 738
                *gmsg << "AUTO mode: The multi bunches will be injected only when the distance between two neighborring bunches " << endl
                      << "is bellow the limitation. The control parameter is set to " << paraMb << endl;
            }
            //////
            else
                throw OpalException("TrackRun::execute()",
                                    "MBMODE name \"" + Attributes::getString(itsAttr[MBMODE]) + "\" unknown.");
gsell's avatar
gsell committed
739 740 741
        }

    }
742
}
gsell's avatar
gsell committed
743

744 745 746 747 748 749 750 751
void TrackRun::setupFieldsolver() {
    fs = FieldSolver::find(Attributes::getString(itsAttr[FIELDSOLVER]));
    fs->initCartesianFields();
    Track::block->bunch->setSolver(fs);
    if (fs->hasPeriodicZ())
        Track::block->bunch->setBCForDCBeam();
    else
        Track::block->bunch->setBCAllOpen();
gsell's avatar
gsell committed
752
}
753

754
double TrackRun::setDistributionParallelT(Beam *beam) {
755 756 757 758 759 760 761 762 763 764

    // If multipacting flag is not set, get distribution(s).
    if (!Attributes::getBool(itsAttr[MULTIPACTING])) {
        /*
         * Distribution(s) can be set via a single distribution or a list
         * (array) of distributions. If an array is defined the first in the
         * list is treated as the primary distribution. All others are added to
         * it to create the full distribution.
         */
        std::vector<std::string> distributionArray
765 766
            = Attributes::getStringArray(itsAttr[DISTRIBUTION]);
        const size_t numberOfDistributions = distributionArray.size();
767

768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792
        if (numberOfDistributions == 0) {
            dist = Distribution::find(defaultDistribution);
        } else {
            dist = Distribution::find(distributionArray.at(0));
            dist->setNumberOfDistributions(numberOfDistributions);

            if (numberOfDistributions > 1) {
                *gmsg << endl
                      << "---------------------------------" << endl
                      << "Found more than one distribution:" << endl << endl;
                *gmsg << "Main Distribution" << endl
                      << "---------------------------------" << endl
                      << distributionArray.at(0) << endl << endl
                      << "Secondary Distribution(s)" << endl
                      << "---------------------------------" << endl;

                for (size_t i = 1; i < numberOfDistributions; ++ i) {
                    Distribution *distribution = Distribution::find(distributionArray.at(i));
                    distribution->setNumberOfDistributions(numberOfDistributions);
                    distrs_m.push_back(distribution);

                    *gmsg << distributionArray.at(i) << endl;
                }
                *gmsg << endl
                      << "---------------------------------" << endl << endl;
793
            }
794
        }
795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826
    }

    /*
     * Initialize distributions.
     */
    size_t numberOfParticles = beam->getNumberOfParticles();
    if (!OPAL->hasBunchAllocated()) {
        if (!OPAL->inRestartRun()) {
            if (!Attributes::getBool(itsAttr[MULTIPACTING])) {
                /*
                 * Here we are not doing a restart or doing a mulitpactor run
                 * and we do not have a bunch already allocated.
                 */
                Track::block->bunch->setDistribution(dist,
                                                     distrs_m,
                                                     numberOfParticles,
                                                     Options::scan);

                /*
                 * If this is an injected beam (rather than an emitted beam), we
                 * make sure it doesn't have any particles at z < 0.
                 */
                Vector_t rMin;
                Vector_t rMax;
                Track::block->bunch->get_bounds(rMin, rMax);

                OPAL->setGlobalPhaseShift(0.5 * dist->GetTEmission());
            }
        } else {
            /*
             * Read in beam from restart file.
             */
827
            dist->DoRestartOpalT(*Track::block->bunch, numberOfParticles, OPAL->getRestartStep(), phaseSpaceSink_m);
828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848
        }
    } else if (OPAL->hasBunchAllocated() && Options::scan) {
        /*
         * We are in scan mode and have already allocated a bunch. So, we need to
         * do a couple more things.
         */
        Track::block->bunch->setDistribution(dist,
                                             distrs_m,
                                             numberOfParticles,
                                             Options::scan);
        Track::block->bunch->resetIfScan();
        Track::block->bunch->LastSection = 1;

        OPAL->setGlobalPhaseShift(0.5 * dist->GetTEmission());
    }

    // Return charge per macroparticle.
    return beam->getCharge() * beam->getCurrent() / beam->getFrequency() / numberOfParticles;

}

849
void TrackRun::findPhasesForMaxEnergy(bool writeToFile) const {
850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903
    if (Options::autoPhase == 0 ||
        OpalData::getInstance()->inRestartRun() ||
        OpalData::getInstance()->hasBunchAllocated()) return;

    std::queue<unsigned long long> localTrackSteps;
    std::queue<double> dtAllTracks;
    std::queue<double> zStop;

    const PartBunch *bunch = Track::block->bunch;
    if (bunch == NULL) {
        bunch = Track::block->slbunch;
    }
    if (bunch == NULL) {
        throw OpalException("findPhasesForMaxEnergy()",
                            "no valid PartBunch object found");
    }

    Vector_t meanP = bunch->get_pmean();
    Vector_t meanR = bunch->get_rmean();
    if (bunch->getTotalNum() == 0) {
        meanP = bunch->get_pmean_Distribution();
        meanR = 0.0;
    }
    AutophaseTracker apTracker(*Track::block->use->fetchLine(),
                               Track::block->reference,
                               bunch->getT(),
                               meanR(2),
                               meanP(2));

    {
        std::vector<unsigned long long>::const_iterator it = Track::block->localTimeSteps.begin();
        std::vector<unsigned long long>::const_iterator end = Track::block->localTimeSteps.end();
        for (; it != end; ++ it) {
            localTrackSteps.push(*it);
        }
    }

    {
        std::vector<double>::const_iterator it = Track::block->dT.begin();
        std::vector<double>::const_iterator end = Track::block->dT.end();
        for (; it != end; ++ it) {
            dtAllTracks.push(*it);
        }
    }

    {
        std::vector<double>::const_iterator it = Track::block->zstop.begin();
        std::vector<double>::const_iterator end = Track::block->zstop.end();
        for (; it != end; ++ it) {
            zStop.push(*it);
        }
    }

    apTracker.execute(dtAllTracks, zStop, localTrackSteps);
904

905 906 907
    if (writeToFile) {
        std::string fileName = OpalData::getInstance()->getInputBasename() + "_phases.in";
        apTracker.save(fileName);
908
    }
909
}
910

911
void TrackRun::executeAutophaseTracker() {
912

913
    Beam *beam = Beam::find(Attributes::getString(itsAttr[BEAM]));
914
    /*double charge =*/ setDistributionParallelT(beam);
915

916
    Track::block->bunch->setdT(Track::block->dT.front());
917 918 919
    Track::block->bunch->dtScInit_m = Track::block->dtScInit;
    Track::block->bunch->deltaTau_m = Track::block->deltaTau;
    Track::block->bunch->setT(Track::block->t0_m);
920

921
    // Track::block->bunch->setCharge(charge);
922

923 924
    double couplingConstant = 1.0 / (4 * pi * epsilon_0);
    Track::block->bunch->setCouplingConstant(couplingConstant);
925 926


927 928 929
    Track::block->bunch->calcBeamParameters();

    findPhasesForMaxEnergy(true);
930

931 932 933 934 935 936 937 938 939 940 941 942 943
}

#ifdef HAVE_AMR_SOLVER
std::vector<std::string> TrackRun::filterString(std::string str) {
  /*
    pid tokens[0]
    il == tokens[2], ih == tokens[3]
    jl == tokens[5], jh == tokens[6]
    kl == tokens[9], kh == tokens[9]
   */

  // charakters to remove from the string
  char chars[] = "Node=;vn_mDmain{[][][]}:,";
944 945

  for (unsigned int i = 0; i < strlen(chars); ++i)
946 947 948
    std::replace(str.begin(), str.end(), chars[i], ' ');

  std::vector<std::string> tokens;
949

950 951 952 953 954 955 956 957 958 959 960 961
  // filter spaces
  std::istringstream iss(str);
  copy(std::istream_iterator<std::string>(iss),
       std::istream_iterator<std::string>(),
       std::back_inserter<std::vector<std::string> >(tokens));
  return tokens;
}

std::pair<Box,unsigned int> TrackRun::getBlGrids(std::string str){

  std::vector<std::string> tokens = filterString(str);

962
  unsigned int theGrid;
963 964 965 966 967 968 969 970 971 972 973 974 975 976
  std::istringstream (tokens[0]) >> theGrid;

  int ilo,ihi,jlo,jhi,klo,khi;

  std::istringstream (tokens[2]) >> ilo;
  std::istringstream (tokens[5]) >> jlo;
  std::istringstream (tokens[8]) >> klo;
  std::istringstream (tokens[3]) >> ihi;
  std::istringstream (tokens[6]) >> jhi;
  std::istringstream (tokens[9]) >> khi;

  Inform m2a("AMR ",INFORM_ALL_NODES);
  /*
  m2a << "Grid " << tokens[0]
977 978
      << " i (" << tokens[2] << " ... " << tokens[3] << ")"
      << " j (" << tokens[5] << " ... " << tokens[6] << ")"
979 980 981 982 983 984 985
      << " k (" << tokens[8] << " ... " << tokens[9] << ")" << " myNode " << Ippl::myNode() << endl;
  */
  IntVect loEnd(ilo,jlo,klo);
  IntVect hiEnd(ihi,jhi,khi);
  Box bx(loEnd,hiEnd);

  return std::pair<Box,unsigned int>(bx,theGrid);
986
}
987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 1138 1139 1140

Amr* TrackRun::setupAMRSolver() {
    Amr* amrptr;
    //	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();

    NDIndex<3> ipplDom = Track::block->bunch->getFieldLayout().getDomain();

    BoxArray lev0_grids(Ippl::getNodes());

    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) {
        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) {
        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];
        }

    // 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);

    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);

    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);

    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);

            //  This allocates 11 spots -- 1 for Q, 3 for v, 3 for E, 3 for B, 1 for ID.
            //  IMPORTANT: THIS ORDERING IS ASSUMED WHEN WE FILL E AT THE PARTICLE LOCATION
            //             IN THE MOVEKICK CALL -- if Evec no longer starts at component 4
            //             then you must change "start_comp_for_e" in Accel_advance.cpp
            /*
              Q      : 0
              Vvec   : 1, 2, 3 the velocity
              Evec   : 4, 5, 6 the electric field at the particle location
              Bvec   : 7, 8, 9 the electric field at the particle location
              id+1   : 10 (we add 1 to make the particle ID > 0)
            */

            // This is the charge
            attr[0] = Track::block->bunch->Q[i];

            // These are the velocity components
            double gamma=sqrt(1+ dot(Track::block->bunch->P[i],Track::block->bunch->P[i]));
            for (unsigned int k=0; k<DIM; k++)
                attr[k+1] = Track::block->bunch->P[i](k) * Physics::c /gamma;

            // These are E and B
            for (unsigned int k=4; k<10; k++)
                attr[k]= 0.0;

            //
            // The Particle stuff in AMR requires ids > 0
            //   (because we flip the sign to make them invalid)
            // So we just make id->id+1 here.
            int particle_id = Track::block->bunch->ID[i] + 1;
            attr[3*DIM + 1] = particle_id;

            amrptr->addOneParticle(particle_id, Ippl::myNode(), x, attr);
        }

    // It is essential that we call this routine since the particles
    //    may not currently be defined on the same processor as the grid
    //    that will hold them in the AMR world.
    amrptr->RedistributeParticles();

    // This part of the call must come after we add the particles
    // since this one calls post_init which does the field solve.
    amrptr->FinalizeInit(strt_time, stop_time);

    amrptr->writePlotFile();

    *gmsg << "A M R Initialization DONE" << endl;

    return amrptr;
}

1141
#endif