Project
Loading...
Searching...
No Matches
PoissonSolver.cxx
Go to the documentation of this file.
1// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3// All rights not expressly granted are reserved.
4//
5// This software is distributed under the terms of the GNU General Public
6// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7//
8// In applying this license CERN does not waive the privileges and immunities
9// granted to it by virtue of its status as an Intergovernmental Organization
10// or submit itself to any jurisdiction.
11
20
23#include "Framework/Logger.h"
24#include <numeric>
25#include <fmt/core.h>
28#include "DataFormatsTPC/Defs.h"
29
30#ifdef WITH_OPENMP
31#include <omp.h>
32#endif
33
34using namespace o2::tpc;
35
36template <typename DataT>
37void PoissonSolver<DataT>::poissonSolver3D(DataContainer& matricesV, const DataContainer& matricesCharge, const int symmetry)
38{
39 using timer = std::chrono::high_resolution_clock;
40 auto start = timer::now();
42 poissonMultiGrid3D(matricesV, matricesCharge, symmetry);
43 } else {
44 poissonMultiGrid3D2D(matricesV, matricesCharge, symmetry);
45 }
46 auto stop = timer::now();
47 std::chrono::duration<float> time = stop - start;
48 const float totalTime = time.count();
49 LOGP(detail, "poissonSolver3D took {}s", totalTime);
50}
51
52template <typename DataT>
54{
55 poissonMultiGrid2D(matricesV, matricesCharge);
56}
57
58template <typename DataT>
59void PoissonSolver<DataT>::poissonMultiGrid2D(DataContainer& matricesV, const DataContainer& matricesCharge, const int iPhi)
60{
62 const DataT gridSpacingR = getSpacingR();
63 const DataT gridSpacingZ = getSpacingZ();
64 const DataT ratioZ = gridSpacingR * gridSpacingR / (gridSpacingZ * gridSpacingZ); // ratio_{Z} = gridSize_{r} / gridSize_{z}
65
66 int nGridRow = 0; // number grid
67 int nGridCol = 0; // number grid
68
69 int nnRow = mParamGrid.NRVertices;
70 while (nnRow >>= 1) {
71 ++nGridRow;
72 }
73
74 int nnCol = mParamGrid.NZVertices;
75 while (nnCol >>= 1) {
76 ++nGridCol;
77 }
78
79 // Check that number of mParamGrid.NRVertices and mParamGrid.NZVertices is suitable for multi grid
80 if (!isPowerOfTwo(mParamGrid.NRVertices - 1)) {
81 LOGP(error, "PoissonMultiGrid2D: PoissonMultiGrid - Error in the number of mParamGrid.NRVertices. Must be 2**M + 1");
82 return;
83 }
84 if (!isPowerOfTwo(mParamGrid.NZVertices - 1)) {
85 LOGP(error, "PoissonMultiGrid2D: PoissonMultiGrid - Error in the number of mParamGrid.NZVertices. Must be 2**N + 1");
86 return;
87 }
88
89 const int nLoop = std::max(nGridRow, nGridCol); // Calculate the number of nLoop for the binary expansion
90
91 LOGP(detail, "{}", fmt::format("PoissonMultiGrid2D: nGridRow={}, nGridCol={}, nLoop={}, nMGCycle={}", nGridRow, nGridCol, nLoop, MGParameters::nMGCycle));
92
93 unsigned int iOne = 1; // index
94 unsigned int jOne = 1; // index
95 int tnRRow = mParamGrid.NRVertices;
96 int tnZColumn = mParamGrid.NZVertices;
97
98 // Vector for storing multi grid array
99 std::vector<Vector> tvArrayV(nLoop); // potential <--> error
100 std::vector<Vector> tvChargeFMG(nLoop); // charge is restricted in full multiGrid
101 std::vector<Vector> tvCharge(nLoop); // charge <--> residue
102 std::vector<Vector> tvResidue(nLoop); // residue calculation
103
104 // Allocate memory for temporary grid
105 for (int count = 1; count <= nLoop; ++count) {
106 tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1;
107 tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1;
108 // if one just address to matrixV
109 const int index = count - 1;
110 tvResidue[index].resize(tnRRow, tnZColumn, 1);
111 tvChargeFMG[index].resize(tnRRow, tnZColumn, 1);
112 tvArrayV[index].resize(tnRRow, tnZColumn, 1);
113 tvCharge[index].resize(tnRRow, tnZColumn, 1);
114
115 if (count == 1) {
116 for (int iphi = iPhi; iphi <= iPhi; ++iphi) {
117 for (int ir = 0; ir < mParamGrid.NRVertices; ++ir) {
118 for (int iz = 0; iz < mParamGrid.NZVertices; ++iz) {
119 tvChargeFMG[index](ir, iz, iphi) = matricesCharge(iz, ir, iphi);
120 tvCharge[index](ir, iz, iphi) = matricesCharge(iz, ir, iphi);
121 tvArrayV[index](ir, iz, iphi) = matricesV(iz, ir, iphi);
122 }
123 }
124 }
125 } else {
126 restrict2D(tvChargeFMG[index], tvChargeFMG[count - 2], tnRRow, tnZColumn, 0);
127 }
128 iOne *= 2;
129 jOne *= 2;
130 }
131
133 if (MGParameters::cycleType == CycleType::FCycle) {
134
135 LOGP(detail, "PoissonMultiGrid2D: Do full cycle");
136 // FMG
137 // 1) Relax on the coarsest grid
138 iOne /= 2;
139 jOne /= 2;
140 tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1;
141 tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1;
142 const DataT h = gridSpacingR * nLoop;
143 const DataT h2 = h * h;
144 const DataT tempRatio = ratioZ * iOne * iOne / (jOne * jOne);
145 const DataT tempFourth = 1 / (2 + 2 * tempRatio);
146
147 std::vector<DataT> coefficient1(tnRRow);
148 std::vector<DataT> coefficient2(tnRRow);
149 calcCoefficients2D(1, tnRRow - 1, h, coefficient1, coefficient2);
150
151 relax2D(tvArrayV[nLoop - 1], tvChargeFMG[nLoop - 1], tnRRow, tnZColumn, h2, tempFourth, tempRatio, coefficient1, coefficient2);
152
153 // Do VCycle from nLoop H to h
154 for (int count = nLoop - 2; count >= 0; --count) {
155 iOne /= 2;
156 jOne /= 2;
157
158 tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1;
159 tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1;
160
161 interp2D(tvArrayV[count], tvArrayV[count + 1], tnRRow, tnZColumn, iPhi);
162
163 // Copy the relax charge to the tvCharge
164 tvCharge[count] = tvChargeFMG[count]; // copy
165
166 // Do V cycle
167 for (int mgCycle = 0; mgCycle < MGParameters::nMGCycle; ++mgCycle) {
168 vCycle2D(count + 1, nLoop, MGParameters::nPre, MGParameters::nPost, gridSpacingR, ratioZ, tvArrayV, tvCharge, tvResidue);
169 }
170 }
171 } else if (MGParameters::cycleType == CycleType::VCycle) {
172 // 2. VCycle
173 LOGP(detail, "PoissonMultiGrid2D: Do V cycle");
174
175 int gridFrom = 1;
176 int gridTo = nLoop;
177
178 // Do MGCycle
179 for (int mgCycle = 0; mgCycle < MGParameters::nMGCycle; ++mgCycle) {
180 vCycle2D(gridFrom, gridTo, MGParameters::nPre, MGParameters::nPost, gridSpacingR, ratioZ, tvArrayV, tvCharge, tvResidue);
181 }
182 } else if (MGParameters::cycleType == CycleType::WCycle) {
183 // 3. W Cycle (TODO:)
184 int gridFrom = 1;
185 int gridTo = nLoop;
186 // Do MGCycle
187 for (Int_t mgCycle = 0; mgCycle < MGParameters::nMGCycle; mgCycle++) {
188 wCycle2D(gridFrom, gridTo, MGParameters::gamma, MGParameters::nPre, MGParameters::nPost, gridSpacingR, ratioZ, tvArrayV, tvCharge, tvResidue);
189 }
190 }
191
192 // fill output
193 for (int iphi = iPhi; iphi <= iPhi; ++iphi) {
194 for (int ir = 0; ir < mParamGrid.NRVertices; ++ir) {
195 for (int iz = 0; iz < mParamGrid.NZVertices; ++iz) {
196 matricesV(iz, ir, iphi) = tvArrayV[0](ir, iz, iphi);
197 }
198 }
199 }
200}
201
202template <typename DataT>
203void PoissonSolver<DataT>::poissonMultiGrid3D2D(DataContainer& matricesV, const DataContainer& matricesCharge, const int symmetry)
204{
205 LOGP(detail, "{}", fmt::format("PoissonMultiGrid3D2D: in Poisson Solver 3D multiGrid semi coarsening mParamGrid.NRVertices={}, cols={}, mParamGrid.NPhiVertices={}", mParamGrid.NZVertices, mParamGrid.NRVertices, mParamGrid.NPhiVertices));
206
207 // Check that the number of mParamGrid.NRVertices and mParamGrid.NZVertices is suitable for a binary expansion
208 if (!isPowerOfTwo((mParamGrid.NRVertices - 1))) {
209 LOGP(error, "PoissonMultiGrid3D2D: Poisson3DMultiGrid - Error in the number of mParamGrid.NRVertices. Must be 2**M + 1");
210 return;
211 }
212 if (!isPowerOfTwo((mParamGrid.NZVertices - 1))) {
213 LOGP(error, "PoissonMultiGrid3D2D: Poisson3DMultiGrid - Error in the number of mParamGrid.NZVertices. Must be 2**N + 1");
214 return;
215 }
216 if (mParamGrid.NPhiVertices <= 3) {
217 LOGP(error, "PoissonMultiGrid3D2D: Poisson3DMultiGrid - Error in the number of mParamGrid.NPhiVertices. Must be larger than 3");
218 return;
219 }
220
221 const DataT gridSpacingR = getSpacingR();
222 const DataT gridSpacingZ = getSpacingZ();
223 const DataT gridSpacingPhi = getSpacingPhi();
224 const DataT ratioPhi = gridSpacingR * gridSpacingR / (gridSpacingPhi * gridSpacingPhi); // ratio_{phi} = gridSize_{r} / gridSize_{phi}
225 const DataT ratioZ = gridSpacingR * gridSpacingR / (gridSpacingZ * gridSpacingZ); // ratio_{Z} = gridSize_{r} / gridSize_{z}
226
227 // Solve Poisson's equation in cylindrical coordinates by multiGrid technique
228 // Allow for different size grid spacing in R and Z directions
229 int nGridRow = 0; // number grid
230 int nGridCol = 0; // number grid
231 int nnRow = mParamGrid.NRVertices;
232 int nnCol = mParamGrid.NZVertices;
233
234 while (nnRow >>= 1) {
235 ++nGridRow;
236 }
237 while (nnCol >>= 1) {
238 ++nGridCol;
239 }
240
241 const int maxVal = std::max(nGridRow, nGridCol); // Calculate the number of nLoop for the binary expansion
242 const size_t nLoop = (maxVal > MGParameters::maxLoop) ? MGParameters::maxLoop : maxVal;
243 unsigned int iOne = 1; // index i in gridSize r (original)
244 unsigned int jOne = 1; // index j in gridSize z (original)
245
246 std::vector<Vector> tvArrayV(nLoop); // potential <--> error
247 std::vector<Vector> tvChargeFMG(nLoop); // charge is restricted in full multiGrid
248 std::vector<Vector> tvCharge(nLoop); // charge <--> residue
249 std::vector<Vector> tvPrevArrayV(nLoop); // error calculation
250 std::vector<Vector> tvResidue(nLoop); // residue calculation
251
252 for (unsigned int count = 1; count <= nLoop; count++) {
253 const unsigned int tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1;
254 const unsigned int tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1;
255 const unsigned int index = count - 1;
256 tvResidue[index].resize(tnRRow, tnZColumn, mParamGrid.NPhiVertices);
257 tvPrevArrayV[index].resize(tnRRow, tnZColumn, mParamGrid.NPhiVertices);
258
259 // memory for the finest grid is from parameters
260 tvChargeFMG[index].resize(tnRRow, tnZColumn, mParamGrid.NPhiVertices);
261 tvArrayV[index].resize(tnRRow, tnZColumn, mParamGrid.NPhiVertices);
262 tvCharge[index].resize(tnRRow, tnZColumn, mParamGrid.NPhiVertices);
263
264 if (count == 1) {
265 for (int iphi = 0; iphi < mParamGrid.NPhiVertices; ++iphi) {
266 for (int ir = 0; ir < mParamGrid.NRVertices; ++ir) {
267 for (int iz = 0; iz < mParamGrid.NZVertices; ++iz) {
268 tvChargeFMG[index](ir, iz, iphi) = matricesCharge(iz, ir, iphi);
269 tvArrayV[index](ir, iz, iphi) = matricesV(iz, ir, iphi);
270 }
271 }
272 }
273 } else {
274 restrict3D(tvChargeFMG[index], tvChargeFMG[count - 2], tnRRow, tnZColumn, mParamGrid.NPhiVertices, mParamGrid.NPhiVertices);
275 restrictBoundary3D(tvArrayV[index], tvArrayV[count - 2], tnRRow, tnZColumn, mParamGrid.NPhiVertices, mParamGrid.NPhiVertices);
276 }
277 iOne *= 2; // doubling
278 jOne *= 2; // doubling
279 }
280
281 std::vector<DataT> coefficient1(mParamGrid.NRVertices); // coefficient1(mParamGrid.NRVertices) for storing (1 + h_{r}/2r_{i}) from central differences in r direction
282 std::vector<DataT> coefficient2(mParamGrid.NRVertices); // coefficient2(mParamGrid.NRVertices) for storing (1 + h_{r}/2r_{i}) from central differences in r direction
283 std::vector<DataT> coefficient3(mParamGrid.NRVertices); // coefficient3(mParamGrid.NRVertices) for storing (1/r_{i}^2) from central differences in phi direction
284 std::vector<DataT> coefficient4(mParamGrid.NRVertices); // coefficient4(mParamGrid.NRVertices) for storing 1/2
285 std::vector<DataT> inverseCoefficient4(mParamGrid.NRVertices); // inverse of coefficient4(mParamGrid.NRVertices)
286
287 // Case full multi grid (FMG)
288 if (MGParameters::cycleType == CycleType::FCycle) {
289 // 1) Relax on the coarsest grid
290 iOne /= 2;
291 jOne /= 2;
292 int tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1;
293 int tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1;
294
295 const DataT h = getSpacingR() * iOne;
296 const DataT h2 = h * h;
297 const DataT iOne2 = iOne * iOne;
298 const DataT tempRatioPhi = ratioPhi * iOne2; // Used tobe divided by ( m_one * m_one ) when m_one was != 1
299 const DataT tempRatioZ = ratioZ * iOne2 / (jOne * jOne);
300
301 calcCoefficients(1, tnRRow - 1, h, tempRatioZ, tempRatioPhi, coefficient1, coefficient2, coefficient3, coefficient4);
302
303 // relax on the coarsest level
304 relax3D(tvArrayV[nLoop - 1], tvChargeFMG[nLoop - 1], tnRRow, tnZColumn, mParamGrid.NPhiVertices, symmetry, h2, tempRatioZ, coefficient1, coefficient2, coefficient3, coefficient4);
305
306 // 2) Do multiGrid v-cycle from coarsest to finest
307 for (int count = nLoop - 2; count >= 0; --count) {
308 // move to finer grid
309 iOne /= 2;
310 jOne /= 2;
311 tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1;
312 tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1;
313
314 // 2) a) Interpolate potential for h -> 2h (coarse -> fine)
315 interp3D(tvArrayV[count], tvArrayV[count + 1], tnRRow, tnZColumn, mParamGrid.NPhiVertices, mParamGrid.NPhiVertices);
316
317 // 2) c) Copy the restricted charge to charge for calculation
318 tvCharge[count] = tvChargeFMG[count]; // copy
319
320 // 2) c) Do V cycle MGParameters::nMGCycle times at most
321 for (int mgCycle = 0; mgCycle < MGParameters::nMGCycle; ++mgCycle) {
322 // Copy the potential to temp array for convergence calculation
323 tvPrevArrayV[count] = tvArrayV[count];
324
325 // 2) c) i) Call V cycle from grid count+1 (current fine level) to nLoop (coarsest)
326 vCycle3D2D(symmetry, count + 1, nLoop, MGParameters::nPre, MGParameters::nPost, ratioZ, ratioPhi, tvArrayV, tvCharge, tvResidue, coefficient1, coefficient2, coefficient3, coefficient4, inverseCoefficient4);
327
328 const DataT convergenceError = getConvergenceError(tvArrayV[count], tvPrevArrayV[count]);
329
331 if (convergenceError <= sConvergenceError) {
332 break;
333 }
334 }
335 }
336 }
337
338 // fill output
339 for (int iphi = 0; iphi < mParamGrid.NPhiVertices; ++iphi) {
340 for (int ir = 0; ir < mParamGrid.NRVertices; ++ir) {
341 for (int iz = 0; iz < mParamGrid.NZVertices; ++iz) {
342 matricesV(iz, ir, iphi) = tvArrayV[0](ir, iz, iphi);
343 }
344 }
345 }
346}
347
348template <typename DataT>
349void PoissonSolver<DataT>::poissonMultiGrid3D(DataContainer& matricesV, const DataContainer& matricesCharge, const int symmetry)
350{
351 const DataT gridSpacingR = getSpacingR();
352 const DataT gridSpacingZ = getSpacingZ();
353 const DataT ratioZ = gridSpacingR * gridSpacingR / (gridSpacingZ * gridSpacingZ); // ratio_{Z} = gridSize_{r} / gridSize_{z}
354
355 LOGP(detail, "{}", fmt::format("PoissonMultiGrid3D: in Poisson Solver 3D multi grid full coarsening mParamGrid.NRVertices={}, cols={}, mParamGrid.NPhiVertices={}", mParamGrid.NRVertices, mParamGrid.NZVertices, mParamGrid.NPhiVertices));
356
357 // Check that the number of mParamGrid.NRVertices and mParamGrid.NZVertices is suitable for a binary expansion
358 if (!isPowerOfTwo((mParamGrid.NRVertices - 1))) {
359 LOGP(error, "PoissonMultiGrid3D: Poisson3DMultiGrid - Error in the number of mParamGrid.NRVertices. Must be 2**M + 1");
360 return;
361 }
362 if (!isPowerOfTwo((mParamGrid.NZVertices - 1))) {
363 LOGP(error, "PoissonMultiGrid3D: Poisson3DMultiGrid - Error in the number of mParamGrid.NZVertices. Must be 2**N + 1");
364 return;
365 }
366 if (mParamGrid.NPhiVertices <= 3) {
367 LOGP(error, "PoissonMultiGrid3D: Poisson3DMultiGrid - Error in the number of mParamGrid.NPhiVertices. Must be larger than 3");
368 return;
369 }
370
371 // Solve Poisson's equation in cylindrical coordinates by multi grid technique
372 // Allow for different size grid spacing in R and Z directions
373 int nGridRow = 0; // number grid
374 int nGridCol = 0; // number grid
375 int nGridPhi = 0;
376
377 int nnRow = mParamGrid.NRVertices;
378 while (nnRow >>= 1) {
379 ++nGridRow;
380 }
381
382 int nnCol = mParamGrid.NZVertices;
383 while (nnCol >>= 1) {
384 ++nGridCol;
385 }
386
387 int nnPhi = mParamGrid.NPhiVertices;
388 while (nnPhi % 2 == 0) {
389 ++nGridPhi;
390 nnPhi /= 2;
391 }
392
393 LOGP(detail, "{}", fmt::format("PoissonMultiGrid3D: nGridRow={}, nGridCol={}, nGridPhi={}", nGridRow, nGridCol, nGridPhi));
394 const int nLoop = std::max({nGridRow, nGridCol, nGridPhi}); // Calculate the number of nLoop for the binary expansion
395
396 // Vector for storing multi grid array
397 unsigned int iOne = 1; // index i in gridSize r (original)
398 unsigned int jOne = 1; // index j in gridSize z (original)
399 unsigned int kOne = 1; // index k in gridSize phi
400 int tnRRow = mParamGrid.NRVertices;
401 int tnZColumn = mParamGrid.NZVertices;
402 int tPhiSlice = mParamGrid.NPhiVertices;
403
404 // 1) Memory allocation for multi grid
405 std::vector<Vector> tvArrayV(nLoop); // potential <--> error
406 std::vector<Vector> tvChargeFMG(nLoop); // charge is restricted in full multiGrid
407 std::vector<Vector> tvCharge(nLoop); // charge <--> residue
408 std::vector<Vector> tvPrevArrayV(nLoop); // error calculation
409 std::vector<Vector> tvResidue(nLoop); // residue calculation
410
411 std::vector<DataT> coefficient1(mParamGrid.NRVertices); // coefficient1(mParamGrid.NRVertices) for storing (1 + h_{r}/2r_{i}) from central differences in r direction
412 std::vector<DataT> coefficient2(mParamGrid.NRVertices); // coefficient2(mParamGrid.NRVertices) for storing (1 + h_{r}/2r_{i}) from central differences in r direction
413 std::vector<DataT> coefficient3(mParamGrid.NRVertices); // coefficient3(mParamGrid.NRVertices) for storing (1/r_{i}^2) from central differences in phi direction
414 std::vector<DataT> coefficient4(mParamGrid.NRVertices); // coefficient4(mParamGrid.NRVertices) for storing 1/2
415 std::vector<DataT> inverseCoefficient4(mParamGrid.NRVertices); // inverse of coefficient4(mParamGrid.NRVertices)
416
417 for (int count = 1; count <= nLoop; ++count) {
418 // tnRRow,tnZColumn in new grid
419 tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1;
420 tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1;
421 tPhiSlice = kOne == 1 ? mParamGrid.NPhiVertices : mParamGrid.NPhiVertices / kOne;
422 tPhiSlice = tPhiSlice < nnPhi ? nnPhi : tPhiSlice;
423
424 // allocate memory for residue
425 const int index = count - 1;
426 tvResidue[index].resize(tnRRow, tnZColumn, tPhiSlice);
427 tvPrevArrayV[index].resize(tnRRow, tnZColumn, tPhiSlice);
428 tvChargeFMG[index].resize(tnRRow, tnZColumn, tPhiSlice);
429 tvArrayV[index].resize(tnRRow, tnZColumn, tPhiSlice);
430 tvCharge[index].resize(tnRRow, tnZColumn, tPhiSlice);
431
432 // memory for the finest grid is from parameters
433 if (count == 1) {
434 for (int iphi = 0; iphi < mParamGrid.NPhiVertices; ++iphi) {
435 for (int ir = 0; ir < mParamGrid.NRVertices; ++ir) {
436 for (int iz = 0; iz < mParamGrid.NZVertices; ++iz) {
437 tvChargeFMG[index](ir, iz, iphi) = matricesCharge(iz, ir, iphi);
438 tvArrayV[index](ir, iz, iphi) = matricesV(iz, ir, iphi);
439 }
440 }
441 }
442 tvCharge[index] = tvChargeFMG[index];
443 }
444 iOne *= 2; // doubling
445 jOne *= 2; // doubling
446 kOne *= 2;
447 }
448
449 // Case full multi grid (FMG)
450 if (MGParameters::cycleType == CycleType::FCycle) {
451 // Restrict the charge to coarser grid
452 iOne = 2;
453 jOne = 2;
454 kOne = 2;
455 int otPhiSlice = mParamGrid.NPhiVertices;
456
457 // 1) Restrict Charge and Boundary to coarser grid
458 for (int count = 2; count <= nLoop; ++count) {
459 // tnRRow,tnZColumn in new grid
460 tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1;
461 tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1;
462 tPhiSlice = kOne == 1 ? mParamGrid.NPhiVertices : mParamGrid.NPhiVertices / kOne;
463 tPhiSlice = tPhiSlice < nnPhi ? nnPhi : tPhiSlice;
464
465 LOGP(detail, "{}", fmt::format("PoissonMultiGrid3D: Restrict3D, tnRRow={}, tnZColumn={}, newPhiSlice={}, oldPhiSlice={}", tnRRow, tnZColumn, tPhiSlice, otPhiSlice));
466 restrict3D(tvChargeFMG[count - 1], tvChargeFMG[count - 2], tnRRow, tnZColumn, tPhiSlice, otPhiSlice);
467 // copy boundary values of V
468 restrictBoundary3D(tvArrayV[count - 1], tvArrayV[count - 2], tnRRow, tnZColumn, tPhiSlice, otPhiSlice);
469 otPhiSlice = tPhiSlice;
470
471 iOne *= 2; // doubling
472 jOne *= 2; // doubling
473 kOne *= 2;
474 }
475
476 // Relax on the coarsest grid
477 // FMG
478 // 2) Relax on the coarsest grid
479 // move to the coarsest + 1
480 iOne /= 2;
481 jOne /= 2;
482 kOne /= 2;
483 tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1;
484 tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1;
485 tPhiSlice = kOne == 1 ? mParamGrid.NPhiVertices : mParamGrid.NPhiVertices / kOne;
486 tPhiSlice = tPhiSlice < nnPhi ? nnPhi : tPhiSlice;
487 otPhiSlice = tPhiSlice;
488
489 const DataT h = gridSpacingR * iOne;
490 const DataT h2 = h * h;
491 const DataT gridSizePhiInv = tPhiSlice * getGridSizePhiInv(); // h_{phi}
492 const DataT tempRatioPhi = h2 * gridSizePhiInv * gridSizePhiInv; // ratio_{phi} = gridSize_{r} / gridSize_{phi}
493 const DataT tempRatioZ = ratioZ * iOne * iOne / (jOne * jOne);
494
495 calcCoefficients(1, tnRRow - 1, h, tempRatioZ, tempRatioPhi, coefficient1, coefficient2, coefficient3, coefficient4);
496
497 // 3) Relax on the coarsest grid
498 relax3D(tvArrayV[nLoop - 1], tvChargeFMG[nLoop - 1], tnRRow, tnZColumn, tPhiSlice, symmetry, h2, tempRatioZ, coefficient1, coefficient2, coefficient3, coefficient4);
499
500 // 4) V Cycle from coarsest to finest
501 for (int count = nLoop - 2; count >= 0; --count) {
502 // move to finer grid
503 std::fill(std::begin(coefficient1), std::end(coefficient1), 0);
504 std::fill(std::begin(coefficient2), std::end(coefficient2), 0);
505 std::fill(std::begin(coefficient3), std::end(coefficient3), 0);
506 std::fill(std::begin(coefficient4), std::end(coefficient4), 0);
507 std::fill(std::begin(inverseCoefficient4), std::end(inverseCoefficient4), 0);
508
509 iOne /= 2;
510 jOne /= 2;
511 kOne /= 2;
512
513 tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1;
514 tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1;
515 tPhiSlice = kOne == 1 ? mParamGrid.NPhiVertices : mParamGrid.NPhiVertices / kOne;
516 tPhiSlice = tPhiSlice < nnPhi ? nnPhi : tPhiSlice;
517
518 // 4) a) interpolate from 2h --> h grid
519 interp3D(tvArrayV[count], tvArrayV[count + 1], tnRRow, tnZColumn, tPhiSlice, otPhiSlice);
520
521 // Copy the relax charge to the tvCharge
522 if (count > 0) {
523 tvCharge[count] = tvChargeFMG[count];
524 }
525
526 using timer = std::chrono::high_resolution_clock;
527 auto start = timer::now();
528 for (int mgCycle = 0; mgCycle < MGParameters::nMGCycle; ++mgCycle) {
529 // copy to store previous potential
530 tvPrevArrayV[count] = tvArrayV[count];
531
532 vCycle3D(symmetry, count + 1, nLoop, MGParameters::nPre, MGParameters::nPost, ratioZ, tvArrayV, tvCharge, tvResidue, coefficient1, coefficient2, coefficient3, coefficient4, inverseCoefficient4);
533
534 // converge error
535 const DataT convergenceError = getConvergenceError(tvArrayV[count], tvPrevArrayV[count]);
536 // if already converge just break move to finer grid
537 if (convergenceError <= sConvergenceError) {
538 LOGP(detail, "Cycle converged. Continue to next cycle...");
539 break;
540 }
541 if (count <= 1 && !(mgCycle % 10)) {
542 auto stop = timer::now();
543 std::chrono::duration<float> time = stop - start;
544 const float totalTime = time.count();
545 const float timePerCycle = totalTime / (mgCycle + 1);
546 const float remaining = timePerCycle * (MGParameters::nMGCycle - mgCycle);
547 LOGP(detail, "Cycle {} out of {} for current V cycle {}. Processed time {}s with {}s per cycle. Max remaining time for current cycle {}s. Convergence {} > {}", mgCycle, MGParameters::nMGCycle, count, time.count(), timePerCycle, remaining, convergenceError, sConvergenceError);
548 }
549 if (mgCycle == (MGParameters::nMGCycle - 1)) {
550 LOGP(warning, "Cycle {} did not convergence! Current convergence error is larger than expected convergence error: {} > {}", mgCycle, convergenceError, sConvergenceError);
551 }
552 }
553 // keep old slice information
554 otPhiSlice = tPhiSlice;
555 }
556 } else if (MGParameters::cycleType == CycleType::VCycle) {
557 // V-cycle
558 int gridFrom = 1;
559 int gridTo = nLoop;
560
561 for (int mgCycle = 0; mgCycle < MGParameters::nMGCycle; ++mgCycle) {
562 // copy to store previous potential
563 tvPrevArrayV[0] = tvArrayV[0];
564
565 // Do V Cycle from the coarsest to finest grid
566 vCycle3D(symmetry, gridFrom, gridTo, MGParameters::nPre, MGParameters::nPost, ratioZ, tvArrayV, tvCharge, tvResidue, coefficient1, coefficient2, coefficient3, coefficient4, inverseCoefficient4);
567
568 // convergence error
569 const DataT convergenceError = getConvergenceError(tvArrayV[0], tvPrevArrayV[0]);
570
571 // if error already achieved then stop mg iteration
572 if (convergenceError <= sConvergenceError) {
573 break;
574 }
575 }
576 }
577
578 // fill output
579 for (int iphi = 0; iphi < mParamGrid.NPhiVertices; ++iphi) {
580 for (int ir = 0; ir < mParamGrid.NRVertices; ++ir) {
581 for (int iz = 0; iz < mParamGrid.NZVertices; ++iz) {
582 matricesV(iz, ir, iphi) = tvArrayV[0](ir, iz, iphi);
583 }
584 }
585 }
586}
587
588template <typename DataT>
589void PoissonSolver<DataT>::wCycle2D(const int gridFrom, const int gridTo, const int gamma, const int nPre, const int nPost, const DataT gridSizeR, const DataT ratio,
590 std::vector<Vector>& tvArrayV, std::vector<Vector>& tvCharge, std::vector<Vector>& tvResidue)
591{
592 unsigned int iOne = 1 << (gridFrom - 1);
593 unsigned int jOne = 1 << (gridFrom - 1);
594
595 int tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1;
596 int tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1;
597
598 std::vector<DataT> coefficient1(mParamGrid.NRVertices);
599 std::vector<DataT> coefficient2(mParamGrid.NZVertices);
600
601 // 1) Go to coarsest level
602 for (int count = gridFrom; count <= gridTo - 2; ++count) {
603 const int index = count - 1;
604 const DataT h = gridSizeR * iOne;
605 const DataT h2 = h * h;
606 const DataT ih2 = 1 / h2;
607 const DataT tempRatio = ratio * iOne * iOne / (jOne * jOne);
608 const DataT tempFourth = 1 / (2 + 2 * tempRatio);
609 const DataT inverseTempFourth = 1 / tempFourth;
610 calcCoefficients2D(1, tnRRow - 1, h, coefficient1, coefficient2);
611 Vector matricesCurrentV = tvArrayV[index];
612 Vector matricesCurrentCharge = tvCharge[index];
613 Vector residue = tvResidue[index];
614
615 // 1) Pre-Smoothing: Gauss-Seidel Relaxation or Jacobi
616 for (int jPre = 1; jPre <= nPre; ++jPre) {
617 relax2D(matricesCurrentV, matricesCurrentCharge, tnRRow, tnZColumn, h2, tempFourth, tempRatio, coefficient1, coefficient2);
618 }
619
620 // 2) Residue calculation
621 residue2D(residue, matricesCurrentV, matricesCurrentCharge, tnRRow, tnZColumn, ih2, inverseTempFourth, tempRatio, coefficient1, coefficient2);
622
623 iOne *= 2;
624 jOne *= 2;
625 tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1;
626 tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1;
627
628 matricesCurrentCharge = tvCharge[count];
629 matricesCurrentV = tvArrayV[count];
630
631 // 3) Restriction
632 restrict2D(matricesCurrentCharge, residue, tnRRow, tnZColumn, 0);
633 }
634
635 // Do V cycle from: gridTo-1 to gridTo gamma times
636 for (int iGamma = 0; iGamma < gamma; ++iGamma) {
637 vCycle2D(gridTo - 1, gridTo, nPre, nPost, gridSizeR, ratio, tvArrayV, tvCharge, tvResidue);
638 }
639
640 // Go to finest grid
641 for (int count = gridTo - 2; count >= gridFrom; --count) {
642 iOne /= 2;
643 jOne /= 2;
644
645 const DataT h = gridSizeR * iOne;
646 const DataT h2 = h * h;
647 const DataT tempRatio = ratio * iOne * iOne / (jOne * jOne);
648 const DataT tempFourth = 1 / (2 + 2 * tempRatio);
649
650 tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1;
651 tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1;
652 const Vector matricesCurrentCharge = tvCharge[count - 1];
653 Vector matricesCurrentV = tvArrayV[count - 1];
654 const Vector matricesCurrentVC = tvArrayV[count];
655
656 // 6) Interpolation/Prolongation
657 addInterp2D(matricesCurrentV, matricesCurrentVC, tnRRow, tnZColumn, 0);
658
659 calcCoefficients2D(1, tnRRow - 1, h, coefficient1, coefficient2);
660
661 // 7) Post-Smoothing: Gauss-Seidel Relaxation
662 for (Int_t jPost = 1; jPost <= nPost; ++jPost) {
663 relax2D(matricesCurrentV, matricesCurrentCharge, tnRRow, tnZColumn, h2, tempFourth, tempRatio, coefficient1, coefficient2);
664 } // end post smoothing
665 }
666}
667
668template <typename DataT>
669void PoissonSolver<DataT>::vCycle2D(const int gridFrom, const int gridTo, const int nPre, const int nPost, const DataT gridSizeR, const DataT ratio, std::vector<Vector>& tvArrayV,
670 std::vector<Vector>& tvCharge, std::vector<Vector>& tvResidue)
671{
672 unsigned int iOne = 1 << (gridFrom - 1);
673 unsigned int jOne = 1 << (gridFrom - 1);
674
675 int tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1;
676 int tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1;
677
678 std::vector<DataT> coefficient1(mParamGrid.NRVertices);
679 std::vector<DataT> coefficient2(mParamGrid.NZVertices);
680
681 // 1) Go to coarsest level
682 for (int count = gridFrom; count <= gridTo - 1; ++count) {
683 const int index = count - 1;
684 const DataT h = gridSizeR * iOne;
685 const DataT h2 = h * h;
686 const DataT ih2 = 1 / h2;
687 const DataT tempRatio = ratio * iOne * iOne / (jOne * jOne);
688 const DataT tempFourth = 1 / (2 + 2 * tempRatio);
689 const DataT inverseTempFourth = 1 / tempFourth;
690 calcCoefficients2D(1, tnRRow - 1, h, coefficient1, coefficient2);
691
692 for (int jPre = 1; jPre <= nPre; ++jPre) {
693 relax2D(tvArrayV[index], tvCharge[index], tnRRow, tnZColumn, h2, tempFourth, tempRatio, coefficient1, coefficient2);
694 }
695
696 // 2) Residue calculation
697 residue2D(tvResidue[index], tvArrayV[index], tvCharge[index], tnRRow, tnZColumn, ih2, inverseTempFourth, tempRatio, coefficient1, coefficient2);
698
699 iOne *= 2;
700 jOne *= 2;
701 tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1;
702 tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1;
703
704 // 3) Restriction
705 restrict2D(tvCharge[count], tvResidue[index], tnRRow, tnZColumn, 0);
706
707 // 4) Zeroing coarser V
708 std::fill(tvArrayV[count].begin(), tvArrayV[count].end(), 0); // is this necessary???
709 }
710
711 // 5) coarsest grid
712 const DataT h = gridSizeR * iOne;
713 const DataT h2 = h * h;
714 const DataT tempRatio = ratio * iOne * iOne / (jOne * jOne);
715 const DataT tempFourth = 1 / (2 + 2 * tempRatio);
716 calcCoefficients2D(1, tnRRow - 1, h, coefficient1, coefficient2);
717
718 relax2D(tvArrayV[gridTo - 1], tvCharge[gridTo - 1], tnRRow, tnZColumn, h2, tempFourth, tempRatio, coefficient1, coefficient2);
719
720 // Go to finest grid
721 for (int count = gridTo - 1; count >= gridFrom; count--) {
722 const int index = count - 1;
723 iOne /= 2;
724 jOne /= 2;
725
726 const DataT hTmp = gridSizeR * iOne;
727 const DataT h2Tmp = hTmp * hTmp;
728 const DataT tempRatioTmp = ratio * iOne * iOne / (jOne * jOne);
729 const DataT tempFourthTmp = 1 / (2 + 2 * tempRatioTmp);
730
731 tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1;
732 tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1;
733
734 // 6) Interpolation/Prolongation
735 addInterp2D(tvArrayV[index], tvArrayV[count], tnRRow, tnZColumn, 0);
736
737 calcCoefficients2D(1, tnRRow - 1, hTmp, coefficient1, coefficient2);
738
739 // 7) Post-Smoothing: Gauss-Seidel Relaxation
740 for (int jPost = 1; jPost <= nPost; ++jPost) {
741 relax2D(tvArrayV[index], tvCharge[index], tnRRow, tnZColumn, h2Tmp, tempFourthTmp, tempRatioTmp, coefficient1, coefficient2);
742 } // end post smoothing
743 }
744}
745
746template <typename DataT>
747void PoissonSolver<DataT>::vCycle3D2D(const int symmetry, const int gridFrom, const int gridTo, const int nPre, const int nPost, const DataT ratioZ, const DataT ratioPhi,
748 std::vector<Vector>& tvArrayV, std::vector<Vector>& tvCharge, std::vector<Vector>& tvResidue, std::vector<DataT>& coefficient1,
749 std::vector<DataT>& coefficient2, std::vector<DataT>& coefficient3, std::vector<DataT>& coefficient4, std::vector<DataT>& inverseCoefficient4) const
750{
751 unsigned int iOne = 1 << (gridFrom - 1);
752 unsigned int jOne = 1 << (gridFrom - 1);
753 int tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1;
754 int tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1;
755
756 for (int count = gridFrom; count <= gridTo - 1; ++count) {
757 const int index = count - 1;
758 const DataT h = getSpacingR() * iOne;
759 const DataT h2 = h * h;
760 const DataT ih2 = 1 / h2;
761 const int iOne2 = iOne * iOne;
762 const DataT tempRatioPhi = ratioPhi * iOne2; // Used tobe divided by ( m_one * m_one ) when m_one was != 1
763 const DataT tempRatioZ = ratioZ * iOne2 / (jOne * jOne);
764
765 calcCoefficients(1, tnRRow - 1, h, tempRatioZ, tempRatioPhi, coefficient1, coefficient2, coefficient3, coefficient4);
766 for (unsigned int i = 1; i < tnRRow - 1; ++i) {
767 inverseCoefficient4[i] = 1 / coefficient4[i];
768 }
769
770 // Info("VCycle3D2D","Before Pre-smoothing");
771 // 1) Pre-Smoothing: Gauss-Seidel Relaxation or Jacobi
772 for (int jPre = 1; jPre <= nPre; ++jPre) {
773 relax3D(tvArrayV[index], tvCharge[index], tnRRow, tnZColumn, mParamGrid.NPhiVertices, symmetry, h2, tempRatioZ, coefficient1, coefficient2, coefficient3, coefficient4);
774 } // end pre smoothing
775
776 // 2) Residue calculation
777 residue3D(tvResidue[index], tvArrayV[index], tvCharge[index], tnRRow, tnZColumn, mParamGrid.NPhiVertices, symmetry, ih2, tempRatioZ, coefficient1, coefficient2, coefficient3, inverseCoefficient4);
778
779 iOne *= 2;
780 jOne *= 2;
781 tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1;
782 tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1;
783
784 // 3) Restriction
785 restrict3D(tvCharge[count], tvResidue[index], tnRRow, tnZColumn, mParamGrid.NPhiVertices, mParamGrid.NPhiVertices);
786
787 // 4) Zeroing coarser V
788 std::fill(tvArrayV[count].begin(), tvArrayV[count].end(), 0);
789 }
790
791 // coarsest grid
792 const DataT hTmp = getSpacingR() * iOne;
793 const DataT h2Tmp = hTmp * hTmp;
794
795 const int iOne2Tmp = iOne * iOne;
796 const DataT tempRatioPhiTmp = ratioPhi * iOne2Tmp;
797 const DataT tempRatioZTmp = ratioZ * iOne2Tmp / (jOne * jOne);
798
799 calcCoefficients(1, tnRRow - 1, hTmp, tempRatioZTmp, tempRatioPhiTmp, coefficient1, coefficient2, coefficient3, coefficient4);
800
801 // 3) Relax on the coarsest grid
802 relax3D(tvArrayV[gridTo - 1], tvCharge[gridTo - 1], tnRRow, tnZColumn, mParamGrid.NPhiVertices, symmetry, h2Tmp, tempRatioZTmp, coefficient1, coefficient2, coefficient3, coefficient4);
803
804 // back to fine
805 for (int count = gridTo - 1; count >= gridFrom; --count) {
806 const int index = count - 1;
807 iOne /= 2;
808 jOne /= 2;
809
810 tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1;
811 tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1;
812
813 const DataT h = getSpacingR() * iOne;
814 const DataT h2 = h * h;
815 const int iOne2 = iOne * iOne;
816 const DataT tempRatioPhi = ratioPhi * iOne2; // Used tobe divided by ( m_one * m_one ) when m_one was != 1
817 const DataT tempRatioZ = ratioZ * iOne2 / (jOne * jOne);
818
819 // 4) Interpolation/Prolongation
820 addInterp3D(tvArrayV[index], tvArrayV[count], tnRRow, tnZColumn, mParamGrid.NPhiVertices, mParamGrid.NPhiVertices);
821
822 calcCoefficients(1, tnRRow - 1, h, tempRatioZ, tempRatioPhi, coefficient1, coefficient2, coefficient3, coefficient4);
823
824 // 5) Post-Smoothing: Gauss-Seidel Relaxation
825 for (int jPost = 1; jPost <= nPost; ++jPost) {
826 relax3D(tvArrayV[index], tvCharge[index], tnRRow, tnZColumn, mParamGrid.NPhiVertices, symmetry, h2, tempRatioZ, coefficient1, coefficient2, coefficient3, coefficient4);
827 } // end post smoothing
828 }
829}
830
831template <typename DataT>
832void PoissonSolver<DataT>::vCycle3D(const int symmetry, const int gridFrom, const int gridTo, const int nPre, const int nPost, const DataT ratioZ, std::vector<Vector>& tvArrayV,
833 std::vector<Vector>& tvCharge, std::vector<Vector>& tvResidue, std::vector<DataT>& coefficient1, std::vector<DataT>& coefficient2, std::vector<DataT>& coefficient3,
834 std::vector<DataT>& coefficient4, std::vector<DataT>& inverseCoefficient4) const
835{
836 const DataT gridSpacingR = getSpacingR();
837
838 unsigned int iOne = 1 << (gridFrom - 1);
839 unsigned int jOne = 1 << (gridFrom - 1);
840 unsigned int kOne = 1 << (gridFrom - 1);
841
842 int nnPhi = mParamGrid.NPhiVertices;
843 while (nnPhi % 2 == 0) {
844 nnPhi /= 2;
845 }
846
847 int tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1;
848 int tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1;
849 int tPhiSlice = kOne == 1 ? mParamGrid.NPhiVertices : mParamGrid.NPhiVertices / kOne;
850 tPhiSlice = tPhiSlice < nnPhi ? nnPhi : tPhiSlice;
851
852 for (int count = gridFrom; count <= gridTo - 1; ++count) {
853 const int index = count - 1;
854 const int otPhiSlice = tPhiSlice;
855 const DataT h = gridSpacingR * iOne;
856 const DataT h2 = h * h;
857 const DataT ih2 = 1 / h2;
858 const DataT tempGridSizePhiInv = tPhiSlice * getGridSizePhiInv(); // phi now is multiGrid
859 const DataT tempRatioPhi = h2 * tempGridSizePhiInv * tempGridSizePhiInv; // ratio_{phi} = gridSize_{r} / gridSize_{phi}
860 const DataT tempRatioZ = ratioZ * iOne * iOne / (jOne * jOne);
861
862 calcCoefficients(1, tnRRow - 1, h, tempRatioZ, tempRatioPhi, coefficient1, coefficient2, coefficient3, coefficient4);
863
864 for (int i = 1; i < tnRRow - 1; ++i) {
865 inverseCoefficient4[i] = 1 / coefficient4[i];
866 }
867
868 // 1) Pre-Smoothing: Gauss-Seidel Relaxation or Jacobi
869 for (int jPre = 1; jPre <= nPre; ++jPre) {
870 relax3D(tvArrayV[index], tvCharge[index], tnRRow, tnZColumn, tPhiSlice, symmetry, h2, tempRatioZ, coefficient1, coefficient2, coefficient3, coefficient4);
871 } // end pre smoothing
872
873 // 2) Residue calculation
874 residue3D(tvResidue[index], tvArrayV[index], tvCharge[index], tnRRow, tnZColumn, tPhiSlice, symmetry, ih2, tempRatioZ, coefficient1, coefficient2, coefficient3, inverseCoefficient4);
875
876 iOne *= 2;
877 jOne *= 2;
878 kOne *= 2;
879 tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1;
880 tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1;
881 tPhiSlice = mParamGrid.NPhiVertices / kOne;
882 tPhiSlice = tPhiSlice < nnPhi ? nnPhi : tPhiSlice;
883
884 // 3) Restriction
885 restrict3D(tvCharge[count], tvResidue[index], tnRRow, tnZColumn, tPhiSlice, otPhiSlice);
886
887 // 4) Zeroing coarser V
888 std::fill(tvArrayV[count].begin(), tvArrayV[count].end(), 0);
889 }
890
891 // coarsest grid
892 const DataT hTmp = gridSpacingR * iOne;
893 const DataT h2Tmp = hTmp * hTmp;
894 const DataT tempGridSizePhiInvTmp = tPhiSlice * getGridSizePhiInv(); // phi now is multiGrid
895 const DataT tempRatioPhiTmp = h2Tmp * tempGridSizePhiInvTmp * tempGridSizePhiInvTmp; // ratio_{phi} = gridSize_{r} / gridSize_{phi}
896 const DataT tempRatioZTmp = ratioZ * iOne * iOne / (jOne * jOne);
897
898 calcCoefficients(1, tnRRow - 1, hTmp, tempRatioZTmp, tempRatioPhiTmp, coefficient1, coefficient2, coefficient3, coefficient4);
899
900 // 3) Relax on the coarsest grid
901 relax3D(tvArrayV[gridTo - 1], tvCharge[gridTo - 1], tnRRow, tnZColumn, tPhiSlice, symmetry, h2Tmp, tempRatioZTmp, coefficient1, coefficient2, coefficient3, coefficient4);
902 // back to fine
903 for (int count = gridTo - 1; count >= gridFrom; --count) {
904 const int index = count - 1;
905 const int otPhiSlice = tPhiSlice;
906 iOne /= 2;
907 jOne /= 2;
908 kOne /= 2;
909
910 tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1;
911 tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1;
912 tPhiSlice = kOne == 1 ? mParamGrid.NPhiVertices : mParamGrid.NPhiVertices / kOne;
913 tPhiSlice = tPhiSlice < nnPhi ? nnPhi : tPhiSlice;
914
915 const DataT h = gridSpacingR * iOne;
916 const DataT h2 = h * h;
917 const DataT tempGridSizePhiInv = tPhiSlice * getGridSizePhiInv();
918 const DataT tempRatioPhi = h2 * tempGridSizePhiInv * tempGridSizePhiInv; // ratio_{phi} = gridSize_{r} / gridSize_{phi}
919 const DataT tempRatioZ = ratioZ * iOne * iOne / (jOne * jOne);
920
921 // 4) Interpolation/Prolongation
922 addInterp3D(tvArrayV[index], tvArrayV[count], tnRRow, tnZColumn, tPhiSlice, otPhiSlice);
923
924 calcCoefficients(1, tnRRow - 1, h, tempRatioZ, tempRatioPhi, coefficient1, coefficient2, coefficient3, coefficient4);
925
926 // 5) Post-Smoothing: Gauss-Seidel Relaxation
927 for (int jPost = 1; jPost <= nPost; ++jPost) {
928 relax3D(tvArrayV[index], tvCharge[index], tnRRow, tnZColumn, tPhiSlice, symmetry, h2, tempRatioZ, coefficient1, coefficient2, coefficient3, coefficient4);
929 }
930 }
931}
932
933template <typename DataT>
934void PoissonSolver<DataT>::residue2D(Vector& residue, const Vector& matricesCurrentV, const Vector& matricesCurrentCharge, const int tnRRow, const int tnZColumn, const DataT ih2, const DataT inverseTempFourth,
935 const DataT tempRatio, std::vector<DataT>& coefficient1, std::vector<DataT>& coefficient2)
936{
937 const int iPhi = 0;
938#pragma omp parallel for num_threads(sNThreads)
939 for (int i = 1; i < tnRRow - 1; ++i) {
940 for (int j = 1; j < tnZColumn - 1; ++j) {
941 residue(i, j, iPhi) = ih2 * (coefficient1[i] * matricesCurrentV(i + 1, j, iPhi) + coefficient2[i] * matricesCurrentV(i - 1, j, iPhi) + tempRatio * (matricesCurrentV(i, j + 1, iPhi) + matricesCurrentV(i, j - 1, iPhi)) - inverseTempFourth * matricesCurrentV(i, j, iPhi)) + matricesCurrentCharge(i, j, iPhi);
942 } // end cols
943 } // end nRRow
944
945 // Boundary points.
946 for (int i = 0; i < tnRRow; ++i) {
947 residue(i, 0, iPhi) = residue(i, tnZColumn - 1, iPhi) = 0.0;
948 }
949
950 for (int j = 0; j < tnZColumn; ++j) {
951 residue(0, j, iPhi) = residue(tnRRow - 1, j, iPhi) = 0.0;
952 }
953}
954
955template <typename DataT>
956void PoissonSolver<DataT>::residue3D(Vector& residue, const Vector& matricesCurrentV, const Vector& matricesCurrentCharge, const int tnRRow, const int tnZColumn, const int tnPhi, const int symmetry,
957 const DataT ih2, const DataT tempRatioZ, const std::vector<DataT>& coefficient1, const std::vector<DataT>& coefficient2, const std::vector<DataT>& coefficient3, const std::vector<DataT>& inverseCoefficient4) const
958{
959#pragma omp parallel for num_threads(sNThreads) // parallising this loop is possible - but using more than 2 cores makes it slower -
960 for (int m = 0; m < tnPhi; ++m) {
961 int mp1 = m + 1;
962 int signPlus = 1;
963 int mm1 = m - 1;
964 int signMinus = 1;
965
966 // Reflection symmetry in phi (e.g. symmetry at sector boundaries, or half sectors, etc.)
967 if (symmetry == 1) {
968 if (mp1 > tnPhi - 1) {
969 mp1 = tnPhi - 2;
970 }
971 if (mm1 < 0) {
972 mm1 = 1;
973 }
974 }
975 // Anti-symmetry in phi
976 else if (symmetry == -1) {
977 if (mp1 > tnPhi - 1) {
978 mp1 = tnPhi - 2;
979 signPlus = -1;
980 }
981 if (mm1 < 0) {
982 mm1 = 1;
983 signMinus = -1;
984 }
985 } else { // No Symmetries in phi, no boundaries, the calculation is continuous across all phi
986 if (mp1 > tnPhi - 1) {
987 mp1 = m + 1 - tnPhi;
988 }
989 if (mm1 < 0) {
990 mm1 = m - 1 + tnPhi;
991 }
992 }
993
994 for (int j = 1; j < tnZColumn - 1; ++j) {
995 for (int i = 1; i < tnRRow - 1; ++i) {
996 residue(i, j, m) = ih2 * (coefficient2[i] * matricesCurrentV(i - 1, j, m) + tempRatioZ * (matricesCurrentV(i, j - 1, m) + matricesCurrentV(i, j + 1, m)) + coefficient1[i] * matricesCurrentV(i + 1, j, m) +
997 coefficient3[i] * (signPlus * matricesCurrentV(i, j, mp1) + signMinus * matricesCurrentV(i, j, mm1)) - inverseCoefficient4[i] * matricesCurrentV(i, j, m)) +
998 matricesCurrentCharge(i, j, m);
999 } // end cols
1000 } // end mParamGrid.NRVertices
1001 }
1002}
1003
1004template <typename DataT>
1005void PoissonSolver<DataT>::interp3D(Vector& matricesCurrentV, const Vector& matricesCurrentVC, const int tnRRow, const int tnZColumn, const int newPhiSlice, const int oldPhiSlice) const
1006{
1007 // Do restrict 2 D for each slice
1008 if (newPhiSlice == 2 * oldPhiSlice) {
1009 for (int m = 0; m < newPhiSlice; m += 2) {
1010 // assuming no symmetry
1011 int mm = m / 2;
1012 int mmPlus = mm + 1;
1013 int mp1 = m + 1;
1014
1015 // round
1016 if (mmPlus > oldPhiSlice - 1) {
1017 mmPlus = mm + 1 - oldPhiSlice;
1018 }
1019 if (mp1 > newPhiSlice - 1) {
1020 mp1 = m + 1 - newPhiSlice;
1021 }
1022
1023 for (int j = 2; j < tnZColumn - 1; j += 2) {
1024 for (int i = 2; i < tnRRow - 1; i += 2) {
1025 const int iHalf = i / 2;
1026 const int jHalf = j / 2;
1027 matricesCurrentV(i, j, m) = matricesCurrentVC(iHalf, jHalf, mm);
1028 // point on corner lines at phi direction
1029 matricesCurrentV(i, j, mp1) = (matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf, mmPlus)) / 2;
1030 }
1031 }
1032
1033 for (int j = 1; j < tnZColumn - 1; j += 2) {
1034 for (int i = 2; i < tnRRow - 1; i += 2) {
1035 const int iHalf = i / 2;
1036 const int jHalf = j / 2;
1037 matricesCurrentV(i, j, m) = (matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf + 1, mm)) / 2;
1038 // point on corner lines at phi direction
1039 matricesCurrentV(i, j, mp1) = (matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf + 1, mm) + matricesCurrentVC(iHalf, jHalf, mmPlus) + matricesCurrentVC(iHalf, jHalf + 1, mmPlus)) / 4;
1040 }
1041 }
1042
1043 for (int j = 2; j < tnZColumn - 1; j += 2) {
1044 for (int i = 1; i < tnRRow - 1; i += 2) {
1045 const int iHalf = i / 2;
1046 const int jHalf = j / 2;
1047 matricesCurrentV(i, j, m) = (matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf + 1, jHalf, mm)) / 2;
1048 // point on line at phi direction
1049 matricesCurrentV(i, j, mp1) = ((matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf, mmPlus)) + (matricesCurrentVC(iHalf + 1, jHalf, mmPlus) + matricesCurrentVC(iHalf + 1, jHalf, mm))) / 4;
1050 }
1051 }
1052
1053 for (int j = 1; j < tnZColumn - 1; j += 2) {
1054 for (int i = 1; i < tnRRow - 1; i += 2) {
1055 const int iHalf = i / 2;
1056 const int jHalf = j / 2;
1057 matricesCurrentV(i, j, m) = ((matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf + 1, mm)) + (matricesCurrentVC(iHalf + 1, jHalf, mm) + matricesCurrentVC(iHalf + 1, jHalf + 1, mm))) / 4;
1058 // point at the center at phi direction
1059 matricesCurrentV(i, j, mp1) = ((matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf + 1, mm) + matricesCurrentVC(iHalf, jHalf, mmPlus) + matricesCurrentVC(iHalf, jHalf + 1, mmPlus)) +
1060 (matricesCurrentVC(iHalf + 1, jHalf, mm) + matricesCurrentVC(iHalf + 1, jHalf + 1, mm) + matricesCurrentVC(iHalf + 1, jHalf, mmPlus) + matricesCurrentVC(iHalf + 1, jHalf + 1, mmPlus))) /
1061 8;
1062 }
1063 }
1064 }
1065
1066 } else {
1067#pragma omp parallel for num_threads(sNThreads) // no change
1068 for (int m = 0; m < newPhiSlice; ++m) {
1069 interp2D(matricesCurrentV, matricesCurrentVC, tnRRow, tnZColumn, m);
1070 }
1071 }
1072}
1073
1074template <typename DataT>
1075void PoissonSolver<DataT>::interp2D(Vector& matricesCurrentV, const Vector& matricesCurrentVC, const int tnRRow, const int tnZColumn, const int iphi) const
1076{
1077 for (int j = 2; j < tnZColumn - 1; j += 2) {
1078 for (int i = 2; i < tnRRow - 1; i += 2) {
1079 const int jHalf = j / 2;
1080 matricesCurrentV(i, j, iphi) = matricesCurrentVC(i / 2, jHalf, iphi);
1081 }
1082 }
1083
1084 for (int j = 1; j < tnZColumn - 1; j += 2) {
1085 for (int i = 2; i < tnRRow - 1; i += 2) {
1086 const int iHalf = i / 2;
1087 const int jHalf = j / 2;
1088 matricesCurrentV(i, j, iphi) = (matricesCurrentVC(iHalf, jHalf, iphi) + matricesCurrentVC(iHalf, jHalf + 1, iphi)) / 2;
1089 }
1090 }
1091
1092 for (int j = 2; j < tnZColumn - 1; j += 2) {
1093 for (int i = 1; i < tnRRow - 1; i += 2) {
1094 const int iHalf = i / 2;
1095 const int jHalf = j / 2;
1096 matricesCurrentV(i, j, iphi) = (matricesCurrentVC(iHalf, jHalf, iphi) + matricesCurrentVC(iHalf + 1, jHalf, iphi)) / 2;
1097 }
1098 }
1099
1100 // only if full
1101 if (MGParameters::gtType == GridTransferType::Full) {
1102 for (int j = 1; j < tnZColumn - 1; j += 2) {
1103 for (int i = 1; i < tnRRow - 1; i += 2) {
1104 const int iHalf = i / 2;
1105 const int jHalf = j / 2;
1106 matricesCurrentV(i, j, iphi) = (matricesCurrentVC(iHalf, jHalf, iphi) + matricesCurrentVC(iHalf, jHalf + 1, iphi) + matricesCurrentVC(iHalf + 1, jHalf, iphi) + matricesCurrentVC(iHalf + 1, jHalf + 1, iphi)) / 4;
1107 }
1108 }
1109 }
1110}
1111
1112template <typename DataT>
1113void PoissonSolver<DataT>::addInterp3D(Vector& matricesCurrentV, const Vector& matricesCurrentVC, const int tnRRow, const int tnZColumn, const int newPhiSlice, const int oldPhiSlice) const
1114{
1115 // Do restrict 2 D for each slice
1116 if (newPhiSlice == 2 * oldPhiSlice) {
1117 for (int m = 0; m < newPhiSlice; m += 2) {
1118 // assuming no symmetry
1119 int mm = m / 2;
1120 int mmPlus = mm + 1;
1121 int mp1 = m + 1;
1122
1123 // round
1124 if (mmPlus > (oldPhiSlice)-1) {
1125 mmPlus = mm + 1 - (oldPhiSlice);
1126 }
1127 if (mp1 > (newPhiSlice)-1) {
1128 mp1 = m + 1 - (newPhiSlice);
1129 }
1130
1131 for (int j = 2; j < tnZColumn - 1; j += 2) {
1132 for (int i = 2; i < tnRRow - 1; i += 2) {
1133 const int iHalf = i / 2;
1134 const int jHalf = j / 2;
1135 matricesCurrentV(i, j, m) += matricesCurrentVC(iHalf, jHalf, mm);
1136 // point on corner lines at phi direction
1137 matricesCurrentV(i, j, mp1) += (matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf, mmPlus)) / 2;
1138 }
1139 }
1140
1141 for (int j = 1; j < tnZColumn - 1; j += 2) {
1142 for (int i = 2; i < tnRRow - 1; i += 2) {
1143 const int iHalf = i / 2;
1144 const int jHalf = j / 2;
1145 matricesCurrentV(i, j, m) += (matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf + 1, mm)) / 2;
1146 // point on corner lines at phi direction
1147 matricesCurrentV(i, j, mp1) += (matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf + 1, mm) + matricesCurrentVC(iHalf, jHalf, mmPlus) + matricesCurrentVC(iHalf, jHalf + 1, mmPlus)) / 4;
1148 }
1149 }
1150
1151 for (int j = 2; j < tnZColumn - 1; j += 2) {
1152 for (int i = 1; i < tnRRow - 1; i += 2) {
1153 const int iHalf = i / 2;
1154 const int jHalf = j / 2;
1155 matricesCurrentV(i, j, m) += (matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf + 1, jHalf, mm)) / 2;
1156 // point on line at phi direction
1157 matricesCurrentV(i, j, mp1) += ((matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf, mmPlus)) + (matricesCurrentVC(iHalf + 1, jHalf, mmPlus) + matricesCurrentVC(iHalf + 1, jHalf, mm))) / 4;
1158 }
1159 }
1160
1161 for (int j = 1; j < tnZColumn - 1; j += 2) {
1162 for (int i = 1; i < tnRRow - 1; i += 2) {
1163 const int iHalf = i / 2;
1164 const int jHalf = j / 2;
1165 matricesCurrentV(i, j, m) += ((matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf + 1, mm)) + (matricesCurrentVC(iHalf + 1, jHalf, mm) + matricesCurrentVC(iHalf + 1, jHalf + 1, mm))) / 4;
1166 // point at the center at phi direction
1167 matricesCurrentV(i, j, mp1) += ((matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf + 1, mm) + matricesCurrentVC(iHalf, jHalf, mmPlus) + matricesCurrentVC(iHalf, jHalf + 1, mmPlus)) +
1168 (matricesCurrentVC(iHalf + 1, jHalf, mm) + matricesCurrentVC(iHalf + 1, jHalf + 1, mm) + matricesCurrentVC(iHalf + 1, jHalf, mmPlus) + matricesCurrentVC(iHalf + 1, jHalf + 1, mmPlus))) /
1169 8;
1170 }
1171 }
1172 }
1173
1174 } else {
1175#pragma omp parallel for num_threads(sNThreads) // no change
1176 for (int m = 0; m < newPhiSlice; m++) {
1177 addInterp2D(matricesCurrentV, matricesCurrentVC, tnRRow, tnZColumn, m);
1178 }
1179 }
1180}
1181
1182template <typename DataT>
1183void PoissonSolver<DataT>::addInterp2D(Vector& matricesCurrentV, const Vector& matricesCurrentVC, const int tnRRow, const int tnZColumn, const int tnPhi) const
1184{
1185 for (int j = 2; j < tnZColumn - 1; j += 2) {
1186 for (int i = 2; i < tnRRow - 1; i += 2) {
1187 matricesCurrentV(i, j, tnPhi) = matricesCurrentV(i, j, tnPhi) + matricesCurrentVC(i / 2, j / 2, tnPhi);
1188 }
1189 }
1190
1191 for (int j = 1; j < tnZColumn - 1; j += 2) {
1192 for (int i = 2; i < tnRRow - 1; i += 2) {
1193 const int iHalf = i / 2;
1194 const int jHalf = j / 2;
1195 matricesCurrentV(i, j, tnPhi) = matricesCurrentV(i, j, tnPhi) + (matricesCurrentVC(iHalf, jHalf, tnPhi) + matricesCurrentVC(iHalf, jHalf + 1, tnPhi)) / 2;
1196 }
1197 }
1198
1199 for (int j = 2; j < tnZColumn - 1; j += 2) {
1200 for (int i = 1; i < tnRRow - 1; i += 2) {
1201 const int iHalf = i / 2;
1202 const int jHalf = j / 2;
1203 matricesCurrentV(i, j, tnPhi) = matricesCurrentV(i, j, tnPhi) + (matricesCurrentVC(iHalf, jHalf, tnPhi) + matricesCurrentVC(iHalf + 1, jHalf, tnPhi)) / 2;
1204 }
1205 }
1206
1207 // only if full
1208 if (MGParameters::gtType == GridTransferType::Full) {
1209 for (int j = 1; j < tnZColumn - 1; j += 2) {
1210 for (int i = 1; i < tnRRow - 1; i += 2) {
1211 const int iHalf = i / 2;
1212 const int jHalf = j / 2;
1213 matricesCurrentV(i, j, tnPhi) = matricesCurrentV(i, j, tnPhi) + (matricesCurrentVC(iHalf, jHalf, tnPhi) + matricesCurrentVC(iHalf, jHalf + 1, tnPhi) + matricesCurrentVC(iHalf + 1, jHalf, tnPhi) + matricesCurrentVC(iHalf + 1, jHalf + 1, tnPhi)) / 4;
1214 }
1215 }
1216 }
1217}
1218
1219template <typename DataT>
1220void PoissonSolver<DataT>::relax3D(Vector& matricesCurrentV, const Vector& matricesCurrentCharge, const int tnRRow, const int tnZColumn, const int iPhi, const int symmetry, const DataT h2,
1221 const DataT tempRatioZ, const std::vector<DataT>& coefficient1, const std::vector<DataT>& coefficient2, const std::vector<DataT>& coefficient3, const std::vector<DataT>& coefficient4) const
1222{
1223 // Gauss-Seidel (Read Black}
1224 if (MGParameters::relaxType == RelaxType::GaussSeidel) {
1225 // for each slice
1226 for (int iPass = 1; iPass <= 2; ++iPass) {
1227 const int msw = (iPass % 2) ? 1 : 2;
1228 for (int m = 0; m < iPhi; ++m) {
1229 const int jsw = ((msw + m) % 2) ? 1 : 2;
1230 int mp1 = m + 1;
1231 int signPlus = 1;
1232 int mm1 = m - 1;
1233 int signMinus = 1;
1234 // Reflection symmetry in phi (e.g. symmetry at sector boundaries, or half sectors, etc.)
1235 if (symmetry == 1) {
1236 if (mp1 > iPhi - 1) {
1237 mp1 = iPhi - 2;
1238 }
1239 if (mm1 < 0) {
1240 mm1 = 1;
1241 }
1242 }
1243 // Anti-symmetry in phi
1244 else if (symmetry == -1) {
1245 if (mp1 > iPhi - 1) {
1246 mp1 = iPhi - 2;
1247 signPlus = -1;
1248 }
1249 if (mm1 < 0) {
1250 mm1 = 1;
1251 signMinus = -1;
1252 }
1253 } else { // No Symmetries in phi, no boundaries, the calculation is continuous across all phi
1254 if (mp1 > iPhi - 1) {
1255 mp1 = m + 1 - iPhi;
1256 }
1257 if (mm1 < 0) {
1258 mm1 = m - 1 + iPhi;
1259 }
1260 }
1261 int isw = jsw;
1262 for (int j = 1; j < tnZColumn - 1; ++j, isw = 3 - isw) {
1263 for (int i = isw; i < tnRRow - 1; i += 2) {
1264 (matricesCurrentV)(i, j, m) = (coefficient2[i] * (matricesCurrentV)(i - 1, j, m) + tempRatioZ * ((matricesCurrentV)(i, j - 1, m) + (matricesCurrentV)(i, j + 1, m)) + coefficient1[i] * (matricesCurrentV)(i + 1, j, m) + coefficient3[i] * (signPlus * (matricesCurrentV)(i, j, mp1) + signMinus * (matricesCurrentV)(i, j, mm1)) + (h2 * (matricesCurrentCharge)(i, j, m))) * coefficient4[i];
1265 } // end cols
1266 } // end mParamGrid.NRVertices
1267 } // end phi
1268 } // end sweep
1269 } else if (MGParameters::relaxType == RelaxType::Jacobi) {
1270 // for each slice
1271 for (int m = 0; m < iPhi; ++m) {
1272 int mp1 = m + 1;
1273 int signPlus = 1;
1274 int mm1 = m - 1;
1275 int signMinus = 1;
1276
1277 // Reflection symmetry in phi (e.g. symmetry at sector boundaries, or half sectors, etc.)
1278 if (symmetry == 1) {
1279 if (mp1 > iPhi - 1) {
1280 mp1 = iPhi - 2;
1281 }
1282 if (mm1 < 0) {
1283 mm1 = 1;
1284 }
1285 }
1286 // Anti-symmetry in phi
1287 else if (symmetry == -1) {
1288 if (mp1 > iPhi - 1) {
1289 mp1 = iPhi - 2;
1290 signPlus = -1;
1291 }
1292 if (mm1 < 0) {
1293 mm1 = 1;
1294 signMinus = -1;
1295 }
1296 } else { // No Symmetries in phi, no boundaries, the calculation is continuous across all phi
1297 if (mp1 > iPhi - 1) {
1298 mp1 = m + 1 - iPhi;
1299 }
1300 if (mm1 < 0) {
1301 mm1 = m - 1 + iPhi;
1302 }
1303 }
1304 // Jacobian
1305 for (int j = 1; j < tnZColumn - 1; ++j) {
1306 for (int i = 1; i < tnRRow - 1; ++i) {
1307 (matricesCurrentV)(i, j, m) = (coefficient2[i] * (matricesCurrentV)(i - 1, j, m) + tempRatioZ * ((matricesCurrentV)(i, j - 1, m) + (matricesCurrentV)(i, j + 1, m)) + coefficient1[i] * (matricesCurrentV)(i + 1, j, m) + coefficient3[i] * (signPlus * (matricesCurrentV)(i, j, mp1) + signMinus * (matricesCurrentV)(i, j, mm1)) + (h2 * (matricesCurrentCharge)(i, j, m))) * coefficient4[i];
1308 } // end cols
1309 } // end mParamGrid.NRVertices
1310 } // end phi
1311 } else {
1312 // Case weighted Jacobi
1313 // TODO
1314 }
1315}
1316
1317template <typename DataT>
1318void PoissonSolver<DataT>::relax2D(Vector& matricesCurrentV, const Vector& matricesCurrentCharge, const int tnRRow, const int tnZColumn, const DataT h2, const DataT tempFourth, const DataT tempRatio,
1319 std::vector<DataT>& coefficient1, std::vector<DataT>& coefficient2)
1320{
1321 // Gauss-Seidel
1322 const int iPhi = 0;
1323 if (MGParameters::relaxType == RelaxType::GaussSeidel) {
1324 int jsw = 1;
1325 for (int iPass = 1; iPass <= 2; ++iPass, jsw = 3 - jsw) {
1326 int isw = jsw;
1327 for (int j = 1; j < tnZColumn - 1; ++j, isw = 3 - isw) {
1328 for (int i = isw; i < tnRRow - 1; i += 2) {
1329 matricesCurrentV(i, j, iPhi) = tempFourth * (coefficient1[i] * matricesCurrentV(i + 1, j, iPhi) + coefficient2[i] * matricesCurrentV(i - 1, j, iPhi) +
1330 tempRatio * (matricesCurrentV(i, j + 1, iPhi) + matricesCurrentV(i, j - 1, iPhi)) + (h2 * matricesCurrentCharge(i, j, iPhi)));
1331 } // end cols
1332 } // end mParamGrid.NRVertices
1333 } // end pass red-black
1334 } else if (MGParameters::relaxType == RelaxType::Jacobi) {
1335 for (int j = 1; j < tnZColumn - 1; ++j) {
1336 for (int i = 1; i < tnRRow - 1; ++i) {
1337 matricesCurrentV(i, j, iPhi) = tempFourth * (coefficient1[i] * matricesCurrentV(i + 1, j, iPhi) + coefficient2[i] * matricesCurrentV(i - 1, j, iPhi) +
1338 tempRatio * (matricesCurrentV(i, j + 1, iPhi) + matricesCurrentV(i, j - 1, iPhi)) + (h2 * matricesCurrentCharge(i, j, iPhi)));
1339 } // end cols
1340 } // end mParamGrid.NRVertices
1341 } else if (MGParameters::relaxType == RelaxType::WeightedJacobi) {
1342 // Weighted Jacobi
1343 // TODO
1344 }
1345}
1346
1347template <typename DataT>
1348void PoissonSolver<DataT>::restrictBoundary3D(Vector& matricesCurrentCharge, const Vector& residue, const int tnRRow, const int tnZColumn, const int newPhiSlice, const int oldPhiSlice) const
1349{
1350 // in case of full 3d and the mParamGrid.NPhiVertices is also coarsening
1351 if (2 * newPhiSlice == oldPhiSlice) {
1352 for (int m = 0, mm = 0; m < newPhiSlice; ++m, mm += 2) {
1353 // for boundary
1354 for (int j = 0, jj = 0; j < tnZColumn; ++j, jj += 2) {
1355 matricesCurrentCharge(0, j, m) = residue(0, jj, mm);
1356 matricesCurrentCharge(tnRRow - 1, j, m) = residue((tnRRow - 1) * 2, jj, mm);
1357 }
1358
1359 // for boundary
1360 for (int i = 0, ii = 0; i < tnRRow; ++i, ii += 2) {
1361 matricesCurrentCharge(i, 0, m) = residue(ii, 0, mm);
1362 matricesCurrentCharge(i, tnZColumn - 1, m) = residue(ii, (tnZColumn - 1) * 2, mm);
1363 }
1364 } // end phis
1365 } else {
1366 for (int m = 0; m < newPhiSlice; ++m) {
1367 restrictBoundary2D(matricesCurrentCharge, residue, tnRRow, tnZColumn, m);
1368 }
1369 }
1370}
1371
1372template <typename DataT>
1373void PoissonSolver<DataT>::restrictBoundary2D(Vector& matricesCurrentCharge, const Vector& residue, const int tnRRow, const int tnZColumn, const int tnPhi) const
1374{
1375 // for boundary
1376 for (int j = 0, jj = 0; j < tnZColumn; ++j, jj += 2) {
1377 matricesCurrentCharge(0, j, tnPhi) = residue(0, jj, tnPhi);
1378 matricesCurrentCharge(tnRRow - 1, j, tnPhi) = residue((tnRRow - 1) * 2, jj, tnPhi);
1379 }
1380
1381 // for boundary
1382 for (int i = 0, ii = 0; i < tnRRow; ++i, ii += 2) {
1383 matricesCurrentCharge(i, 0, tnPhi) = residue(ii, 0, tnPhi);
1384 matricesCurrentCharge(i, tnZColumn - 1, tnPhi) = residue(ii, (tnZColumn - 1) * 2, tnPhi);
1385 }
1386}
1387
1388template <typename DataT>
1389void PoissonSolver<DataT>::restrict3D(Vector& matricesCurrentCharge, const Vector& residue, const int tnRRow, const int tnZColumn, const int newPhiSlice, const int oldPhiSlice) const
1390{
1391 if (2 * newPhiSlice == oldPhiSlice) {
1392 int mm = 0;
1393 for (int m = 0; m < newPhiSlice; m++, mm += 2) {
1394 // assuming no symmetry
1395 int mp1 = mm + 1;
1396 int mm1 = mm - 1;
1397
1398 if (mp1 > (oldPhiSlice)-1) {
1399 mp1 = mm + 1 - (oldPhiSlice);
1400 }
1401 if (mm1 < 0) {
1402 mm1 = mm - 1 + (oldPhiSlice);
1403 }
1404
1405 for (int i = 1, ii = 2; i < tnRRow - 1; ++i, ii += 2) {
1406 for (int j = 1, jj = 2; j < tnZColumn - 1; ++j, jj += 2) {
1407
1408 // at the same plane
1409 const int iip1 = ii + 1;
1410 const int iim1 = ii - 1;
1411 const int jjp1 = jj + 1;
1412 const int jjm1 = jj - 1;
1413 const DataT s1 = residue(iip1, jj, mm) + residue(iim1, jj, mm) + residue(ii, jjp1, mm) + residue(ii, jjm1, mm) + residue(ii, jj, mp1) + residue(ii, jj, mm1);
1414
1415 const DataT s2 = (residue(iip1, jjp1, mm) + residue(iip1, jjm1, mm) + residue(iip1, jj, mp1) + residue(iip1, jj, mm1)) +
1416 (residue(iim1, jjm1, mm) + residue(iim1, jjp1, mm) + residue(iim1, jj, mp1) + residue(iim1, jj, mm1)) +
1417 residue(ii, jjm1, mp1) + residue(ii, jjp1, mm1) + residue(ii, jjm1, mm1) + residue(ii, jjp1, mp1);
1418
1419 const DataT s3 = (residue(iip1, jjp1, mp1) + residue(iip1, jjm1, mp1) + residue(iip1, jjp1, mm1) + residue(iip1, jjm1, mm1)) +
1420 (residue(iim1, jjm1, mm1) + residue(iim1, jjp1, mm1) + residue(iim1, jjm1, mp1) + residue(iim1, jjp1, mp1));
1421
1422 matricesCurrentCharge(i, j, m) = residue(ii, jj, mm) / 8 + s1 / 16 + s2 / 32 + s3 / 64;
1423 } // end cols
1424 } // end mParamGrid.NRVertices
1425
1426 // for boundary
1427 for (int j = 0, jj = 0; j < tnZColumn; ++j, jj += 2) {
1428 matricesCurrentCharge(0, j, m) = residue(0, jj, mm);
1429 matricesCurrentCharge(tnRRow - 1, j, m) = residue((tnRRow - 1) * 2, jj, mm);
1430 }
1431
1432 // for boundary
1433 for (int i = 0, ii = 0; i < tnRRow; ++i, ii += 2) {
1434 matricesCurrentCharge(i, 0, m) = residue(ii, 0, mm);
1435 matricesCurrentCharge(i, tnZColumn - 1, m) = residue(ii, (tnZColumn - 1) * 2, mm);
1436 }
1437 } // end phis
1438
1439 } else {
1440 for (int m = 0; m < newPhiSlice; ++m) {
1441 restrict2D(matricesCurrentCharge, residue, tnRRow, tnZColumn, m);
1442 }
1443 }
1444}
1445
1446template <typename DataT>
1447void PoissonSolver<DataT>::restrict2D(Vector& matricesCurrentCharge, const Vector& residue, const int tnRRow, const int tnZColumn, const int iphi) const
1448{
1449 for (int i = 1, ii = 2; i < tnRRow - 1; ++i, ii += 2) {
1450 for (int j = 1, jj = 2; j < tnZColumn - 1; ++j, jj += 2) {
1451 const int iip1 = ii + 1;
1452 const int iim1 = ii - 1;
1453 const int jjp1 = jj + 1;
1454 const int jjm1 = jj - 1;
1455 if (MGParameters::gtType == GridTransferType::Half) {
1456 // half
1457 matricesCurrentCharge(i, j, iphi) = residue(ii, jj, iphi) / 2 + (residue(iip1, jj, iphi) + residue(iim1, jj, iphi) + residue(ii, jjp1, iphi) + residue(ii, jjm1, iphi)) / 8;
1458 } else if (MGParameters::gtType == GridTransferType::Full) {
1459 matricesCurrentCharge(i, j, iphi) = residue(ii, jj, iphi) / 4 + (residue(iip1, jj, iphi) + residue(iim1, jj, iphi) + residue(ii, jjp1, iphi) + residue(ii, jjm1, iphi)) / 8 +
1460 (residue(iip1, jjp1, iphi) + residue(iim1, jjp1, iphi) + residue(iip1, jjm1, iphi) + residue(iim1, jjm1, iphi)) / 16;
1461 }
1462 } // end cols
1463 } // end mParamGrid.NRVertices
1464 // boundary
1465 // for boundary
1466 for (int j = 0, jj = 0; j < tnZColumn; ++j, jj += 2) {
1467 matricesCurrentCharge(0, j, iphi) = residue(0, jj, iphi);
1468 matricesCurrentCharge(tnRRow - 1, j, iphi) = residue((tnRRow - 1) * 2, jj, iphi);
1469 }
1470 // for boundary
1471 for (int i = 0, ii = 0; i < tnRRow; ++i, ii += 2) {
1472 matricesCurrentCharge(i, 0, iphi) = residue(ii, 0, iphi);
1473 matricesCurrentCharge(i, tnZColumn - 1, iphi) = residue(ii, (tnZColumn - 1) * 2, iphi);
1474 }
1475}
1476
1477template <typename DataT>
1478DataT PoissonSolver<DataT>::getConvergenceError(const Vector& matricesCurrentV, Vector& prevArrayV) const
1479{
1480 std::vector<DataT> errorArr(prevArrayV.getNphi());
1481
1482 // subtract the two matrices
1483 std::transform(prevArrayV.begin(), prevArrayV.end(), matricesCurrentV.begin(), prevArrayV.begin(), std::minus<DataT>());
1484
1485#pragma omp parallel for num_threads(sNThreads) // parallising this loop is possible - but using more than 2 cores makes it slower -
1486 for (unsigned int m = 0; m < prevArrayV.getNphi(); ++m) {
1487 // square each entry in the vector and sum them up
1488 const auto phiStep = prevArrayV.getNr() * prevArrayV.getNz(); // number of points in one phi slice
1489 const auto start = prevArrayV.begin() + m * phiStep;
1490 const auto end = start + phiStep;
1491 errorArr[m] = std::inner_product(start, end, start, DataT(0)); // inner product "Sum (matrix[a]*matrix[a])"
1492 }
1493 // return largest error
1494 return *std::max_element(std::begin(errorArr), std::end(errorArr));
1495}
1496
1497template <typename DataT>
1498void PoissonSolver<DataT>::calcCoefficients(unsigned int from, unsigned int to, const DataT h, const DataT tempRatioZ, const DataT tempRatioPhi, std::vector<DataT>& coefficient1, std::vector<DataT>& coefficient2, std::vector<DataT>& coefficient3, std::vector<DataT>& coefficient4) const
1499{
1500 for (unsigned int i = from; i < to; ++i) {
1501 const DataT radiusInv = 1 / (TPCParameters<DataT>::IFCRADIUS + i * h);
1502 const DataT hRadiusTmp = h * radiusInv / 2;
1503 coefficient1[i] = 1 + hRadiusTmp;
1504 coefficient2[i] = 1 - hRadiusTmp;
1505 coefficient3[i] = tempRatioPhi * radiusInv * radiusInv;
1506 coefficient4[i] = 1 / (2 * (1 + tempRatioZ + coefficient3[i]));
1507 }
1508}
1509
1510template <typename DataT>
1511void PoissonSolver<DataT>::calcCoefficients2D(unsigned int from, unsigned int to, const DataT h, std::vector<DataT>& coefficient1, std::vector<DataT>& coefficient2) const
1512{
1513 for (int i = from; i < to; ++i) {
1514 DataT radiusInvHalf = h / (2 * (TPCParameters<DataT>::IFCRADIUS + i * h));
1515 coefficient1[i] = 1 + radiusInvHalf;
1516 coefficient2[i] = 1 - radiusInvHalf;
1517 }
1518}
1519
1520template <typename DataT>
1522{
1523 return MGParameters::normalizeGridToOneSector ? (INVTWOPI * SECTORSPERSIDE) : INVTWOPI;
1524}
1525
1526template class o2::tpc::PoissonSolver<double>;
1527template class o2::tpc::PoissonSolver<float>;
This class provides a simple method to store values on a large 3-Dim grid with ROOT io functionality.
int16_t time
Definition RawEventData.h:4
int32_t i
This class provides implementation of Poisson equation solver by MultiGrid Method Original version of...
uint32_t j
Definition RawData.h:0
this is a simple 3D-matrix class with the possibility to resize the dimensions
Class for time synchronization of RawReader instances.
void poissonSolver3D(DataContainer &matricesV, const DataContainer &matricesCharge, const int symmetry)
static DataT getConvergenceError()
void poissonSolver2D(DataContainer &matricesV, const DataContainer &matricesCharge)
const GLfloat * m
Definition glcorearb.h:4066
GLint GLsizei count
Definition glcorearb.h:399
GLuint GLuint end
Definition glcorearb.h:469
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat s1
Definition glcorearb.h:5034
GLuint index
Definition glcorearb.h:781
GLuint start
Definition glcorearb.h:469
constexpr float mm
Definition SpecsV2.h:28
Global TPC definitions and constants.
Definition SimTraits.h:167
constexpr unsigned char SECTORSPERSIDE
Definition Defs.h:40
Enum< T >::Iterator begin(Enum< T >)
Definition Defs.h:173
const int iz
Definition TrackUtils.h:69
static int gamma
number of iteration at coarsest level !TODO SET TO REASONABLE VALUE!
static int nMGCycle
number of multi grid cycle (V type)
static CycleType cycleType
cycleType follow CycleType
static int maxLoop
the number of tree-deep of multi grid
static bool normalizeGridToOneSector
the grid in phi direction is squashed from 2 Pi to (2 Pi / SECTORSPERSIDE). This can used to get the ...
static bool isFull3D
< Parameters choice for MultiGrid algorithm
static GridTransferType gtType
gtType grid transfer type follow GridTransferType
static RelaxType relaxType
relaxType follow RelaxType
static int nPre
number of iteration for pre smoothing
static int nPost
number of iteration for post smoothing
o2::InteractionRecord ir(0, 0)