OPAL (Object Oriented Parallel Accelerator Library)  2024.1
OPAL
GreenWakeFunction.cpp
Go to the documentation of this file.
1 //
2 // Class GreenWakeFunction
3 //
4 // Copyright (c) 2008 - 2020, Paul Scherrer Institut, Villigen PSI, Switzerland
5 // All rights reserved
6 //
7 // This file is part of OPAL.
8 //
9 // OPAL is free software: you can redistribute it and/or modify
10 // it under the terms of the GNU General Public License as published by
11 // the Free Software Foundation, either version 3 of the License, or
12 // (at your option) any later version.
13 //
14 // You should have received a copy of the GNU General Public License
15 // along with OPAL. If not, see <https://www.gnu.org/licenses/>.
16 //
18 
21 
22 #include "gsl/gsl_fft_real.h"
23 #include "gsl/gsl_fft_halfcomplex.h"
24 
25 #include <fstream>
26 #include <istream>
27 #include <iostream> // Needed for stream I/O
28 #include <iomanip> // Needed for I/O manipulators
29 
30 //IFF: TEST
31 //#define ENABLE_WAKE_DEBUG
32 //#define ENABLE_WAKE_DUMP
33 //#define ENABLE_WAKE_TESTS_FFT_OUT
34 
35 
36 const std::map<WakeDirection, std::string> GreenWakeFunction::wakeDirectiontoString_s = {
37  {WakeDirection::TRANSVERSAL, "TRANSVERSAL"},
38  {WakeDirection::LONGITUDINAL, "LONGITUDINAL"},
39 };
40 
57  std::vector<Filter *> filters,
58  int NBIN,
59  double Z0,
60  double radius,
61  double sigma,
62  int acMode,
63  double tau,
64  WakeDirection direction,
65  bool constLength,
66  std::string fname):
67  WakeFunction(name, NBIN),
68  lineDensity_m(),
69  //~ FftWField_m(0),
70  NBin_m(NBIN),
71  Z0_m(Z0),
72  radius_m(radius),
73  sigma_m(sigma),
74  acMode_m(acMode),
75  tau_m(tau),
76  direction_m(direction),
77  constLength_m(constLength),
78  filename_m(fname),
79  filters_m(filters.begin(), filters.end()) {
80 #ifdef ENABLE_WAKE_DEBUG
81  *gmsg << "* ************* W A K E ************************************************************ " << endl;
82  *gmsg << "* Entered GreenWakeFunction::GreenWakeFunction " << '\n';
83  *gmsg << "* ********************************************************************************** " << endl;
84 #endif
85 }
86 
88  //~ if (FftWField_m != 0) {
89  //~ delete[] FftWField_m;
90  //~ }
91 }
92 
93 
104 std::pair<int, int> GreenWakeFunction::distrIndices(int vectLen) {
105 
106  std::pair<int, int> dist;
107 
108  //IFF: properly distribute remainder
109  int rem = vectLen - (vectLen / Ippl::getNodes()) * Ippl::getNodes();
110  int tmp = (rem > Ippl::myNode()) ? 1 : 0;
111  int locBunchRange = vectLen / Ippl::getNodes() + tmp;
112 
113  dist.first = locBunchRange * Ippl::myNode() + (1 - tmp) * rem;
114  dist.second = dist.first + locBunchRange - 1;
115 
116  return dist;
117 }
118 
120 
121  Vector_t rmin, rmax;
122  double charge = bunch->getChargePerParticle();
123  // CKR: was bunch->Q[1] changed it;
124  //FIXME: why 1? bunch,getTotalCharge()
125  // or bunch->getChargePerParticle()?
126  double K = 0; // constant to normalize the lineDensity_m to 1
127  double spacing, mindist;
128  std::vector<double> OutEnergy(NBin_m);
129 
130  bunch->calcBeamParameters();
131  bunch->get_bounds(rmin, rmax);
132  //FIXME IFF: do we have unitless r's here? is that what we want?
133 
134  mindist = rmin(2);
135  switch(direction_m) {
137  spacing = std::abs(rmax(2) - rmin(2));
138  break; //FIXME: Kann mann das Spacing immer ändern?
139  }
141  spacing = rmax(0) * rmax(0) + rmax(1) * rmax(1);
142  break;
143  }
144  default: {
145  throw GeneralClassicException("GreenWakeFunction::apply",
146  "Invalid direction specified");
147  }
148  }
149  PAssert(NBin_m > 0);
150  spacing /= (NBin_m - 1); //FIXME: why -1? CKR: because grid spacings = grid points - 1
151 
152  // Calculate the Wakefield if needed
153  if (FftWField_m.empty()) {
154  FftWField_m.resize(2*NBin_m-1);
155  if (!filename_m.empty()) {
156  setWakeFromFile(NBin_m, spacing);
157  } else {
158  CalcWakeFFT(spacing);
159  }
160  } else if (!constLength_m) {
161  CalcWakeFFT(spacing);
162  }
163 
164  // Calculate the line density of the particle bunch
165  std::pair<double, double> meshInfo;
166  bunch->calcLineDensity(nBins_m, lineDensity_m, meshInfo);
167 
168 #ifdef ENABLE_WAKE_DEBUG
169  *gmsg << "* ************* W A K E ************************************************************ " << endl;
170  *gmsg << "* GreenWakeFunction::apply lineDensity_m.size() = " << lineDensity_m.size() << endl;
171  *gmsg << "* ********************************************************************************** " << endl;
172 #endif
173 
174  // smooth the line density of the particle bunch
175  for (std::vector<Filter *>::const_iterator fit = filters_m.begin(); fit != filters_m.end(); ++fit) {
176  (*fit)->apply(lineDensity_m);
177  }
178 
179  for (unsigned int i = 0; i < lineDensity_m.size(); i++) {
180  K += lineDensity_m[i];
181  }
182  K = 1 / K;
183 
184  // compute the kick due to the wakefield
185  compEnergy(K, charge, lineDensity_m, OutEnergy.data());
186 
187  // Add the right OutEnergy[i] to all the particles
188  //FIXME: can we specify LONG AND TRANS?
189  switch(direction_m) {
191  for (unsigned int i = 0; i < bunch->getLocalNum(); i++) {
192 
193  //FIXME: Stimmt das????????? (von den einheiten)
194  // calculate bin containing particle
195  int idx = (int)(floor((bunch->R[i](2) - mindist) / spacing));
196  //IFF: should be ok
197  if (idx == NBin_m) idx--;
198  PAssert(idx >= 0 && idx < NBin_m);
199  double dE = OutEnergy[idx];
200  bunch->Ef[i](2) += dE;
201 
202  }
203  break;
204  }
206  for (unsigned int i = 0; i < bunch->getLocalNum(); i++) {
207 
208  // calculate bin containing particle
209  int idx = (int)(floor((bunch->R[i](2) - mindist) / spacing));
210  //IFF: should be ok
211  if (idx == NBin_m) idx--;
212  PAssert(idx >= 0 && idx < NBin_m);
213  double dE = OutEnergy[idx];
214 
215  // ACHTUNG spacing auch in transversal richtung
216  double dist = sqrt(bunch->R[i](0) * bunch->R[i](0) + bunch->R[i](1) * bunch->R[i](1));
217  PAssert(dist > 0);
218  bunch->Ef[i](0) += dE * bunch->R[i](0) / dist;
219  bunch->Ef[i](1) += dE * bunch->R[i](1) / dist;
220 
221  }
222  break;
223  }
224  default: {
225  throw GeneralClassicException("GreenWakeFunction::apply",
226  "Invalid direction specified");
227  }
228  }
229 
230 #ifdef ENABLE_WAKE_DUMP
231  ofstream f2("OutEnergy.dat");
232  f2 << "# Energy of the Wake calculated in Opal\n"
233  << "# Z0 = " << Z0_m << "\n"
234  << "# radius = " << radius << "\n"
235  << "# sigma = " << sigma << "\n"
236  << "# c = " << c << "\n"
237  << "# acMode = " << acMode << "\n"
238  << "# tau = " << tau << "\n"
239  << "# direction = " << direction << "\n"
240  << "# spacing = " << spacing << "\n"
241  << "# Lbunch = " << NBin_m << "\n";
242  for (int i = 0; i < NBin_m; i++) {
243  f2 << i + 1 << " " << OutEnergy[i] << "\n";
244  }
245  f2.flush();
246  f2.close();
247 #endif
248 }
249 
259 void GreenWakeFunction::compEnergy(const double K,
260  const double charge,
261  const double* lambda,
262  double* OutEnergy) {
263  int N = 2 * NBin_m - 1;
264  // Allocate Space for the zero padded lambda and its Fourier Transformed
265  std::vector<double> pLambda(N);
266 
267  gsl_fft_halfcomplex_wavetable *hc;
268  gsl_fft_real_wavetable *real = gsl_fft_real_wavetable_alloc(N);
269  gsl_fft_real_workspace *work = gsl_fft_real_workspace_alloc(N);
270 
271  // fill the arrays with data
272  for (int i = 0; i < NBin_m; i ++) {
273  pLambda[N-i-1] = 0.0;
274  pLambda[i] = lambda[i];
275  }
276 
277  //FFT of the lambda
278  gsl_fft_real_transform(pLambda.data(), 1, N, real, work);
279  gsl_fft_real_wavetable_free(real);
280 
281  // convolution -> multiplication in Fourier space
282  pLambda[0] *= FftWField_m[0];
283  for (int i = 1; i < N; i += 2) {
284  double temp = pLambda[i];
285  pLambda[i] = FftWField_m[i] * pLambda[i] - FftWField_m[i+1] * pLambda[i+1];
286  pLambda[i+1] = FftWField_m[i] * pLambda[i+1] + FftWField_m[i+1] * temp;
287  }
288 
289  // inverse transform to get c, the convolution of a and b;
290  hc = gsl_fft_halfcomplex_wavetable_alloc(N);
291 
292  gsl_fft_halfcomplex_inverse(pLambda.data(), 1, N, hc, work);
293 
294  // Write the result to the output:
295  for (int i = 0; i < NBin_m; i ++) {
296  OutEnergy[i] = -charge * K * pLambda[i] / (2.0 * NBin_m) * N; // CKR: I doubt that the multiplication with N is correct,
297  // put it here to get the same result as with FFTW
298  // My suspicion: S. Pauli has forgotten that if you
299  // do an fft followed by an inverse fft you'll get
300  // N times your original data
301 
302  }
303 
304 
305  gsl_fft_halfcomplex_wavetable_free(hc);
306  gsl_fft_real_workspace_free(work);
307 }
308 
319 void GreenWakeFunction::compEnergy(const double K,
320  const double charge,
321  std::vector<double> lambda,
322  double* OutEnergy) {
323  int N = 2 * NBin_m - 1;
324  // Allocate Space for the zero padded lambda and its Fourier Transformed
325  std::vector<double> pLambda(N);
326 
327  gsl_fft_halfcomplex_wavetable *hc;
328  gsl_fft_real_wavetable *real = gsl_fft_real_wavetable_alloc(N);
329  gsl_fft_real_workspace *work = gsl_fft_real_workspace_alloc(N);
330 
331  // fill the arrays with data
332  for (int i = 0; i < NBin_m; i ++) {
333  pLambda[N-i-1] = 0.0;
334  pLambda[i] = lambda[i];
335  }
336 
337  //FFT of the lambda
338  gsl_fft_real_transform(pLambda.data(), 1, N, real, work);
339  gsl_fft_real_wavetable_free(real);
340 
341 
342  // Convolution -> just a multiplication in Fourier space
343  pLambda[0] *= FftWField_m[0];
344  for (int i = 1; i < N; i += 2) {
345  double temp = pLambda[i];
346  pLambda[i] = FftWField_m[i] * pLambda[i] - FftWField_m[i+1] * pLambda[i+1];
347  pLambda[i+1] = FftWField_m[i] * pLambda[i+1] + FftWField_m[i+1] * temp;
348  }
349 
350  // IFFT
351  hc = gsl_fft_halfcomplex_wavetable_alloc(N);
352 
353  gsl_fft_halfcomplex_inverse(pLambda.data(), 1, N, hc, work);
354 
355  // Write the result to the output:
356  for (int i = 0; i < NBin_m; i ++) {
357  OutEnergy[i] = -charge * K * pLambda[i] / (2.0 * NBin_m) * N; // CKR: I doubt that the multiplication with N is correct,
358  // put it here to get the same result as with FFTW
359  // My suspicion: S. Pauli has forgotten that if you
360  // do an fft followed by an inverse fft you'll get
361  // N times your original data
362 
363  }
364 
365 
366  gsl_fft_halfcomplex_wavetable_free(hc);
367  gsl_fft_real_workspace_free(work);
368 }
369 
376 void GreenWakeFunction::CalcWakeFFT(double spacing) {
377  // Set integration properties
378  double a = 1, b = 1000000;
379  unsigned int N = 1000000;
380  int M = 2 * NBin_m - 1;
381 
382  gsl_fft_real_wavetable *real = gsl_fft_real_wavetable_alloc(M);
383  gsl_fft_real_workspace *work = gsl_fft_real_workspace_alloc(M);
384 
385  const std::pair<int, int> myDist = distrIndices(NBin_m);
386  const int lowIndex = myDist.first;
387  const int hiIndex = myDist.second;
388 
389  for (int i = 0; i < M; i ++) {
390  FftWField_m[i] = 0.0;
391  }
392 
397  // if (Ippl::myNode() != Ippl::getNodes()-1) {
398  for (int i = lowIndex; i <= hiIndex; i ++) {
399  Wake w(i * spacing, Z0_m, radius_m, sigma_m, acMode_m, tau_m, direction_m);
400  FftWField_m[i] = simpson(w, a, b, N);
401  }
402  // } else {
403  // //IFF: changed to <= with new distr
404  // for (int i = lowIndex; i <= hiIndex; i ++) {
405  // Wake w(i*spacing, Z0_m, radius_m, sigma_m, acMode_m, tau_m, direction_m);
406  // FftWField[i] = simpson(w,a,b,N);
407  // }
408  // }
409 
413  reduce(&(FftWField_m[0]), &(FftWField_m[0]) + NBin_m, &(FftWField_m[0]), OpAddAssign());
414 
415 #ifdef ENABLE_WAKE_TESTS_FFT_OUT
416  std::vector<double> wf(2*NBin_m-1);
417  for (int i = 0; i < 2 * NBin_m - 1; ++ i) {
418  wf[i] = FftWField_m[i];
419  }
420 #endif
421  // calculate the FFT of the Wakefield
422  gsl_fft_real_transform(FftWField_m.data(), 1, M, real, work);
423 
424 
425 #ifdef ENABLE_WAKE_TESTS_FFT_OUT
426  ofstream f2("FFTwake.dat");
427  f2 << "# FFT of the Wake calculated in Opal" << "\n"
428  << "# Z0 = " << Z0_m << "\n"
429  << "# radius = " << radius_m << "\n"
430  << "# sigma = " << sigma_m << "\n"
431  << "# mode = " << acMode_m << "\n"
432  << "# tau = " << tau_m << "\n"
433  << "# direction = " << direction_m << "\n"
434  << "# spacing = " << spacing << "\n"
435  << "# Lbunch = " << NBin_m << "\n";
436 
437  f2 << "0\t" << FftWField_m[0] << "\t0.0\t" << wf[0] << "\n";
438  for (int i = 1; i < M; i += 2) {
439  f2 << (i + 1) / 2 << "\t"
440  << FftWField_m[i] << "\t"
441  << FftWField_m[i + 1] << "\t"
442  << wf[(i+1)/2] << "\n";
443  }
444 
445 
446  f2.flush();
447  f2.close();
448 #endif
449  gsl_fft_real_wavetable_free(real);
450  gsl_fft_real_workspace_free(work);
451 }
452 
456 void GreenWakeFunction::setWakeFromFile(int NBin_m, double spacing) {
457  Inform msg("readSDDS ");
458  std::string name;
459  char temp[256];
460  int Np;
461  gsl_fft_real_wavetable *real;
462  gsl_fft_real_workspace *work;
463  std::ifstream fs;
464 
465  fs.open(filename_m.c_str());
466 
467  if (fs.fail()) {
468  throw GeneralClassicException("GreenWakeFunction::setWakeFromFile",
469  "Open file operation failed, please check if '"
470  + filename_m + "' really exists.");
471  }
472 
473  fs >> name;
474  msg << " SSDS1 read = " << name << endl;
475  if (name.compare("SDDS1") != 0) {
476  throw GeneralClassicException("GreenWakeFunction::setWakeFromFile",
477  " No SDDS1 File. A SDDS1 file should start with a SDDS1 String. Check file '"
478  + filename_m + "' ");
479  }
480 
481  for (int i = 0; i < 6; i++) {
482  fs.getline(temp, 256);
483  msg << "line " << i << " : " << temp << endl;
484  }
485 
486  fs >> Np;
487  msg << " header read" << endl;
488  if (Np <= 0) {
489  throw GeneralClassicException("GreenWakeFunction::setWakeFromFile",
490  " The particle number should be bigger than zero! Please check the first line of file '"
491  + filename_m + "'.");
492  }
493 
494  msg << " Np = " << Np << endl;
495  std::vector<double> wake(Np);
496  std::vector<double> dist(Np);
497 
498  // read the wakefunction
499  for (int i = 0; i < Np; i ++) {
500  if (!fs.eof()) {
501  fs >> dist[i] >> wake[i];
502  }
503  if (fs.eof()) {
504  throw GeneralClassicException("GreenWakeFunction::setWakeFromFile",
505  " End of file reached before the whole wakefield is imported, please check file '"
506  + filename_m + "'.");
507  }
508  }
509  // if needed interpolate the wake in a way that the wake form the file fits to the wake needs in the code (??)
510 
511  FftWField_m.resize(NBin_m);
512 
513  for (int i = 0; i < NBin_m; i ++) {
514  int j = 0;
515  while(dist[j] < i * spacing) {
516  j++;
517  }
518  -- j; // i * spacing should probably between dist[j] and dist[j+1]
519  // linear interpolation
520  FftWField_m[i] = wake[j] + ((wake[j+1] - wake[j]) / (dist[j+1] - dist[j]) * (i * spacing - dist[j]));
521  }
522 
523  real = gsl_fft_real_wavetable_alloc(NBin_m);
524  work = gsl_fft_real_workspace_alloc(NBin_m);
525 
526  gsl_fft_real_transform(FftWField_m.data(), 1, NBin_m, real, work);
527 
528  gsl_fft_real_wavetable_free(real);
529  gsl_fft_real_workspace_free(work);
530 }
531 
534 }
535 
537  return wakeDirectiontoString_s.at(direction);
538 }
const unsigned int nBins_m
Definition: WakeFunction.h:52
void setWakeFromFile(int NBin, double spacing)
reads in the wakefield from file
Tps< T > sqrt(const Tps< T > &x)
Square root.
Definition: TpsMath.h:91
constexpr double c
The velocity of light in m/s.
Definition: Physics.h:45
std::string filename_m
filename of the wakefield
bool constLength_m
true if the length of the particle bunch is considered as constant
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
double simpson(F &f, double a, double b, unsigned int N)
Simpson-Integration from the function f from a to b with N steps.
static int myNode()
Definition: IpplInfo.cpp:691
static std::string getWakeDirectionString(const WakeDirection &direction)
ParticleAttrib< Vector_t > Ef
void apply(PartBunchBase< double, 3 > *bunch) override
void get_bounds(Vector_t &rmin, Vector_t &rmax) const
void calcBeamParameters()
GreenWakeFunction(const std::string &name, std::vector< Filter * > filters, int NBIN, double Z0, double radius, double sigma, int acMode, double tau, WakeDirection direction, bool constLength, std::string fname)
clearpage the user may choose between constant or variable radius This model includes fringe fields begin
Definition: multipole_t.tex:6
Inform & endl(Inform &inf)
Definition: Inform.cpp:42
static int getNodes()
Definition: IpplInfo.cpp:670
FRONT * fs
Definition: hypervolume.cpp:59
size_t getLocalNum() const
Definition: Inform.h:42
int NBin_m
divides the particle bunch in NBin slices
int acMode_m
conductivity either 1=&quot;AC&quot; or 2=&quot;DC&quot;
std::vector< double > lineDensity_m
save the line Density of the particle bunch
static const std::map< WakeDirection, std::string > wakeDirectiontoString_s
double getChargePerParticle() const
get the macro particle charge
const std::string name
std::vector< double > FftWField_m
FFT of the zero padded wakefield.
std::vector< Filter * > filters_m
ParticlePos_t & R
WakeType
Definition: WakeFunction.h:28
and that you know you can do these things To protect your we need to make restrictions that forbid anyone to deny you these rights or to ask you to surrender the rights These restrictions translate to certain responsibilities for you if you distribute copies of the or if you modify it For if you distribute copies of such a whether gratis or for a you must give the recipients all the rights that you have You must make sure that receive or can get the source code And you must show them these terms so they know their rights We protect your rights with two distribute and or modify the software for each author s protection and we want to make certain that everyone understands that there is no warranty for this free software If the software is modified by someone else and passed we want its recipients to know that what they have is not the so that any problems introduced by others will not reflect on the original authors reputations any free program is threatened constantly by software patents We wish to avoid the danger that redistributors of a free program will individually obtain patent in effect making the program proprietary To prevent we have made it clear that any patent must be licensed for everyone s free use or not licensed at all The precise terms and conditions for distribution and modification follow GNU GENERAL PUBLIC LICENSE TERMS AND CONDITIONS FOR DISTRIBUTION AND MODIFICATION This License applies to any program or other work which contains a notice placed by the copyright holder saying it may be distributed under the terms of this General Public License The refers to any such program or work
Definition: LICENSE:43
double Z0_m
impedance
double radius_m
radius
void calcLineDensity(unsigned int nBins, std::vector< double > &lineDensity, std::pair< double, double > &meshInfo)
calculates the 1d line density (not normalized) and append it to a file.
#define PAssert(c)
Definition: PAssert.h:102
double sigma_m
material constant
double tau_m
material constant
void compEnergy(const double K, const double charge, const double *lambda, double *OutEnergy)
just a Testfunction! Calculate the energy of the Wakefunction with the lambda
PETE_TUTree< FnFloor, typename T::PETE_Expr_t > floor(const PETE_Expr< T > &l)
Definition: PETE.h:733
WakeDirection direction_m
direction either &quot;Longitudinal&quot; - &quot;Transversal&quot;
WakeDirection
FLieGenerator< T, N > real(const FLieGenerator< std::complex< T >, N > &)
Take real part of a complex generator.
Inform * gmsg
Definition: Main.cpp:70
void CalcWakeFFT(double spacing)
Calculate the FFT of the Wakefunction.
end
Definition: multipole_t.tex:9
virtual WakeType getType() const override
std::pair< int, int > distrIndices(int vectLen)
given a vector of length N, distribute the indexes among the available processors ...
bool reduce(Communicate &, InputIterator, InputIterator, OutputIterator, const ReduceOp &, bool *IncludeVal=0)
Definition: GlobalComm.hpp:55