// Copyright CERN and copyright holders of ALICE O2. This software is

// distributed under the terms of the GNU General Public License v3 (GPL

// Version 3), copied verbatim in the file "COPYING".

//

// See http://alice-o2.web.cern.ch/license for full licensing information.

//

// In applying this license CERN does not waive the privileges and immunities

// granted to it by virtue of its status as an Intergovernmental Organization

// or submit itself to any jurisdiction.

/// \file AliTPCPoissonSolver.cxx

/// \brief This class provides implementation of Poisson Eq

/// solver by MultiGrid Method

///

///

///

/// \author Rifki Sadikin <rifki.sadikin@cern.ch>, Indonesian Institute of Sciences

/// \date Nov 20, 2017

#include <TMath.h>

#include "AliTPCPoissonSolver.h"

/// \cond CLASSIMP

ClassImp(AliTPCPoissonSolver);

/// \endcond

const Double_t AliTPCPoissonSolver::fgkTPCZ0 = 249.7; ///< nominal gating grid position

const Double_t AliTPCPoissonSolver::fgkIFCRadius = 83.5; ///< radius which renders the "18 rod manifold" best -> compare calc. of Jim Thomas

const Double_t AliTPCPoissonSolver::fgkOFCRadius = 254.5; ///< Mean Radius of the Outer Field Cage (252.55 min, 256.45 max) (cm)

const Double_t AliTPCPoissonSolver::fgkZOffSet = 0.2; ///< Offset from CE: calculate all distortions closer to CE as if at this point

const Double_t AliTPCPoissonSolver::fgkCathodeV = -100000.0; ///< Cathode Voltage (volts)

const Double_t AliTPCPoissonSolver::fgkGG = -70.0; ///< Gating Grid voltage (volts)

const Double_t AliTPCPoissonSolver::fgkdvdE = 0.0024; ///< [cm/V] drift velocity dependency on the E field (from Magboltz for NeCO2N2 at standard environment)

const Double_t AliTPCPoissonSolver::fgkEM = -1.602176487e-19 / 9.10938215e-31; ///< charge/mass in [C/kg]

const Double_t AliTPCPoissonSolver::fgke0 = 8.854187817e-12; ///< vacuum permittivity [A·s/(V·m)]

Double_t AliTPCPoissonSolver::fgExactErr = 1e-4;

Double_t AliTPCPoissonSolver::fgConvergenceError = 1e-3;

/// constructor

///

AliTPCPoissonSolver::AliTPCPoissonSolver()

: TNamed("poisson solver", "solver"),

fErrorConvergenceNorm2{new TVectorD(fMgParameters.nMGCycle)},

fErrorConvergenceNormInf{new TVectorD(fMgParameters.nMGCycle)},

fError{new TVectorD(fMgParameters.nMGCycle)}

{

// default strategy

}

/// Constructor

/// \param name name of the object

/// \param title title of the object

AliTPCPoissonSolver::AliTPCPoissonSolver(const char* name, const char* title)

: TNamed(name, title),

fErrorConvergenceNorm2{new TVectorD(fMgParameters.nMGCycle)},

fErrorConvergenceNormInf{new TVectorD(fMgParameters.nMGCycle)},

fError{new TVectorD(fMgParameters.nMGCycle)}

{

/// constructor

}

/// destructor

AliTPCPoissonSolver::~AliTPCPoissonSolver()

{

/// virtual destructor

delete[] fExactSolution;

delete fErrorConvergenceNorm2;

delete fErrorConvergenceNormInf;

delete fError;

}

/// Provides poisson solver in 2D

///

/// Based on the strategy (relaxation, multi grid or FFT)

///

/// \param matrixV TMatrixD& potential in matrix

/// \param matrixCharge TMatrixD& charge density in matrix (side effect

/// \param nRRow Int_t number of nRRow in the grid

/// \param nZColumn Int_t number of nZColumn in the grid

/// \param maxIteration Int_t maximum iteration for relaxation method

///

/// \return A fixed number that has nothing to do with what the function does

void AliTPCPoissonSolver::PoissonSolver2D(TMatrixD& matrixV, TMatrixD& matrixCharge, Int_t nRRow, Int_t nZColumn,

Int_t maxIteration)

{

switch (fStrategy) {

case kMultiGrid:

PoissonMultiGrid2D(matrixV, matrixCharge, nRRow, nZColumn);

break;

default:

PoissonRelaxation2D(matrixV, matrixCharge, nRRow, nZColumn, maxIteration);

}

}

/// Provides poisson solver in Cylindrical 3D (TPC geometry)

///

/// Strategy based on parameter settings (fStrategy and fMgParameters)provided

/// * Cascaded multi grid with S.O.R

/// * Geometric MultiGrid

/// * Cycles: V, W, Full

/// * Relaxation: Jacobi, Weighted-Jacobi, Gauss-Seidel

/// * Grid transfer operators: Full, Half

/// * Spectral Methods (TODO)

///

/// \param matricesV TMatrixD** potential in 3D matrix

/// \param matricesCharge TMatrixD** charge density in 3D matrix (side effect)

/// \param nRRow Int_t number of nRRow in the r direction of TPC

/// \param nZColumn Int_t number of nZColumn in z direction of TPC

/// \param phiSlice Int_t number of phiSlice in phi direction of T{C

/// \param maxIteration Int_t maximum iteration for relaxation method

/// \param symmetry Int_t symmetry or not

///

/// \pre Charge density distribution in **matricesCharge** is known and boundary values for **matricesV** are set

/// \post Numerical solution for potential distribution is calculated and stored in each rod at **matricesV**

void AliTPCPoissonSolver::PoissonSolver3D(TMatrixD** matricesV, TMatrixD** matricesCharge,

Int_t nRRow, Int_t nZColumn, Int_t phiSlice, Int_t maxIteration,

Int_t symmetry)

{

switch (fStrategy) {

case kMultiGrid:

if (fMgParameters.isFull3D) {

PoissonMultiGrid3D(matricesV, matricesCharge, nRRow, nZColumn, phiSlice, symmetry);

} else {

PoissonMultiGrid3D2D(matricesV, matricesCharge, nRRow, nZColumn, phiSlice, symmetry);

}

break;

default:

PoissonRelaxation3D(matricesV, matricesCharge, nRRow, nZColumn, phiSlice, maxIteration, symmetry);

}

}

