36template <
typename DataT>
39 using timer = std::chrono::high_resolution_clock;
40 auto start = timer::now();
42 poissonMultiGrid3D(matricesV, matricesCharge, symmetry);
44 poissonMultiGrid3D2D(matricesV, matricesCharge, symmetry);
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);
52template <
typename DataT>
55 poissonMultiGrid2D(matricesV, matricesCharge);
58template <
typename DataT>
62 const DataT gridSpacingR = getSpacingR();
63 const DataT gridSpacingZ = getSpacingZ();
64 const DataT ratioZ = gridSpacingR * gridSpacingR / (gridSpacingZ * gridSpacingZ);
69 int nnRow = mParamGrid.NRVertices;
74 int nnCol = mParamGrid.NZVertices;
80 if (!isPowerOfTwo(mParamGrid.NRVertices - 1)) {
81 LOGP(error,
"PoissonMultiGrid2D: PoissonMultiGrid - Error in the number of mParamGrid.NRVertices. Must be 2**M + 1");
84 if (!isPowerOfTwo(mParamGrid.NZVertices - 1)) {
85 LOGP(error,
"PoissonMultiGrid2D: PoissonMultiGrid - Error in the number of mParamGrid.NZVertices. Must be 2**N + 1");
89 const int nLoop = std::max(nGridRow, nGridCol);
91 LOGP(detail,
"{}", fmt::format(
"PoissonMultiGrid2D: nGridRow={}, nGridCol={}, nLoop={}, nMGCycle={}", nGridRow, nGridCol, nLoop,
MGParameters::nMGCycle));
93 unsigned int iOne = 1;
94 unsigned int jOne = 1;
95 int tnRRow = mParamGrid.NRVertices;
96 int tnZColumn = mParamGrid.NZVertices;
99 std::vector<Vector> tvArrayV(nLoop);
100 std::vector<Vector> tvChargeFMG(nLoop);
101 std::vector<Vector> tvCharge(nLoop);
102 std::vector<Vector> tvResidue(nLoop);
106 tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1;
107 tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 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);
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);
126 restrict2D(tvChargeFMG[
index], tvChargeFMG[
count - 2], tnRRow, tnZColumn, 0);
135 LOGP(detail,
"PoissonMultiGrid2D: Do full cycle");
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;
144 const DataT tempRatio = ratioZ * iOne * iOne / (jOne * jOne);
145 const DataT tempFourth = 1 / (2 + 2 * tempRatio);
147 std::vector<DataT> coefficient1(tnRRow);
148 std::vector<DataT> coefficient2(tnRRow);
149 calcCoefficients2D(1, tnRRow - 1,
h, coefficient1, coefficient2);
151 relax2D(tvArrayV[nLoop - 1], tvChargeFMG[nLoop - 1], tnRRow, tnZColumn, h2, tempFourth, tempRatio, coefficient1, coefficient2);
158 tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1;
159 tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1;
161 interp2D(tvArrayV[
count], tvArrayV[
count + 1], tnRRow, tnZColumn, iPhi);
173 LOGP(detail,
"PoissonMultiGrid2D: Do V cycle");
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);
202template <
typename DataT>
205 LOGP(detail,
"{}", fmt::format(
"PoissonMultiGrid3D2D: in Poisson Solver 3D multiGrid semi coarsening mParamGrid.NRVertices={}, cols={}, mParamGrid.NPhiVertices={}", mParamGrid.NZVertices, mParamGrid.NRVertices, mParamGrid.NPhiVertices));
208 if (!isPowerOfTwo((mParamGrid.NRVertices - 1))) {
209 LOGP(error,
"PoissonMultiGrid3D2D: Poisson3DMultiGrid - Error in the number of mParamGrid.NRVertices. Must be 2**M + 1");
212 if (!isPowerOfTwo((mParamGrid.NZVertices - 1))) {
213 LOGP(error,
"PoissonMultiGrid3D2D: Poisson3DMultiGrid - Error in the number of mParamGrid.NZVertices. Must be 2**N + 1");
216 if (mParamGrid.NPhiVertices <= 3) {
217 LOGP(error,
"PoissonMultiGrid3D2D: Poisson3DMultiGrid - Error in the number of mParamGrid.NPhiVertices. Must be larger than 3");
221 const DataT gridSpacingR = getSpacingR();
222 const DataT gridSpacingZ = getSpacingZ();
223 const DataT gridSpacingPhi = getSpacingPhi();
224 const DataT ratioPhi = gridSpacingR * gridSpacingR / (gridSpacingPhi * gridSpacingPhi);
225 const DataT ratioZ = gridSpacingR * gridSpacingR / (gridSpacingZ * gridSpacingZ);
231 int nnRow = mParamGrid.NRVertices;
232 int nnCol = mParamGrid.NZVertices;
234 while (nnRow >>= 1) {
237 while (nnCol >>= 1) {
241 const int maxVal = std::max(nGridRow, nGridCol);
243 unsigned int iOne = 1;
244 unsigned int jOne = 1;
246 std::vector<Vector> tvArrayV(nLoop);
247 std::vector<Vector> tvChargeFMG(nLoop);
248 std::vector<Vector> tvCharge(nLoop);
249 std::vector<Vector> tvPrevArrayV(nLoop);
250 std::vector<Vector> tvResidue(nLoop);
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;
256 tvResidue[
index].resize(tnRRow, tnZColumn, mParamGrid.NPhiVertices);
257 tvPrevArrayV[
index].resize(tnRRow, tnZColumn, mParamGrid.NPhiVertices);
260 tvChargeFMG[
index].resize(tnRRow, tnZColumn, mParamGrid.NPhiVertices);
261 tvArrayV[
index].resize(tnRRow, tnZColumn, mParamGrid.NPhiVertices);
262 tvCharge[
index].resize(tnRRow, tnZColumn, mParamGrid.NPhiVertices);
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);
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);
281 std::vector<DataT> coefficient1(mParamGrid.NRVertices);
282 std::vector<DataT> coefficient2(mParamGrid.NRVertices);
283 std::vector<DataT> coefficient3(mParamGrid.NRVertices);
284 std::vector<DataT> coefficient4(mParamGrid.NRVertices);
285 std::vector<DataT> inverseCoefficient4(mParamGrid.NRVertices);
292 int tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1;
293 int tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1;
295 const DataT h = getSpacingR() * iOne;
297 const DataT iOne2 = iOne * iOne;
298 const DataT tempRatioPhi = ratioPhi * iOne2;
299 const DataT tempRatioZ = ratioZ * iOne2 / (jOne * jOne);
301 calcCoefficients(1, tnRRow - 1,
h, tempRatioZ, tempRatioPhi, coefficient1, coefficient2, coefficient3, coefficient4);
304 relax3D(tvArrayV[nLoop - 1], tvChargeFMG[nLoop - 1], tnRRow, tnZColumn, mParamGrid.NPhiVertices, symmetry, h2, tempRatioZ, coefficient1, coefficient2, coefficient3, coefficient4);
311 tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1;
312 tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1;
315 interp3D(tvArrayV[
count], tvArrayV[
count + 1], tnRRow, tnZColumn, mParamGrid.NPhiVertices, mParamGrid.NPhiVertices);
326 vCycle3D2D(symmetry,
count + 1, nLoop,
MGParameters::nPre,
MGParameters::nPost, ratioZ, ratioPhi, tvArrayV, tvCharge, tvResidue, coefficient1, coefficient2, coefficient3, coefficient4, inverseCoefficient4);
328 const DataT convergenceError = getConvergenceError(tvArrayV[
count], tvPrevArrayV[
count]);
331 if (convergenceError <= sConvergenceError) {
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);
348template <
typename DataT>
351 const DataT gridSpacingR = getSpacingR();
352 const DataT gridSpacingZ = getSpacingZ();
353 const DataT ratioZ = gridSpacingR * gridSpacingR / (gridSpacingZ * gridSpacingZ);
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));
358 if (!isPowerOfTwo((mParamGrid.NRVertices - 1))) {
359 LOGP(error,
"PoissonMultiGrid3D: Poisson3DMultiGrid - Error in the number of mParamGrid.NRVertices. Must be 2**M + 1");
362 if (!isPowerOfTwo((mParamGrid.NZVertices - 1))) {
363 LOGP(error,
"PoissonMultiGrid3D: Poisson3DMultiGrid - Error in the number of mParamGrid.NZVertices. Must be 2**N + 1");
366 if (mParamGrid.NPhiVertices <= 3) {
367 LOGP(error,
"PoissonMultiGrid3D: Poisson3DMultiGrid - Error in the number of mParamGrid.NPhiVertices. Must be larger than 3");
377 int nnRow = mParamGrid.NRVertices;
378 while (nnRow >>= 1) {
382 int nnCol = mParamGrid.NZVertices;
383 while (nnCol >>= 1) {
387 int nnPhi = mParamGrid.NPhiVertices;
388 while (nnPhi % 2 == 0) {
393 LOGP(detail,
"{}", fmt::format(
"PoissonMultiGrid3D: nGridRow={}, nGridCol={}, nGridPhi={}", nGridRow, nGridCol, nGridPhi));
394 const int nLoop = std::max({nGridRow, nGridCol, nGridPhi});
397 unsigned int iOne = 1;
398 unsigned int jOne = 1;
399 unsigned int kOne = 1;
400 int tnRRow = mParamGrid.NRVertices;
401 int tnZColumn = mParamGrid.NZVertices;
402 int tPhiSlice = mParamGrid.NPhiVertices;
405 std::vector<Vector> tvArrayV(nLoop);
406 std::vector<Vector> tvChargeFMG(nLoop);
407 std::vector<Vector> tvCharge(nLoop);
408 std::vector<Vector> tvPrevArrayV(nLoop);
409 std::vector<Vector> tvResidue(nLoop);
411 std::vector<DataT> coefficient1(mParamGrid.NRVertices);
412 std::vector<DataT> coefficient2(mParamGrid.NRVertices);
413 std::vector<DataT> coefficient3(mParamGrid.NRVertices);
414 std::vector<DataT> coefficient4(mParamGrid.NRVertices);
415 std::vector<DataT> inverseCoefficient4(mParamGrid.NRVertices);
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;
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);
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);
455 int otPhiSlice = mParamGrid.NPhiVertices;
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;
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);
468 restrictBoundary3D(tvArrayV[
count - 1], tvArrayV[
count - 2], tnRRow, tnZColumn, tPhiSlice, otPhiSlice);
469 otPhiSlice = tPhiSlice;
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;
489 const DataT h = gridSpacingR * iOne;
491 const DataT gridSizePhiInv = tPhiSlice * getGridSizePhiInv();
492 const DataT tempRatioPhi = h2 * gridSizePhiInv * gridSizePhiInv;
493 const DataT tempRatioZ = ratioZ * iOne * iOne / (jOne * jOne);
495 calcCoefficients(1, tnRRow - 1,
h, tempRatioZ, tempRatioPhi, coefficient1, coefficient2, coefficient3, coefficient4);
498 relax3D(tvArrayV[nLoop - 1], tvChargeFMG[nLoop - 1], tnRRow, tnZColumn, tPhiSlice, symmetry, h2, tempRatioZ, coefficient1, coefficient2, coefficient3, coefficient4);
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);
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;
519 interp3D(tvArrayV[
count], tvArrayV[
count + 1], tnRRow, tnZColumn, tPhiSlice, otPhiSlice);
526 using timer = std::chrono::high_resolution_clock;
527 auto start = timer::now();
532 vCycle3D(symmetry,
count + 1, nLoop,
MGParameters::nPre,
MGParameters::nPost, ratioZ, tvArrayV, tvCharge, tvResidue, coefficient1, coefficient2, coefficient3, coefficient4, inverseCoefficient4);
535 const DataT convergenceError = getConvergenceError(tvArrayV[
count], tvPrevArrayV[
count]);
537 if (convergenceError <= sConvergenceError) {
538 LOGP(detail,
"Cycle converged. Continue to next cycle...");
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);
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);
550 LOGP(warning,
"Cycle {} did not convergence! Current convergence error is larger than expected convergence error: {} > {}", mgCycle, convergenceError, sConvergenceError);
554 otPhiSlice = tPhiSlice;
563 tvPrevArrayV[0] = tvArrayV[0];
566 vCycle3D(symmetry, gridFrom, gridTo,
MGParameters::nPre,
MGParameters::nPost, ratioZ, tvArrayV, tvCharge, tvResidue, coefficient1, coefficient2, coefficient3, coefficient4, inverseCoefficient4);
569 const DataT convergenceError = getConvergenceError(tvArrayV[0], tvPrevArrayV[0]);
572 if (convergenceError <= sConvergenceError) {
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);
588template <
typename DataT>
590 std::vector<Vector>& tvArrayV, std::vector<Vector>& tvCharge, std::vector<Vector>& tvResidue)
592 unsigned int iOne = 1 << (gridFrom - 1);
593 unsigned int jOne = 1 << (gridFrom - 1);
595 int tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1;
596 int tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1;
598 std::vector<DataT> coefficient1(mParamGrid.NRVertices);
599 std::vector<DataT> coefficient2(mParamGrid.NZVertices);
604 const DataT h = gridSizeR * iOne;
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);
616 for (
int jPre = 1; jPre <= nPre; ++jPre) {
617 relax2D(matricesCurrentV, matricesCurrentCharge, tnRRow, tnZColumn, h2, tempFourth, tempRatio, coefficient1, coefficient2);
621 residue2D(residue, matricesCurrentV, matricesCurrentCharge, tnRRow, tnZColumn, ih2, inverseTempFourth, tempRatio, coefficient1, coefficient2);
625 tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1;
626 tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1;
628 matricesCurrentCharge = tvCharge[
count];
629 matricesCurrentV = tvArrayV[
count];
632 restrict2D(matricesCurrentCharge, residue, tnRRow, tnZColumn, 0);
636 for (
int iGamma = 0; iGamma < gamma; ++iGamma) {
637 vCycle2D(gridTo - 1, gridTo, nPre, nPost, gridSizeR, ratio, tvArrayV, tvCharge, tvResidue);
645 const DataT h = gridSizeR * iOne;
647 const DataT tempRatio = ratio * iOne * iOne / (jOne * jOne);
648 const DataT tempFourth = 1 / (2 + 2 * tempRatio);
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];
657 addInterp2D(matricesCurrentV, matricesCurrentVC, tnRRow, tnZColumn, 0);
659 calcCoefficients2D(1, tnRRow - 1,
h, coefficient1, coefficient2);
662 for (Int_t jPost = 1; jPost <= nPost; ++jPost) {
663 relax2D(matricesCurrentV, matricesCurrentCharge, tnRRow, tnZColumn, h2, tempFourth, tempRatio, coefficient1, coefficient2);
668template <
typename DataT>
670 std::vector<Vector>& tvCharge, std::vector<Vector>& tvResidue)
672 unsigned int iOne = 1 << (gridFrom - 1);
673 unsigned int jOne = 1 << (gridFrom - 1);
675 int tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1;
676 int tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1;
678 std::vector<DataT> coefficient1(mParamGrid.NRVertices);
679 std::vector<DataT> coefficient2(mParamGrid.NZVertices);
684 const DataT h = gridSizeR * iOne;
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);
692 for (
int jPre = 1; jPre <= nPre; ++jPre) {
693 relax2D(tvArrayV[
index], tvCharge[
index], tnRRow, tnZColumn, h2, tempFourth, tempRatio, coefficient1, coefficient2);
697 residue2D(tvResidue[
index], tvArrayV[
index], tvCharge[
index], tnRRow, tnZColumn, ih2, inverseTempFourth, tempRatio, coefficient1, coefficient2);
701 tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1;
702 tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1;
705 restrict2D(tvCharge[
count], tvResidue[
index], tnRRow, tnZColumn, 0);
712 const DataT h = gridSizeR * iOne;
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);
718 relax2D(tvArrayV[gridTo - 1], tvCharge[gridTo - 1], tnRRow, tnZColumn, h2, tempFourth, tempRatio, coefficient1, coefficient2);
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);
731 tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1;
732 tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1;
735 addInterp2D(tvArrayV[
index], tvArrayV[
count], tnRRow, tnZColumn, 0);
737 calcCoefficients2D(1, tnRRow - 1, hTmp, coefficient1, coefficient2);
740 for (
int jPost = 1; jPost <= nPost; ++jPost) {
741 relax2D(tvArrayV[
index], tvCharge[
index], tnRRow, tnZColumn, h2Tmp, tempFourthTmp, tempRatioTmp, coefficient1, coefficient2);
746template <
typename DataT>
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
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;
758 const DataT h = getSpacingR() * iOne;
760 const DataT ih2 = 1 / h2;
761 const int iOne2 = iOne * iOne;
762 const DataT tempRatioPhi = ratioPhi * iOne2;
763 const DataT tempRatioZ = ratioZ * iOne2 / (jOne * jOne);
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];
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);
777 residue3D(tvResidue[
index], tvArrayV[
index], tvCharge[
index], tnRRow, tnZColumn, mParamGrid.NPhiVertices, symmetry, ih2, tempRatioZ, coefficient1, coefficient2, coefficient3, inverseCoefficient4);
781 tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1;
782 tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1;
785 restrict3D(tvCharge[
count], tvResidue[
index], tnRRow, tnZColumn, mParamGrid.NPhiVertices, mParamGrid.NPhiVertices);
792 const DataT hTmp = getSpacingR() * iOne;
793 const DataT h2Tmp = hTmp * hTmp;
795 const int iOne2Tmp = iOne * iOne;
796 const DataT tempRatioPhiTmp = ratioPhi * iOne2Tmp;
797 const DataT tempRatioZTmp = ratioZ * iOne2Tmp / (jOne * jOne);
799 calcCoefficients(1, tnRRow - 1, hTmp, tempRatioZTmp, tempRatioPhiTmp, coefficient1, coefficient2, coefficient3, coefficient4);
802 relax3D(tvArrayV[gridTo - 1], tvCharge[gridTo - 1], tnRRow, tnZColumn, mParamGrid.NPhiVertices, symmetry, h2Tmp, tempRatioZTmp, coefficient1, coefficient2, coefficient3, coefficient4);
810 tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1;
811 tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1;
813 const DataT h = getSpacingR() * iOne;
815 const int iOne2 = iOne * iOne;
816 const DataT tempRatioPhi = ratioPhi * iOne2;
817 const DataT tempRatioZ = ratioZ * iOne2 / (jOne * jOne);
820 addInterp3D(tvArrayV[
index], tvArrayV[
count], tnRRow, tnZColumn, mParamGrid.NPhiVertices, mParamGrid.NPhiVertices);
822 calcCoefficients(1, tnRRow - 1,
h, tempRatioZ, tempRatioPhi, coefficient1, coefficient2, coefficient3, coefficient4);
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);
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
836 const DataT gridSpacingR = getSpacingR();
838 unsigned int iOne = 1 << (gridFrom - 1);
839 unsigned int jOne = 1 << (gridFrom - 1);
840 unsigned int kOne = 1 << (gridFrom - 1);
842 int nnPhi = mParamGrid.NPhiVertices;
843 while (nnPhi % 2 == 0) {
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;
854 const int otPhiSlice = tPhiSlice;
855 const DataT h = gridSpacingR * iOne;
857 const DataT ih2 = 1 / h2;
858 const DataT tempGridSizePhiInv = tPhiSlice * getGridSizePhiInv();
859 const DataT tempRatioPhi = h2 * tempGridSizePhiInv * tempGridSizePhiInv;
860 const DataT tempRatioZ = ratioZ * iOne * iOne / (jOne * jOne);
862 calcCoefficients(1, tnRRow - 1,
h, tempRatioZ, tempRatioPhi, coefficient1, coefficient2, coefficient3, coefficient4);
864 for (
int i = 1;
i < tnRRow - 1; ++
i) {
865 inverseCoefficient4[
i] = 1 / coefficient4[
i];
869 for (
int jPre = 1; jPre <= nPre; ++jPre) {
870 relax3D(tvArrayV[
index], tvCharge[
index], tnRRow, tnZColumn, tPhiSlice, symmetry, h2, tempRatioZ, coefficient1, coefficient2, coefficient3, coefficient4);
874 residue3D(tvResidue[
index], tvArrayV[
index], tvCharge[
index], tnRRow, tnZColumn, tPhiSlice, symmetry, ih2, tempRatioZ, coefficient1, coefficient2, coefficient3, inverseCoefficient4);
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;
885 restrict3D(tvCharge[
count], tvResidue[
index], tnRRow, tnZColumn, tPhiSlice, otPhiSlice);
892 const DataT hTmp = gridSpacingR * iOne;
893 const DataT h2Tmp = hTmp * hTmp;
894 const DataT tempGridSizePhiInvTmp = tPhiSlice * getGridSizePhiInv();
895 const DataT tempRatioPhiTmp = h2Tmp * tempGridSizePhiInvTmp * tempGridSizePhiInvTmp;
896 const DataT tempRatioZTmp = ratioZ * iOne * iOne / (jOne * jOne);
898 calcCoefficients(1, tnRRow - 1, hTmp, tempRatioZTmp, tempRatioPhiTmp, coefficient1, coefficient2, coefficient3, coefficient4);
901 relax3D(tvArrayV[gridTo - 1], tvCharge[gridTo - 1], tnRRow, tnZColumn, tPhiSlice, symmetry, h2Tmp, tempRatioZTmp, coefficient1, coefficient2, coefficient3, coefficient4);
905 const int otPhiSlice = tPhiSlice;
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;
915 const DataT h = gridSpacingR * iOne;
917 const DataT tempGridSizePhiInv = tPhiSlice * getGridSizePhiInv();
918 const DataT tempRatioPhi = h2 * tempGridSizePhiInv * tempGridSizePhiInv;
919 const DataT tempRatioZ = ratioZ * iOne * iOne / (jOne * jOne);
922 addInterp3D(tvArrayV[
index], tvArrayV[
count], tnRRow, tnZColumn, tPhiSlice, otPhiSlice);
924 calcCoefficients(1, tnRRow - 1,
h, tempRatioZ, tempRatioPhi, coefficient1, coefficient2, coefficient3, coefficient4);
927 for (
int jPost = 1; jPost <= nPost; ++jPost) {
928 relax3D(tvArrayV[
index], tvCharge[
index], tnRRow, tnZColumn, tPhiSlice, symmetry, h2, tempRatioZ, coefficient1, coefficient2, coefficient3, coefficient4);
933template <
typename DataT>
935 const DataT tempRatio, std::vector<DataT>& coefficient1, std::vector<DataT>& coefficient2)
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);
946 for (
int i = 0;
i < tnRRow; ++
i) {
947 residue(
i, 0, iPhi) = residue(
i, tnZColumn - 1, iPhi) = 0.0;
950 for (
int j = 0;
j < tnZColumn; ++
j) {
951 residue(0,
j, iPhi) = residue(tnRRow - 1,
j, iPhi) = 0.0;
955template <
typename DataT>
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
959#pragma omp parallel for num_threads(sNThreads)
960 for (
int m = 0;
m < tnPhi; ++
m) {
968 if (mp1 > tnPhi - 1) {
976 else if (symmetry == -1) {
977 if (mp1 > tnPhi - 1) {
986 if (mp1 > tnPhi - 1) {
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);
1004template <
typename DataT>
1008 if (newPhiSlice == 2 * oldPhiSlice) {
1009 for (
int m = 0;
m < newPhiSlice;
m += 2) {
1012 int mmPlus =
mm + 1;
1016 if (mmPlus > oldPhiSlice - 1) {
1017 mmPlus =
mm + 1 - oldPhiSlice;
1019 if (mp1 > newPhiSlice - 1) {
1020 mp1 =
m + 1 - newPhiSlice;
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);
1029 matricesCurrentV(
i,
j, mp1) = (matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf, mmPlus)) / 2;
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;
1039 matricesCurrentV(
i,
j, mp1) = (matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf + 1, mm) + matricesCurrentVC(iHalf, jHalf, mmPlus) + matricesCurrentVC(iHalf, jHalf + 1, mmPlus)) / 4;
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;
1049 matricesCurrentV(
i,
j, mp1) = ((matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf, mmPlus)) + (matricesCurrentVC(iHalf + 1, jHalf, mmPlus) + matricesCurrentVC(iHalf + 1, jHalf, mm))) / 4;
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;
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))) /
1067#pragma omp parallel for num_threads(sNThreads)
1068 for (
int m = 0;
m < newPhiSlice; ++
m) {
1069 interp2D(matricesCurrentV, matricesCurrentVC, tnRRow, tnZColumn,
m);
1074template <
typename DataT>
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);
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;
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;
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;
1112template <
typename DataT>
1116 if (newPhiSlice == 2 * oldPhiSlice) {
1117 for (
int m = 0;
m < newPhiSlice;
m += 2) {
1120 int mmPlus =
mm + 1;
1124 if (mmPlus > (oldPhiSlice)-1) {
1125 mmPlus =
mm + 1 - (oldPhiSlice);
1127 if (mp1 > (newPhiSlice)-1) {
1128 mp1 =
m + 1 - (newPhiSlice);
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);
1137 matricesCurrentV(
i,
j, mp1) += (matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf, mmPlus)) / 2;
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;
1147 matricesCurrentV(
i,
j, mp1) += (matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf + 1, mm) + matricesCurrentVC(iHalf, jHalf, mmPlus) + matricesCurrentVC(iHalf, jHalf + 1, mmPlus)) / 4;
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;
1157 matricesCurrentV(
i,
j, mp1) += ((matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf, mmPlus)) + (matricesCurrentVC(iHalf + 1, jHalf, mmPlus) + matricesCurrentVC(iHalf + 1, jHalf, mm))) / 4;
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;
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))) /
1175#pragma omp parallel for num_threads(sNThreads)
1176 for (
int m = 0;
m < newPhiSlice;
m++) {
1177 addInterp2D(matricesCurrentV, matricesCurrentVC, tnRRow, tnZColumn,
m);
1182template <
typename DataT>
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);
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;
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;
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;
1219template <
typename DataT>
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
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;
1235 if (symmetry == 1) {
1236 if (mp1 > iPhi - 1) {
1244 else if (symmetry == -1) {
1245 if (mp1 > iPhi - 1) {
1254 if (mp1 > iPhi - 1) {
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];
1271 for (
int m = 0;
m < iPhi; ++
m) {
1278 if (symmetry == 1) {
1279 if (mp1 > iPhi - 1) {
1287 else if (symmetry == -1) {
1288 if (mp1 > iPhi - 1) {
1297 if (mp1 > iPhi - 1) {
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];
1317template <
typename DataT>
1319 std::vector<DataT>& coefficient1, std::vector<DataT>& coefficient2)
1325 for (
int iPass = 1; iPass <= 2; ++iPass, jsw = 3 - 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)));
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)));
1347template <
typename DataT>
1351 if (2 * newPhiSlice == oldPhiSlice) {
1352 for (
int m = 0, mm = 0;
m < newPhiSlice; ++
m,
mm += 2) {
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);
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);
1366 for (
int m = 0;
m < newPhiSlice; ++
m) {
1367 restrictBoundary2D(matricesCurrentCharge, residue, tnRRow, tnZColumn,
m);
1372template <
typename DataT>
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);
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);
1388template <
typename DataT>
1391 if (2 * newPhiSlice == oldPhiSlice) {
1393 for (
int m = 0;
m < newPhiSlice;
m++,
mm += 2) {
1398 if (mp1 > (oldPhiSlice)-1) {
1399 mp1 =
mm + 1 - (oldPhiSlice);
1402 mm1 =
mm - 1 + (oldPhiSlice);
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) {
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);
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);
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));
1422 matricesCurrentCharge(
i,
j,
m) = residue(ii, jj, mm) / 8 +
s1 / 16 + s2 / 32 + s3 / 64;
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);
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);
1440 for (
int m = 0;
m < newPhiSlice; ++
m) {
1441 restrict2D(matricesCurrentCharge, residue, tnRRow, tnZColumn,
m);
1446template <
typename DataT>
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;
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;
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;
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);
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);
1477template <
typename DataT>
1480 std::vector<DataT> errorArr(prevArrayV.getNphi());
1483 std::transform(prevArrayV.begin(), prevArrayV.end(), matricesCurrentV.begin(), prevArrayV.begin(), std::minus<DataT>());
1485#pragma omp parallel for num_threads(sNThreads)
1486 for (
unsigned int m = 0;
m < prevArrayV.getNphi(); ++
m) {
1488 const auto phiStep = prevArrayV.getNr() * prevArrayV.getNz();
1489 const auto start = prevArrayV.begin() +
m * phiStep;
1494 return *std::max_element(std::begin(errorArr), std::end(errorArr));
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
1500 for (
unsigned int i = from;
i < to; ++
i) {
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]));
1510template <
typename DataT>
1513 for (
int i = from;
i < to; ++
i) {
1515 coefficient1[
i] = 1 + radiusInvHalf;
1516 coefficient2[
i] = 1 - radiusInvHalf;
1520template <
typename DataT>
This class provides a simple method to store values on a large 3-Dim grid with ROOT io functionality.
This class provides implementation of Poisson equation solver by MultiGrid Method Original version of...
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)
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat s1
Global TPC definitions and constants.
constexpr unsigned char SECTORSPERSIDE
Enum< T >::Iterator begin(Enum< T >)
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)