Commit 051beb73 authored by Kaman Tülin's avatar Kaman Tülin

Set meanR and the quaternion of the rotation of the particle bunch in METHOD

parent 386c405f
......@@ -690,11 +690,6 @@ double PartBunch::getMaxdEBins() {
void PartBunch::computeSelfFields(int binNumber) {
IpplTimings::startTimer(selfFieldTimer_m);
// Need to set mean R to shift geometry from file to account for beam movement
// No rotations necessary, i.e. quaternion stays default unit-quatertion
// TODO: Include transversal offsets here!
globalMeanR_m = Vector_t(0.0, 0.0, get_sPos());
/// Set initial charge density to zero. Create image charge
/// potential field.
rho_m = 0.0;
......@@ -947,11 +942,6 @@ void PartBunch::computeSelfFields() {
rho_m = 0.0;
eg_m = Vector_t(0.0);
// Need to set mean R to shift geometry from file to account for beam movement
// No rotations necessary, i.e. quaternion stays default unit-quatertion
// TODO: Include transversal offsets here!
globalMeanR_m = Vector_t(0.0, 0.0, get_sPos());
if(fs_m->hasValidSolver()) {
//use the mesh that is already set
//if (fs_m->getFieldSolverType() == "SAAMG")
......@@ -1216,15 +1206,10 @@ void PartBunch::computeSelfFields_cycl(double gamma) {
* -) Fixed an error where gamma was not taken into account correctly in direction of movement (y in cyclotron)
* -) Comment: There is no account for image charges in the cyclotron tracker (yet?)!
*/
void PartBunch::computeSelfFields_cycl(double gamma,
Vector_t const meanR,
Quaternion_t const quaternion) {
void PartBunch::computeSelfFields_cycl(double gamma) {
IpplTimings::startTimer(selfFieldTimer_m);
globalMeanR_m = meanR;
globalToLocalQuaternion_m = quaternion;
/// set initial charge density to zero.
rho_m = 0.0;
......@@ -1411,12 +1396,9 @@ void PartBunch::computeSelfFields_cycl(double gamma,
* -) Fixed an error where gamma was not taken into account correctly in direction of movement (y in cyclotron)
* -) Comment: There is no account for image charges in the cyclotron tracker (yet?)!
*/
void PartBunch::computeSelfFields_cycl(int bin, Vector_t const meanR, Quaternion_t const quaternion) {
void PartBunch::computeSelfFields_cycl(int bin) {
IpplTimings::startTimer(selfFieldTimer_m);
globalMeanR_m = meanR;
globalToLocalQuaternion_m = quaternion;
/// set initial charge dentsity to zero.
rho_m = 0.0;
......
......@@ -277,17 +277,8 @@ public:
/** /brief used for self fields with binned distribution */
void computeSelfFields(int b);
//void computeSelfFields_cycl(double gamma);
//void computeSelfFields_cycl(int b);
// Replaced computeSelfFields_cycl() with versions that have meanR and the quaternion of the
// rotation of the particle bunch in order to take into account the rotation
// when finding the boundary conditions for the fieldsolver. -DW
void computeSelfFields_cycl(double gamma, Vector_t const meanR=Vector_t(0.0, 0.0, 0.0),
Quaternion_t const quaternion=Quaternion_t(1.0, 0.0, 0.0, 0.0));
void computeSelfFields_cycl(int b, Vector_t const meanR=Vector_t(0.0, 0.0, 0.0),
Quaternion_t const quaternion=Quaternion_t(1.0, 0.0, 0.0, 0.0));
void computeSelfFields_cycl(double gamma);
void computeSelfFields_cycl(int b);
void resetInterpolationCache(bool clearCache = false);
......@@ -410,8 +401,11 @@ public:
inline void setNumBunch(int n);
inline int getNumBunch() const;
inline void setGlobalMeanR(Vector_t globalMeanR) {globalMeanR_m = globalMeanR;}
inline Vector_t getGlobalMeanR() {return globalMeanR_m;}
inline Vektor<double, 4> getGlobalToLocalQuaternion() {return globalToLocalQuaternion_m;}
inline void setGlobalToLocalQuaternion(Quaternion_t globalToLocalQuaternion) {
globalToLocalQuaternion_m = globalToLocalQuaternion;}
inline Quaternion_t getGlobalToLocalQuaternion() {return globalToLocalQuaternion_m;}
void setSteptoLastInj(int n);
int getSteptoLastInj();
......
......@@ -2187,7 +2187,9 @@ void ParallelCyclotronTracker::Tracker_RK4() {
repartition();
itsBunch->computeSelfFields_cycl(temp_meangamma, 0.001 * meanR, quaternionToYAxis);
itsBunch->setGlobalMeanR(0.001 * meanR);
itsBunch->setGlobalToLocalQuaternion(quaternionToYAxis);
itsBunch->computeSelfFields_cycl(temp_meangamma);
IpplTimings::startTimer(TransformTimer_m);
......@@ -3070,7 +3072,9 @@ void ParallelCyclotronTracker::Tracker_Generic() {
for(int b = 0; b < itsBunch->getLastemittedBin(); b++) {
itsBunch->setBinCharge(b, itsBunch->getChargePerParticle());
itsBunch->computeSelfFields_cycl(b, 0.001 * meanR, quaternionToYAxis);
itsBunch->setGlobalMeanR(0.001 * meanR);
itsBunch->setGlobalToLocalQuaternion(quaternionToYAxis);
itsBunch->computeSelfFields_cycl(b);
}
itsBunch->Q = itsBunch->getChargePerParticle();
......@@ -3108,7 +3112,9 @@ void ParallelCyclotronTracker::Tracker_Generic() {
repartition();
itsBunch->computeSelfFields_cycl(temp_meangamma, 0.001 * meanR, quaternionToYAxis);
itsBunch->setGlobalMeanR(0.001 * meanR);
itsBunch->setGlobalToLocalQuaternion(quaternionToYAxis);
itsBunch->computeSelfFields_cycl(temp_meangamma);
IpplTimings::startTimer(TransformTimer_m);
......
......@@ -2246,12 +2246,14 @@ void ParallelTTracker::computeSpaceChargeFields() {
binNumber < itsBunch->GetNumberOfEnergyBins(); ++binNumber) {
itsBunch->setBinCharge(binNumber);
itsBunch->setGlobalMeanR(itsBunch->get_centroid());
itsBunch->computeSelfFields(binNumber);
itsBunch->Q = Q_back;
}
} else {
itsBunch->setGlobalMeanR(itsBunch->get_centroid());
itsBunch->computeSelfFields();
}
......@@ -3217,4 +3219,4 @@ Vector_t ParallelTTracker::calcMeanP() const {
}
reduce(meanP, meanP, OpAddAssign());
return meanP / Vector_t(itsBunch->getTotalNum());
}
\ No newline at end of file
}
......@@ -25,27 +25,25 @@
#define MIN2(a,b) (((a) < (b)) ? (a) : (b))
#define MAX2(a,b) (((a) > (b)) ? (a) : (b))
ArbitraryDomain::ArbitraryDomain(
BoundaryGeometry * bgeom,
Vector_t nr,
Vector_t hr,
std::string interpl) {
bgeom_m = bgeom;
minCoords_m = bgeom->getmincoords();
maxCoords_m = bgeom->getmaxcoords();
setNr(nr);
setHr(hr);
startId = 0;
if(interpl == "CONSTANT")
interpolationMethod = CONSTANT;
else if(interpl == "LINEAR")
interpolationMethod = LINEAR;
else if(interpl == "QUADRATIC")
interpolationMethod = QUADRATIC;
ArbitraryDomain::ArbitraryDomain( BoundaryGeometry * bgeom,
Vector_t nr,
Vector_t hr,
std::string interpl){
bgeom_m = bgeom;
minCoords_m = bgeom->getmincoords();
maxCoords_m = bgeom->getmaxcoords();
setNr(nr);
setHr(hr);
startId = 0;
if (interpl == "CONSTANT")
interpolationMethod = CONSTANT;
else if(interpl == "LINEAR")
interpolationMethod = LINEAR;
else if(interpl == "QUADRATIC")
interpolationMethod = QUADRATIC;
}
ArbitraryDomain::~ArbitraryDomain() {
......@@ -57,7 +55,7 @@ void ArbitraryDomain::Compute(Vector_t hr) {
setHr(hr);
}
/*
void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId) {
setHr(hr);
......@@ -224,16 +222,18 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId) {
}
}
}
*/
void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId, Vector_t globalMeanR, Vektor<double, 4> globalToLocalQuaternion){
void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId){
setHr(hr);
globalMeanR_m = getCentroid();
globalToLocalQuaternion_m = globalToLocalQuaternion;
localToGlobalQuaternion_m[0] = globalToLocalQuaternion[0];
globalMeanR_m = getGlobalMeanR();
globalToLocalQuaternion_m = getGlobalToLocalQuaternion();
localToGlobalQuaternion_m[0] = globalToLocalQuaternion_m[0];
for (int i=1; i<4; i++)
localToGlobalQuaternion_m[i] = -globalToLocalQuaternion[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;
......@@ -612,7 +612,7 @@ inline void ArbitraryDomain::crossProduct(double A[], double B[], double C[]) {
C[2] = A[0] * B[1] - A[1] * B[0];
}
inline void ArbitraryDomain::rotateWithQuaternion(Vector_t & v, Vektor<double, 4> const quaternion) {
inline void ArbitraryDomain::rotateWithQuaternion(Vector_t & v, Quaternion_t const quaternion) {
// rotates a Vector_t (3 elements) using a quaternion.
// Flip direction of rotation by quaternionVectorcomponent *= -1
......@@ -625,7 +625,7 @@ inline void ArbitraryDomain::rotateWithQuaternion(Vector_t & v, Vektor<double, 4
+ 2.0 * quaternionScalarComponent * cross(quaternionVectorComponent, v);
}
inline void ArbitraryDomain::rotateXAxisWithQuaternion(Vector_t & v, Vektor<double, 4> const quaternion) {
inline void ArbitraryDomain::rotateXAxisWithQuaternion(Vector_t & v, Quaternion_t const quaternion) {
// rotates the positive xaxis using a quaternion.
v(0) = quaternion(0) * quaternion(0)
......@@ -637,7 +637,7 @@ inline void ArbitraryDomain::rotateXAxisWithQuaternion(Vector_t & v, Vektor<doub
v(2) = 2.0 * (quaternion(1) * quaternion(3) - quaternion(0) * quaternion(2));
}
inline void ArbitraryDomain::rotateYAxisWithQuaternion(Vector_t & v, Vektor<double, 4> const quaternion) {
inline void ArbitraryDomain::rotateYAxisWithQuaternion(Vector_t & v, Quaternion_t const quaternion) {
// rotates the positive yaxis using a quaternion.
v(0) = 2.0 * (quaternion(1) * quaternion(2) - quaternion(0) * quaternion(3));
......@@ -650,7 +650,7 @@ inline void ArbitraryDomain::rotateYAxisWithQuaternion(Vector_t & v, Vektor<doub
v(2) = 2.0 * (quaternion(2) * quaternion(3) + quaternion(0) * quaternion(1));
}
inline void ArbitraryDomain::rotateZAxisWithQuaternion(Vector_t & v, Vektor<double, 4> const quaternion) {
inline void ArbitraryDomain::rotateZAxisWithQuaternion(Vector_t & v, Quaternion_t const quaternion) {
// rotates the positive zaxis using a quaternion.
v(0) = 2.0 * (quaternion(1) * quaternion(3) + quaternion(0) * quaternion(2));
v(1) = 2.0 * (quaternion(2) * quaternion(3) - quaternion(0) * quaternion(1));
......
......@@ -20,7 +20,7 @@ class ArbitraryDomain : public IrregularDomain {
public:
ArbitraryDomain(BoundaryGeometry *bgeom, Vector_t nr, Vector_t hr, std::string interpl);
ArbitraryDomain(BoundaryGeometry *bgeom, Vector_t nr, Vector_t hr, Vector_t globalMeanR, Vektor<double, 4> globalToLocalQuaternion, std::string interpl);
ArbitraryDomain(BoundaryGeometry *bgeom, Vector_t nr, Vector_t hr, Vector_t globalMeanR, Quaternion_t globalToLocalQuaternion, std::string interpl);
~ArbitraryDomain();
......@@ -41,9 +41,8 @@ public:
int getNumXY(int idz);
/// calculates intersection
void Compute(Vector_t hr);
void Compute(Vector_t hr, NDIndex<3> localId);
// calculates intersection with rotated and shifted geometry
void Compute(Vector_t hr, NDIndex<3> localId, Vector_t globalMeanR, Vektor<double, 4> globalToLocalQuaternion);
void Compute(Vector_t hr, NDIndex<3> localId);
int getStartId() {return startId;}
......@@ -83,10 +82,9 @@ private:
// meanR to shift from global to local frame
Vector_t globalMeanR_m;
Vektor<double, 4> globalToLocalQuaternion_m;
Vektor<double, 4> localToGlobalQuaternion_m;
Quaternion_t globalToLocalQuaternion_m;
Quaternion_t localToGlobalQuaternion_m;
Vector_t globalMinR_m;
int startId;
// Here we store the number of nodes in a xy layer for a given z coordinate
......@@ -127,11 +125,11 @@ private:
void QuadraticInterpolation(int idx, int idy, int idz, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor);
// Rotate positive axes with quaternion -DW
inline void rotateWithQuaternion(Vector_t &v, Vektor<double, 4> const quaternion);
inline void rotateWithQuaternion(Vector_t &v, Quaternion_t const quaternion);
inline void rotateXAxisWithQuaternion(Vector_t &v, Vektor<double, 4> const quaternion);
inline void rotateYAxisWithQuaternion(Vector_t &v, Vektor<double, 4> const quaternion);
inline void rotateZAxisWithQuaternion(Vector_t &v, Vektor<double, 4> const quaternion);
inline void rotateXAxisWithQuaternion(Vector_t &v, Quaternion_t const quaternion);
inline void rotateYAxisWithQuaternion(Vector_t &v, Quaternion_t const quaternion);
inline void rotateZAxisWithQuaternion(Vector_t &v, Quaternion_t const quaternion);
};
#endif //#ifdef HAVE_SAAMG_SOLVER
......
......@@ -133,10 +133,6 @@ void BoxCornerDomain::Compute(Vector_t hr){
void BoxCornerDomain::Compute(Vector_t hr, NDIndex<3> localId){
}
void BoxCornerDomain::Compute(Vector_t hr, NDIndex<3> localId, Vector_t globalMeanR, Vektor<double, 4> globalToLocalQuaternion){
}
void BoxCornerDomain::getBoundaryStencil(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor) {
// determine which interpolation method we use for points near the boundary
......
......@@ -101,8 +101,6 @@ public:
void Compute(Vector_t hr);
void Compute(Vector_t hr, NDIndex<3> localId);
void Compute(Vector_t hr, NDIndex<3> localId, Vector_t globalMeanR, Vektor<double, 4> globalToLocalQuaternion);
double getXRangeMin() { return -A_m; }
double getXRangeMax() { return A_m; }
......
......@@ -191,10 +191,6 @@ void EllipticDomain::Compute(Vector_t hr, NDIndex<3> localId){
}
}
void EllipticDomain::Compute(Vector_t hr, NDIndex<3> localId, Vector_t globalMeanR, Vektor<double, 4> globalToLocalQuaternion){
}
void EllipticDomain::getBoundaryStencil(int x, int y, int z, double &W, double &E, double &S, double &N, double &F, double &B, double &C, double &scaleFactor) {
scaleFactor = 1.0;
......
......@@ -42,7 +42,6 @@ public:
/// calculates intersection
void Compute(Vector_t hr);
void Compute(Vector_t hr, NDIndex<3> localId);
void Compute(Vector_t hr, NDIndex<3> localId, Vector_t globalMeanR, Vektor<double, 4> globalToLocalQuaternion);
double getXRangeMin() { return -SemiMajor; }
double getXRangeMax() { return SemiMajor; }
......
......@@ -24,7 +24,6 @@ public:
*/
virtual void Compute(Vector_t hr) = 0;
virtual void Compute(Vector_t hr, NDIndex<3> localId) = 0;
virtual void Compute(Vector_t hr, NDIndex<3> localId, Vector_t globalMeanR, Vektor<double, 4> globalToLocalQuaternion) = 0;
/** method to get the number of gridpoints in a given z plane
* \param z coordinate of the z plane
* \return int number of grid nodes in the given z plane
......@@ -88,11 +87,12 @@ public:
double getMinZ() { return zMin_m; }
double getMaxZ() { return zMax_m; }
void setCentroid(Vector_t rmean) { rMean_m = rmean;}
Vector_t getCentroid() { return rMean_m; }
void setGlobalMeanR(Vector_t rmean) { rMean_m = rmean;}
Vector_t getGlobalMeanR() { return rMean_m; }
void setOrigin(Vector_t rmin) { rMin_m = rmin;}
Vector_t getOrigin() { return rMin_m; }
void setGlobalToLocalQuaternion(Quaternion_t globalToLocalQuaternion){
globalToLocalQuaternion_m = globalToLocalQuaternion;}
Quaternion_t getGlobalToLocalQuaternion() { return globalToLocalQuaternion_m;}
virtual double getXRangeMin() = 0;
virtual double getXRangeMax() = 0;
......@@ -118,8 +118,7 @@ protected:
/// mean position of bunch (m)
Vector_t rMean_m;
/// min position of the bunch
Vector_t rMin_m;
Quaternion_t globalToLocalQuaternion_m;
};
#endif //#ifdef HAVE_SAAMG_SOLVER
......
......@@ -256,14 +256,14 @@ void MGPoissonSolver::computePotential(Field_t &rho, Vector_t hr) {
nr_m[2] = orig_nr_m[2];
//bp->setMinMaxZ(itsBunch_m->get_origin()[2], itsBunch_m->get_maxExtend()[2]);
bp->setCentroid(itsBunch_m->get_centroid());
bp->setOrigin(itsBunch_m->get_origin());
bp->setGlobalMeanR(itsBunch_m->getGlobalMeanR());
bp->setGlobalToLocalQuaternion(itsBunch_m->getGlobalToLocalQuaternion());
bp->setNr(nr_m);
NDIndex<3> localId = layout_m->getLocalNDIndex();
IpplTimings::startTimer(FunctionTimer1_m);
if ( bp->getType() == "Geometric" ) {
bp->Compute(hr, localId, itsBunch_m->getGlobalMeanR(), itsBunch_m->getGlobalToLocalQuaternion());
bp->Compute(hr, localId);
} else if ( bp->getType() == "Elliptic" )
bp->Compute(hr);
IpplTimings::stopTimer(FunctionTimer1_m);
......
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