/// Solve Poisson's Equation by Relaxation Technique in 2D (assuming cylindrical symmetry)

///

/// Solve Poisson's equation in a cylindrical coordinate system. The matrixV matrix must be filled with the

/// boundary conditions on the first and last nRRow, and the first and last nZColumn. The remainder of the

/// array can be blank or contain a preliminary guess at the solution. The Charge density matrix contains

/// the enclosed spacecharge density at each point. The charge density matrix can be full of zero's if

/// you wish to solve Laplace equation however it should not contain random numbers or you will get

/// random numbers back as a solution.

/// Poisson's equation is solved by iteratively relaxing the matrix to the final solution. In order to

/// speed up the convergence to the best solution, this algorithm does a binary expansion of the solution

/// space. First it solves the problem on a very sparse grid by skipping nRRow and nZColumn in the original

/// matrix. Then it doubles the number of points and solves the problem again. Then it doubles the

/// number of points and solves the problem again. This happens several times until the maximum number

/// of points has been included in the array.

///

/// NOTE: In order for this algorithm to work, the number of nRRow and nZColumn must be a power of 2 plus one.

/// So nRRow == 2**M + 1 and nZColumn == 2**N + 1. The number of nRRow and nZColumn can be different.

///

/// Method for relaxation: S.O.R Weighted Jacobi

///

/// \param matrixV TMatrixD& potential in matrix

/// \param matrixCharge TMatrixD& charge density in matrix (side effect

/// \param nRRow Int_t number of nRRow in the grid

/// \param nZColumn Int_t number of nZColumn in the grid

/// \param maxIteration Int_t maximum iteration for relaxation method

///

/// \return A fixed number that has nothing to do with what the function does

///

///

/// Original code by Jim Thomas (STAR TPC Collaboration)

void AliTPCPoissonSolver::PoissonRelaxation2D(TMatrixD& matrixV, TMatrixD& matrixCharge, Int_t nRRow, Int_t nZColumn,

Int_t maxIteration)

{

const Float_t gridSizeR = (AliTPCPoissonSolver::fgkOFCRadius - AliTPCPoissonSolver::fgkIFCRadius) / (nRRow - 1);

const Float_t gridSizeZ = AliTPCPoissonSolver::fgkTPCZ0 / (nZColumn - 1);

const Float_t ratio = gridSizeR * gridSizeR / (gridSizeZ * gridSizeZ);

TMatrixD arrayEr(nRRow, nZColumn);

TMatrixD arrayEz(nRRow, nZColumn);

//Check that number of nRRow and nZColumn is suitable for a binary expansion

if (!IsPowerOfTwo(nRRow - 1)) {

Error("PoissonRelaxation2D", "PoissonRelaxation - Error in the number of nRRow. Must be 2**M - 1");

return;

}

if (!IsPowerOfTwo(nZColumn - 1)) {

Error("PoissonRelaxation2D", "PoissonRelaxation - Error in the number of nZColumn. Must be 2**N - 1");

return;

}

// Solve Poisson's equation in cylindrical coordinates by relaxation technique

// Allow for different size grid spacing in R and Z directions

// Use a binary expansion of the size of the matrix to speed up the solution of the problem

Int_t iOne = (nRRow - 1) / 4;

Int_t jOne = (nZColumn - 1) / 4;

// Coarse until nLoop

Int_t nLoop = 1 + (int)(0.5 + TMath::Log2((double)TMath::Max(iOne, jOne)));

// Loop while the matrix expands & the resolution increases.

for (Int_t count = 0; count < nLoop; count++) {

Float_t tempGridSizeR = gridSizeR * iOne;

Float_t tempRatio = ratio * iOne * iOne / (jOne * jOne);

Float_t tempFourth = 1.0 / (2.0 + 2.0 * tempRatio);

// Do this the standard C++ way to avoid gcc extensions for Float_t coefficient1[nRRow]

std::vector<float> coefficient1(nRRow);

std::vector<float> coefficient2(nRRow);

for (Int_t i = iOne; i < nRRow - 1; i += iOne) {

Float_t radius = AliTPCPoissonSolver::fgkIFCRadius + i * gridSizeR;

coefficient1[i] = 1.0 + tempGridSizeR / (2 * radius);

coefficient2[i] = 1.0 - tempGridSizeR / (2 * radius);

}

TMatrixD sumChargeDensity(nRRow, nZColumn);

// average charge at the coarse point

for (Int_t i = iOne; i < nRRow - 1; i += iOne) {

Float_t radius = AliTPCPoissonSolver::fgkIFCRadius + iOne * gridSizeR;

for (Int_t j = jOne; j < nZColumn - 1; j += jOne) {

if (iOne == 1 && jOne == 1) {

sumChargeDensity(i, j) = matrixCharge(i, j);

} else {

// Add up all enclosed charge density contributions within 1/2 unit in all directions

Float_t weight = 0.0;

Float_t sum = 0.0;

sumChargeDensity(i, j) = 0.0;

for (Int_t ii = i - iOne / 2; ii <= i + iOne / 2; ii++) {

for (Int_t jj = j - jOne / 2; jj <= j + jOne / 2; jj++) {

if (ii == i - iOne / 2 || ii == i + iOne / 2 || jj == j - jOne / 2 || jj == j + jOne / 2) {

weight = 0.5;

} else {

weight = 1.0;

}

sumChargeDensity(i, j) += matrixCharge(ii, jj) * weight * radius;

sum += weight * radius;

}

}

sumChargeDensity(i, j) /= sum;

}

sumChargeDensity(i, j) *= tempGridSizeR * tempGridSizeR; // just saving a step later on

}

}

// Iterate on the current level

for (Int_t k = 1; k <= maxIteration; k++) {

// Solve Poisson's Equation

// Over-relaxation index, must be >= 1 but < 2. Arrange for it to evolve from 2 => 1

// as iteration increase.

Float_t overRelax = 1.0 + TMath::Sqrt(TMath::Cos((k * TMath::PiOver2()) / maxIteration));

Float_t overRelaxM1 = overRelax - 1.0;

Float_t overRelaxTemp4, overRelaxCoefficient5;

overRelaxTemp4 = overRelax * tempFourth;

overRelaxCoefficient5 = overRelaxM1 / overRelaxTemp4;

for (Int_t i = iOne; i < nRRow - 1; i += iOne) {

for (Int_t j = jOne; j < nZColumn - 1; j += jOne) {

// S.O.R

//

matrixV(i, j) = (coefficient2[i] * matrixV(i - iOne, j) + tempRatio * (matrixV(i, j - jOne) + matrixV(i, j + jOne)) - overRelaxCoefficient5 * matrixV(i, j) + coefficient1[i] * matrixV(i + iOne, j) + sumChargeDensity(i, j)) * overRelaxTemp4;

}

}

// if already at maxIteration

// TODO: stop when it converged

if (k == maxIteration) {

// After full solution is achieved, copy low resolution solution into higher res array

// Interpolate solution

for (Int_t i = iOne; i < nRRow - 1; i += iOne) {

for (Int_t j = jOne; j < nZColumn - 1; j += jOne) {

if (iOne > 1) {

matrixV(i + iOne / 2, j) = (matrixV(i + iOne, j) + matrixV(i, j)) / 2;

if (i == iOne) {

matrixV(i - iOne / 2, j) = (matrixV(0, j) + matrixV(iOne, j)) / 2;

}

}

if (jOne > 1) {

matrixV(i, j + jOne / 2) = (matrixV(i, j + jOne) + matrixV(i, j)) / 2;

if (j == jOne) {

matrixV(i, j - jOne / 2) = (matrixV(i, 0) + matrixV(i, jOne)) / 2;

}

}

if (iOne > 1 && jOne > 1) {

matrixV(i + iOne / 2, j + jOne / 2) = (matrixV(i + iOne, j + jOne) + matrixV(i, j)) / 2;

if (i == iOne) {

matrixV(i - iOne / 2, j - jOne / 2) = (matrixV(0, j - jOne) + matrixV(iOne, j)) / 2;

}

if (j == jOne) {

matrixV(i - iOne / 2, j - jOne / 2) = (matrixV(i - iOne, 0) + matrixV(i, jOne)) / 2;

}

// Note that this leaves a point at the upper left and lower right corners uninitialized.

// -> Not a big deal.

}

}

}

}

}

iOne = iOne / 2;

if (iOne < 1) {

iOne = 1;

}

jOne = jOne / 2;

if (jOne < 1) {

jOne = 1;

}

sumChargeDensity.Clear();

}

}

