Commit 7ac870e1 authored by winklehner_d's avatar winklehner_d

Added spiral inflector option to cyclotron. Changed some things in SAAMG solver. -DW

parent 7de72700
......@@ -58,6 +58,7 @@ Cyclotron::Cyclotron(const Cyclotron &right):
phiinit_m(right.phiinit_m),
zinit_m(right.zinit_m),
pzinit_m(right.pzinit_m),
spiral_flag_m(right.spiral_flag_m),
type_m(right.type_m),
harm_m(right.harm_m),
bscale_m(right.bscale_m),
......@@ -129,6 +130,14 @@ double Cyclotron::getPZinit() const {
return pzinit_m;
}
void Cyclotron::setSpiralFlag(bool spiral_flag) {
spiral_flag_m = spiral_flag;
}
bool Cyclotron::getSpiralFlag() const {
return spiral_flag_m;
}
void Cyclotron::setFieldMapFN(std::string f) {
fmapfn_m = f;
}
......@@ -1351,4 +1360,4 @@ void Cyclotron::getFieldFromFile_BandRF(const double &scaleFactor) {
void Cyclotron::getDimensions(double &zBegin, double &zEnd) const
{ }
#undef CHECK_CYC_FSCANF_EOF
\ No newline at end of file
#undef CHECK_CYC_FSCANF_EOF
......@@ -198,6 +198,9 @@ public:
void setFMHighE(double e);
virtual double getFMHighE() const;
void setSpiralFlag(bool spiral_flag);
virtual bool getSpiralFlag() const;
virtual bool apply(const size_t &id, const double &t, double E[], double B[]);
virtual bool apply(const size_t &id, const double &t, Vector_t &E, Vector_t &B);
......@@ -232,6 +235,8 @@ private:
double zinit_m;
double pzinit_m;
bool spiral_flag_m;
std::string type_m; /* what type of field we use */
double harm_m;
......@@ -287,4 +292,4 @@ private:
};
#endif // CLASSIC_Cyclotron_HH
\ No newline at end of file
#endif // CLASSIC_Cyclotron_HH
......@@ -373,7 +373,13 @@ double PartBunch::getLaserEnergy() const {
return 0.0;
}
/// \brief Return the fieldsolver type if we have a fieldsolver
std::string PartBunch::getFieldSolverType() const {
if(fs_m)
return fs_m->getFieldSolverType();
else
return "";
}
/** \brief After each Schottky scan we delete all the particles.
......@@ -1052,7 +1058,6 @@ void PartBunch::computeSelfFields() {
if(fs_m->getFieldSolverType() == "SAAMG")
resizeMesh();
INFOMSG(level3 << "mesh size" << hr_m << endl);
//scatter charges onto grid
this->Q *= this->dt;
this->Q.scatter(this->rho_m, this->R, IntrplCIC_t());
......@@ -1321,6 +1326,7 @@ void PartBunch::computeSelfFields_cycl(double gamma) {
/// mesh the whole domain
if(fs_m->getFieldSolverType() == "SAAMG")
resizeMesh();
/// scatter particles charge onto grid.
this->Q.scatter(this->rho_m, this->R, IntrplCIC_t());
......
......@@ -385,6 +385,8 @@ public:
bool hasFieldSolver();
std::string getFieldSolverType() const;
void setLPath(double s);
double getLPath() const;
......
......@@ -354,6 +354,28 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) {
Cyclotron *elptr = dynamic_cast<Cyclotron *>(cycl.clone());
myElements.push_back(elptr);
// Is this a Spiral Inflector Simulation? If yes, we'll give the user some
// useful information
spiral_flag = elptr->getSpiralFlag();
if(spiral_flag) {
*gmsg << endl << "* This is a Spiral Inflector Simulation! This means the following:" << endl;
*gmsg << "* 1.) It is up to the user to provide appropriate geometry, electric and magnetic fields!" << endl;
*gmsg << "* (Use BANDRF type cyclotron and use RFMAPFN to load both magnetic" << endl;
*gmsg << "* and electric fields, setting SUPERPOSE to an array of TRUE values.)" << endl;
*gmsg << "* 2.) It is strongly recommended to use the SAAMG fieldsolver," << endl;
*gmsg << "* FFT does not give the correct results (boundaty conditions are missing)." << endl;
*gmsg << "* 3.) The whole geometry will be meshed and used for the fieldsolve." << endl;
*gmsg << "* There will be no transformations of the bunch into a local frame und consequently," << endl;
*gmsg << "* the problem will be treated non-relativistically!" << endl;
*gmsg << "* (This is not an issue for spiral inflectors as they are typically < 100 keV/amu.)" << endl;
*gmsg << endl << "* Note: For now, multi-bunch mode (MBM) needs to be de-activated for spiral inflector" << endl;
*gmsg << "* and space charge needs to be solved every time-step. numBunch_m and scSolveFreq are reset." << endl;
numBunch_m = 1;
}
// Fresh run (no restart):
if(!OpalData::getInstance()->inRestartRun()) {
......@@ -1295,7 +1317,7 @@ void ParallelCyclotronTracker::Tracker_LF() {
// 2.after each revolution
// 3.existing bunches is less than the specified bunches
// 4.FORCE mode, or AUTO mode with flagTransition = true
// Note: restart from 1 < BunchCount < numBunch_m must be avoided.
// Note: restart from 1 < BnchCount < numBunch_m must be avoided.
*gmsg << "step " << step_m << ", inject a new bunch... ... ..." << endl;
BunchCount_m++;
......@@ -2787,7 +2809,14 @@ void ParallelCyclotronTracker::Tracker_Generic() {
// Read in some control parameters
const int SinglePartDumpFreq = Options::sptDumpFreq;
const int resetBinFreq = Options::rebinFreq;
const int scSolveFreq = Options::scSolveFreq;
int scSolveFreq;
if (spiral_flag)
scSolveFreq = 1;
else
scSolveFreq = Options::scSolveFreq;
const int boundpDestroyFreq = Options::boundpDestroyFreq;
const bool doDumpAfterEachTurn = Options::psDumpEachTurn;
......@@ -2815,7 +2844,7 @@ void ParallelCyclotronTracker::Tracker_Generic() {
int stepsNextCheck = step_m + stepsPerTurn; // Steps to next check for transition
//int stepToLastInj = itsBunch->getSteptoLastInj(); // TODO: Do we need this? -DW
// Record how many bunches have already been injected. ONLY FOR MPM
// Record how many bunches have already been injected. ONLY FOR MBM
BunchCount_m = itsBunch->getNumBunch();
BinCount_m = BunchCount_m;
......@@ -3096,40 +3125,61 @@ void ParallelCyclotronTracker::Tracker_Generic() {
} else {
// --- Single bunch mode --- //
double temp_meangamma = sqrt(1.0 + dot(PreviousMeanP, PreviousMeanP));
// If we are doing a spiral inflector simulation and are using the SAAMG solver
// we don't rotate or translate the bunch and gamma is 1.0 (non-relativistic).
if (spiral_flag and itsBunch->getFieldSolverType() == "SAAMG") {
Quaternion_t quaternionToYAxis;
itsBunch->R *= Vector_t(0.001); // mm --> m
getQuaternionTwoVectors(PreviousMeanP, yaxis, quaternionToYAxis);
IpplTimings::stopTimer(TransformTimer_m);
globalToLocal(itsBunch->R, quaternionToYAxis, meanR);
itsBunch->setGlobalMeanR(Vector_t(0.0, 0.0, 0.0));
itsBunch->setGlobalToLocalQuaternion(Quaternion_t(1.0, 0.0, 0.0, 0.0));
itsBunch->R *= Vector_t(0.001); // mm --> m
itsBunch->computeSelfFields_cycl(1.0);
if((step_m + 1) % boundpDestroyFreq == 0)
itsBunch->boundp_destroy();
else
itsBunch->boundp();
IpplTimings::startTimer(TransformTimer_m);
IpplTimings::stopTimer(TransformTimer_m);
itsBunch->R *= Vector_t(1000.0); // m --> mm
repartition();
} else {
double temp_meangamma = sqrt(1.0 + dot(PreviousMeanP, PreviousMeanP));
itsBunch->setGlobalMeanR(0.001 * meanR);
itsBunch->setGlobalToLocalQuaternion(quaternionToYAxis);
itsBunch->computeSelfFields_cycl(temp_meangamma);
Quaternion_t quaternionToYAxis;
IpplTimings::startTimer(TransformTimer_m);
getQuaternionTwoVectors(PreviousMeanP, yaxis, quaternionToYAxis);
//scale coordinates back
itsBunch->R *= Vector_t(1000.0); // m --> mm
globalToLocal(itsBunch->R, quaternionToYAxis, meanR);
// Transform coordinates back to global
localToGlobal(itsBunch->R, quaternionToYAxis, meanR);
itsBunch->R *= Vector_t(0.001); // mm --> m
// Transform self field back to global frame (rotate only)
localToGlobal(itsBunch->Ef, quaternionToYAxis);
localToGlobal(itsBunch->Bf, quaternionToYAxis);
if((step_m + 1) % boundpDestroyFreq == 0)
itsBunch->boundp_destroy();
else
itsBunch->boundp();
IpplTimings::stopTimer(TransformTimer_m);
repartition();
itsBunch->setGlobalMeanR(0.001 * meanR);
itsBunch->setGlobalToLocalQuaternion(quaternionToYAxis);
itsBunch->computeSelfFields_cycl(temp_meangamma);
IpplTimings::startTimer(TransformTimer_m);
//scale coordinates back
itsBunch->R *= Vector_t(1000.0); // m --> mm
// Transform coordinates back to global
localToGlobal(itsBunch->R, quaternionToYAxis, meanR);
// Transform self field back to global frame (rotate only)
localToGlobal(itsBunch->Ef, quaternionToYAxis);
localToGlobal(itsBunch->Bf, quaternionToYAxis);
}
}
IpplTimings::stopTimer(TransformTimer_m);
......
......@@ -227,6 +227,8 @@ private:
double referencePsi;
double referencePhi;
bool spiral_flag = false;
Vector_t PreviousMeanP;
bool previousH5Local;
......@@ -501,4 +503,4 @@ inline
double ParallelCyclotronTracker::calculateAngle2(double x, double y)
{ return atan2(y,x); }
#endif // OPAL_ParallelCyclotronTracker_HH
\ No newline at end of file
#endif // OPAL_ParallelCyclotronTracker_HH
......@@ -23,7 +23,6 @@
#include "Structure/BoundaryGeometry.h"
#include "Physics/Physics.h"
// Class OpalCyclotron
// ------------------------------------------------------------------------
......@@ -109,6 +108,10 @@ OpalCyclotron::OpalCyclotron():
itsAttr[FMHIGHE] = Attributes::makeReal
("FMHIGHE","Maximal Energy [MeV] the fieldmap can accept. Used in GAUSSMATCHED", -1.0);
itsAttr[SPIRAL] = Attributes::makeBool
("SPIRAL","Flag whether or not this is a spiral inflector simulation", false);
registerStringAttribute("FMAPFN");
registerStringAttribute("GEOMETRY");
registerStringAttribute("RFMAPFN");
......@@ -191,6 +194,8 @@ void OpalCyclotron::update() {
double fmLowE = Attributes::getReal(itsAttr[FMLOWE]);
double fmHighE = Attributes::getReal(itsAttr[FMHIGHE]);
bool spiral_flag = Attributes::getBool(itsAttr[SPIRAL]);
cycl->setFieldMapFN(fmapfm);
cycl->setSymmetry(symmetry);
......@@ -218,8 +223,7 @@ void OpalCyclotron::update() {
cycl->setFMLowE(fmLowE);
cycl->setFMHighE(fmHighE);
cycl->setSpiralFlag(spiral_flag);
std::vector<std::string> fm_str = Attributes::getStringArray(itsAttr[RFMAPFN]);
std::vector<double> scale_str = Attributes::getRealArray(itsAttr[ESCALE]);
......
......@@ -45,10 +45,10 @@ public:
RFMAPFN, // The filename of the fieldmap
BSCALE, // A scalar to scale the B-field
ESCALE, // A scalar to scale the RF field
TCR1, //trim coil r1 (mm)
TCR2, //trim coil r2 (mm)
MBTC, //max bfield of trim coil (kG)
SLPTC, //slope of the rising edge
TCR1, // trim coil r1 (mm)
TCR2, // trim coil r2 (mm)
MBTC, // max bfield of trim coil (kG)
SLPTC, // slope of the rising edge
RFPHI, // the initial phase of RF field
SUPERPOSE, // whether the electric field map are superposed or not
MINZ, // minimal vertical extend of the machine
......@@ -57,6 +57,7 @@ public:
MAXR, // maximal radial extend of the machine
FMLOWE, // minimal energy of the field map
FMHIGHE, // maximal energy of the field map
SPIRAL, // flag whether or not this is a spiral inflector simulation
SIZE
};
......
......@@ -54,7 +54,7 @@ ArbitraryDomain::~ArbitraryDomain() {
}
void ArbitraryDomain::Compute(Vector_t hr){
setHr(hr);
globalMeanR_m = getGlobalMeanR();
......@@ -191,7 +191,7 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId){
globalToLocalQuaternion_m = getGlobalToLocalQuaternion();
localToGlobalQuaternion_m[0] = globalToLocalQuaternion_m[0];
for (int i=1; i<4; i++)
localToGlobalQuaternion_m[i] = -globalToLocalQuaternion_m[i];
localToGlobalQuaternion_m[i] = -globalToLocalQuaternion_m[i];
int zGhostOffsetLeft = (localId[2].first()== 0) ? 0 : 1;
int zGhostOffsetRight = (localId[2].last() == nr[2] - 1) ? 0 : 1;
......@@ -209,31 +209,54 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId){
IntersectLoZ.clear();
IntersectHiZ.clear();
//calculate intersection
// Calculate intersection
Vector_t P, saveP, dir, I;
//Reference Point
// Get minimum coordinates of boundary geometry
Vector_t geomMinCoords = bgeom_m->getmincoords();
// Reference Point that is inside (allowed particle region)
// This is dangerous, who knows if this is really an "inside" point? -DW
// TODO: Change to more generalized approach
Vector_t P0 = geomCentroid_m;
for (int idz = localId[2].first()-zGhostOffsetLeft; idz <= localId[2].last()+zGhostOffsetRight; idz++) {
saveP[2] = (idz - (nr[2]-1)/2.0)*hr[2];
for (int idy = localId[1].first()-yGhostOffsetLeft; idy <= localId[1].last()+yGhostOffsetRight; idy++) {
saveP[1] = (idy - (nr[1]-1)/2.0)*hr[1];
//saveP[2] = (idz - (nr[2]-1)/2.0)*hr[2];
saveP[2] = geomMinCoords[2] + (idz + 0.5) * hr[2];
for (int idy = localId[1].first()-yGhostOffsetLeft; idy <= localId[1].last()+yGhostOffsetRight; idy++) {
//saveP[1] = (idy - (nr[1]-1)/2.0)*hr[1];
saveP[1] = geomMinCoords[1] + (idy + 0.5) * hr[1];
for (int idx = localId[0].first()-xGhostOffsetLeft; idx <= localId[0].last()+xGhostOffsetRight; idx++) {
saveP[0] = (idx - (nr[0]-1)/2.0)*hr[0];
// We cannot assume that the geometry is symmetric about the xy, xz, and yz planes!
// In my spiral inflector simulation, this is not the case for z direction for
// example. -DW
// saveP[0] = (idx - (nr[0]-1)/2.0)*hr[0];
saveP[0] = geomMinCoords[0] + (idx + 0.5) * hr[0];
P = saveP;
rotateWithQuaternion(P, localToGlobalQuaternion_m);
P += geomCentroid_m;
//P += geomCentroid_m; //sorry, this doesn't make sense. -DW
P += globalMeanR_m;
if (bgeom_m->fastIsInside(P0, P) % 2 == 0) {
P0 = P;
std::tuple<int, int, int> pos(idx, idy, idz);
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
*gmsg << "INSIDE" << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
#endif
// Commenting this out seems to reduce the number of INSIDE points that fail
// the intersection test. It was meant to make the algorithm faster, but
// is not strictly necessary... -DW
//P0 = P;
std::tuple<int, int, int> pos(idx, idy, idz);
rotateZAxisWithQuaternion(dir, localToGlobalQuaternion_m);
if (bgeom_m->intersectRayBoundary(P, dir, I)) {
I -= geomCentroid_m;
rotateWithQuaternion(I, globalToLocalQuaternion_m);
//I -= geomCentroid_m;
I -= globalMeanR_m;
rotateWithQuaternion(I, globalToLocalQuaternion_m);
IntersectHiZ.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[2]));
} else {
#ifdef DEBUG_INTERSECT_RAY_BOUNDARY
......@@ -242,7 +265,8 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId){
}
if (bgeom_m->intersectRayBoundary(P, -dir, I)) {
I -= geomCentroid_m;
//I -= geomCentroid_m;
I -= globalMeanR_m;
rotateWithQuaternion(I, globalToLocalQuaternion_m);
IntersectLoZ.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[2]));
} else {
......@@ -252,8 +276,10 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId){
}
rotateYAxisWithQuaternion(dir, localToGlobalQuaternion_m);
if (bgeom_m->intersectRayBoundary(P, dir, I)) {
I -= geomCentroid_m;
//I -= geomCentroid_m;
I -= globalMeanR_m;
rotateWithQuaternion(I, globalToLocalQuaternion_m);
IntersectHiY.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[1]));
} else {
......@@ -263,7 +289,8 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId){
}
if (bgeom_m->intersectRayBoundary(P, -dir, I)) {
I -= geomCentroid_m;
//I -= geomCentroid_m;
I -= globalMeanR_m;
rotateWithQuaternion(I, globalToLocalQuaternion_m);
IntersectLoY.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[1]));
} else {
......@@ -271,10 +298,12 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId){
*gmsg << "ydir=-1" << -dir << " x,y,z= " << idx << "," << idy << "," << idz << " P=" << P <<" I=" << I << endl;
#endif
}
rotateXAxisWithQuaternion(dir, localToGlobalQuaternion_m);
if (bgeom_m->intersectRayBoundary(P, dir, I)) {
I -= geomCentroid_m;
//I -= geomCentroid_m;
I -= globalMeanR_m;
rotateWithQuaternion(I, globalToLocalQuaternion_m);
IntersectHiX.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[0]));
} else {
......@@ -284,7 +313,8 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId){
}
if (bgeom_m->intersectRayBoundary(P, -dir, I)){
I -= geomCentroid_m;
//I -= geomCentroid_m;
I -= globalMeanR_m;
rotateWithQuaternion(I, globalToLocalQuaternion_m);
IntersectLoX.insert(std::pair< std::tuple<int, int, int>, double >(pos, I[0]));
} else {
......@@ -300,6 +330,7 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId){
}
}
}
//number of ghost nodes to the left
int numGhostNodesLeft = 0;
if(localId[2].first() != 0) {
......
......@@ -261,6 +261,7 @@ void MGPoissonSolver::computePotential(Field_t &rho, Vector_t hr) {
bp->setGlobalMeanR(itsBunch_m->getGlobalMeanR());
bp->setGlobalToLocalQuaternion(itsBunch_m->getGlobalToLocalQuaternion());
bp->setNr(nr_m);
NDIndex<3> localId = layout_m->getLocalNDIndex();
IpplTimings::startTimer(FunctionTimer1_m);
......@@ -277,6 +278,7 @@ void MGPoissonSolver::computePotential(Field_t &rho, Vector_t hr) {
if (Teuchos::is_null(RHS))
RHS = rcp(new Epetra_Vector(*Map));
RHS->PutScalar(0.0);
// get charge densities from IPPL field and store in Epetra vector (RHS)
IpplTimings::startTimer(FunctionTimer3_m);
int id = 0;
......
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