OPAL (Object Oriented Parallel Accelerator Library) 2022.1
OPAL
FFTBoxPoissonSolver.cpp
Go to the documentation of this file.
1// -*- C++ -*-
2/***************************************************************************
3 *
4 *
5 * FFTBoxPoissonSolver.cc
6 *
7 *
8 ***************************************************************************/
9
10// include files
12
14#include "Physics/Physics.h"
15#include "Utility/IpplTimings.h"
16#include <fstream>
18// a little helper class to specialize the action of the Green's function
19// calculation. This should be specialized for each dimension
20// to the proper action for computing the Green's function. The first
21// template parameter is the full type of the Field to compute, and the second
22// is the dimension of the data, which should be specialized.
23
24template<unsigned int Dim>
26
27template<>
29 template<class T, class FT, class FT2>
30 static void calculate(Vektor<T, 3> &hrsq, FT &grn, FT2 *grnI) {
31 grn = grnI[0] * hrsq[0] + grnI[1] * hrsq[1] + grnI[2] * hrsq[2];
32 grn = 1.0 / sqrt(grn);
33 grn[0][0][0] = grn[0][0][1];
34 }
35};
36
37FFTBoxPoissonSolver::FFTBoxPoissonSolver(Mesh_t *mesh, FieldLayout_t *fl, std::string greensFunction, double boxSize):
38 mesh_m(mesh),
39 layout_m(fl),
40 mesh2_m(0),
41 layout2_m(0),
42 greensFunction_m(greensFunction),
43 a_m(boxSize)
44{
45 int i;
47
48 // For efficiency in the FFT's, we can use a parallel decomposition
49 // which can be serial in the first dimension.
50 // e_dim_tag decomp[3];
51 // // e_dim_tag decomp2[3];
52 // for(int d = 0; d < 3; ++d) {
53 // decomp[d] = layout_m->getRequestedDistribution(d);
54 // //decomp2[d] = layout_m->getRequestedDistribution(d);
55 // }
56
57 // The FFT's require double-sized field sizes in order to (more closely
58 // do not understand this ...)
59 // simulate an isolated system. The FFT of the charge density field, rho,
60 // would otherwise mimic periodic boundary conditions, i.e. as if there were
61 // several beams set a periodic distance apart. The double-sized fields
62 // alleviate this problem.
63 for(i = 0; i < 3; i++) {
65 nr_m[i] = domain_m[i].length();
66 //domain2_m[i] = Index( 2*nr_m[i] + 1);
67 }
68
69 // create double sized mesh and layout objects for the use in the FFT's
70 //mesh2_m = new Mesh_t(domain2_m);
71 //layout2_m = new FieldLayout_t(*mesh2_m, decomp);
72 //rho2_m.initialize(*mesh2_m, *layout2_m);
75
76 bool sineTransformDims[3];
77 for(int d = 0; d < 3; ++d)
78 sineTransformDims[d] = true;
79
80 sine_m = new SINE_t(grntr_m.getDomain(), sineTransformDims, false);
81
82 // these are fields that are used for calculating the Green's function.
83 // they eliminate some calculation at each time-step.
84 for(i = 0; i < 3; ++i) {
87 domain_m[i] * domain_m[i],
88 (nr_m[i] - domain_m[i]) *
89 (nr_m[i] - domain_m[i]));
90 }
91
93
94 if(greensFunction_m == std::string("INTEGRATED")) {
99
104 } else {
107 }
108}
109
110FFTBoxPoissonSolver::FFTBoxPoissonSolver(PartBunch &beam, std::string greensFunction):
111 mesh_m(&beam.getMesh()),
112 layout_m(&beam.getFieldLayout()),
113 mesh2_m(0),
114 layout2_m(0),
115 greensFunction_m(greensFunction)
116{
117 int i;
119
120 // For efficiency in the FFT's, we can use a parallel decomposition
121 // which can be serial in the first dimension.
122 // e_dim_tag decomp[3];
123 // for(int d = 0; d < 3; ++d) {
124 // decomp[d] = layout_m->getRequestedDistribution(d);
125 // }
126
127 // The FFT's require double-sized field sizes in order to (more closely
128 // do not understand this ...)
129 // simulate an isolated system. The FFT of the charge density field, rho,
130 // would otherwise mimic periodic boundary conditions, i.e. as if there were
131 // several beams set a periodic distance apart. The double-sized fields
132 // alleviate this problem.
133 for(i = 0; i < 3; i++) {
134 hr_m[i] = mesh_m->get_meshSpacing(i);
135 nr_m[i] = domain_m[i].length();
136 }
137
138 //rho_m.initialize(*mesh_m, *layout_m);
141
142 bool sineTransformDims[3];
143 for(int d = 0; d < 3; ++d)
144 sineTransformDims[d] = true;
145
146 // create the FFT object
147 sine_m = new SINE_t(grntr_m.getDomain(), sineTransformDims, false);
148
149 // these are fields that are used for calculating the Green's function.
150 // they eliminate some calculation at each time-step.
151 for(i = 0; i < 3; ++i) {
154 domain_m[i] * domain_m[i],
155 (nr_m[i] - domain_m[i]) *
156 (nr_m[i] - domain_m[i]));
157 }
158 GreensFunctionTimer_m = IpplTimings::getTimer("SF: GreenFTotal");
159
160 if(greensFunction_m == std::string("INTEGRATED")) {
165
170 } else {
173 }
174}
175
177// destructor
179 // delete the FFT object
180 delete sine_m;
181}
182
184// given a charge-density field rho and a set of mesh spacings hr,
185// compute the electric potential from the image charge by solving
186// the Poisson's equation
187
189 // use grid of complex doubled in both dimensions
190 // and store rho in lower left quadrant of doubled grid
191 //rho2_m = 0.0;
192
193 //rho2_m[domain_m] = rho[domain_m];
194
195 // needed in greens function
196 hr_m = hr;
197 // FFT double-sized charge density
198 // we do a backward transformation so that we dont have to account for the normalization factor
199 // that is used in the forward transformation of the IPPL FFT
200 sine_m->transform(-1, rho);
201
202 // must be called if the mesh size has changed
203 // have to check if we can do G with h = (1,1,1)
204 // and rescale later
205
206 // Do image charge.
207 // The minus sign is due to image charge.
208 // Convolute transformed charge density with shifted green's function.
212
213 // Multiply transformed charge density and
214 // transformed Green's function.
215 rho = - rho * grntr_m;
216
217 // Inverse FFT to find image charge potential, rho2_m equals the electrostatic potential.
218 sine_m->transform(+1, rho);
219
220 // Re-use rho to store image potential. Flip z coordinate since this is a mirror image.
221 Index I = nr_m[0];
222 Index J = nr_m[1];
223 Index K = nr_m[2];
224 rho[I][J][K] = rho[I][J][nr_m[2] - K - 1];
225}
226
227
229// given a charge-density field rho and a set of mesh spacings hr,
230// compute the electric field and put in eg by solving the Poisson's equation
231
233 // use grid of complex doubled in both dimensions
234 // and store rho in lower left quadrant of doubled grid
235 //rho2_m = 0.0;
236
237 //rho2_m[domain_m] = rho[domain_m];
238
239 // needed in greens function
240 hr_m = hr;
241 // FFT double-sized charge density
242 // we do a backward transformation so that we dont have to account for the normalization factor
243 // that is used in the forward transformation of the IPPL FFT
244 sine_m->transform(-1, rho);
245 // must be called if the mesh size has changed
246 // have to check if we can do G with h = (1,1,1)
247 // and rescale later
249 if(greensFunction_m == std::string("INTEGRATED"))
251 else
254 // multiply transformed charge density
255 // and transformed Green function
256 // Don't divide by (2*nx_m)*(2*ny_m), as Ryne does;
257 // this normalization is done in POOMA's fft routine.
258 rho *= grntr_m;
259
260 // inverse FFT, rho2_m equals to the electrostatic potential
261 sine_m->transform(+1, rho);
262 // end convolution!
263
264 // back to physical grid
265 // reuse the charge density field to store the electrostatic potential
266 //rho[domain_m] = rho2_m[domain_m];
267}
269// calculate the FFT of the Green's function for the given field
271
272 Vector_t hrsq(hr_m * hr_m);
276 // Green's function calculation complete at this point.
277 // The next step is to FFT it.
278
279 // we do a backward transformation so that we dont have to account for the normalization factor
280 // that is used in the forward transformation of the IPPL FFT
281
283 sine_m->transform(-1, grntr_m);
285}
286
306
307 tmpgreen = 0.0;
308 double tmpgrn, r;
310
312
316 for(int k = idx[2].first(); k < std::min(nr_m[2] + 2, idx[2].last() + 3); k++) {
317 for(int j = idx[1].first(); j < std::min(nr_m[1] + 2, idx[1].last() + 3); j++) {
318 for(int i = idx[0].first(); i < std::min(nr_m[0] + 2, idx[0].last() + 3); i++) {
319
320 Vector_t vv = Vector_t(0.0);
321 vv(0) = i * hr_m[0] - hr_m[0] / 2;
322 vv(1) = j * hr_m[1] - hr_m[1] / 2;
323 vv(2) = k * hr_m[2] - hr_m[2] / 2;
324
325 r = sqrt(vv(0) * vv(0) + vv(1) * vv(1) + vv(2) * vv(2));
326 tmpgrn = -vv(2) * vv(2) * atan(vv(0) * vv(1) / (vv(2) * r)) / 2;
327 tmpgrn += -vv(1) * vv(1) * atan(vv(0) * vv(2) / (vv(1) * r)) / 2;
328 tmpgrn += -vv(0) * vv(0) * atan(vv(1) * vv(2) / (vv(0) * r)) / 2;
329 tmpgrn += vv(1) * vv(2) * log(vv(0) + r);
330 tmpgrn += vv(0) * vv(2) * log(vv(1) + r);
331 tmpgrn += vv(0) * vv(1) * log(vv(2) + r);
332
333 tmpgreen[i][j][k] = tmpgrn / (hr_m[0] * hr_m[1] * hr_m[2]);
334
335 }
336 }
337 }
339
341
342 grntr_m = 0.0;
343
344 Index I = nr_m[0] + 1;
345 Index J = nr_m[1] + 1;
346 Index K = nr_m[2] + 1;
347
348 // the actual integration
349 grntr_m[I][J][K] = tmpgreen[I+1][J+1][K+1];
350 grntr_m[I][J][K] += tmpgreen[I+1][J][K];
351 grntr_m[I][J][K] += tmpgreen[I][J+1][K];
352 grntr_m[I][J][K] += tmpgreen[I][J][K+1];
353 grntr_m[I][J][K] -= tmpgreen[I+1][J+1][K];
354 grntr_m[I][J][K] -= tmpgreen[I+1][J][K+1];
355 grntr_m[I][J][K] -= tmpgreen[I][J+1][K+1];
356 grntr_m[I][J][K] -= tmpgreen[I][J][K];
357
360
361 //assign seems to have problems when we need values that are on another CPU, i.e. [I+1]
362 /*assign(rho2_m[I][J][K] ,
363 tmpgreen[I+1][J+1][K+1] - tmpgreen[I][J+1][K+1] - tmpgreen[I+1][J][K+1] + tmpgreen[I][J][K+1] - tmpgreen[I+1][J+1][K] +
364 tmpgreen[I][J+1][K] + tmpgreen[I+1][J][K] - tmpgreen[I][J][K]);*/
365
366 grntr_m[0][0][0] = grntr_m[0][0][1];
367
368 Index IE(nr_m[0] + 1, 2 * nr_m[0] - 1);
369 Index JE(nr_m[1] + 1, 2 * nr_m[1] - 1);
370 Index KE(nr_m[2] + 1, 2 * nr_m[2] - 1);
371
372 grntr_m[IE][JE][KE] = grntr_m[2*nr_m[0] - IE][2*nr_m[1] - JE][2*nr_m[2] - KE];
373 grntr_m[IE][J ][K ] = grntr_m[2*nr_m[0] - IE][J][K];
374 grntr_m[I ][JE][K ] = grntr_m[I][2*nr_m[1] - JE][K];
375 grntr_m[I ][J ][KE] = grntr_m[I][J][2*nr_m[2] - KE];
376 grntr_m[IE][JE][K ] = grntr_m[2*nr_m[0] - IE][2*nr_m[1] - JE][K];
377 grntr_m[IE][J ][KE] = grntr_m[2*nr_m[0] - IE][J][2*nr_m[2] - KE];
378 grntr_m[I ][JE][KE] = grntr_m[I][2*nr_m[1] - JE][2*nr_m[2] - KE];
380
382 sine_m->transform(-1, grntr_m);
384}
385
386
388
389 tmpgreen = 0.0;
390 double tmpgrn, r;
392 Field_t grn2, ggrn2;
394 grn2 = 0.0;
395 ggrn2.initialize(*mesh2_m, *layout2_m);
396 ggrn2 = 0.0;
397
399
400 for(int k = idx[2].first(); k < std::min(nr_m[2] + 2, idx[2].last() + 3); k++) {
401 for(int j = idx[1].first(); j < std::min(nr_m[1] + 2, idx[1].last() + 3); j++) {
402 for(int i = idx[0].first(); i < std::min(nr_m[0] + 2, idx[0].last() + 3); i++) {
403
404 Vector_t vv = Vector_t(0.0);
405 vv(0) = i * hr_m[0] - hr_m[0] / 2;
406 vv(1) = j * hr_m[1] - hr_m[1] / 2;
407 vv(2) = k * hr_m[2] - hr_m[2] / 2 + zshift;
408
409 r = sqrt(vv(0) * vv(0) + vv(1) * vv(1) + vv(2) * vv(2));
410 tmpgrn = -vv(2) * vv(2) * atan(vv(0) * vv(1) / (vv(2) * r)) / 2;
411 tmpgrn += -vv(1) * vv(1) * atan(vv(0) * vv(2) / (vv(1) * r)) / 2;
412 tmpgrn += -vv(0) * vv(0) * atan(vv(1) * vv(2) / (vv(0) * r)) / 2;
413 tmpgrn += vv(1) * vv(2) * log(vv(0) + r);
414 tmpgrn += vv(0) * vv(2) * log(vv(1) + r);
415 tmpgrn += vv(0) * vv(1) * log(vv(2) + r);
416
417 tmpgreen[i][j][k] = tmpgrn / (hr_m[0] * hr_m[1] * hr_m[2]);
418
419 }
420 }
421 }
423
425
426 for(int k = idx[2].first(); k < std::min(nr_m[2] + 2, idx[2].last() + 3); k++) {
427 for(int j = idx[1].first(); j < std::min(nr_m[1] + 2, idx[1].last() + 3); j++) {
428 for(int i = idx[0].first(); i < std::min(nr_m[0] + 2, idx[0].last() + 3); i++) {
429
430 Vector_t vv = Vector_t(0.0);
431 vv(0) = i * hr_m[0] - hr_m[0] / 2;
432 vv(1) = j * hr_m[1] - hr_m[1] / 2;
433 vv(2) = zshift - hr_m[2] * (nr_m[2] - k) - hr_m[2] / 2;
434 //vv(2) = zshift-k*hr_m[2]-hr_m[2]/2;
435
436 r = sqrt(vv(0) * vv(0) + vv(1) * vv(1) + vv(2) * vv(2));
437 tmpgrn = -vv(2) * vv(2) * atan(vv(0) * vv(1) / (vv(2) * r)) / 2;
438 tmpgrn += -vv(1) * vv(1) * atan(vv(0) * vv(2) / (vv(1) * r)) / 2;
439 tmpgrn += -vv(0) * vv(0) * atan(vv(1) * vv(2) / (vv(0) * r)) / 2;
440 tmpgrn += vv(1) * vv(2) * log(vv(0) + r);
441 tmpgrn += vv(0) * vv(2) * log(vv(1) + r);
442 tmpgrn += vv(0) * vv(1) * log(vv(2) + r);
443
444 grn2[i][j][k] = tmpgrn / (hr_m[0] * hr_m[1] * hr_m[2]);
445
446 }
447 }
448 }
450
452
453 grntr_m = 0.0;
454
455 Index I = nr_m[0] + 1;
456 Index J = nr_m[1] + 1;
457 Index K = nr_m[2] + 1;
458
459 // the actual integration
460 grntr_m[I][J][K] = tmpgreen[I+1][J+1][K+1];
461 grntr_m[I][J][K] += tmpgreen[I+1][J][K];
462 grntr_m[I][J][K] += tmpgreen[I][J+1][K];
463 grntr_m[I][J][K] += tmpgreen[I][J][K+1];
464 grntr_m[I][J][K] -= tmpgreen[I+1][J+1][K];
465 grntr_m[I][J][K] -= tmpgreen[I+1][J][K+1];
466 grntr_m[I][J][K] -= tmpgreen[I][J+1][K+1];
467 grntr_m[I][J][K] -= tmpgreen[I][J][K];
468 tmpgreen = 0.0;
469 //later replace ggrn2 with tmpgreen
470 ggrn2[I][J][K] = grn2[I+1][J+1][K+1];
471 ggrn2[I][J][K] += grn2[I+1][J][K];
472 ggrn2[I][J][K] += grn2[I][J+1][K];
473 ggrn2[I][J][K] += grn2[I][J][K+1];
474 ggrn2[I][J][K] -= grn2[I+1][J+1][K];
475 ggrn2[I][J][K] -= grn2[I+1][J][K+1];
476 ggrn2[I][J][K] -= grn2[I][J+1][K+1];
477 ggrn2[I][J][K] -= grn2[I][J][K];
478
479 Index IE(nr_m[0] + 1, 2 * nr_m[0] - 1);
480 Index JE(nr_m[1] + 1, 2 * nr_m[1] - 1);
481 Index KE(nr_m[2], 2 * nr_m[2] - 1);
482
495 grntr_m[IE][J ][K ] = grntr_m[2*nr_m[0] - IE][J][K];
496 grntr_m[I ][JE][K ] = grntr_m[I][2*nr_m[1] - JE][K];
497 grntr_m[IE][JE][K ] = grntr_m[2*nr_m[0] - IE][2*nr_m[1] - JE][K];
498
499 grntr_m[I ][J ][KE] = ggrn2[I][J][KE - nr_m[2]];
500 grntr_m[IE][J ][KE] = ggrn2[2*nr_m[0] - IE][J][KE - nr_m[2]];
501 grntr_m[I ][JE][KE] = ggrn2[I][2*nr_m[1] - JE][KE - nr_m[2]];
502 grntr_m[IE][JE][KE] = ggrn2[2*nr_m[0] - IE][2*nr_m[1] - JE][KE - nr_m[2]];
503
505
507 sine_m->transform(-1, grntr_m);
509
510}
511
513 os << "* ************* F F T P o i s s o n S o l v e r ************************************ " << endl;
514 os << "* h " << hr_m << '\n';
515 os << "* ********************************************************************************** " << endl;
516 return os;
517}
518
519/***************************************************************************
520 * $RCSfile: FFTBoxPoissonSolver.cc,v $ $Author: adelmann $
521 * $Revision: 1.6 $ $Date: 2001/08/16 09:36:08 $
522 ***************************************************************************/
Tps< T > log(const Tps< T > &x)
Natural logarithm.
Definition: TpsMath.h:182
Tps< T > sqrt(const Tps< T > &x)
Square root.
Definition: TpsMath.h:91
T::PETE_Expr_t::PETE_Return_t min(const PETE_Expr< T > &expr, NDIndex< D > &loc)
Definition: ReductionLoc.h:76
PETE_TUTree< FnArcTan, typename T::PETE_Expr_t > atan(const PETE_Expr< T > &l)
Definition: PETE.h:727
PETE_TTTree< OpWhere, typename Cond_t::PETE_Expr_t, typename True_t::PETE_Expr_t, PETE_Scalar< Vektor< T, Dim > > > where(const PETE_Expr< Cond_t > &c, const PETE_Expr< True_t > &t, const Vektor< T, Dim > &f)
Inform & endl(Inform &inf)
Definition: Inform.cpp:42
bool lt(double x, double y)
static void calculate(Vektor< T, 3 > &hrsq, FT &grn, FT2 *grnI)
IpplTimings::TimerRef ShIntGreensFunctionTimer4_m
IpplTimings::TimerRef GreensFunctionTimer_m
FieldLayout_t * layout_m
void computePotential(Field_t &rho, Vector_t hr, double zshift)
Inform & print(Inform &os) const
IpplTimings::TimerRef ShIntGreensFunctionTimer2_m
FFT< SineTransform, 3, double > SINE_t
IpplTimings::TimerRef IntGreensFunctionTimer4_m
FieldLayout_t * layout2_m
IpplTimings::TimerRef IntGreensFunctionTimer3_m
void shiftedIntGreensFunction(double zshift)
compute the shifted integrated Green function as described in Three-dimensional quasistatic model for...
void integratedGreensFunction()
compute the integrated Green function as described in Three-dimensional quasistatic model for high br...
IpplTimings::TimerRef GreensFunctionTimer4_m
IpplTimings::TimerRef IntGreensFunctionTimer2_m
FFTBoxPoissonSolver(PartBunch &bunch, std::string greensFuntion)
IpplTimings::TimerRef ShIntGreensFunctionTimer3_m
IpplTimings::TimerRef GreensFunctionTimer1_m
IpplTimings::TimerRef IntGreensFunctionTimer1_m
IpplTimings::TimerRef ShIntGreensFunctionTimer1_m
Definition: Vektor.h:32
void initialize(Layout_t &)
const NDIndex< Dim > & getDomain() const
Definition: BareField.h:152
NDIndex< Dim > getLocalNDIndex()
const NDIndex< Dim > & getDomain() const
Definition: FieldLayout.h:325
MFLOAT get_meshSpacing(unsigned d) const
Definition: Index.h:237
Definition: Inform.h:42
static TimerRef getTimer(const char *nm)
Definition: IpplTimings.h:182
static void stopTimer(TimerRef t)
Definition: IpplTimings.h:192
static void startTimer(TimerRef t)
Definition: IpplTimings.h:187
Vektor< double, 3 > Vector_t
Definition: Vektor.h:6