/// Solve Poisson's Equation by MultiGrid Technique in 2D (assuming cylindrical symmetry)

///

/// NOTE: In order for this algorithm to work, the number of nRRow and nZColumn must be a power of 2 plus one.

/// So nRRow == 2**M + 1 and nZColumn == 2**N + 1. The number of nRRow and nZColumn can be different.

///

/// \param matrixV TMatrixD& potential in matrix

/// \param matrixCharge TMatrixD& charge density in matrix (side effect

/// \param nRRow Int_t number of nRRow

/// \param nZColumn Int_t number of nZColumn

/// \param maxIteration Int_t maximum iteration for relaxation method

///

/// \return A fixed number that has nothing to do with what the function does

void AliTPCPoissonSolver::PoissonMultiGrid2D(TMatrixD& matrixV, TMatrixD& matrixCharge, Int_t nRRow, Int_t nZColumn)

{

/// Geometry of TPC -- should be use AliTPCParams instead

const Float_t gridSizeR = (AliTPCPoissonSolver::fgkOFCRadius - AliTPCPoissonSolver::fgkIFCRadius) / (nRRow - 1);

const Float_t gridSizeZ = AliTPCPoissonSolver::fgkTPCZ0 / (nZColumn - 1);

const Float_t ratio = gridSizeR * gridSizeR / (gridSizeZ * gridSizeZ);

Int_t nGridRow = 0; // number grid

Int_t nGridCol = 0; // number grid

Int_t nnRow;

Int_t nnCol;

nnRow = nRRow;

while (nnRow >>= 1) {

nGridRow++;

}

nnCol = nZColumn;

while (nnCol >>= 1) {

nGridCol++;

}

//Check that number of nRRow and nZColumn is suitable for multi grid

if (!IsPowerOfTwo(nRRow - 1)) {

Error("PoissonMultiGrid2D", "PoissonMultiGrid - Error in the number of nRRow. Must be 2**M - 1");

return;

}

if (!IsPowerOfTwo(nZColumn - 1)) {

Error("PoissonMultiGrid2D", "PoissonMultiGrid - Error in the number of nZColumn. Must be 2**N - 1");

return;

}

Int_t nLoop = TMath::Max(nGridRow, nGridCol); // Calculate the number of nLoop for the binary expansion

Info("PoissonMultiGrid2D", "%s", Form("nGridRow=%d, nGridCol=%d, nLoop=%d, nMGCycle=%d", nGridRow, nGridCol, nLoop, fMgParameters.nMGCycle));

Float_t h, h2, radius;

Int_t iOne = 1; // in/dex

Int_t jOne = 1; // index

Int_t tnRRow = nRRow, tnZColumn = nZColumn;

Int_t count;

Float_t tempRatio, tempFourth;

// Vector for storing multi grid array

std::vector<TMatrixD*> tvChargeFMG(nLoop);

std::vector<TMatrixD*> tvArrayV(nLoop);

std::vector<TMatrixD*> tvCharge(nLoop);

std::vector<TMatrixD*> tvResidue(nLoop);

// Allocate memory for temporary grid

for (count = 1; count <= nLoop; count++) {

tnRRow = iOne == 1 ? nRRow : nRRow / iOne + 1;

tnZColumn = jOne == 1 ? nZColumn : nZColumn / jOne + 1;

// if one just address to matrixV

tvResidue[count - 1] = new TMatrixD(tnRRow, tnZColumn);

if (count == 1) {

tvChargeFMG[count - 1] = &matrixCharge;

tvArrayV[count - 1] = &matrixV;

tvCharge[count - 1] = &matrixCharge;

} else {

tvArrayV[count - 1] = new TMatrixD(tnRRow, tnZColumn);

tvCharge[count - 1] = new TMatrixD(tnRRow, tnZColumn);

tvChargeFMG[count - 1] = new TMatrixD(tnRRow, tnZColumn);

Restrict2D(*tvChargeFMG[count - 1], *tvChargeFMG[count - 2], tnRRow, tnZColumn);

}

iOne = 2 * iOne;

jOne = 2 * jOne;

}

/// full multi grid

if (fMgParameters.cycleType == kFCycle) {

Info("PoissonMultiGrid2D", "Do full cycle");

// FMG

// 1) Relax on the coarsest grid

iOne = iOne / 2;

jOne = jOne / 2;

tnRRow = iOne == 1 ? nRRow : nRRow / iOne + 1;

tnZColumn = jOne == 1 ? nZColumn : nZColumn / jOne + 1;

h = gridSizeR * count;

h2 = h * h;

tempRatio = ratio * iOne * iOne / (jOne * jOne);

tempFourth = 1.0 / (2.0 + 2.0 * tempRatio);

std::vector<float> coefficient1(tnRRow);

std::vector<float> coefficient2(tnRRow);

for (Int_t i = 1; i < tnRRow - 1; i++) {

radius = AliTPCPoissonSolver::fgkIFCRadius + i * h;

coefficient1[i] = 1.0 + h / (2 * radius);

coefficient2[i] = 1.0 - h / (2 * radius);

}

Relax2D(*tvArrayV[nLoop - 1], *tvChargeFMG[nLoop - 1], tnRRow, tnZColumn, h2, tempFourth, tempRatio, coefficient1,

coefficient2);

// Do VCycle from nLoop H to h

for (count = nLoop - 2; count >= 0; count--) {

iOne = iOne / 2;

jOne = jOne / 2;

tnRRow = iOne == 1 ? nRRow : nRRow / iOne + 1;

tnZColumn = jOne == 1 ? nZColumn : nZColumn / jOne + 1;

Interp2D(*tvArrayV[count], *tvArrayV[count + 1], tnRRow, tnZColumn);

// Copy the relax charge to the tvCharge

*tvCharge[count] = *tvChargeFMG[count]; //copy

//tvCharge[count]->Print();

// Do V cycle

for (Int_t mgCycle = 0; mgCycle < fMgParameters.nMGCycle; mgCycle++) {

VCycle2D(nRRow, nZColumn, count + 1, nLoop, fMgParameters.nPre, fMgParameters.nPost, gridSizeR, ratio, tvArrayV,

tvCharge, tvResidue);

}

}

} else if (fMgParameters.cycleType == kVCycle) {

// 2. VCycle

Info("PoissonMultiGrid2D", "Do V cycle");

Int_t gridFrom = 1;

Int_t gridTo = nLoop;

// Do MGCycle

for (Int_t mgCycle = 0; mgCycle < fMgParameters.nMGCycle; mgCycle++) {

VCycle2D(nRRow, nZColumn, gridFrom, gridTo, fMgParameters.nPre, fMgParameters.nPost, gridSizeR, ratio, tvArrayV,

tvCharge, tvResidue);

}

} else if (fMgParameters.cycleType == kWCycle) {

// 3. W Cycle (TODO:)

Int_t gridFrom = 1;

//nLoop = nLoop >= 4 ? 4 : nLoop;

Int_t gridTo = nLoop;

//Int_t gamma = 1;

// Do MGCycle

for (Int_t mgCycle = 0; mgCycle < fMgParameters.nMGCycle; mgCycle++) {

WCycle2D(nRRow, nZColumn, gridFrom, gridTo, fMgParameters.gamma, fMgParameters.nPre, fMgParameters.nPost,

gridSizeR, ratio, tvArrayV, tvCharge, tvResidue);

}

}

// Deallocate memory

for (count = nLoop; count >= 1; count--) {

// if one just address to matrixV

if (count > 1) {

delete tvArrayV[count - 1];

delete tvCharge[count - 1];

delete tvChargeFMG[count - 1];

}

delete tvResidue[count - 1];

}

}

