Commit 97f26824 authored by kraus's avatar kraus

standalone Autophase Tracker

parent e45cd295
......@@ -897,7 +897,7 @@ double RFCavity::getAutoPhaseEstimate(const double &E0, const double &t0, const
}
}
const int prevPrecision = Ippl::Info->precision(8);
INFOMSG("estimated phi= " << tmp_phi << " rad, "
INFOMSG(level2 << "estimated phase= " << tmp_phi << " rad, "
<< "Ekin= " << E[N - 1] << " MeV" << setprecision(prevPrecision) << endl);
return tmp_phi;
......@@ -952,7 +952,7 @@ double RFCavity::getAutoPhaseEstimate(const double &E0, const double &t0, const
}
const int prevPrecision = Ippl::Info->precision(8);
INFOMSG("estimated phi= " << tmp_phi << " rad, "
INFOMSG(level2 << "estimated phase= " << tmp_phi << " rad, "
<< "Ekin= " << E[N - 1] << " MeV" << setprecision(prevPrecision) << endl);
return phi;
......
......@@ -583,43 +583,8 @@ double TravelingWave::getAutoPhaseEstimate(const double &E0, const double &t0, c
}
const int prevPrecision = Ippl::Info->precision(8);
INFOMSG("estimated phi= " << tmp_phi << " rad, "
INFOMSG(level2 << "estimated phase= " << tmp_phi << " rad, "
<< "Ekin= " << E[N3 - 1] << " MeV" << setprecision(prevPrecision) << endl);
// std::stringstream fn;
// fn << getName() << ".dat";
// ofstream ckrtest(fn.str().c_str());
// ckrtest << setprecision(12);
// for (int i = 0; i < N1; ++ i) {
// ckrtest << startField_m + F[i].first << "\t"
// << t[i] << "\t"
// << scale_m * cos(frequency_m * t[i] + phi) * F[i].second << "\t"
// << E[i] << "\t"
// << cos(frequency_m * t[i] + phi) << "\t"
// << i << endl;
// }
// for (int i = N1; i < N3 - N1 + 1; ++ i) {
// int I = (i - N1) % N2 + N1;
// int J = (i - N1 + N4) % N2 + N1;
// double Z = startField_m + F[I].first + floor((i - N1) / N2) * Dz;
// ckrtest << Z << "\t"
// << t[i] << "\t"
// << scaleCore_m * (cos(frequency_m * t[i] + phaseC1 + phi) * F[I].second +
// cos(frequency_m * t[i] + phaseC2 + phi) * F[J].second) << "\t"
// << E[i] << "\t"
// << frequency_m * t[i] + phi << "\t"
// << I << endl;
// }
// for (int i = N3 - N1 + 1; i < N3; ++ i) {
// int I = i - N3 - 1 + 2 * N1 + N2;
// double Z = startField_m + F[I].first + ((NumCells_m - 1) * Mode_m - 1) * Dz;
// ckrtest << Z << "\t"
// << t[i] << "\t"
// << scale_m * cos(frequency_m * t[i] + phaseE + phi) * F[I].second << "\t"
// << E[i] << "\t"
// << frequency_m * t[i] + phi << "\t"
// << I << endl;
// }
// ckrtest.close();
return tmp_phi;
}
phi = tmp_phi - floor(tmp_phi / Physics::two_pi + 0.5) * Physics::two_pi;
......@@ -657,7 +622,7 @@ double TravelingWave::getAutoPhaseEstimate(const double &E0, const double &t0, c
const int prevPrecision = Ippl::Info->precision(8);
INFOMSG("estimated phi= " << tmp_phi << " rad, "
INFOMSG(level2 << "estimated phase= " << tmp_phi << " rad, "
<< "Ekin= " << E[N3 - 1] << " MeV" << setprecision(prevPrecision) << endl);
return phi;
......
......@@ -4,6 +4,9 @@
#include "Utilities/OpalOptions.h"
#include "Utilities/Timer.h"
#include <iostream>
#include <fstream>
#define DTFRACTION 2
extern Inform *gmsg;
......@@ -121,7 +124,7 @@ void AutophaseTracker::execute(const std::queue<double> &dtAllTracks,
<< "step " << step << ", "
<< "t = " << itsBunch_m.getT() << " [s], "
<< "E = " << getEnergyMeV(refP) << " [MeV]\n\n"
<< "start phase scan ... \n");
<< "start phase scan ... " << endl;);
double guess = guessCavityPhase(next);
optimizeCavityPhase(next,
......@@ -129,7 +132,7 @@ void AutophaseTracker::execute(const std::queue<double> &dtAllTracks,
step,
newDt);
INFOMSG(level2 << "\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n");
INFOMSG(level2 << "\n%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << endl);
track(endNext,
step,
......@@ -147,9 +150,10 @@ void AutophaseTracker::execute(const std::queue<double> &dtAllTracks,
maxStepsAutoPhasing.pop();
}
printCavityPhases();
sendCavityPhases();
*gmsg << level1 << "\n\nfinished autophasing" << clock.time() << "\n\n" << endl;
*gmsg << level1 << "\n\nfinished autophasing " << clock.time() << "\n\n" << endl;
}
void AutophaseTracker::track(double uptoSPos,
......@@ -199,8 +203,9 @@ void AutophaseTracker::track(double uptoSPos,
maxStepsAutoPhasing.front() = step;
if(step % 5000 == 0) {
INFOMSG("step = " << step << ", spos = " << refR(2) << " [m], t= " << itsBunch_m.getT() << " [s], "
<< "E = " << getEnergyMeV(refP) << " [MeV] " << endl);
INFOMSG(level2 << "step = " << step << ", spos = " << refR(2) << " [m], "
<< "t= " << itsBunch_m.getT() << " [s], " << "E = " << getEnergyMeV(refP) << " [MeV] "
<< endl);
}
}
......@@ -255,12 +260,12 @@ double AutophaseTracker::guessCavityPhase(const std::shared_ptr<Component> &cavi
orig_phi = element->getPhasem();
apVeto = element->getAutophaseVeto();
if(apVeto) {
INFOMSG(" ----> APVETO -----> "
INFOMSG(level2 << " ----> APVETO -----> "
<< element->getName() << endl);
return orig_phi;
}
INFOMSG(cavity->getName() << ", "
<< "phi = " << orig_phi << endl);
INFOMSG(level2 << cavity->getName() << ", "
<< "phase = " << orig_phi << " rad" << endl);
Phimax = element->getAutoPhaseEstimate(getEnergyMeV(refP),
itsBunch_m.getT() + tErr,
......@@ -271,12 +276,12 @@ double AutophaseTracker::guessCavityPhase(const std::shared_ptr<Component> &cavi
orig_phi = element->getPhasem();
apVeto = element->getAutophaseVeto();
if(apVeto) {
INFOMSG(" ----> APVETO -----> "
INFOMSG(level2 << " ----> APVETO -----> "
<< element->getName() << endl);
return orig_phi;
}
INFOMSG(cavity->getName() << ", "
<< "phi = " << orig_phi << endl);
INFOMSG(level2 << cavity->getName() << ", "
<< "phase = " << orig_phi << " rad" << endl);
Phimax = element->getAutoPhaseEstimate(getEnergyMeV(refP),
itsBunch_m.getT() + tErr,
......@@ -469,4 +474,71 @@ void AutophaseTracker::receiveCavityPhases() {
}
delete mess;
}
void AutophaseTracker::printCavityPhases() {
std::shared_ptr<Component> cavity = NULL;
const double RADDEG = 180.0 / Physics::pi;
const double globalTimeShift = OpalData::getInstance()->getGlobalPhaseShift();
*gmsg << level1 << "\n-------------------------------------------------------------------------------------\n";
while (true) {
std::shared_ptr<Component> next = getNextCavity(cavity);
if (next == NULL) break;
std::string name = next->getName();
double frequency;
double phase;
if (next->getType() == ElementBase::TRAVELINGWAVE) {
phase = static_cast<TravelingWave *>(next.get())->getPhasem();
frequency = static_cast<TravelingWave *>(next.get())->getFrequencym();
} else {
phase = static_cast<RFCavity *>(next.get())->getPhasem();
frequency = static_cast<RFCavity *>(next.get())->getFrequencym();
}
phase -= globalTimeShift * frequency;
*gmsg << (cavity == NULL? "": "\n")
<< name
<< ": phi = phi_nom + phi_maxE + global phase shift = " << phase * RADDEG << " degree, "
<< "(global phase shift = " << -globalTimeShift *frequency *RADDEG << " degree) \n";
cavity = next;
}
*gmsg << "-------------------------------------------------------------------------------------\n"
<< endl;
}
void AutophaseTracker::save(const std::string &fname) {
std::ofstream out(fname);
std::shared_ptr<Component> cavity = NULL;
const double globalTimeShift = OpalData::getInstance()->getGlobalPhaseShift();
while (true) {
std::shared_ptr<Component> next = getNextCavity(cavity);
if (next == NULL) break;
std::string name = next->getName();
double frequency;
double phase;
if (next->getType() == ElementBase::TRAVELINGWAVE) {
phase = static_cast<TravelingWave *>(next.get())->getPhasem();
frequency = static_cast<TravelingWave *>(next.get())->getFrequencym();
} else {
phase = static_cast<RFCavity *>(next.get())->getPhasem();
frequency = static_cast<RFCavity *>(next.get())->getFrequencym();
}
phase -= globalTimeShift * frequency;
out << name << "_lag = " << phase << ";\n";
cavity = next;
}
out.close();
}
\ No newline at end of file
......@@ -61,6 +61,8 @@ public:
const std::queue<double> &maxZ,
const std::queue<unsigned long long> &maxTrackSteps);
void save(const std::string &fname);
virtual void visitBeamline(const Beamline &bl);
AP_VISITELEMENT(AlignWrapper)
......@@ -112,6 +114,7 @@ private:
double getEndCavity(const std::shared_ptr<Component> &);
void sendCavityPhases();
void receiveCavityPhases();
void printCavityPhases();
OpalBeamline itsOpalBeamline_m;
FieldList::iterator currentAPCavity_m;
......
......@@ -190,7 +190,7 @@ SeyNum_m(0.0),
timeIntegrationTimer1Loop1_m(IpplTimings::getTimer("TIntegration1Loop1")),
timeIntegrationTimer1Loop2_m(IpplTimings::getTimer("TIntegration1Loop2")),
timeIntegrationTimer2Loop1_m(IpplTimings::getTimer("TIntegration2Loop1")),
timeIntegrationTimer2Loop2_m(IpplTimings::getTimer("TIntegration2Loop2"))
timeIntegrationTimer2Loop2_m(IpplTimings::getTimer("TIntegration2Loop2"))
{
for (std::vector<unsigned long long>::const_iterator it = maxSteps.begin(); it != maxSteps.end(); ++ it) {
......@@ -338,42 +338,6 @@ void ParallelTTracker::updateRFElement(std::string elName, double maxPhase) {
}
}
void ParallelTTracker::printRFPhases() {
/**
All RF-Elements gets updated, where the phiShift is the
global phase shift in units of seconds.
*/
FieldList &cl = cavities_m;
const double RADDEG = 180.0 / Physics::pi;
const double globalTimeShift = OpalData::getInstance()->getGlobalPhaseShift();
*gmsg << level1 << "\n-------------------------------------------------------------------------------------\n";
for (FieldList::iterator it = cl.begin(); it != cl.end(); ++it) {
std::shared_ptr<Component> element(it->getElement());
std::string name = element->getName();
double frequency;
double phase;
if (element->getType() == ElementBase::TRAVELINGWAVE) {
phase = static_cast<TravelingWave *>(element.get())->getPhasem();
frequency = static_cast<TravelingWave *>(element.get())->getFrequencym();
} else {
phase = static_cast<RFCavity *>(element.get())->getPhasem();
frequency = static_cast<RFCavity *>(element.get())->getFrequencym();
}
*gmsg << (it == cl.begin()? "": "\n")
<< name
<< ": phi = phi_nom + phi_maxE + global phase shift = " << phase * RADDEG << " degree, "
<< "(global phase shift = " << -globalTimeShift *frequency *RADDEG << " degree) \n";
}
*gmsg << "-------------------------------------------------------------------------------------\n"
<< endl;
}
void ParallelTTracker::handleAutoPhasing() {
typedef std::vector<MaxPhasesT>::iterator iterator_t;
......@@ -388,8 +352,6 @@ void ParallelTTracker::handleAutoPhasing() {
for(; it < end; ++ it) {
updateRFElement((*it).first, (*it).second);
}
printRFPhases();
}
void ParallelTTracker::execute() {
......@@ -473,11 +435,11 @@ void ParallelTTracker::executeDefaultTracker() {
//get number of elements in the bunch
numDeviceElements = itsBunch->getLocalNum();
//allocate memory on device
//allocate memory on device
r_ptr = dksbase.allocateMemory<Vector_t>(numDeviceElements, ierr);
p_ptr = dksbase.allocateMemory<Vector_t>(numDeviceElements, ierr);
x_ptr = dksbase.allocateMemory<Vector_t>(numDeviceElements, ierr);
lastSec_ptr = dksbase.allocateMemory<long>(numDeviceElements, ierr);
dt_ptr = dksbase.allocateMemory<double>(numDeviceElements, ierr);
orient_ptr = dksbase.allocateMemory<Vector_t>(itsOpalBeamline_m.sections_m.size(), ierr);
......@@ -491,7 +453,7 @@ void ParallelTTracker::executeDefaultTracker() {
//write orientations to device
dksbase.writeData<Vector_t>(orient_ptr, orientation, nsec);
//page lock itsBunch->X, itsBunch->R, itsBunch-P
dksbase.registerHostMemory(&itsBunch->R[0], numDeviceElements);
dksbase.registerHostMemory(&itsBunch->P[0], numDeviceElements);
......@@ -1605,16 +1567,16 @@ void ParallelTTracker::timeIntegration1(BorisPusher & pusher) {
dksbase.writeDataAsync<Vector_t>(x_ptr, &itsBunch->X[0], itsBunch->getLocalNum(), stream2);
}
//write LastSection to device
dksbase.writeDataAsync<long>(lastSec_ptr, &itsBunch->LastSection[0], itsBunch->getLocalNum(),
dksbase.writeDataAsync<long>(lastSec_ptr, &itsBunch->LastSection[0], itsBunch->getLocalNum(),
stream2);
//calc push
dksbase.callParallelTTrackerPushTransform(x_ptr, p_ptr, lastSec_ptr, orient_ptr,
itsBunch->getLocalNum(),
dksbase.callParallelTTrackerPushTransform(x_ptr, p_ptr, lastSec_ptr, orient_ptr,
itsBunch->getLocalNum(),
itsOpalBeamline_m.sections_m.size(),
NULL, itsBunch->getdT(), Physics::c, false, stream2);
//read R from device
dksbase.readDataAsync<Vector_t>(x_ptr, &itsBunch->X[0], itsBunch->getLocalNum(), stream2);
//sync and wait till all tasks and reads are complete
dksbase.syncDevice();
#else
......@@ -1732,7 +1694,7 @@ IpplTimings::startTimer(timeIntegrationTimer2Loop1_m);
//wrote dt to device
dksbase.writeDataAsync<double>(dt_ptr, &itsBunch->dt[0], itsBunch->getLocalNum(), stream1);
//calc push
dksbase.callParallelTTrackerPush(r_ptr, p_ptr, itsBunch->getLocalNum(), dt_ptr,
dksbase.callParallelTTrackerPush(r_ptr, p_ptr, itsBunch->getLocalNum(), dt_ptr,
itsBunch->getdT(), Physics::c, true, stream1);
//read R from device
dksbase.readDataAsync<Vector_t>(r_ptr, &itsBunch->R[0], itsBunch->getLocalNum(), stream1);
......@@ -1742,11 +1704,11 @@ IpplTimings::startTimer(timeIntegrationTimer2Loop1_m);
dksbase.writeDataAsync<Vector_t>(x_ptr, &itsBunch->X[0], itsBunch->getLocalNum(), stream2);
}
//write LastSection to device
dksbase.writeDataAsync<long>(lastSec_ptr, &itsBunch->LastSection[0],
dksbase.writeDataAsync<long>(lastSec_ptr, &itsBunch->LastSection[0],
itsBunch->getLocalNum(), stream2);
//calc push
dksbase.callParallelTTrackerPushTransform(x_ptr, p_ptr, lastSec_ptr, orient_ptr,
itsBunch->getLocalNum(),
dksbase.callParallelTTrackerPushTransform(x_ptr, p_ptr, lastSec_ptr, orient_ptr,
itsBunch->getLocalNum(),
itsOpalBeamline_m.sections_m.size(),
dt_ptr, itsBunch->getdT(), Physics::c, true, stream2);
//read R from device
......@@ -2100,7 +2062,7 @@ void ParallelTTracker::computeExternalFields() {
//assign new device memory size
numDeviceElements = itsBunch->getLocalNum();
//register new pagelocked memory
dksbase.registerHostMemory(&itsBunch->R[0], numDeviceElements);
dksbase.registerHostMemory(&itsBunch->P[0], numDeviceElements);
......@@ -2108,7 +2070,7 @@ void ParallelTTracker::computeExternalFields() {
dksbase.registerHostMemory(&itsBunch->dt[0], numDeviceElements);
dksbase.registerHostMemory(&itsBunch->LastSection[0], numDeviceElements);
//allocate memory on device
//allocate memory on device
r_ptr = dksbase.allocateMemory<Vector_t>(numDeviceElements, ierr);
p_ptr = dksbase.allocateMemory<Vector_t>(numDeviceElements, ierr);
x_ptr = dksbase.allocateMemory<Vector_t>(numDeviceElements, ierr);
......@@ -2828,4 +2790,4 @@ void ParallelTTracker::writePhaseSpace(const long long step, const double &sposR
// mode:c
// c-basic-offset: 4
// indent-tabs-mode:nil
// End:
// End:
\ No newline at end of file
......@@ -151,10 +151,8 @@ int main(int argc, char *argv[]) {
// DTA
std::cout.precision(16);
std::cout.setf(std::ios::scientific, std::ios::floatfield);
std::cout.setf(std::ios::showpos);
std::cerr.precision(16);
std::cerr.setf(std::ios::scientific, std::ios::floatfield);
std::cerr.setf(std::ios::showpos);
// /DTA
// Set global truncation orders.
......
......@@ -66,8 +66,6 @@ SurfacePhysics::SurfacePhysics():
itsAttr[NPART] = Attributes::makeReal("NPART", "Number of particles in bunch");
std::cout << "Debug ==> Number of particles: " << Attributes::getReal(itsAttr[NPART]) << std::endl;
SurfacePhysics *defSurfacePhysics = clone("UNNAMED_SURFACEPHYSICS");
defSurfacePhysics->builtin = true;
......@@ -154,4 +152,4 @@ void SurfacePhysics::print(std::ostream &os) const {
<< "* SIGMA " << Attributes::getReal(itsAttr[SIGMA]) << '\n'
<< "* TAU " << Attributes::getReal(itsAttr[TAU]) << '\n';
os << "* ********************************************************************************** " << std::endl;
}
}
\ No newline at end of file
......@@ -173,6 +173,10 @@ void TrackRun::execute() {
} else if(method == "CYCLOTRON-T") {
setupCyclotronTracker();
} else if(method == "AUTOPHASE") {
executeAutophaseTracker();
return;
} else {
throw OpalException("TrackRun::execute()",
"Method name \"" + method + "\" unknown.");
......@@ -826,7 +830,7 @@ double TrackRun::setDistributionParallelT(Beam *beam) {
}
void TrackRun::findPhasesForMaxEnergy() const {
void TrackRun::findPhasesForMaxEnergy(bool writeToFile) const {
if (Options::autoPhase == 0 ||
OpalData::getInstance()->inRestartRun() ||
OpalData::getInstance()->hasBunchAllocated()) return;
......@@ -881,61 +885,33 @@ void TrackRun::findPhasesForMaxEnergy() const {
}
apTracker.execute(dtAllTracks, zStop, localTrackSteps);
}
ParallelTTracker *TrackRun::setupForAutophase() {
Inform m("setupForAutophase() ");
Beam *beam = Beam::find(Attributes::getString(itsAttr[BEAM]));
DataSink *ds = NULL;
fs = FieldSolver::find(Attributes::getString(itsAttr[FIELDSOLVER]));
fs->initCartesianFields();
Track::block->bunch->setSolver(fs);
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));
if (writeToFile) {
std::string fileName = OpalData::getInstance()->getInputBasename() + "_phases.in";
apTracker.save(fileName);
}
}
double charge = beam->getCharge() * beam->getCurrent() / beam->getFrequency();
void TrackRun::executeAutophaseTracker() {
charge /= beam->getNumberOfParticles();
Beam *beam = Beam::find(Attributes::getString(itsAttr[BEAM]));
double charge = setDistributionParallelT(beam);
Track::block->bunch->setdT(Track::block->dT.front());
Track::block->bunch->dtScInit_m = Track::block->dtScInit;
Track::block->bunch->deltaTau_m = Track::block->deltaTau;
Track::block->bunch->setT(Track::block->t0_m);
Track::block->bunch->resetIfScan();
Track::block->bunch->setCharge(charge);
Track::block->bunch->LastSection = 1;
double couplingConstant = 1.0 / (4 * pi * epsilon_0);
Track::block->bunch->setCouplingConstant(couplingConstant);
Track::block->bunch->setCharge(charge);
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->calcBeamParameters();
findPhasesForMaxEnergy(true);
double coefE = 1.0 / (4 * pi * epsilon_0);
Track::block->bunch->setCouplingConstant(coefE);
size_t N = 1;
#ifdef HAVE_AMR_SOLVER
Amr* null_amrptr = 0;
return 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, N, null_amrptr);
#else
return 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, N);
#endif
}
#ifdef HAVE_AMR_SOLVER
......
......@@ -77,8 +77,8 @@ private:
void setupFieldsolver();
double setDistributionParallelT(Beam *beam);
void findPhasesForMaxEnergy() const;
ParallelTTracker *setupForAutophase();
void findPhasesForMaxEnergy(bool writeToFile = false) const;
void executeAutophaseTracker();
// Pointer to tracking algorithm.
Tracker *itsTracker;
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment