OPAL (Object Oriented Parallel Accelerator Library)  2.2.0
OPAL
PartBunch.cpp
Go to the documentation of this file.
1 // ------------------------------------------------------------------------
2 // $RCSfile: PartBunch.cpp,v $
3 // ------------------------------------------------------------------------
4 // $Revision: 1.1.1.1.2.1 $
5 // ------------------------------------------------------------------------
6 // Copyright: see Copyright.readme
7 // ------------------------------------------------------------------------
8 //
9 // Class PartBunch
10 // Interface to a particle bunch.
11 // Can be used to avoid use of a template in user code.
12 //
13 // ------------------------------------------------------------------------
14 // Class category: Algorithms
15 // ------------------------------------------------------------------------
16 //
17 // $Date: 2004/11/12 18:57:53 $
18 // $Author: adelmann $
19 //
20 // ------------------------------------------------------------------------
21 
22 #include "Algorithms/PartBunch.h"
23 #include "FixedAlgebra/FMatrix.h"
24 #include "FixedAlgebra/FVector.h"
25 #include <iostream>
26 #include <cfloat>
27 #include <fstream>
28 #include <iomanip>
29 #include <memory>
30 #include <utility>
31 
32 
33 #include "Distribution/Distribution.h" // OPAL file
34 #include "Structure/FieldSolver.h" // OPAL file
35 #include "Utilities/GeneralClassicException.h"
36 
37 #include "Algorithms/ListElem.h"
38 
39 #include <gsl/gsl_rng.h>
40 #include <gsl/gsl_histogram.h>
41 #include <gsl/gsl_cdf.h>
42 #include <gsl/gsl_randist.h>
43 #include <gsl/gsl_sf_erf.h>
44 #include <gsl/gsl_qrng.h>
45 
46 #include <boost/format.hpp>
47 
48 //#define DBG_SCALARFIELD
49 //#define FIELDSTDOUT
50 
51 // Class PartBunch
52 // ------------------------------------------------------------------------
53 
54 PartBunch::PartBunch(const PartData *ref): // Layout is set using setSolver()
55  PartBunchBase<double, 3>(new PartBunch::pbase_t(new Layout_t()), ref),
56  interpolationCacheSet_m(false)
57 {
58 
59 }
60 
61 
62 PartBunch::~PartBunch() {
63 
64 }
65 
66 // PartBunch::pbase_t* PartBunch::clone() {
67 // return new pbase_t(new Layout_t());
68 // }
69 
70 
71 void PartBunch::initialize(FieldLayout_t *fLayout) {
72  Layout_t* layout = static_cast<Layout_t*>(&getLayout());
73  layout->getLayout().changeDomain(*fLayout);
74 }
75 
76 void PartBunch::runTests() {
77 
78  Vector_t ll(-0.005);
79  Vector_t ur(0.005);
80 
81  setBCAllPeriodic();
82 
83  NDIndex<3> domain = getFieldLayout().getDomain();
84  for (unsigned int i = 0; i < Dimension; i++)
85  nr_m[i] = domain[i].length();
86 
87  for (int i = 0; i < 3; i++)
88  hr_m[i] = (ur[i] - ll[i]) / nr_m[i];
89 
90  getMesh().set_meshSpacing(&(hr_m[0]));
91  getMesh().set_origin(ll);
92 
93  rho_m.initialize(getMesh(),
94  getFieldLayout(),
95  GuardCellSizes<Dimension>(1),
96  bc_m);
97  eg_m.initialize(getMesh(),
98  getFieldLayout(),
99  GuardCellSizes<Dimension>(1),
100  vbc_m);
101 
102  fs_m->solver_m->test(this);
103 }
104 
105 
106 void PartBunch::do_binaryRepart() {
107  get_bounds(rmin_m, rmax_m);
108 
109  pbase_t* underlyingPbase =
110  dynamic_cast<pbase_t*>(pbase.get());
111 
112  BinaryRepartition(*underlyingPbase);
113  update();
114  get_bounds(rmin_m, rmax_m);
115  boundp();
116 }
117 
118 
119 void PartBunch::computeSelfFields(int binNumber) {
120  IpplTimings::startTimer(selfFieldTimer_m);
121 
124  rho_m = 0.0;
125  Field_t imagePotential = rho_m;
126 
128  eg_m = Vector_t(0.0);
129 
130  if(fs_m->hasValidSolver()) {
132  if(fs_m->getFieldSolverType() == "SAAMG")
133  resizeMesh();
134 
136  this->Q *= this->dt;
137  if(!interpolationCacheSet_m) {
138  if(interpolationCache_m.size() < getLocalNum()) {
139  interpolationCache_m.create(getLocalNum() - interpolationCache_m.size());
140  } else {
141  interpolationCache_m.destroy(interpolationCache_m.size() - getLocalNum(),
142  getLocalNum(),
143  true);
144  }
145  interpolationCacheSet_m = true;
146  this->Q.scatter(this->rho_m, this->R, IntrplCIC_t(), interpolationCache_m);
147  } else {
148  this->Q.scatter(this->rho_m, IntrplCIC_t(), interpolationCache_m);
149  }
150 
151  this->Q /= this->dt;
152  this->rho_m /= getdT();
153 
155  double scaleFactor = 1;
156  // double scaleFactor = Physics::c * getdT();
157  double gammaz = getBinGamma(binNumber);
158 
161  double tmp2 = 1 / hr_m[0] * 1 / hr_m[1] * 1 / hr_m[2] / (scaleFactor * scaleFactor * scaleFactor) / gammaz;
162  rho_m *= tmp2;
163 
166  Vector_t hr_scaled = hr_m * Vector_t(scaleFactor);
167  hr_scaled[2] *= gammaz;
168 
171  imagePotential = rho_m;
172 
173  fs_m->solver_m->computePotential(rho_m, hr_scaled);
174 
176  rho_m *= hr_scaled[0] * hr_scaled[1] * hr_scaled[2];
177 
180  rho_m *= getCouplingConstant();
181 
184  eg_m = -Grad(rho_m, eg_m);
185 
188  eg_m *= Vector_t(gammaz / (scaleFactor), gammaz / (scaleFactor), 1.0 / (scaleFactor * gammaz));
189 
190  // If desired write E-field and potential to terminal
191 #ifdef FIELDSTDOUT
192  // Immediate debug output:
193  // Output potential and e-field along the x-, y-, and z-axes
194  int mx = (int)nr_m[0];
195  int mx2 = (int)nr_m[0] / 2;
196  int my = (int)nr_m[1];
197  int my2 = (int)nr_m[1] / 2;
198  int mz = (int)nr_m[2];
199  int mz2 = (int)nr_m[2] / 2;
200 
201  for (int i=0; i<mx; i++ )
202  *gmsg << "Bin " << binNumber
203  << ", Self Field along x axis E = " << eg_m[i][my2][mz2]
204  << ", Pot = " << rho_m[i][my2][mz2] << endl;
205 
206  for (int i=0; i<my; i++ )
207  *gmsg << "Bin " << binNumber
208  << ", Self Field along y axis E = " << eg_m[mx2][i][mz2]
209  << ", Pot = " << rho_m[mx2][i][mz2] << endl;
210 
211  for (int i=0; i<mz; i++ )
212  *gmsg << "Bin " << binNumber
213  << ", Self Field along z axis E = " << eg_m[mx2][my2][i]
214  << ", Pot = " << rho_m[mx2][my2][i] << endl;
215 #endif
216 
221  Eftmp.gather(eg_m, IntrplCIC_t(), interpolationCache_m);
222  //Eftmp.gather(eg_m, this->R, IntrplCIC_t());
223 
231  double betaC = sqrt(gammaz * gammaz - 1.0) / gammaz / Physics::c;
232 
233  Bf(0) = Bf(0) - betaC * Eftmp(1);
234  Bf(1) = Bf(1) + betaC * Eftmp(0);
235 
236  Ef += Eftmp;
237 
240 
242  NDIndex<3> domain = getFieldLayout().getDomain();
243  Vector_t origin = rho_m.get_mesh().get_origin();
244  double hz = rho_m.get_mesh().get_meshSpacing(2);
245  double zshift = -(2 * origin(2) + (domain[2].first() + domain[2].last() + 1) * hz) * gammaz * scaleFactor;
246 
249  fs_m->solver_m->computePotential(imagePotential, hr_scaled, zshift);
250 
252  imagePotential *= hr_scaled[0] * hr_scaled[1] * hr_scaled[2];
253 
256  imagePotential *= getCouplingConstant();
257 
258  const int dumpFreq = 100;
259 #ifdef DBG_SCALARFIELD
260  VField_t tmp_eg = eg_m;
261 
262 
263  if (Ippl::getNodes() == 1 && (fieldDBGStep_m + 1) % dumpFreq == 0) {
264 #else
265  VField_t tmp_eg;
266 
267  if (false) {
268 #endif
269  INFOMSG(level1 << "*** START DUMPING SCALAR FIELD ***" << endl);
270 
271 
272  std::string SfileName = OpalData::getInstance()->getInputBasename();
273  boost::format phi_fn("data/%1%-phi_scalar-%|2$05|.dat");
274  phi_fn % SfileName % (fieldDBGStep_m / dumpFreq);
275 
276  std::ofstream fstr2(phi_fn.str());
277  fstr2.precision(9);
278 
279  NDIndex<3> myidx = getFieldLayout().getLocalNDIndex();
280  Vector_t origin = rho_m.get_mesh().get_origin();
281  Vector_t spacing(rho_m.get_mesh().get_meshSpacing(0),
282  rho_m.get_mesh().get_meshSpacing(1),
283  rho_m.get_mesh().get_meshSpacing(2));
284  Vector_t rmin, rmax;
285  get_bounds(rmin, rmax);
286 
287  INFOMSG(level1
288  << (rmin(0) - origin(0)) / spacing(0) << "\t"
289  << (rmin(1) - origin(1)) / spacing(1) << "\t"
290  << (rmin(2) - origin(2)) / spacing(2) << "\t"
291  << rmin(2) << endl);
292  for (int x = myidx[0].first(); x <= myidx[0].last(); x++) {
293  for (int y = myidx[1].first(); y <= myidx[1].last(); y++) {
294  for (int z = myidx[2].first(); z <= myidx[2].last(); z++) {
295  fstr2 << std::setw(5) << x + 1
296  << std::setw(5) << y + 1
297  << std::setw(5) << z + 1
298  << std::setw(17) << origin(2) + z * spacing(2)
299  << std::setw(17) << rho_m[x][y][z].get()
300  << std::setw(17) << imagePotential[x][y][z].get() << std::endl;
301  }
302  }
303  }
304  fstr2.close();
305 
306  INFOMSG(level1 << "*** FINISHED DUMPING SCALAR FIELD ***" << endl);
307  }
308 
312  eg_m = -Grad(imagePotential, eg_m);
313 
316  eg_m *= Vector_t(gammaz / (scaleFactor), gammaz / (scaleFactor), 1.0 / (scaleFactor * gammaz));
317 
318  // If desired write E-field and potential to terminal
319 #ifdef FIELDSTDOUT
320  // Immediate debug output:
321  // Output potential and e-field along the x-, y-, and z-axes
322  //int mx = (int)nr_m[0];
323  //int mx2 = (int)nr_m[0] / 2;
324  //int my = (int)nr_m[1];
325  //int my2 = (int)nr_m[1] / 2;
326  //int mz = (int)nr_m[2];
327  //int mz2 = (int)nr_m[2] / 2;
328 
329  for (int i=0; i<mx; i++ )
330  *gmsg << "Bin " << binNumber
331  << ", Image Field along x axis E = " << eg_m[i][my2][mz2]
332  << ", Pot = " << rho_m[i][my2][mz2] << endl;
333 
334  for (int i=0; i<my; i++ )
335  *gmsg << "Bin " << binNumber
336  << ", Image Field along y axis E = " << eg_m[mx2][i][mz2]
337  << ", Pot = " << rho_m[mx2][i][mz2] << endl;
338 
339  for (int i=0; i<mz; i++ )
340  *gmsg << "Bin " << binNumber
341  << ", Image Field along z axis E = " << eg_m[mx2][my2][i]
342  << ", Pot = " << rho_m[mx2][my2][i] << endl;
343 #endif
344 
345 #ifdef DBG_SCALARFIELD
346  if (Ippl::getNodes() == 1 && (fieldDBGStep_m + 1) % dumpFreq == 0) {
347 #else
348  if (false) {
349 #endif
350  INFOMSG(level1 << "*** START DUMPING E FIELD ***" << endl);
351 
352  std::string SfileName = OpalData::getInstance()->getInputBasename();
353  boost::format phi_fn("data/%1%-e_field-%|2$05|.dat");
354  phi_fn % SfileName % (fieldDBGStep_m / dumpFreq);
355 
356  std::ofstream fstr2(phi_fn.str());
357  fstr2.precision(9);
358 
359  Vector_t origin = eg_m.get_mesh().get_origin();
360  Vector_t spacing(eg_m.get_mesh().get_meshSpacing(0),
361  eg_m.get_mesh().get_meshSpacing(1),
362  eg_m.get_mesh().get_meshSpacing(2));
363 
364  NDIndex<3> myidxx = getFieldLayout().getLocalNDIndex();
365  for (int x = myidxx[0].first(); x <= myidxx[0].last(); x++) {
366  for (int y = myidxx[1].first(); y <= myidxx[1].last(); y++) {
367  for (int z = myidxx[2].first(); z <= myidxx[2].last(); z++) {
368  Vector_t ef = eg_m[x][y][z].get() + tmp_eg[x][y][z].get();
369  fstr2 << std::setw(5) << x + 1
370  << std::setw(5) << y + 1
371  << std::setw(5) << z + 1
372  << std::setw(17) << origin(2) + z * spacing(2)
373  << std::setw(17) << ef(0)
374  << std::setw(17) << ef(1)
375  << std::setw(17) << ef(2) << std::endl;
376  }
377  }
378  }
379 
380  fstr2.close();
381 
382  //MPI_File_write_shared(file, (char*)oss.str().c_str(), oss.str().length(), MPI_CHAR, &status);
383  //MPI_File_close(&file);
384 
385  INFOMSG(level1 << "*** FINISHED DUMPING E FIELD ***" << endl);
386  }
387  fieldDBGStep_m++;
388 
393  Eftmp.gather(eg_m, IntrplCIC_t(), interpolationCache_m);
394  //Eftmp.gather(eg_m, this->R, IntrplCIC_t());
395 
404  Bf(0) = Bf(0) + betaC * Eftmp(1);
405  Bf(1) = Bf(1) - betaC * Eftmp(0);
406 
407  Ef += Eftmp;
408 
409  }
410  IpplTimings::stopTimer(selfFieldTimer_m);
411 }
412 
413 void PartBunch::resizeMesh() {
414  double xmin = fs_m->solver_m->getXRangeMin();
415  double xmax = fs_m->solver_m->getXRangeMax();
416  double ymin = fs_m->solver_m->getYRangeMin();
417  double ymax = fs_m->solver_m->getYRangeMax();
418  double zmin = fs_m->solver_m->getZRangeMin();
419  double zmax = fs_m->solver_m->getZRangeMax();
420 
421  if(xmin > rmin_m[0] || xmax < rmax_m[0] ||
422  ymin > rmin_m[1] || ymax < rmax_m[1]) {
423 
424  for (unsigned int n = 0; n < getLocalNum(); n++) {
425 
426  if(R[n](0) < xmin || R[n](0) > xmax ||
427  R[n](1) < ymin || R[n](1) > ymax) {
428 
429  // delete the particle
430  INFOMSG(level2 << "destroyed particle with id=" << ID[n] << endl;);
431  destroy(1, n);
432  }
433 
434  }
435 
436  update();
437  boundp();
438  get_bounds(rmin_m, rmax_m);
439  }
440  Vector_t mymin = Vector_t(xmin, ymin , zmin);
441  Vector_t mymax = Vector_t(xmax, ymax , zmax);
442 
443  for (int i = 0; i < 3; i++)
444  hr_m[i] = (mymax[i] - mymin[i])/nr_m[i];
445 
446  getMesh().set_meshSpacing(&(hr_m[0]));
447  getMesh().set_origin(mymin);
448 
449  rho_m.initialize(getMesh(),
450  getFieldLayout(),
451  GuardCellSizes<Dimension>(1),
452  bc_m);
453  eg_m.initialize(getMesh(),
454  getFieldLayout(),
455  GuardCellSizes<Dimension>(1),
456  vbc_m);
457 
458  update();
459 
460 // setGridIsFixed();
461 }
462 
463 void PartBunch::computeSelfFields() {
464  IpplTimings::startTimer(selfFieldTimer_m);
465  rho_m = 0.0;
466  eg_m = Vector_t(0.0);
467 
468  if(fs_m->hasValidSolver()) {
469  //mesh the whole domain
470  if(fs_m->getFieldSolverType() == "SAAMG")
471  resizeMesh();
472 
473  //scatter charges onto grid
474  this->Q *= this->dt;
475  this->Q.scatter(this->rho_m, this->R, IntrplCIC_t());
476  this->Q /= this->dt;
477  this->rho_m /= getdT();
478 
479  //calculating mesh-scale factor
480  double gammaz = sum(this->P)[2] / getTotalNum();
481  gammaz *= gammaz;
482  gammaz = sqrt(gammaz + 1.0);
483  double scaleFactor = 1;
484  // double scaleFactor = Physics::c * getdT();
485  //and get meshspacings in real units [m]
486  Vector_t hr_scaled = hr_m * Vector_t(scaleFactor);
487  hr_scaled[2] *= gammaz;
488 
489  //double tmp2 = 1/hr_m[0] * 1/hr_m[1] * 1/hr_m[2] / (scaleFactor*scaleFactor*scaleFactor) / gammaz;
490  double tmp2 = 1 / hr_scaled[0] * 1 / hr_scaled[1] * 1 / hr_scaled[2];
491  //divide charge by a 'grid-cube' volume to get [C/m^3]
492  rho_m *= tmp2;
493 
494 #ifdef DBG_SCALARFIELD
495  INFOMSG("*** START DUMPING SCALAR FIELD ***" << endl);
496  std::ofstream fstr1;
497  fstr1.precision(9);
498 
499  std::string SfileName = OpalData::getInstance()->getInputBasename();
500 
501  std::string rho_fn = std::string("data/") + SfileName + std::string("-rho_scalar-") + std::to_string(fieldDBGStep_m);
502  fstr1.open(rho_fn.c_str(), std::ios::out);
503  NDIndex<3> myidx1 = getFieldLayout().getLocalNDIndex();
504  for (int x = myidx1[0].first(); x <= myidx1[0].last(); x++) {
505  for (int y = myidx1[1].first(); y <= myidx1[1].last(); y++) {
506  for (int z = myidx1[2].first(); z <= myidx1[2].last(); z++) {
507  fstr1 << x + 1 << " " << y + 1 << " " << z + 1 << " " << rho_m[x][y][z].get() << std::endl;
508  }
509  }
510  }
511  fstr1.close();
512  INFOMSG("*** FINISHED DUMPING SCALAR FIELD ***" << endl);
513 #endif
514 
515  // charge density is in rho_m
516  fs_m->solver_m->computePotential(rho_m, hr_scaled);
517 
518  //do the multiplication of the grid-cube volume coming
519  //from the discretization of the convolution integral.
520  //this is only necessary for the FFT solver
521  //FIXME: later move this scaling into FFTPoissonSolver
522  if(fs_m->getFieldSolverType() == "FFT" || fs_m->getFieldSolverType() == "FFTBOX")
523  rho_m *= hr_scaled[0] * hr_scaled[1] * hr_scaled[2];
524 
525  // the scalar potential is given back in rho_m in units
526  // [C/m] = [F*V/m] and must be divided by
527  // 4*pi*\epsilon_0 [F/m] resulting in [V]
528  rho_m *= getCouplingConstant();
529 
530  //write out rho
531 #ifdef DBG_SCALARFIELD
532  INFOMSG("*** START DUMPING SCALAR FIELD ***" << endl);
533 
534  std::ofstream fstr2;
535  fstr2.precision(9);
536 
537  std::string phi_fn = std::string("data/") + SfileName + std::string("-phi_scalar-") + std::to_string(fieldDBGStep_m);
538  fstr2.open(phi_fn.c_str(), std::ios::out);
539  NDIndex<3> myidx = getFieldLayout().getLocalNDIndex();
540  for (int x = myidx[0].first(); x <= myidx[0].last(); x++) {
541  for (int y = myidx[1].first(); y <= myidx[1].last(); y++) {
542  for (int z = myidx[2].first(); z <= myidx[2].last(); z++) {
543  fstr2 << x + 1 << " " << y + 1 << " " << z + 1 << " " << rho_m[x][y][z].get() << std::endl;
544  }
545  }
546  }
547  fstr2.close();
548 
549  INFOMSG("*** FINISHED DUMPING SCALAR FIELD ***" << endl);
550 #endif
551 
552  // IPPL Grad divides by hr_m [m] resulting in
553  // [V/m] for the electric field
554  eg_m = -Grad(rho_m, eg_m);
555 
556  eg_m *= Vector_t(gammaz / (scaleFactor), gammaz / (scaleFactor), 1.0 / (scaleFactor * gammaz));
557 
558  //write out e field
559 #ifdef FIELDSTDOUT
560  // Immediate debug output:
561  // Output potential and e-field along the x-, y-, and z-axes
562  int mx = (int)nr_m[0];
563  int mx2 = (int)nr_m[0] / 2;
564  int my = (int)nr_m[1];
565  int my2 = (int)nr_m[1] / 2;
566  int mz = (int)nr_m[2];
567  int mz2 = (int)nr_m[2] / 2;
568 
569  for (int i=0; i<mx; i++ )
570  *gmsg << "Field along x axis Ex = " << eg_m[i][my2][mz2] << " Pot = " << rho_m[i][my2][mz2] << endl;
571 
572  for (int i=0; i<my; i++ )
573  *gmsg << "Field along y axis Ey = " << eg_m[mx2][i][mz2] << " Pot = " << rho_m[mx2][i][mz2] << endl;
574 
575  for (int i=0; i<mz; i++ )
576  *gmsg << "Field along z axis Ez = " << eg_m[mx2][my2][i] << " Pot = " << rho_m[mx2][my2][i] << endl;
577 #endif
578 
579 #ifdef DBG_SCALARFIELD
580  INFOMSG("*** START DUMPING E FIELD ***" << endl);
581  //ostringstream oss;
582  //MPI_File file;
583  //MPI_Status status;
584  //MPI_Info fileinfo;
585  //MPI_File_open(Ippl::getComm(), "rho_scalar", MPI_MODE_WRONLY | MPI_MODE_CREATE, fileinfo, &file);
586  std::ofstream fstr;
587  fstr.precision(9);
588 
589  std::string e_field = std::string("data/") + SfileName + std::string("-e_field-") + std::to_string(fieldDBGStep_m);
590  fstr.open(e_field.c_str(), std::ios::out);
591  NDIndex<3> myidxx = getFieldLayout().getLocalNDIndex();
592  for (int x = myidxx[0].first(); x <= myidxx[0].last(); x++) {
593  for (int y = myidxx[1].first(); y <= myidxx[1].last(); y++) {
594  for (int z = myidxx[2].first(); z <= myidxx[2].last(); z++) {
595  fstr << x + 1 << " " << y + 1 << " " << z + 1 << " " << eg_m[x][y][z].get() << std::endl;
596  }
597  }
598  }
599 
600  fstr.close();
601  fieldDBGStep_m++;
602 
603  //MPI_File_write_shared(file, (char*)oss.str().c_str(), oss.str().length(), MPI_CHAR, &status);
604  //MPI_File_close(&file);
605 
606  INFOMSG("*** FINISHED DUMPING E FIELD ***" << endl);
607 #endif
608 
609  // interpolate electric field at particle positions. We reuse the
610  // cached information about where the particles are relative to the
611  // field, since the particles have not moved since this the most recent
612  // scatter operation.
613  Ef.gather(eg_m, this->R, IntrplCIC_t());
614 
622  double betaC = sqrt(gammaz * gammaz - 1.0) / gammaz / Physics::c;
623 
624  Bf(0) = Bf(0) - betaC * Ef(1);
625  Bf(1) = Bf(1) + betaC * Ef(0);
626  }
627  IpplTimings::stopTimer(selfFieldTimer_m);
628 }
629 
643 void PartBunch::computeSelfFields_cycl(double gamma) {
644 
645  IpplTimings::startTimer(selfFieldTimer_m);
646 
648  rho_m = 0.0;
649 
651  eg_m = Vector_t(0.0);
652 
653  if(fs_m->hasValidSolver()) {
655  if(fs_m->getFieldSolverType() == "SAAMG")
656  resizeMesh();
657 
659  this->Q.scatter(this->rho_m, this->R, IntrplCIC_t());
660 
663  Vector_t hr_scaled = hr_m ;
664  hr_scaled[1] *= gamma;
665 
667  double tmp2 = 1.0 / (hr_scaled[0] * hr_scaled[1] * hr_scaled[2]);
668  rho_m *= tmp2;
669 
670  // If debug flag is set, dump scalar field (charge density 'rho') into file under ./data/
671 #ifdef DBG_SCALARFIELD
672  INFOMSG("*** START DUMPING SCALAR FIELD ***" << endl);
673  std::ofstream fstr1;
674  fstr1.precision(9);
675 
676  std::ostringstream istr;
677  istr << fieldDBGStep_m;
678 
679  std::string SfileName = OpalData::getInstance()->getInputBasename();
680 
681  std::string rho_fn = std::string("data/") + SfileName + std::string("-rho_scalar-") + std::string(istr.str());
682  fstr1.open(rho_fn.c_str(), std::ios::out);
683  NDIndex<3> myidx1 = getFieldLayout().getLocalNDIndex();
684  for (int x = myidx1[0].first(); x <= myidx1[0].last(); x++) {
685  for (int y = myidx1[1].first(); y <= myidx1[1].last(); y++) {
686  for (int z = myidx1[2].first(); z <= myidx1[2].last(); z++) {
687  fstr1 << x + 1 << " " << y + 1 << " " << z + 1 << " " << rho_m[x][y][z].get() << std::endl;
688  }
689  }
690  }
691  fstr1.close();
692  INFOMSG("*** FINISHED DUMPING SCALAR FIELD ***" << endl);
693 #endif
694 
697  fs_m->solver_m->computePotential(rho_m, hr_scaled);
698 
699  //do the multiplication of the grid-cube volume coming
700  //from the discretization of the convolution integral.
701  //this is only necessary for the FFT solver
702  //TODO FIXME: later move this scaling into FFTPoissonSolver
703  if(fs_m->getFieldSolverType() == "FFT" || fs_m->getFieldSolverType() == "FFTBOX")
704  rho_m *= hr_scaled[0] * hr_scaled[1] * hr_scaled[2];
705 
707  rho_m *= getCouplingConstant();
708 
709  // If debug flag is set, dump scalar field (potential 'phi') into file under ./data/
710 #ifdef DBG_SCALARFIELD
711  INFOMSG("*** START DUMPING SCALAR FIELD ***" << endl);
712 
713  std::ofstream fstr2;
714  fstr2.precision(9);
715 
716  std::string phi_fn = std::string("data/") + SfileName + std::string("-phi_scalar-") + std::string(istr.str());
717  fstr2.open(phi_fn.c_str(), std::ios::out);
718  NDIndex<3> myidx = getFieldLayout().getLocalNDIndex();
719  for (int x = myidx[0].first(); x <= myidx[0].last(); x++) {
720  for (int y = myidx[1].first(); y <= myidx[1].last(); y++) {
721  for (int z = myidx[2].first(); z <= myidx[2].last(); z++) {
722  fstr2 << x + 1 << " " << y + 1 << " " << z + 1 << " " << rho_m[x][y][z].get() << std::endl;
723  }
724  }
725  }
726  fstr2.close();
727 
728  INFOMSG("*** FINISHED DUMPING SCALAR FIELD ***" << endl);
729 #endif
730 
732  eg_m = -Grad(rho_m, eg_m);
733 
738  eg_m *= Vector_t(gamma, 1.0 / gamma, gamma);
739 
740 #ifdef FIELDSTDOUT
741  // Immediate debug output:
742  // Output potential and e-field along the x-, y-, and z-axes
743  int mx = (int)nr_m[0];
744  int mx2 = (int)nr_m[0] / 2;
745  int my = (int)nr_m[1];
746  int my2 = (int)nr_m[1] / 2;
747  int mz = (int)nr_m[2];
748  int mz2 = (int)nr_m[2] / 2;
749 
750  for (int i=0; i<mx; i++ )
751  *gmsg << "Field along x axis Ex = " << eg_m[i][my2][mz2] << " Pot = " << rho_m[i][my2][mz2] << endl;
752 
753  for (int i=0; i<my; i++ )
754  *gmsg << "Field along y axis Ey = " << eg_m[mx2][i][mz2] << " Pot = " << rho_m[mx2][i][mz2] << endl;
755 
756  for (int i=0; i<mz; i++ )
757  *gmsg << "Field along z axis Ez = " << eg_m[mx2][my2][i] << " Pot = " << rho_m[mx2][my2][i] << endl;
758 #endif
759 
760 #ifdef DBG_SCALARFIELD
761  // If debug flag is set, dump vector field (electric field) into file under ./data/
762  INFOMSG("*** START DUMPING E FIELD ***" << endl);
763  //ostringstream oss;
764  //MPI_File file;
765  //MPI_Status status;
766  //MPI_Info fileinfo;
767  //MPI_File_open(Ippl::getComm(), "rho_scalar", MPI_MODE_WRONLY | MPI_MODE_CREATE, fileinfo, &file);
768  std::ofstream fstr;
769  fstr.precision(9);
770 
771  std::string e_field = std::string("data/") + SfileName + std::string("-e_field-") + std::string(istr.str());
772  fstr.open(e_field.c_str(), std::ios::out);
773  NDIndex<3> myidxx = getFieldLayout().getLocalNDIndex();
774  for (int x = myidxx[0].first(); x <= myidxx[0].last(); x++) {
775  for (int y = myidxx[1].first(); y <= myidxx[1].last(); y++) {
776  for (int z = myidxx[2].first(); z <= myidxx[2].last(); z++) {
777  fstr << x + 1 << " " << y + 1 << " " << z + 1 << " " << eg_m[x][y][z].get() << std::endl;
778  }
779  }
780  }
781 
782  fstr.close();
783  fieldDBGStep_m++;
784 
785  //MPI_File_write_shared(file, (char*)oss.str().c_str(), oss.str().length(), MPI_CHAR, &status);
786  //MPI_File_close(&file);
787 
788  INFOMSG("*** FINISHED DUMPING E FIELD ***" << endl);
789 #endif
790 
792  Ef.gather(eg_m, this->R, IntrplCIC_t());
793 
795  // Relativistic E&M says gamma*v/c^2 = gamma*beta/c = sqrt(gamma*gamma-1)/c
796  // but because we already transformed E_trans into the moving frame we have to
797  // add 1/gamma so we are using the E_trans from the rest frame -DW
798  double betaC = sqrt(gamma * gamma - 1.0) / gamma / Physics::c;
799 
801  Bf(0) = betaC * Ef(2);
802  Bf(2) = -betaC * Ef(0);
803 
804  }
805 
806  /*
807  *gmsg << "gamma =" << gamma << endl;
808  *gmsg << "dx,dy,dz =(" << hr_m[0] << ", " << hr_m[1] << ", " << hr_m[2] << ") [m] " << endl;
809  *gmsg << "max of bunch is (" << rmax_m(0) << ", " << rmax_m(1) << ", " << rmax_m(2) << ") [m] " << endl;
810  *gmsg << "min of bunch is (" << rmin_m(0) << ", " << rmin_m(1) << ", " << rmin_m(2) << ") [m] " << endl;
811  */
812 
813  IpplTimings::stopTimer(selfFieldTimer_m);
814 }
815 
833 void PartBunch::computeSelfFields_cycl(int bin) {
834  IpplTimings::startTimer(selfFieldTimer_m);
835 
837  rho_m = 0.0;
838 
840  eg_m = Vector_t(0.0);
841 
843  double gamma = getBinGamma(bin);
844 
845  if(fs_m->hasValidSolver()) {
847  if(fs_m->getFieldSolverType() == "SAAMG")
848  resizeMesh();
849 
851  this->Q.scatter(this->rho_m, this->R, IntrplCIC_t());
852 
855  Vector_t hr_scaled = hr_m ;
856  hr_scaled[1] *= gamma;
857 
859  double tmp2 = 1.0 / (hr_scaled[0] * hr_scaled[1] * hr_scaled[2]);
860  rho_m *= tmp2;
861 
862  // If debug flag is set, dump scalar field (charge density 'rho') into file under ./data/
863 #ifdef DBG_SCALARFIELD
864  INFOMSG("*** START DUMPING SCALAR FIELD ***" << endl);
865  std::ofstream fstr1;
866  fstr1.precision(9);
867 
868  std::ostringstream istr;
869  istr << fieldDBGStep_m;
870 
871  std::string SfileName = OpalData::getInstance()->getInputBasename();
872 
873  std::string rho_fn = std::string("data/") + SfileName + std::string("-rho_scalar-") + std::string(istr.str());
874  fstr1.open(rho_fn.c_str(), std::ios::out);
875  NDIndex<3> myidx1 = getFieldLayout().getLocalNDIndex();
876  for (int x = myidx1[0].first(); x <= myidx1[0].last(); x++) {
877  for (int y = myidx1[1].first(); y <= myidx1[1].last(); y++) {
878  for (int z = myidx1[2].first(); z <= myidx1[2].last(); z++) {
879  fstr1 << x + 1 << " " << y + 1 << " " << z + 1 << " " << rho_m[x][y][z].get() << std::endl;
880  }
881  }
882  }
883  fstr1.close();
884  INFOMSG("*** FINISHED DUMPING SCALAR FIELD ***" << endl);
885 #endif
886 
889  fs_m->solver_m->computePotential(rho_m, hr_scaled);
890 
891  // Do the multiplication of the grid-cube volume coming from the discretization of the convolution integral.
892  // This is only necessary for the FFT solver. FIXME: later move this scaling into FFTPoissonSolver
893  if(fs_m->getFieldSolverType() == "FFT" || fs_m->getFieldSolverType() == "FFTBOX")
894  rho_m *= hr_scaled[0] * hr_scaled[1] * hr_scaled[2];
895 
897  rho_m *= getCouplingConstant();
898 
899  // If debug flag is set, dump scalar field (potential 'phi') into file under ./data/
900 #ifdef DBG_SCALARFIELD
901  INFOMSG("*** START DUMPING SCALAR FIELD ***" << endl);
902 
903  std::ofstream fstr2;
904  fstr2.precision(9);
905 
906  std::string phi_fn = std::string("data/") + SfileName + std::string("-phi_scalar-") + std::string(istr.str());
907  fstr2.open(phi_fn.c_str(), std::ios::out);
908  NDIndex<3> myidx = getFieldLayout().getLocalNDIndex();
909  for (int x = myidx[0].first(); x <= myidx[0].last(); x++) {
910  for (int y = myidx[1].first(); y <= myidx[1].last(); y++) {
911  for (int z = myidx[2].first(); z <= myidx[2].last(); z++) {
912  fstr2 << x + 1 << " " << y + 1 << " " << z + 1 << " " << rho_m[x][y][z].get() << std::endl;
913  }
914  }
915  }
916  fstr2.close();
917 
918  INFOMSG("*** FINISHED DUMPING SCALAR FIELD ***" << endl);
919 #endif
920 
922  eg_m = -Grad(rho_m, eg_m);
923 
928  eg_m *= Vector_t(gamma, 1.0 / gamma, gamma);
929 
930 #ifdef FIELDSTDOUT
931  // Immediate debug output:
932  // Output potential and e-field along the x-, y-, and z-axes
933  int mx = (int)nr_m[0];
934  int mx2 = (int)nr_m[0] / 2;
935  int my = (int)nr_m[1];
936  int my2 = (int)nr_m[1] / 2;
937  int mz = (int)nr_m[2];
938  int mz2 = (int)nr_m[2] / 2;
939 
940  for (int i=0; i<mx; i++ )
941  *gmsg << "Bin " << bin
942  << ", Field along x axis Ex = " << eg_m[i][my2][mz2]
943  << ", Pot = " << rho_m[i][my2][mz2] << endl;
944 
945  for (int i=0; i<my; i++ )
946  *gmsg << "Bin " << bin
947  << ", Field along y axis Ey = " << eg_m[mx2][i][mz2]
948  << ", Pot = " << rho_m[mx2][i][mz2] << endl;
949 
950  for (int i=0; i<mz; i++ )
951  *gmsg << "Bin " << bin
952  << ", Field along z axis Ez = " << eg_m[mx2][my2][i]
953  << ", Pot = " << rho_m[mx2][my2][i] << endl;
954 #endif
955 
956  // If debug flag is set, dump vector field (electric field) into file under ./data/
957 #ifdef DBG_SCALARFIELD
958  INFOMSG("*** START DUMPING E FIELD ***" << endl);
959  //ostringstream oss;
960  //MPI_File file;
961  //MPI_Status status;
962  //MPI_Info fileinfo;
963  //MPI_File_open(Ippl::getComm(), "rho_scalar", MPI_MODE_WRONLY | MPI_MODE_CREATE, fileinfo, &file);
964  std::ofstream fstr;
965  fstr.precision(9);
966 
967  std::string e_field = std::string("data/") + SfileName + std::string("-e_field-") + std::string(istr.str());
968  fstr.open(e_field.c_str(), std::ios::out);
969  NDIndex<3> myidxx = getFieldLayout().getLocalNDIndex();
970  for (int x = myidxx[0].first(); x <= myidxx[0].last(); x++) {
971  for (int y = myidxx[1].first(); y <= myidxx[1].last(); y++) {
972  for (int z = myidxx[2].first(); z <= myidxx[2].last(); z++) {
973  fstr << x + 1 << " " << y + 1 << " " << z + 1 << " " << eg_m[x][y][z].get() << std::endl;
974  }
975  }
976  }
977 
978  fstr.close();
979  fieldDBGStep_m++;
980 
981  //MPI_File_write_shared(file, (char*)oss.str().c_str(), oss.str().length(), MPI_CHAR, &status);
982  //MPI_File_close(&file);
983 
984  INFOMSG("*** FINISHED DUMPING E FIELD ***" << endl);
985 #endif
986 
988  Eftmp.gather(eg_m, this->R, IntrplCIC_t());
989 
991  double betaC = sqrt(gamma * gamma - 1.0) / gamma / Physics::c;
992 
994  Bf(0) = Bf(0) + betaC * Eftmp(2);
995  Bf(2) = Bf(2) - betaC * Eftmp(0);
996 
997  Ef += Eftmp;
998  }
999 
1000  /*
1001  *gmsg << "gamma =" << gamma << endl;
1002  *gmsg << "dx,dy,dz =(" << hr_m[0] << ", " << hr_m[1] << ", " << hr_m[2] << ") [m] " << endl;
1003  *gmsg << "max of bunch is (" << rmax_m(0) << ", " << rmax_m(1) << ", " << rmax_m(2) << ") [m] " << endl;
1004  *gmsg << "min of bunch is (" << rmin_m(0) << ", " << rmin_m(1) << ", " << rmin_m(2) << ") [m] " << endl;
1005  */
1006 
1007 
1008  IpplTimings::stopTimer(selfFieldTimer_m);
1009 }
1010 
1011 
1012 // void PartBunch::setMesh(Mesh_t* mesh) {
1013 // Layout_t* layout = static_cast<Layout_t*>(&getLayout());
1014 // // layout->getLayout().setMesh(mesh);
1015 // }
1016 
1017 
1018 // void PartBunch::setFieldLayout(FieldLayout_t* fLayout) {
1019 // Layout_t* layout = static_cast<Layout_t*>(&getLayout());
1020 // // layout->getLayout().setFieldLayout(fLayout);
1021 // // layout->rebuild_neighbor_data();
1022 // layout->getLayout().changeDomain(*fLayout);
1023 // }
1024 
1025 
1026 FieldLayout_t &PartBunch::getFieldLayout() {
1027  Layout_t* layout = static_cast<Layout_t*>(&getLayout());
1028  return dynamic_cast<FieldLayout_t &>(layout->getLayout().getFieldLayout());
1029 }
1030 
1031 void PartBunch::setBCAllPeriodic() {
1032  for (int i = 0; i < 2 * 3; ++i) {
1033 
1034  if (Ippl::getNodes()>1) {
1035  bc_m[i] = new ParallelInterpolationFace<double, Dimension, Mesh_t, Center_t>(i);
1036  //std periodic boundary conditions for gradient computations etc.
1037  vbc_m[i] = new ParallelPeriodicFace<Vector_t, Dimension, Mesh_t, Center_t>(i);
1038  }
1039  else {
1040  bc_m[i] = new InterpolationFace<double, Dimension, Mesh_t, Center_t>(i);
1041  //std periodic boundary conditions for gradient computations etc.
1042  vbc_m[i] = new PeriodicFace<Vector_t, Dimension, Mesh_t, Center_t>(i);
1043  }
1044  getBConds()[i] = ParticlePeriodicBCond;
1045  }
1046  dcBeam_m=true;
1047  INFOMSG(level3 << "BC set P3M, all periodic" << endl);
1048 }
1049 
1050 void PartBunch::setBCAllOpen() {
1051  for (int i = 0; i < 2 * 3; ++i) {
1052  bc_m[i] = new ZeroFace<double, 3, Mesh_t, Center_t>(i);
1053  vbc_m[i] = new ZeroFace<Vector_t, 3, Mesh_t, Center_t>(i);
1054  getBConds()[i] = ParticleNoBCond;
1055  }
1056  dcBeam_m=false;
1057  INFOMSG(level3 << "BC set for normal Beam" << endl);
1058 }
1059 
1060 void PartBunch::setBCForDCBeam() {
1061  for (int i = 0; i < 2 * 3; ++ i) {
1062  if (i >= 4) {
1063  if (Ippl::getNodes() > 1) {
1064  bc_m[i] = new ParallelPeriodicFace<double, 3, Mesh_t, Center_t>(i);
1065  vbc_m[i] = new ParallelPeriodicFace<Vector_t, 3, Mesh_t, Center_t>(i);
1066  } else {
1067  bc_m[i] = new PeriodicFace<double, 3, Mesh_t, Center_t>(i);
1068  vbc_m[i] = new PeriodicFace<Vector_t, 3, Mesh_t, Center_t>(i);
1069  }
1070 
1071  getBConds()[i] = ParticlePeriodicBCond;
1072  } else {
1073  bc_m[i] = new ZeroFace<double, 3, Mesh_t, Center_t>(i);
1074  vbc_m[i] = new ZeroFace<Vector_t, 3, Mesh_t, Center_t>(i);
1075  getBConds()[i] = ParticleNoBCond;
1076  }
1077  }
1078  dcBeam_m=true;
1079  INFOMSG(level3 << "BC set for DC-Beam, longitudinal periodic" << endl);
1080 }
1081 
1082 
1083 void PartBunch::updateDomainLength(Vektor<int, 3>& grid) {
1084  NDIndex<3> domain = getFieldLayout().getDomain();
1085  for (unsigned int i = 0; i < Dimension; i++)
1086  grid[i] = domain[i].length();
1087 }
1088 
1089 
1090 void PartBunch::updateFields(const Vector_t& hr, const Vector_t& origin) {
1091  getMesh().set_meshSpacing(&(hr_m[0]));
1092  getMesh().set_origin(origin);
1093  rho_m.initialize(getMesh(),
1094  getFieldLayout(),
1095  GuardCellSizes<Dimension>(1),
1096  bc_m);
1097  eg_m.initialize(getMesh(),
1098  getFieldLayout(),
1099  GuardCellSizes<Dimension>(1),
1100  vbc_m);
1101 }
1102 
1103 
1138 inline
1139 PartBunch::VectorPair_t PartBunch::getEExtrema() {
1140  const Vector_t maxE = max(eg_m);
1141  // const double maxL = max(dot(eg_m,eg_m));
1142  const Vector_t minE = min(eg_m);
1143  // INFOMSG("MaxE= " << maxE << " MinE= " << minE << endl);
1144  return VectorPair_t(maxE, minE);
1145 }
1146 
1147 
1148 inline
1149 void PartBunch::resetInterpolationCache(bool clearCache) {
1150  interpolationCacheSet_m = false;
1151  if(clearCache) {
1152  interpolationCache_m.destroy(interpolationCache_m.size(), 0, true);
1153  }
1154 }
1155 
1156 void PartBunch::swap(unsigned int i, unsigned int j) {
1157 
1158  // FIXME
1159  PartBunchBase<double, 3>::swap(i, j);
1160 
1161  if (interpolationCacheSet_m)
1162  std::swap(interpolationCache_m[i], interpolationCache_m[j]);
1163 }
1164 
1165 
1166 Inform &PartBunch::print(Inform &os) {
1167  return PartBunchBase<double, 3>::print(os);
1168 }
PartBunch::setBCForDCBeam
void setBCForDCBeam()
Definition: PartBunch.cpp:1060
PartBunchBase::swap
virtual void swap(unsigned int i, unsigned int j)
Definition: PartBunchBase.hpp:2130
PartBunch::setBCAllPeriodic
void setBCAllPeriodic()
Definition: PartBunch.cpp:1031
PartBunchBase< double, 3 >::rmin_m
Vector_t rmin_m
minimal extend of particles
Definition: PartBunchBase.h:590
IpplInfo::getNodes
static int getNodes()
Definition: IpplInfo.cpp:773
PartBunchBase< double, 3 >::P
ParticleAttrib< Vector_t > P
Definition: PartBunchBase.h:484
PartBunchBase< double, 3 >::selfFieldTimer_m
IpplTimings::TimerRef selfFieldTimer_m
timer for selfField calculation
Definition: PartBunchBase.h:537
RegionLayout::changeDomain
void changeDomain(FieldLayout< Dim > &)
Definition: RegionLayout.hpp:486
PartBunchBase< double, 3 >::destroy
void destroy(size_t M, size_t I, bool doNow=false)
PartBunchBase< double, 3 >::Ef
ParticleAttrib< Vector_t > Ef
Definition: PartBunchBase.h:488
PartBunchBase< double, 3 >::Eftmp
ParticleAttrib< Vector_t > Eftmp
Definition: PartBunchBase.h:489
ParticlePeriodicBCond
T ParticlePeriodicBCond(const T t, const T minval, const T maxval)
Definition: ParticleBConds.h:38
ParticleAttrib::create
virtual void create(size_t)
Definition: ParticleAttrib.hpp:74
ParticleSpatialLayout::getLayout
RegionLayout< T, Dim, Mesh > & getLayout()
Definition: ParticleSpatialLayout.h:120
PartBunch::computeSelfFields
void computeSelfFields()
Definition: PartBunch.cpp:463
ParticleAttrib::scatter
void scatter(Field< T, Dim, M, C > &f, const ParticleAttrib< Vektor< PT, Dim > > &pp, const IntOp &intop) const
Definition: ParticleAttrib.h:203
PoissonSolver::getXRangeMin
virtual double getXRangeMin(unsigned short level=0)=0
PoissonSolver::getXRangeMax
virtual double getXRangeMax(unsigned short level=0)=0
max
T::PETE_Expr_t::PETE_Return_t max(const PETE_Expr< T > &expr, NDIndex< D > &loc)
Definition: ReductionLoc.h:123
PartBunchBase< double, 3 >::Q
ParticleAttrib< double > Q
Definition: PartBunchBase.h:485
PoissonSolver::getZRangeMax
virtual double getZRangeMax(unsigned short level=0)=0
ParticleAttrib::size
size_t size(void) const
Definition: ParticleAttrib.h:179
PartBunchBase< double, 3 >::nr_m
Vektor< int, 3 > nr_m
meshsize of cartesian mesh
Definition: PartBunchBase.h:623
PartBunch::swap
void swap(unsigned int i, unsigned int j)
Definition: PartBunch.cpp:1156
PartBunch::print
Inform & print(Inform &os)
Definition: PartBunch.cpp:1166
PartData
Particle reference data.
Definition: PartData.h:38
PartBunchBase< double, 3 >::getCouplingConstant
double getCouplingConstant() const
gmsg
Inform * gmsg
Definition: Main.cpp:21
PartBunch::getLayout
ParticleLayout< double, 3 > & getLayout()
Definition: PartBunch.h:128
FieldSolver::getFieldSolverType
std::string getFieldSolverType()
Definition: FieldSolver.h:90
ParticleNoBCond
T ParticleNoBCond(const T t, const T, const T)
Definition: ParticleBConds.h:32
PartBunch::interpolationCache_m
ParticleAttrib< CacheDataCIC< double, 3U > > interpolationCache_m
Definition: PartBunch.h:125
sum
T::PETE_Expr_t::PETE_Return_t sum(const PETE_Expr< T > &expr)
Definition: PETE.h:1213
PartBunchBase< double, 3 >::getTotalNum
size_t getTotalNum() const
Grad
Field< Vektor< T, 1U >, 1U, Cartesian< 1U, MFLOAT >, Cell > & Grad(Field< T, 1U, Cartesian< 1U, MFLOAT >, Vert > &x, Field< Vektor< T, 1U >, 1U, Cartesian< 1U, MFLOAT >, Cell > &r)
Definition: Cartesian.hpp:2709
PartBunch::vbc_m
BConds< Vector_t, 3, Mesh_t, Center_t > vbc_m
Definition: PartBunch.h:120
level2
Inform & level2(Inform &inf)
Definition: Inform.cpp:46
OpalData::getInstance
static OpalData * getInstance()
Definition: OpalData.cpp:209
PartBunchBase< double, 3 >::pbase
std::shared_ptr< AbstractParticle< double, Dim > > pbase
Definition: PartBunchBase.h:685
PartBunchBase< double, 3 >::VectorPair_t
std::pair< Vector_t, Vector_t > VectorPair_t
Definition: PartBunchBase.h:37
FieldSolver::hasValidSolver
bool hasValidSolver()
Definition: FieldSolver.cpp:543
PartBunchBase
Class: DataSink.
Definition: OpalData.h:29
PartBunch::rho_m
Field_t rho_m
scalar potential
Definition: PartBunch.h:102
PartBunch::runTests
void runTests()
Definition: PartBunch.cpp:76
BinaryRepartition
void BinaryRepartition(FieldLayout< Dim > &layout, BareField< double, Dim > &weights)
Definition: BinaryBalancer.h:68
FieldSolver::solver_m
PoissonSolver * solver_m
the actual solver, should be a base object
Definition: FieldSolver.h:114
PartBunch::updateDomainLength
void updateDomainLength(Vektor< int, 3 > &grid)
Definition: PartBunch.cpp:1083
Physics::c
constexpr double c
The velocity of light in m/s.
Definition: Physics.h:52
PartBunch::interpolationCacheSet_m
bool interpolationCacheSet_m
Definition: PartBunch.h:123
IpplTimings::startTimer
static void startTimer(TimerRef t)
Definition: IpplTimings.h:187
INFOMSG
#define INFOMSG(msg)
Definition: IpplInfo.h:397
PartBunchBase< double, 3 >::rmax_m
Vector_t rmax_m
maximal extend of particles
Definition: PartBunchBase.h:588
Vector_t
Vektor< double, 3 > Vector_t
Definition: Vektor.h:6
UniformCartesian::set_origin
void set_origin(const Vektor< MFLOAT, Dim > &o)
Definition: UniformCartesian.hpp:431
PoissonSolver::computePotential
virtual void computePotential(Field_t &rho, Vector_t hr)=0
PartBunch::getEExtrema
VectorPair_t getEExtrema()
Definition: PartBunch.cpp:1139
PartBunchBase< double, 3 >::fs_m
FieldSolver * fs_m
stores the used field solver
Definition: PartBunchBase.h:626
PartBunchBase< double, 3 >::getLocalNum
size_t getLocalNum() const
ParticleAttrib::destroy
virtual void destroy(size_t M, size_t I, bool optDestroy=true)
Definition: ParticleAttrib.hpp:98
PartBunchBase::print
Inform & print(Inform &os)
Definition: PartBunchBase.hpp:1905
PoissonSolver::getYRangeMin
virtual double getYRangeMin(unsigned short level=0)=0
sqrt
Tps< T > sqrt(const Tps< T > &x)
Square root.
Definition: TpsMath.h:91
PartBunchBase< double, 3 >::dt
ParticleAttrib< double > dt
Definition: PartBunchBase.h:493
PartBunch::do_binaryRepart
void do_binaryRepart()
Definition: PartBunch.cpp:106
PartBunch::initialize
void initialize(FieldLayout_t *fLayout)
Definition: PartBunch.cpp:71
PoissonSolver::getYRangeMax
virtual double getYRangeMax(unsigned short level=0)=0
PartBunch::eg_m
VField_t eg_m
vector field on the grid
Definition: PartBunch.h:105
PartBunchBase< double, 3 >::hr_m
Vector_t hr_m
meshspacing of cartesian mesh
Definition: PartBunchBase.h:621
PartBunch::getMesh
const Mesh_t & getMesh() const
Definition: PartBunch.h:146
PartBunchBase< double, 3 >::Dimension
static const unsigned Dimension
Definition: PartBunchBase.h:39
PartBunchBase< double, 3 >::getBinGamma
double getBinGamma(int bin)
Get gamma of one bin.
PartBunch::setBCAllOpen
void setBCAllOpen()
Definition: PartBunch.cpp:1050
PartBunch
Particle Bunch.
Definition: PartBunch.h:30
PoissonSolver::getZRangeMin
virtual double getZRangeMin(unsigned short level=0)=0
RegionLayout::getFieldLayout
FieldLayout< Dim > & getFieldLayout()
Definition: RegionLayout.h:137
PartBunch::resizeMesh
void resizeMesh()
resize mesh to geometry specified
Definition: PartBunch.cpp:413
PartBunchBase< double, 3 >::Bf
ParticleAttrib< Vector_t > Bf
Definition: PartBunchBase.h:491
Field::initialize
void initialize(Layout_t &)
PartBunch::updateFields
void updateFields(const Vector_t &hr, const Vector_t &origin)
Definition: PartBunch.cpp:1090
level3
Inform & level3(Inform &inf)
Definition: Inform.cpp:47
UniformCartesian::set_meshSpacing
void set_meshSpacing(MFLOAT *const del)
Definition: UniformCartesian.hpp:447
PartBunch::resetInterpolationCache
void resetInterpolationCache(bool clearCache=false)
Definition: PartBunch.cpp:1149
Field::get_mesh
Mesh_t & get_mesh() const
Definition: Field.h:113
PartBunchBase< double, 3 >::get_bounds
void get_bounds(Vector_t &rmin, Vector_t &rmax)
PartBunchBase< double, 3 >::getBConds
ParticleBConds< Position_t, Dimension > & getBConds()
Definition: PartBunchBase.h:432
IpplTimings::stopTimer
static void stopTimer(TimerRef t)
Definition: IpplTimings.h:192
OpalData::getInputBasename
std::string getInputBasename()
get input file name without extension
Definition: OpalData.cpp:717
Inform
Definition: Inform.h:41
PartBunch::bc_m
BConds< double, 3, Mesh_t, Center_t > bc_m
for defining the boundary conditions
Definition: PartBunch.h:119
level1
Inform & level1(Inform &inf)
Definition: Inform.cpp:45
FieldLayout::getLocalNDIndex
NDIndex< Dim > getLocalNDIndex()
Definition: FieldLayout.hpp:1579
min
T::PETE_Expr_t::PETE_Return_t min(const PETE_Expr< T > &expr, NDIndex< D > &loc)
Definition: ReductionLoc.h:95
FieldLayout::getDomain
const NDIndex< Dim > & getDomain() const
Definition: FieldLayout.h:325
IntrplCIC_t
IntCIC IntrplCIC_t
Definition: PBunchDefs.h:24
PoissonSolver::test
virtual void test(PartBunchBase< double, 3 > *bunch)=0
endl
Inform & endl(Inform &inf)
Definition: Inform.cpp:42
PartBunch::PartBunch
PartBunch()=delete
PartBunch::computeSelfFields_cycl
void computeSelfFields_cycl(double gamma)
Calculates the self electric field from the charge density distribution for use in cyclotrons...
Definition: PartBunch.cpp:643
PartBunch::getFieldLayout
FieldLayout_t & getFieldLayout()
Definition: PartBunch.cpp:1026
Generated on Sat Jan 4 2020 14:47:51 for OPAL (Object Oriented Parallel Accelerator Library) by   doxygen 1.8.5