/// 3D - Solve Poisson's Equation in 3D by Relaxation Technique

///

/// NOTE: In order for this algorithm to work, the number of nRRow and nZColumn must be a power of 2 plus one.

/// The number of nRRow and Z Column can be different.

///

/// R Row == 2**M + 1

/// Z Column == 2**N + 1

/// Phi Slice == Arbitrary but greater than 3

///

/// DeltaPhi in Radians

///

/// SYMMETRY = 0 if no phi symmetries, and no phi boundary conditions

/// = 1 if we have reflection symmetry at the boundaries (eg. sector symmetry or half sector symmetries).

///

/// \param matricesV TMatrixD** potential in 3D matrix

/// \param matricesCharge TMatrixD** charge density in 3D matrix (side effect)

/// \param nRRow Int_t number of nRRow in the r direction of TPC

/// \param nZColumn Int_t number of nZColumn in z direction of TPC

/// \param phiSlice Int_t number of phiSlice in phi direction of T{C

/// \param maxIteration Int_t maximum iteration for relaxation method

/// \param symmetry Int_t symmetry or not

///

void AliTPCPoissonSolver::PoissonRelaxation3D(TMatrixD** matricesV, TMatrixD** matricesCharge,

Int_t nRRow, Int_t nZColumn, Int_t phiSlice, Int_t maxIteration,

Int_t symmetry)

