OPAL (Object Oriented Parallel Accelerator Library)  2024.1
OPAL
FM1DElectroStatic_fast.cpp
Go to the documentation of this file.
2 #include "Fields/Fieldmap.hpp"
3 #include "Physics/Physics.h"
4 #include "Physics/Units.h"
6 #include "Utilities/Util.h"
7 
8 #include "gsl/gsl_fft_real.h"
9 
10 #include <fstream>
11 #include <ios>
12 
14  _Fieldmap(filename),
15  accuracy_m(0)
16 {
17 
19  onAxisField_m = nullptr;
20 
21  std::ifstream fieldFile(Filename_m.c_str());
22  if (fieldFile.good()) {
23 
24  bool parsingPassed = readFileHeader(fieldFile);
25  parsingPassed = checkFileData(fieldFile, parsingPassed);
26  fieldFile.close();
27 
28  if (!parsingPassed) {
30  zEnd_m = zBegin_m - 1.0e-3;
31  } else
33 
36 
37  } else {
39  zBegin_m = 0.0;
40  zEnd_m = -1.0e-3;
41  }
42 }
43 
45  freeMap();
46 }
47 
49 {
50  return FM1DElectroStatic_fast(new _FM1DElectroStatic_fast(filename));
51 }
52 
54  if (onAxisField_m == nullptr) {
55 
56  std::ifstream fieldFile(Filename_m.c_str());
57  stripFileHeader(fieldFile);
58 
59  onAxisField_m = new double[numberOfGridPoints_m];
60  double maxEz = readFileData(fieldFile, onAxisField_m);
61  fieldFile.close();
62 
63  std::vector<double> fourierCoefs
65  normalizeField(maxEz, fourierCoefs);
66 
67  double *onAxisFieldP = new double[numberOfGridPoints_m];
68  double *onAxisFieldPP = new double[numberOfGridPoints_m];
69  double *onAxisFieldPPP = new double[numberOfGridPoints_m];
70  computeFieldDerivatives(fourierCoefs, onAxisFieldP,
71  onAxisFieldPP, onAxisFieldPPP);
72  computeInterpolationVectors(onAxisFieldP, onAxisFieldPP,
73  onAxisFieldPPP);
74 
75  prepareForMapCheck(fourierCoefs);
76 
77  delete [] onAxisFieldP;
78  delete [] onAxisFieldPP;
79  delete [] onAxisFieldPPP;
80 
81  INFOMSG(level3 << typeset_msg("read in fieldmap '" + Filename_m + "'", "info")
82  << endl);
83  }
84 }
85 
87  if (onAxisField_m != nullptr) {
88  delete [] onAxisField_m;
89  onAxisField_m = nullptr;
90 
91  gsl_spline_free(onAxisFieldInterpolants_m);
92  gsl_spline_free(onAxisFieldPInterpolants_m);
93  gsl_spline_free(onAxisFieldPPInterpolants_m);
94  gsl_spline_free(onAxisFieldPPPInterpolants_m);
95  gsl_interp_accel_free(onAxisFieldAccel_m);
96  gsl_interp_accel_free(onAxisFieldPAccel_m);
97  gsl_interp_accel_free(onAxisFieldPPAccel_m);
98  gsl_interp_accel_free(onAxisFieldPPPAccel_m);
99  }
100 }
101 
103  Vector_t &B) const {
104 
105  std::vector<double> fieldComponents;
106  computeFieldOnAxis(R(2) - zBegin_m, fieldComponents);
107  computeFieldOffAxis(R, E, B, fieldComponents);
108 
109  return false;
110 
111 }
112 
114  Vector_t &E,
115  Vector_t &/*B*/,
116  const DiffDirection &/*dir*/) const {
117 
118  E(2) += gsl_spline_eval(onAxisFieldPInterpolants_m, R(2) - zBegin_m,
120 
121  return false;
122 
123 }
124 
125 void _FM1DElectroStatic_fast::getFieldDimensions(double &zBegin, double &zEnd) const {
126  zBegin = zBegin_m;
127  zEnd = zEnd_m;
128 }
129 void _FM1DElectroStatic_fast::getFieldDimensions(double &/*xIni*/, double &/*xFinal*/,
130  double &/*yIni*/, double &/*yFinal*/,
131  double &/*zIni*/, double &/*zFinal*/) const {}
132 
134 { }
135 
137  (*msg) << Filename_m
138  << " (1D electrotostatic); zini= "
139  << zBegin_m << " m; zfinal= "
140  << zEnd_m << " m;" << endl;
141 }
142 
144  return 0.0;
145 }
146 
148 { }
149 
150 bool _FM1DElectroStatic_fast::checkFileData(std::ifstream &fieldFile,
151  bool parsingPassed) {
152 
153  double tempDouble;
154  for (unsigned int dataIndex = 0; dataIndex < numberOfGridPoints_m; ++ dataIndex)
155  parsingPassed = parsingPassed
156  && interpretLine<double>(fieldFile, tempDouble);
157 
158  return parsingPassed && interpreteEOF(fieldFile);
159 
160 }
161 
162 void _FM1DElectroStatic_fast::computeFieldDerivatives(std::vector<double> fourierCoefs,
163  double onAxisFieldP[],
164  double onAxisFieldPP[],
165  double onAxisFieldPPP[]) {
166 
167  for (unsigned int zStepIndex = 0; zStepIndex < numberOfGridPoints_m; ++ zStepIndex) {
168 
169  double z = deltaZ_m * zStepIndex;
170  double kz = Physics::two_pi * z / length_m + Physics::pi;
171  onAxisFieldP[zStepIndex] = 0.0;
172  onAxisFieldPP[zStepIndex] = 0.0;
173  onAxisFieldPPP[zStepIndex] = 0.0;
174 
175  int coefIndex = 1;
176  for (unsigned int n = 1; n < accuracy_m; ++ n) {
177 
178  double kn = n * Physics::two_pi / length_m;
179  double coskzn = cos(kz * n);
180  double sinkzn = sin(kz * n);
181 
182  onAxisFieldP[zStepIndex] += kn * (-fourierCoefs.at(coefIndex) * sinkzn
183  - fourierCoefs.at(coefIndex + 1) * coskzn);
184 
185  double derivCoeff = pow(kn, 2.0);
186  onAxisFieldPP[zStepIndex] += derivCoeff * (-fourierCoefs.at(coefIndex) * coskzn
187  + fourierCoefs.at(coefIndex + 1) * sinkzn);
188  derivCoeff *= kn;
189  onAxisFieldPPP[zStepIndex] += derivCoeff * (fourierCoefs.at(coefIndex) * sinkzn
190  + fourierCoefs.at(coefIndex + 1) * coskzn);
191 
192  coefIndex += 2;
193  }
194  }
195 }
196 
198  std::vector<double> fieldComponents) const {
199 
200  double radiusSq = pow(R(0), 2.0) + pow(R(1), 2.0);
201  double transverseEFactor = -fieldComponents.at(1) / 2.0
202  + radiusSq * fieldComponents.at(3) / 16.0;
203 
204  E(0) += R(0) * transverseEFactor;
205  E(1) += R(1) * transverseEFactor;
206  E(2) += fieldComponents.at(0) - fieldComponents.at(2) * radiusSq / 4.0;
207 
208 }
209 
211  std::vector<double> &fieldComponents) const {
212 
213  fieldComponents.push_back(gsl_spline_eval(onAxisFieldInterpolants_m,
214  z, onAxisFieldAccel_m));
215  fieldComponents.push_back(gsl_spline_eval(onAxisFieldPInterpolants_m,
216  z, onAxisFieldPAccel_m));
217  fieldComponents.push_back(gsl_spline_eval(onAxisFieldPPInterpolants_m,
219  fieldComponents.push_back(gsl_spline_eval(onAxisFieldPPPInterpolants_m,
221 }
222 
223 std::vector<double> _FM1DElectroStatic_fast::computeFourierCoefficients(double fieldData[]) {
224 
225  const unsigned int totalSize = 2 * numberOfGridPoints_m - 1;
226  gsl_fft_real_wavetable *waveTable = gsl_fft_real_wavetable_alloc(totalSize);
227  gsl_fft_real_workspace *workSpace = gsl_fft_real_workspace_alloc(totalSize);
228 
229  // Reflect field about minimum z value to ensure that it is periodic.
230  double *fieldDataReflected = new double[totalSize];
231  for (unsigned int dataIndex = 0; dataIndex < numberOfGridPoints_m; ++ dataIndex) {
232  fieldDataReflected[numberOfGridPoints_m - 1 + dataIndex] =
233  fieldData[dataIndex];
234  fieldDataReflected[numberOfGridPoints_m - 1 - dataIndex] =
235  fieldData[dataIndex];
236  }
237 
238  gsl_fft_real_transform(fieldDataReflected, 1,
239  totalSize,
240  waveTable, workSpace);
241 
242  std::vector<double> fourierCoefs;
243  fourierCoefs.push_back(fieldDataReflected[0] / totalSize);
244  for (unsigned int coefIndex = 1; coefIndex + 1 < 2 * accuracy_m; ++ coefIndex)
245  fourierCoefs.push_back(2.0 * fieldDataReflected[coefIndex] / totalSize);
246 
247  delete [] fieldDataReflected;
248  gsl_fft_real_workspace_free(workSpace);
249  gsl_fft_real_wavetable_free(waveTable);
250 
251  return fourierCoefs;
252 
253 }
254 
256  double onAxisFieldPP[],
257  double onAxisFieldPPP[]) {
258 
259  onAxisFieldInterpolants_m = gsl_spline_alloc(gsl_interp_cspline,
261  onAxisFieldPInterpolants_m = gsl_spline_alloc(gsl_interp_cspline,
263  onAxisFieldPPInterpolants_m = gsl_spline_alloc(gsl_interp_cspline,
265  onAxisFieldPPPInterpolants_m = gsl_spline_alloc(gsl_interp_cspline,
267 
268  double *z = new double[numberOfGridPoints_m];
269  for (unsigned int zStepIndex = 0; zStepIndex < numberOfGridPoints_m; ++ zStepIndex)
270  z[zStepIndex] = deltaZ_m * zStepIndex;
271 
272  gsl_spline_init(onAxisFieldInterpolants_m, z,
273  onAxisField_m, numberOfGridPoints_m);
274  gsl_spline_init(onAxisFieldPInterpolants_m, z,
275  onAxisFieldP, numberOfGridPoints_m);
276  gsl_spline_init(onAxisFieldPPInterpolants_m, z,
277  onAxisFieldPP, numberOfGridPoints_m);
278  gsl_spline_init(onAxisFieldPPPInterpolants_m, z,
279  onAxisFieldPPP, numberOfGridPoints_m);
280 
281  onAxisFieldAccel_m = gsl_interp_accel_alloc();
282  onAxisFieldPAccel_m = gsl_interp_accel_alloc();
283  onAxisFieldPPAccel_m = gsl_interp_accel_alloc();
284  onAxisFieldPPPAccel_m = gsl_interp_accel_alloc();
285 
286  delete [] z;
287 }
288 
290 
291  // Convert to m.
293  rEnd_m *= Units::cm2m;
295  zEnd_m *= Units::cm2m;
296 
297 }
298 
299 void _FM1DElectroStatic_fast::normalizeField(double maxEz, std::vector<double> &fourierCoefs) {
300 
301  for (unsigned int dataIndex = 0; dataIndex < numberOfGridPoints_m; ++ dataIndex)
302  onAxisField_m[dataIndex] /= (maxEz * Units::Vpm2MVpm);
303 
304  for (auto fourierIt = fourierCoefs.begin(); fourierIt < fourierCoefs.end(); ++ fourierIt)
305  *fourierIt /= (maxEz * Units::Vpm2MVpm);
306 
307 }
308 
309 double _FM1DElectroStatic_fast::readFileData(std::ifstream &fieldFile,
310  double fieldData[]) {
311 
312  double maxEz = 0.0;
313  for (unsigned int dataIndex = 0; dataIndex < numberOfGridPoints_m; ++ dataIndex) {
314  interpretLine<double>(fieldFile, fieldData[dataIndex]);
315  if (std::abs(fieldData[dataIndex]) > maxEz)
316  maxEz = std::abs(fieldData[dataIndex]);
317  }
318 
319  if (!normalize_m)
320  maxEz = 1.0;
321 
322  return maxEz;
323 }
324 
325 bool _FM1DElectroStatic_fast::readFileHeader(std::ifstream &fieldFile) {
326 
327  std::string tempString;
328  int tempInt;
329 
330  bool parsingPassed = true;
331  try {
332  parsingPassed = interpretLine<std::string, unsigned int>(fieldFile,
333  tempString,
334  accuracy_m);
335  } catch (GeneralClassicException &e) {
336  parsingPassed = interpretLine<std::string, unsigned int, std::string>(fieldFile,
337  tempString,
338  accuracy_m,
339  tempString);
340 
341  tempString = Util::toUpper(tempString);
342  if (tempString != "TRUE" &&
343  tempString != "FALSE")
344  throw GeneralClassicException("_FM1DElectroStatic_fast::readFileHeader",
345  "The third string on the first line of 1D field "
346  "maps has to be either TRUE or FALSE");
347 
348  normalize_m = (tempString == "TRUE");
349  }
350  parsingPassed = parsingPassed &&
351  interpretLine<double, double, unsigned int>(fieldFile,
352  zBegin_m,
353  zEnd_m,
355  parsingPassed = parsingPassed &&
356  interpretLine<double, double, int>(fieldFile,
357  rBegin_m,
358  rEnd_m,
359  tempInt);
360 
362  if (accuracy_m > numberOfGridPoints_m)
364 
365  return parsingPassed;
366 }
367 
368 void _FM1DElectroStatic_fast::stripFileHeader(std::ifstream &fieldFile) {
369 
370  std::string tempString;
371 
372  getLine(fieldFile, tempString);
373  getLine(fieldFile, tempString);
374  getLine(fieldFile, tempString);
375 }
376 
377 void _FM1DElectroStatic_fast::prepareForMapCheck(std::vector<double> &fourierCoefs) {
378  std::vector<double> zSampling(numberOfGridPoints_m);
379  for (unsigned int zStepIndex = 0; zStepIndex < numberOfGridPoints_m; ++ zStepIndex)
380  zSampling[zStepIndex] = deltaZ_m * zStepIndex;
381 
383  length_m,
384  zSampling,
385  fourierCoefs,
388 }
virtual bool getFieldDerivative(const Vector_t &R, Vector_t &E, Vector_t &B, const DiffDirection &dir) const
void computeInterpolationVectors(double onAxisFieldP[], double onAxisFieldPP[], double onAxisFieldPPP[])
MapType Type
Definition: Fieldmap.h:115
void computeFieldDerivatives(std::vector< double > fourierCoefs, double onAxisFieldP[], double onAxisFieldPP[], double onAxisFieldPPP[])
bool readFileHeader(std::ifstream &fieldFile)
gsl_spline * onAxisFieldPPInterpolants_m
On axis field first derivative interpolation structure.
DiffDirection
Definition: Fieldmap.h:55
double zBegin_m
Maximum radius of field.
virtual bool getFieldstrength(const Vector_t &R, Vector_t &E, Vector_t &B) const
constexpr double two_pi
The value of .
Definition: Physics.h:33
void noFieldmapWarning()
Definition: Fieldmap.cpp:618
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
constexpr double cm2m
Definition: Units.h:35
_FM1DElectroStatic_fast(const std::string &filename)
gsl_spline * onAxisFieldInterpolants_m
On axis field data.
void computeFieldOffAxis(const Vector_t &R, Vector_t &E, Vector_t &B, std::vector< double > fieldComponents) const
bool checkFileData(std::ifstream &fieldFile, bool parsingPassed)
static std::string typeset_msg(const std::string &msg, const std::string &title)
Definition: Fieldmap.cpp:649
constexpr double pi
The value of .
Definition: Physics.h:30
std::string toUpper(const std::string &str)
Definition: Util.cpp:147
constexpr double Vpm2MVpm
Definition: Units.h:125
#define INFOMSG(msg)
Definition: IpplInfo.h:348
Inform & endl(Inform &inf)
Definition: Inform.cpp:42
gsl_interp_accel * onAxisFieldPPAccel_m
bool interpreteEOF(std::ifstream &in)
Definition: Fieldmap.cpp:555
double length_m
Longitudinal end of field.
virtual void getInfo(Inform *)
virtual double getFrequency() const
bool normalize_m
Definition: Fieldmap.h:121
virtual void setFrequency(double freq)
double rEnd_m
Minimum radius of field.
static FM1DElectroStatic_fast create(const std::string &filename)
Definition: Inform.h:42
unsigned int accuracy_m
Field grid point spacing.
std::vector< double > computeFourierCoefficients(double fieldData[])
void prepareForMapCheck(std::vector< double > &fourierCoefs)
double deltaZ_m
Number of grid points in field input file.
void disableFieldmapWarning()
Definition: Fieldmap.cpp:610
Tps< T > cos(const Tps< T > &x)
Cosine.
Definition: TpsMath.h:129
double readFileData(std::ifstream &fieldFile, double fieldData[])
void getLine(std::ifstream &in, std::string &buffer)
Definition: Fieldmap.h:122
unsigned int numberOfGridPoints_m
Field length.
void stripFileHeader(std::ifstream &fieldFile)
gsl_interp_accel * onAxisFieldPPPAccel_m
Inform & level3(Inform &inf)
Definition: Inform.cpp:47
constexpr double e
The value of .
Definition: Physics.h:39
std::string Filename_m
Definition: Fieldmap.h:118
virtual void getFieldDimensions(double &zBegin, double &zEnd) const
gsl_spline * onAxisFieldPPPInterpolants_m
On axis field second derivative interpolation structure.
gsl_spline * onAxisFieldPInterpolants_m
On axis field interpolation structure.
std::shared_ptr< _FM1DElectroStatic_fast > FM1DElectroStatic_fast
Definition: Definitions.h:36
Tps< T > pow(const Tps< T > &x, int y)
Integer power.
Definition: TpsMath.h:76
Tps< T > sin(const Tps< T > &x)
Sine.
Definition: TpsMath.h:111
void checkMap(unsigned int accuracy, std::pair< double, double > fieldDimensions, double deltaZ, const std::vector< double > &fourierCoefficients, gsl_spline *splineCoefficients, gsl_interp_accel *splineAccelerator)
Definition: Fieldmap.cpp:448
double zEnd_m
Longitudinal start of field.
void computeFieldOnAxis(double z, std::vector< double > &fieldComponents) const
gsl_interp_accel * onAxisFieldPAccel_m
void normalizeField(double maxEz, std::vector< double > &fourierCoefs)
gsl_interp_accel * onAxisFieldAccel_m
On axis field third derivative interpolation structure.