Commit 845b83b4 authored by Daniel Winklehner's avatar Daniel Winklehner
Browse files

Brought back saving in local frame for cyclotron tracker. Reverted the way...

Brought back saving in local frame for cyclotron tracker. Reverted the way mean E is calculated in PartBunch. Hopefully, restarts in local as well as in global frame now work in the cyclotron tracker.
parent c8b49f25
......@@ -47,6 +47,7 @@
#endif
//#define DBG_SCALARFIELD
#define FIELDSTDOUT
using Physics::pi;
......@@ -976,7 +977,7 @@ void PartBunch::computeSelfFields() {
eg_m *= Vector_t(gammaz / (scaleFactor), gammaz / (scaleFactor), 1.0 / (scaleFactor * gammaz));
//write out e field
#ifdef DBG_SCALARFIELD
#ifdef FIELDSTDOUT
// Immediate debug output:
// Output potential and e-field along the x-, y-, and z-axes
int m1 = (int)nr_m[0]-1;
......@@ -990,7 +991,8 @@ void PartBunch::computeSelfFields() {
for (int i=0; i<m1; i++ )
*gmsg << "Field along z axis E = " << eg_m[m2][m2][i] << " Pot = " << rho_m[m2][m2][i] << endl;
#endif
#ifdef DBG_SCALARFIELD
INFOMSG("*** START DUMPING E FIELD ***" << endl);
//ostringstream oss;
//MPI_File file;
......@@ -1236,7 +1238,7 @@ void PartBunch::computeSelfFields_cycl(double gamma,
/// only hr_scaled is! -DW
eg_m *= Vector_t(gamma, 1.0 / gamma, gamma);
#ifdef DBG_SCALARFIELD
#ifdef FIELDSTDOUT
// Immediate debug output:
// Output potential and e-field along the x-, y-, and z-axes
int m1 = (int)nr_m[0]-1;
......@@ -1250,7 +1252,9 @@ void PartBunch::computeSelfFields_cycl(double gamma,
for (int i=0; i<m1; i++ )
*gmsg << "Field along z axis E = " << eg_m[m2][m2][i] << " Pot = " << rho_m[m2][m2][i] << endl;
#endif
#ifdef DBG_SCALARFIELD
// If debug flag is set, dump vector field (electric field) into file under ./data/
INFOMSG("*** START DUMPING E FIELD ***" << endl);
//ostringstream oss;
......@@ -2337,24 +2341,24 @@ void PartBunch::calcEMean() {
const double totalNp = static_cast<double>(this->getTotalNum());
const double locNp = static_cast<double>(this->getLocalNum());
Vector_t meanP_temp = Vector_t(0.0);
//Vector_t meanP_temp = Vector_t(0.0);
eKin_m = 0.0;
for(unsigned int k = 0; k < locNp; k++)
meanP_temp += P[k];
//eKin_m += sqrt(dot(P[k], P[k]) + 1.0);
//meanP_temp += P[k];
eKin_m += sqrt(dot(P[k], P[k]) + 1.0);
//eKin_m -= locNp;
//eKin_m *= getM() * 1.0e-6;
eKin_m -= locNp;
eKin_m *= getM() * 1.0e-6;
reduce(meanP_temp, meanP_temp, OpAddAssign());
//reduce(eKin_m, eKin_m, OpAddAssign());
//reduce(meanP_temp, meanP_temp, OpAddAssign());
reduce(eKin_m, eKin_m, OpAddAssign());
meanP_temp /= totalNp;
//eKin_m /= totalNp;
//meanP_temp /= totalNp;
eKin_m /= totalNp;
eKin_m = (sqrt(dot(meanP_temp, meanP_temp) + 1.0) - 1.0) * getM() * 1.0e-6;
//eKin_m = (sqrt(dot(meanP_temp, meanP_temp) + 1.0) - 1.0) * getM() * 1.0e-6;
}
......
......@@ -152,7 +152,6 @@ public:
/// set multipacting flag
virtual void setMpacflg(bool mpacflg) {};
virtual void setPr(double x) { } ;
virtual void setPt(double x) { } ;
virtual void setPz(double x) { } ;
......@@ -160,7 +159,9 @@ public:
virtual void setTheta(double x) { } ;
virtual void setZ(double x) { } ;
virtual void setBeGa(double x) { } ;
virtual void setPhi(double x) { } ;
virtual void setPsi(double x) { } ;
virtual void setPreviousH5Local(bool x) { } ;
// standing wave structures
FieldList cavities_m;
......
This diff is collapsed.
......@@ -177,6 +177,9 @@ public:
inline void setTheta(double x) {referenceTheta = x;}
inline void setZ(double x) {referenceZ = x;}
inline void setBeGa(double x) {bega = x;}
inline void setPhi(double x) {referencePhi = x;}
inline void setPsi(double x) {referencePsi = x;}
inline void setPreviousH5Local(bool x) {previousH5Local = x;}
void bgf_main_collision_test();
void initializeBoundaryGeometry();
......@@ -221,8 +224,13 @@ private:
double referencePz = 0.0;
double referencePtot;
double referencePsi;
double referencePhi;
Vector_t PreviousMeanP;
bool previousH5Local;
double sinRefTheta_m;
double cosRefTheta_m;
......@@ -364,23 +372,43 @@ private:
// global reference frame to the local reference frame.
// phi is the angle of the bunch measured counter-clockwise from the positive x-axis.
void globalToLocal(ParticleAttrib<Vector_t> & vectorArray, double phi, Vector_t const translationToGlobal = 0);
void globalToLocal(ParticleAttrib<Vector_t> & vectorArray,
double phi, Vector_t const translationToGlobal = 0);
// Transform the x- and y-parts of a particle attribute (position, momentum, fields) from the
// local reference frame to the global reference frame.
void localToGlobal(ParticleAttrib<Vector_t> & vectorArray, double phi, Vector_t const translationToGlobal = 0);
void localToGlobal(ParticleAttrib<Vector_t> & vectorArray,
double phi, Vector_t const translationToGlobal = 0);
// Overloaded version of globalToLocal using a quaternion for 3D rotation
inline void globalToLocal(ParticleAttrib<Vector_t> & vectorArray, Quaternion_t const quaternion, Vector_t const meanR = Vector_t(0.0));
inline void globalToLocal(ParticleAttrib<Vector_t> & vectorArray,
Quaternion_t const quaternion,
Vector_t const meanR = Vector_t(0.0));
// Overloaded version of localToGlobal using a quaternion for 3D rotation
inline void localToGlobal(ParticleAttrib<Vector_t> & vectorArray, Quaternion_t const quaternion, Vector_t const meanR = Vector_t(0.0));
inline void localToGlobal(ParticleAttrib<Vector_t> & vectorArray,
Quaternion_t const quaternion,
Vector_t const meanR = Vector_t(0.0));
// Overloaded version of globalToLocal using phi and theta for pseudo 3D rotation
inline void globalToLocal(ParticleAttrib<Vector_t> & particleVectors, double const phi, double const psi, Vector_t const meanR = Vector_t(0.0));
inline void globalToLocal(ParticleAttrib<Vector_t> & particleVectors,
double const phi, double const psi,
Vector_t const meanR = Vector_t(0.0));
// Overloaded version of localToGlobal using phi and theta for pseudo 3D rotation
inline void localToGlobal(ParticleAttrib<Vector_t> & particleVectors, double const phi, double const psi, Vector_t const meanR = Vector_t(0.0));
inline void localToGlobal(ParticleAttrib<Vector_t> & particleVectors,
double const phi, double const psi,
Vector_t const meanR = Vector_t(0.0));
// Overloaded version of globalToLocal using phi and theta for pseudo 3D rotation, single vector
inline void globalToLocal(Vector_t & myVector,
double const phi, double const psi,
Vector_t const meanR = Vector_t(0.0));
// Overloaded version of localToGlobal using phi and theta for pseudo 3D rotation, single vector
inline void localToGlobal(Vector_t & myVector,
double const phi, double const psi,
Vector_t const meanR = Vector_t(0.0));
// Rotate the particles by an angle and axis defined in the quaternion (4-vector w,x,y,z)
inline void rotateWithQuaternion(ParticleAttrib<Vector_t> & vectorArray, Quaternion_t const quaternion);
......@@ -397,6 +425,8 @@ private:
// Some rotations
inline void rotateAroundZ(ParticleAttrib<Vector_t> & particleVectors, double const phi);
inline void rotateAroundX(ParticleAttrib<Vector_t> & particleVectors, double const psi);
inline void rotateAroundZ(Vector_t & myVector, double const phi);
inline void rotateAroundX(Vector_t & myVector, double const psi);
// Push particles for time h.
// Apply effects of RF Gap Crossings.
......@@ -424,7 +454,9 @@ private:
void singleParticleDump();
void bunchDumpPhaseSpaceStatData();
void bunchDumpStatData();
void bunchDumpPhaseSpaceData();
void evaluateSpaceChargeField();
......
......@@ -806,6 +806,7 @@ void Distribution::DoRestartOpalCycl(PartBunch &beam, size_t Np, int restartStep
beam.create(localN);
if(strcmp(opalFlavour, "opal-t") == 0) {
*gmsg << "Restart from hdf5 file generated by OPAL-t" << endl;
// force the initial time to zero
......@@ -854,8 +855,34 @@ void Distribution::DoRestartOpalCycl(PartBunch &beam, size_t Np, int restartStep
beam.ID[n] = larray[n];
} else {
*gmsg << "Restart from hdf5 file generated by OPAL-cycl" << endl;
rc = H5ReadStepAttribFloat64(H5file, "AZIMUTH", &phi_m);
if(rc != H5_SUCCESS)
ERRORMSG("H5 rc= " << rc << " in " << __FILE__ << " @ line " << __LINE__ << endl);
rc = H5ReadStepAttribFloat64(H5file, "ELEVATION", &psi_m);
if(rc != H5_SUCCESS)
ERRORMSG("H5 rc= " << rc << " in " << __FILE__ << " @ line " << __LINE__ << endl);
h5_int64_t localDump = 0;
previousH5Local_m = false;
rc = H5ReadStepAttribInt64(H5file, "LOCAL", &localDump);
if(rc != H5_SUCCESS) {
ERRORMSG("H5 rc= " << rc << " in " << __FILE__ << " @ line " << __LINE__ << endl);
throw OpalException("Error during restart for Cyclotron Tracker:",
"You are trying to restart from a legacy file that doesn't contain\
information on local/global frame. We are working on legacy support, but for now you have to use\
OPAL 1.3.0!");
} else {
if (localDump > 0) previousH5Local_m = true;
}
double actualT;
rc = H5ReadStepAttribFloat64(H5file, "TIME", &actualT);
if(rc != H5_SUCCESS)
......
......@@ -185,6 +185,10 @@ public:
double GetTheta() {return referenceTheta_m;}
double GetZ() {return referenceZ_m;}
double GetPhi() {return phi_m;}
double GetPsi() {return psi_m;}
bool GetPreviousH5Local() {return previousH5Local_m;}
private:
Distribution(const std::string &name, Distribution *parent);
......@@ -391,6 +395,7 @@ private:
double sigmaFall_m;
double cutoff_m;
// Cyclotron stuff
double referencePr_m;
double referencePt_m;
double referencePz_m = 0.0;
......@@ -398,6 +403,11 @@ private:
double referenceTheta_m;
double referenceZ_m = 0.0;
double bega_m;
// Cyclotron for restart in local mode
double phi_m;
double psi_m;
bool previousH5Local_m;
};
......
......@@ -256,7 +256,7 @@ void ArbitraryDomain::Compute(Vector_t hr, NDIndex<3> localId, Vector_t globalMe
//Reference Point inside the boundary for IsoDar Geo
Vector_t P0 = Vector_t(0,0,bgeom_m->getmincoords()[2]+hr[2]);
//Reference Point where the boundary geometry is cylinder
P0 = Vector_t(0,0,0);
//P0 = Vector_t(0,0,0); // Uncomment for cylinder Benchmarking -DW
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++) {
......
......@@ -930,7 +930,7 @@ void DataSink::writePhaseSpace(PartBunch &beam, Vector_t FDext[], double sposHea
int DataSink::writePhaseSpace_cycl(PartBunch &beam, Vector_t FDext[], double meanEnergy,
double refPr, double refPt, double refPz,
double refR, double refTheta, double refZ,
double azimuth, double elevation) {
double azimuth, double elevation, int localFrame) {
if (!doHDF5_m) return -1;
//if (beam.getLocalNum() == 0) return -1; //TEMP for testing -DW
......@@ -1125,6 +1125,10 @@ int DataSink::writePhaseSpace_cycl(PartBunch &beam, Vector_t FDext[], double mea
if(rc != H5_SUCCESS)
ERRORMSG("H5 rc= " << rc << " in " << __FILE__ << " @ line " << __LINE__ << endl);
rc = H5WriteStepAttribInt64(H5file_m, "LOCAL", (h5_int64_t *)&localFrame, 1);
if(rc != H5_SUCCESS)
ERRORMSG("H5 rc= " << rc << " in " << __FILE__ << " @ line " << __LINE__ << endl);
/// Write particle mass and charge per particle. (Consider making these file attributes.)
double mpart = 1.0e-9 * beam.getM();
double Q = beam.getCharge();
......
......@@ -191,7 +191,7 @@ public:
int writePhaseSpace_cycl(PartBunch &beam, Vector_t FDext[], double E,
double refPr, double refPt, double refPz,
double refR, double refTheta, double refZ,
double azimuth, double elevation);
double azimuth, double elevation, int localFrame);
/** \brief Dumps Phase Space for Envelope trakcer to H5 file.
*
......
......@@ -650,6 +650,11 @@ void TrackRun::execute() {
itsTracker->setR(dist->GetR());
itsTracker->setTheta(dist->GetTheta());
itsTracker->setZ(dist->GetZ());
// The following is for restarts in local frame
itsTracker->setPhi(dist->GetPhi());
itsTracker->setPsi(dist->GetPsi());
itsTracker->setPreviousH5Local(dist->GetPreviousH5Local());
}
// statistical data are calculated (rms, eps etc.)
......
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