{

const Float_t gridSizeR = (AliTPCPoissonSolver::fgkOFCRadius - AliTPCPoissonSolver::fgkIFCRadius) / (nRRow - 1);

const Float_t gridSizePhi = TMath::TwoPi() / phiSlice;

const Float_t gridSizeZ = AliTPCPoissonSolver::fgkTPCZ0 / (nZColumn - 1);

const Float_t ratioPhi = gridSizeR * gridSizeR / (gridSizePhi * gridSizePhi);

const Float_t ratioZ = gridSizeR * gridSizeR / (gridSizeZ * gridSizeZ);

Info("PoissonRelaxation3D", "%s", Form("in Poisson Solver 3D relaxation nRRow=%d, cols=%d, phiSlice=%d \n", nRRow, nZColumn, phiSlice));

// Check that the number of nRRow and nZColumn is suitable for a binary expansion

if (!IsPowerOfTwo((nRRow - 1))) {

Error("PoissonRelaxation3D", "Poisson3DRelaxation - Error in the number of nRRow. Must be 2**M - 1");

return;

}

if (!IsPowerOfTwo((nZColumn - 1))) {

Error("PoissonRelaxation3D", "Poisson3DRelaxation - Error in the number of nZColumn. Must be 2**N - 1");

return;

}

if (phiSlice <= 3) {

Error("PoissonRelaxation3D", "Poisson3DRelaxation - Error in the number of phiSlice. Must be larger than 3");

return;

}

if (phiSlice > 1000) {

Error("PoissonRelaxation3D", "Poisson3D phiSlice > 1000 is not allowed (nor wise) ");

return;

}

// Solve Poisson's equation in cylindrical coordinates by relaxation technique

// Allow for different size grid spacing in R and Z directions

// Use a binary expansion of the matrix to speed up the solution of the problem

Int_t nLoop, mPlus, mMinus, signPlus, signMinus;

Int_t iOne = (nRRow - 1) / 4;

Int_t jOne = (nZColumn - 1) / 4;

nLoop = TMath::Max(iOne, jOne); // Calculate the number of nLoop for the binary expansion

nLoop = 1 + (int)(0.5 + TMath::Log2((double)nLoop)); // Solve for N in 2**N

TMatrixD* matricesSumChargeDensity[1000]; // Create temporary arrays to store low resolution charge arrays

std::vector<float> coefficient1(

nRRow); // Do this the standard C++ way to avoid gcc extensions for Float_t coefficient1[nRRow]

std::vector<float> coefficient2(

nRRow); // Do this the standard C++ way to avoid gcc extensions for Float_t coefficient1[nRRow]

std::vector<float> coefficient3(

nRRow); // Do this the standard C++ way to avoid gcc extensions for Float_t coefficient1[nRRow]

std::vector<float> coefficient4(

nRRow); // Do this the standard C++ way to avoid gcc extensions for Float_t coefficient1[nRRow]

std::vector<float> overRelaxCoefficient4(nRRow); // Do this the standard C++ way to avoid gcc extensions

std::vector<float> overRelaxCoefficient5(nRRow); // Do this the standard C++ way to avoid gcc extensions

for (Int_t i = 0; i < phiSlice; i++) {

matricesSumChargeDensity[i] = new TMatrixD(nRRow, nZColumn);

}

///// Test of Convergence

TMatrixD* prevArrayV[phiSlice];

for (Int_t m = 0; m < phiSlice; m++) {

prevArrayV[m] = new TMatrixD(nRRow, nZColumn);

}

/////

// START the master loop and do the binary expansion

for (Int_t count = 0; count < nLoop; count++) {

Float_t tempGridSizeR = gridSizeR * iOne;

Float_t tempRatioPhi = ratioPhi * iOne * iOne;

Float_t tempRatioZ = ratioZ * iOne * iOne / (jOne * jOne);

for (Int_t i = iOne; i < nRRow - 1; i += iOne) {

Float_t radius = AliTPCPoissonSolver::fgkIFCRadius + i * gridSizeR;

coefficient1[i] = 1.0 + tempGridSizeR / (2 * radius);

coefficient2[i] = 1.0 - tempGridSizeR / (2 * radius);

coefficient3[i] = tempRatioPhi / (radius * radius);

coefficient4[i] = 0.5 / (1.0 + tempRatioZ + coefficient3[i]);

}

for (Int_t m = 0; m < phiSlice; m++) {

TMatrixD& matrixCharge = *matricesCharge[m];

TMatrixD& sumChargeDensity = *matricesSumChargeDensity[m];

for (Int_t i = iOne; i < nRRow - 1; i += iOne) {

Float_t radius = AliTPCPoissonSolver::fgkIFCRadius + i * gridSizeR;

for (Int_t j = jOne; j < nZColumn - 1; j += jOne) {

if (iOne == 1 && jOne == 1) {

sumChargeDensity(i, j) = matrixCharge(i, j);

} else { // Add up all enclosed charge density contributions within 1/2 unit in all directions

Float_t weight = 0.0;

Float_t sum = 0.0;

sumChargeDensity(i, j) = 0.0;

for (Int_t ii = i - iOne / 2; ii <= i + iOne / 2; ii++) {

for (Int_t jj = j - jOne / 2; jj <= j + jOne / 2; jj++) {

if (ii == i - iOne / 2 || ii == i + iOne / 2 || jj == j - jOne / 2 || jj == j + jOne / 2) {

weight = 0.5;

} else {

weight = 1.0;

}

sumChargeDensity(i, j) += matrixCharge(ii, jj) * weight * radius;

sum += weight * radius;

}

}

sumChargeDensity(i, j) /= sum;

}

sumChargeDensity(i, j) *= tempGridSizeR * tempGridSizeR; // just saving a step later on

}

}

}

for (Int_t k = 1; k <= maxIteration; k++) {

if (count == nLoop - 1) {

//// Test of Convergence

for (Int_t m = 0; m < phiSlice; m++) {

(*prevArrayV[m]) = (*matricesV[m]);

}

////

}

Float_t overRelax = 1.0 + TMath::Sqrt(TMath::Cos((k * TMath::PiOver2()) / maxIteration));

Float_t overRelaxM1 = overRelax - 1.0;

for (Int_t i = iOne; i < nRRow - 1; i += iOne) {

overRelaxCoefficient4[i] = overRelax * coefficient4[i];

overRelaxCoefficient5[i] = overRelaxM1 / overRelaxCoefficient4[i];

}

for (Int_t m = 0; m < phiSlice; m++) {

mPlus = m + 1;

signPlus = 1;

mMinus = m - 1;

signMinus = 1;

if (symmetry == 1) { // Reflection symmetry in phi (e.g. symmetry at sector boundaries, or half sectors, etc.)

if (mPlus > phiSlice - 1) {

mPlus = phiSlice - 2;

}

if (mMinus < 0) {

mMinus = 1;

}

} else if (symmetry == -1) { // Anti-symmetry in phi

if (mPlus > phiSlice - 1) {

mPlus = phiSlice - 2;

signPlus = -1;

}

if (mMinus < 0) {

mMinus = 1;

signMinus = -1;

}

} else { // No Symmetries in phi, no boundaries, the calculation is continuous across all phi

if (mPlus > phiSlice - 1) {

mPlus = m + 1 - phiSlice;

}

if (mMinus < 0) {

mMinus = m - 1 + phiSlice;

}

}

TMatrixD& matrixV = *matricesV[m];

TMatrixD& matrixVP = *matricesV[mPlus];

TMatrixD& matrixVM = *matricesV[mMinus];

TMatrixD& sumChargeDensity = *matricesSumChargeDensity[m];

Double_t* matrixVFast = matrixV.GetMatrixArray();

Double_t* matrixVPFast = matrixVP.GetMatrixArray();

Double_t* matrixVMFast = matrixVM.GetMatrixArray();

Double_t* sumChargeDensityFast = sumChargeDensity.GetMatrixArray();

if (fStrategy == kRelaxation) {

// slow implementation

for (Int_t i = iOne; i < nRRow - 1; i += iOne) {

for (Int_t j = jOne; j < nZColumn - 1; j += jOne) {

matrixV(i, j) = (coefficient2[i] * matrixV(i - iOne, j) + tempRatioZ * (matrixV(i, j - jOne) + matrixV(i, j + jOne)) - overRelaxCoefficient5[i] * matrixV(i, j) + coefficient1[i] * matrixV(i + iOne, j) + coefficient3[i] * (signPlus * matrixVP(i, j) + signMinus * matrixVM(i, j)) + sumChargeDensity(i, j)) * overRelaxCoefficient4[i];

// Note: over-relax the solution at each step. This speeds up the convergence.

}

}

} else {

for (Int_t i = iOne; i < nRRow - 1; i += iOne) {

Double_t* matrixVFastI = &(matrixVFast[i * nZColumn]);

Double_t* matrixVPFastI = &(matrixVPFast[i * nZColumn]);

Double_t* matrixVMFastI = &(matrixVMFast[i * nZColumn]);

Double_t* sumChargeDensityFastI = &(sumChargeDensityFast[i * nZColumn]);

for (Int_t j = jOne; j < nZColumn - 1; j += jOne) {

Double_t /*resSlow*/ resFast;

resFast = (coefficient2[i] * matrixVFastI[j - nZColumn * iOne] + tempRatioZ * (matrixVFastI[j - jOne] + matrixVFastI[j + jOne]) - overRelaxCoefficient5[i] * matrixVFastI[j] + coefficient1[i] * matrixVFastI[j + nZColumn * iOne] + coefficient3[i] * (signPlus * matrixVPFastI[j] + signMinus * matrixVMFastI[j]) + sumChargeDensityFastI[j]) * overRelaxCoefficient4[i];

matrixVFastI[j] = resFast;

// Note: over-relax the solution at each step. This speeds up the convergence.

} // end j

} //end i

} // end phi

// After full solution is achieved, copy low resolution solution into higher res array

if (k == maxIteration) {

for (Int_t i = iOne; i < nRRow - 1; i += iOne) {

for (Int_t j = jOne; j < nZColumn - 1; j += jOne) {

if (iOne > 1) {

matrixV(i + iOne / 2, j) = (matrixV(i + iOne, j) + matrixV(i, j)) / 2;

if (i == iOne) {

matrixV(i - iOne / 2, j) = (matrixV(0, j) + matrixV(iOne, j)) / 2;

}

}

if (jOne > 1) {

matrixV(i, j + jOne / 2) = (matrixV(i, j + jOne) + matrixV(i, j)) / 2;

if (j == jOne) {

matrixV(i, j - jOne / 2) = (matrixV(i, 0) + matrixV(i, jOne)) / 2;

}

}

if (iOne > 1 && jOne > 1) {

matrixV(i + iOne / 2, j + jOne / 2) = (matrixV(i + iOne, j + jOne) + matrixV(i, j)) / 2;

if (i == iOne) {

matrixV(i - iOne / 2, j - jOne / 2) = (matrixV(0, j - jOne) + matrixV(iOne, j)) / 2;

}

if (j == jOne) {

matrixV(i - iOne / 2, j - jOne / 2) = (matrixV(i - iOne, 0) + matrixV(i, jOne)) / 2;

}

// Note that this leaves a point at the upper left and lower right corners uninitialized. Not a big deal.

}

}

}

}

}

if (count == nLoop - 1) {

(*fErrorConvergenceNormInf)(k - 1) = GetConvergenceError(matricesV, prevArrayV, phiSlice);

(*fError)(k - 1) = GetExactError(matricesV, prevArrayV, phiSlice);

// if error already achieved then stop mg iteration

fIterations = k - 1;

if ((*fErrorConvergenceNormInf)(k - 1) <= fgConvergenceError) {

Info("PoissonRelaxation3D", "%s", Form("Exact Err: %f, Iteration : %d", (*fError)(k - 1), k - 1));

break;

}

if (k == maxIteration) {

Info("PoissonRelaxation3D", "%s", Form("Exact Err: %f, Iteration : %d", (*fError)(k - 1), k - 1));

}

}

}

iOne = iOne / 2;

if (iOne < 1) {

iOne = 1;

}

jOne = jOne / 2;

if (jOne < 1) {

jOne = 1;

}

}

for (Int_t k = 0; k < phiSlice; k++) {

matricesSumChargeDensity[k]->Delete();

}

for (Int_t m = 0; m < phiSlice; m++) {

delete prevArrayV[m];

}

}

/// 3D - Solve Poisson's Equation in 3D by MultiGrid with constant phi slices

///

/// NOTE: In order for this algorithm to work, the number of nRRow and nZColumn must be a power of 2 plus one.

/// The number of nRRow and Z Column can be different.

///

/// R Row == 2**M + 1

/// Z Column == 2**N + 1

/// Phi Slice == Arbitrary but greater than 3

///

/// Solving: \f$ \nabla^{2}V(r,\phi,z) = - f(r,\phi,z) \f$

///

/// Algorithm for MultiGrid Full Cycle (FMG)

/// - Relax on the coarsest grid

/// - Do from coarsest to finest

/// - Interpolate potential from coarse -> fine

/// - Do V-Cycle to the current coarse level to the coarsest

/// - Stop if converged

///

/// DeltaPhi in Radians

/// \param matricesV TMatrixD** potential in 3D matrix \f$ V(r,\phi,z) \f$

/// \param matricesCharge TMatrixD** charge density in 3D matrix (side effect) \f$ - f(r,\phi,z) \f$

/// \param nRRow Int_t number of nRRow in the r direction of TPC

/// \param nZColumn Int_t number of nZColumn in z direction of TPC

/// \param phiSlice Int_t number of phiSlice in phi direction of T{C

/// \param maxIteration Int_t maximum iteration for relaxation method (NOT USED)

/// \param symmetry Int_t symmetry (TODO for symmetry = 1)

//

/// SYMMETRY = 0 if no phi symmetries, and no phi boundary condition

/// = 1 if we have reflection symmetry at the boundaries (eg. sector symmetry or half sector symmetries).

///

void AliTPCPoissonSolver::PoissonMultiGrid3D2D(TMatrixD** matricesV, TMatrixD** matricesCharge, Int_t nRRow,

Int_t nZColumn, Int_t phiSlice, Int_t symmetry)

{

const Float_t gridSizeR =

(AliTPCPoissonSolver::fgkOFCRadius - AliTPCPoissonSolver::fgkIFCRadius) / (nRRow - 1); // h_{r}

const Float_t gridSizePhi = TMath::TwoPi() / phiSlice; // h_{phi}

const Float_t gridSizeZ = AliTPCPoissonSolver::fgkTPCZ0 / (nZColumn - 1); // h_{z}

const Float_t ratioPhi =

gridSizeR * gridSizeR / (gridSizePhi * gridSizePhi); // ratio_{phi} = gridSize_{r} / gridSize_{phi}

const Float_t ratioZ = gridSizeR * gridSizeR / (gridSizeZ * gridSizeZ); // ratio_{Z} = gridSize_{r} / gridSize_{z}

// error tolerate

//const Float_t ERR = 1e-8;

Double_t convergenceError;

Info("PoissonMultiGrid3D2D", "%s", Form("in Poisson Solver 3D multiGrid semi coarsening nRRow=%d, cols=%d, phiSlice=%d \n", nRRow, nZColumn, phiSlice));

// Check that the number of nRRow and nZColumn is suitable for a binary expansion

if (!IsPowerOfTwo((nRRow - 1))) {

Error("PoissonMultiGrid3D2D", "Poisson3DMultiGrid - Error in the number of nRRow. Must be 2**M + 1");

return;

}

if (!IsPowerOfTwo((nZColumn - 1))) {

Error("PoissonMultiGrid3D2D", "Poisson3DMultiGrid - Error in the number of nZColumn. Must be 2**N - 1");

return;

}

if (phiSlice <= 3) {

Error("PoissonMultiGrid3D2D", "Poisson3DMultiGrid - Error in the number of phiSlice. Must be larger than 3");

return;

}

if (phiSlice > 1000) {

Error("PoissonMultiGrid3D2D", "Poisson3D phiSlice > 1000 is not allowed (nor wise) ");

return;

}

// Solve Poisson's equation in cylindrical coordinates by multiGrid technique

// Allow for different size grid spacing in R and Z directions

Int_t nGridRow = 0; // number grid

Int_t nGridCol = 0; // number grid

Int_t nnRow;

Int_t nnCol;

nnRow = nRRow;

while (nnRow >>= 1) {

nGridRow++;

}

nnCol = nZColumn;

while (nnCol >>= 1) {

nGridCol++;

}

Int_t nLoop = TMath::Max(nGridRow, nGridCol); // Calculate the number of nLoop for the binary expansion

nLoop = (nLoop > fMgParameters.maxLoop) ? fMgParameters.maxLoop : nLoop;

Int_t count;

Int_t iOne = 1; // index i in gridSize r (original)

Int_t jOne = 1; // index j in gridSize z (original)

Int_t tnRRow = nRRow, tnZColumn = nZColumn;

std::vector<TMatrixD**> tvChargeFMG(nLoop); // charge is restricted in full multiGrid

std::vector<TMatrixD**> tvArrayV(nLoop); // potential <--> error

std::vector<TMatrixD**> tvCharge(nLoop); // charge <--> residue

std::vector<TMatrixD**> tvResidue(nLoop); // residue calculation

std::vector<TMatrixD**> tvPrevArrayV(nLoop); // error calculation

for (count = 1; count <= nLoop; count++) {

tnRRow = iOne == 1 ? nRRow : nRRow / iOne + 1;

tnZColumn = jOne == 1 ? nZColumn : nZColumn / jOne + 1;

tvResidue[count - 1] = new TMatrixD*[phiSlice];

tvPrevArrayV[count - 1] = new TMatrixD*[phiSlice];

for (Int_t k = 0; k < phiSlice; k++) {

tvResidue[count - 1][k] = new TMatrixD(tnRRow, tnZColumn);

tvPrevArrayV[count - 1][k] = new TMatrixD(tnRRow, tnZColumn);

}

// memory for the finest grid is from parameters

if (count == 1) {

tvChargeFMG[count - 1] = matricesCharge;

tvArrayV[count - 1] = matricesV;

tvCharge[count - 1] = matricesCharge;

} else {

// allocate for coarser grid

tvChargeFMG[count - 1] = new TMatrixD*[phiSlice];

tvArrayV[count - 1] = new TMatrixD*[phiSlice];

tvCharge[count - 1] = new TMatrixD*[phiSlice];

for (Int_t k = 0; k < phiSlice; k++) {

tvArrayV[count - 1][k] = new TMatrixD(tnRRow, tnZColumn);

tvCharge[count - 1][k] = new TMatrixD(tnRRow, tnZColumn);

tvChargeFMG[count - 1][k] = new TMatrixD(tnRRow, tnZColumn);

}

Restrict3D(tvChargeFMG[count - 1], tvChargeFMG[count - 2], tnRRow, tnZColumn, phiSlice, phiSlice);

RestrictBoundary3D(tvArrayV[count - 1], tvArrayV[count - 2], tnRRow, tnZColumn, phiSlice, phiSlice);

}

iOne = 2 * iOne; // doubling

jOne = 2 * jOne; // doubling

}

Float_t h, h2, radius;

Float_t tempRatioPhi, tempRatioZ;

std::vector<float> coefficient1(

nRRow); // coefficient1(nRRow) for storing (1 + h_{r}/2r_{i}) from central differences in r direction

std::vector<float> coefficient2(

nRRow); // coefficient2(nRRow) for storing (1 + h_{r}/2r_{i}) from central differences in r direction

std::vector<float> coefficient3(

nRRow); // coefficient3(nRRow) for storing (1/r_{i}^2) from central differences in phi direction

std::vector<float> coefficient4(nRRow); // coefficient4(nRRow) for storing 1/2

std::vector<float> inverseCoefficient4(nRRow); // inverse of coefficient4(nRRow)

// Case full multi grid (FMG)

if (fMgParameters.cycleType == kFCycle) {

// 1) Relax on the coarsest grid

iOne = iOne / 2;

jOne = jOne / 2;

tnRRow = iOne == 1 ? nRRow : nRRow / iOne + 1;

tnZColumn = jOne == 1 ? nZColumn : nZColumn / jOne + 1;

h = gridSizeR * iOne;

h2 = h * h;

tempRatioPhi = ratioPhi * iOne * iOne; // Used tobe divided by ( m_one * m_one ) when m_one was != 1

tempRatioZ = ratioZ * iOne * iOne / (jOne * jOne);

for (Int_t i = 1; i < tnRRow - 1; i++) {

radius = AliTPCPoissonSolver::fgkIFCRadius + i * h;

coefficient1[i] = 1.0 + h / (2 * radius);

coefficient2[i] = 1.0 - h / (2 * radius);

coefficient3[i] = tempRatioPhi / (radius * radius);

coefficient4[i] = 0.5 / (1.0 + tempRatioZ + coefficient3[i]);

}

// relax on the coarsest level

Relax3D(tvArrayV[nLoop - 1], tvChargeFMG[nLoop - 1], tnRRow, tnZColumn, phiSlice, symmetry, h2, tempRatioZ,

coefficient1,

coefficient2, coefficient3, coefficient4);

// 2) Do multiGrid v-cycle from coarsest to finest

for (count = nLoop - 2; count >= 0; count--) {

// move to finer grid

iOne = iOne / 2;

jOne = jOne / 2;

tnRRow = iOne == 1 ? nRRow : nRRow / iOne + 1;

tnZColumn = jOne == 1 ? nZColumn : nZColumn / jOne + 1;

// 2) a) Interpolate potential for h -> 2h (coarse -> fine)

Interp3D(tvArrayV[count], tvArrayV[count + 1], tnRRow, tnZColumn, phiSlice, phiSlice);

// 2) c) Copy the restricted charge to charge for calculation

for (Int_t m = 0; m < phiSlice; m++) {

*tvCharge[count][m] = *tvChargeFMG[count][m]; //copy

}

// 2) c) Do V cycle fMgParameters.nMGCycle times at most

for (Int_t mgCycle = 0; mgCycle < fMgParameters.nMGCycle; mgCycle++) {

// Copy the potential to temp array for convergence calculation

for (Int_t m = 0; m < phiSlice; m++) {

*tvPrevArrayV[count][m] = *tvArrayV[count][m]; //copy

}

// 2) c) i) Call V cycle from grid count+1 (current fine level) to nLoop (coarsest)

VCycle3D2D(nRRow, nZColumn, phiSlice, symmetry, count + 1, nLoop, fMgParameters.nPre, fMgParameters.nPost,

gridSizeR, ratioZ, ratioPhi, tvArrayV, tvCharge, tvResidue, coefficient1, coefficient2, coefficient3,

coefficient4, inverseCoefficient4);

convergenceError = GetConvergenceError(tvArrayV[count], tvPrevArrayV[count], phiSlice);

if (count == 0) {

(*fErrorConvergenceNormInf)(mgCycle) = convergenceError;

(*fError)(mgCycle) = GetExactError(matricesV, tvPrevArrayV[count], phiSlice);

}

/// if already converge just break move to finer grid

if (convergenceError <= fgConvergenceError) {

fIterations = mgCycle + 1;

break;

}

}

}

} // Case V multi grid (VMG)

else if (fMgParameters.cycleType == kVCycle) {

Int_t gridFrom = 1;

Int_t gridTo = nLoop;

// do v cycle fMgParameters.nMGCycle from the coarsest to finest

for (Int_t mgCycle = 0; mgCycle < fMgParameters.nMGCycle; mgCycle++) {

// copy to store previous potential

for (Int_t m = 0; m < phiSlice; m++) {

*tvPrevArrayV[0][m] = *tvArrayV[0][m]; //copy

}

// Do V Cycle for constant phiSlice

VCycle3D2D(nRRow, nZColumn, phiSlice, symmetry, gridFrom, gridTo, fMgParameters.nPre, fMgParameters.nPost,

gridSizeR, ratioZ, ratioPhi, tvArrayV, tvCharge, tvResidue, coefficient1, coefficient2, coefficient3,

coefficient4, inverseCoefficient4);

// convergence error

convergenceError = GetConvergenceError(tvArrayV[0], tvPrevArrayV[0], phiSlice);

(*fErrorConvergenceNormInf)(mgCycle) = convergenceError;

(*fError)(mgCycle) = GetExactError(matricesV, tvPrevArrayV[0], phiSlice);

// if error already achieved then stop mg iteration

if (convergenceError <= fgConvergenceError) {

fIterations = mgCycle + 1;

break;

}

}

}

// Deallocate memory

for (count = 1; count <= nLoop; count++) {

delete[] tvResidue[count - 1];

delete[] tvPrevArrayV[count - 1];

if (count > 1) {

delete[] tvChargeFMG[count - 1];

delete[] tvArrayV[count - 1];

delete[] tvCharge[count - 1];

}

}

}