From f5fda7f5beae5da3b2007e303e3bcfd7cd8d8fc4 Mon Sep 17 00:00:00 2001 From: Matthias Kleiner Date: Thu, 19 Aug 2021 15:06:26 +0200 Subject: [PATCH 1/4] [SC] Adding possibility to use float for calculations - setting default values for omega tau - adding function to calculate the drift path for an electron - adding simple debug function - poisson solver: cleaning + using more than 1000 phi bins --- .../include/TPCSpaceCharge/RegularGrid3D.h | 26 +- .../include/TPCSpaceCharge/SpaceCharge.h | 29 +- .../include/TPCSpaceCharge/Vector.h | 3 +- .../macro/calculateDistortionsCorrections.C | 48 +- .../TPC/spacecharge/src/PoissonSolver.cxx | 430 +++++++++--------- Detectors/TPC/spacecharge/src/SpaceCharge.cxx | 186 +++++++- .../spacecharge/src/TPCSpacechargeLinkDef.h | 15 +- Detectors/TPC/spacecharge/src/TriCubic.cxx | 1 + 8 files changed, 465 insertions(+), 273 deletions(-) diff --git a/Detectors/TPC/spacecharge/include/TPCSpaceCharge/RegularGrid3D.h b/Detectors/TPC/spacecharge/include/TPCSpaceCharge/RegularGrid3D.h index c3c6c2420a998..213d96b1f2092 100644 --- a/Detectors/TPC/spacecharge/include/TPCSpaceCharge/RegularGrid3D.h +++ b/Detectors/TPC/spacecharge/include/TPCSpaceCharge/RegularGrid3D.h @@ -137,19 +137,19 @@ struct RegularGrid3D { private: inline static auto& mParamGrid = ParameterSpaceCharge::Instance(); - static constexpr unsigned int FDIM = 3; ///< dimensions of the grid (only 3 supported) - static constexpr unsigned int FZ = 0; ///< index for x coordinate - static constexpr unsigned int FR = 1; ///< index for y coordinate - static constexpr unsigned int FPHI = 2; ///< index for z coordinate - const Vector mMin{}; ///< min vertices positions of the grid - const Vector mMax{}; ///< max vertices positions of the grid - const Vector mSpacing{}; ///< spacing of the grid - const Vector mInvSpacing{}; ///< inverse spacing of grid - const Vector sMaxIndex{{static_cast(mParamGrid.NZVertices - 1.), static_cast(mParamGrid.NRVertices - 1), static_cast(mParamGrid.NPhiVertices - 1)}}; ///< max index which is on the grid in all dimensions - const Vector sNdim{{static_cast(mParamGrid.NZVertices), static_cast(mParamGrid.NRVertices), static_cast(mParamGrid.NPhiVertices)}}; ///< number of vertices for each dimension - std::vector mZVertices{}; ///< positions of vertices in x direction - std::vector mRVertices{}; ///< positions of vertices in y direction - std::vector mPhiVertices{}; ///< positions of vertices in z direction + static constexpr unsigned int FDIM = 3; ///< dimensions of the grid (only 3 supported) + static constexpr unsigned int FZ = 0; ///< index for x coordinate + static constexpr unsigned int FR = 1; ///< index for y coordinate + static constexpr unsigned int FPHI = 2; ///< index for z coordinate + const Vector mMin{}; ///< min vertices positions of the grid + const Vector mMax{}; ///< max vertices positions of the grid + const Vector mSpacing{}; ///< spacing of the grid + const Vector mInvSpacing{}; ///< inverse spacing of grid + const Vector sMaxIndex{{static_cast(mParamGrid.NZVertices - 1.), static_cast(mParamGrid.NRVertices - 1), static_cast(mParamGrid.NPhiVertices - 1)}}; ///< max index which is on the grid in all dimensions + const Vector sNdim{{static_cast(mParamGrid.NZVertices), static_cast(mParamGrid.NRVertices), static_cast(mParamGrid.NPhiVertices)}}; ///< number of vertices for each dimension + std::vector mZVertices{}; ///< positions of vertices in x direction + std::vector mRVertices{}; ///< positions of vertices in y direction + std::vector mPhiVertices{}; ///< positions of vertices in z direction void initLists(); diff --git a/Detectors/TPC/spacecharge/include/TPCSpaceCharge/SpaceCharge.h b/Detectors/TPC/spacecharge/include/TPCSpaceCharge/SpaceCharge.h index 4263db3f11aae..56962333b7a31 100644 --- a/Detectors/TPC/spacecharge/include/TPCSpaceCharge/SpaceCharge.h +++ b/Detectors/TPC/spacecharge/include/TPCSpaceCharge/SpaceCharge.h @@ -28,6 +28,7 @@ class TH3; class TH3D; +class TH3F; namespace o2 { @@ -50,11 +51,15 @@ class SpaceCharge using DataContainer = DataContainer3D; using GridProp = GridProperties; using TriCubic = TriCubicInterpolator; + using TH3DataT = std::conditional_t::value, TH3D, TH3F>; // datatype for TH3 (TH3F for DataT==float and TH3D for DataT==double) public: /// default constructor /// grid granularity has to set before constructing an object using the static function setGrid(nZVertices, nRVertices, nPhiVertices)! - SpaceCharge() = default; + /// \param omegaTau \omega \tau value + /// \param t1 value for t1 + /// \param t2 value for t2 + SpaceCharge(const DataT omegaTau = -0.35f, const DataT t1 = 1, const DataT t2 = 1) { setOmegaTauT1T2(omegaTau, t1, t2); }; /// \param nZVertices number of vertices of the grid in z direction /// \param nRVertices number of vertices of the grid in z direction @@ -182,7 +187,7 @@ class SpaceCharge /// \param z global z coordinate /// \param r global r coordinate /// \param phi global phi coordinate - DataT getChargeCyl(const DataT z, const DataT r, const DataT phi, const Side side) const; + DataT getDensityCyl(const DataT z, const DataT r, const DataT phi, const Side side) const; /// get the potential for given coordinate /// \param z global z coordinate @@ -538,7 +543,7 @@ class SpaceCharge /// \param omegaTau \omega \tau value /// \param t1 value for t1 see: ??? /// \param t2 value for t2 see: ??? - void setOmegaTauT1T2(const DataT omegaTau, const DataT t1, const DataT t2) + void setOmegaTauT1T2(const DataT omegaTau = -0.35f, const DataT t1 = 1, const DataT t2 = 1) { const DataT wt0 = t2 * omegaTau; mC0 = 1 / (1 + wt0 * wt0); @@ -567,7 +572,7 @@ class SpaceCharge static void setNThreads(const int nThreads) { sNThreads = nThreads; - o2::tpc::TriCubicInterpolator::setNThreads(nThreads); + TriCubic::setNThreads(nThreads); } /// set which kind of numerical integration is used for calcution of the integrals int Er/Ez dz, int Ephi/Ez dz, int Ez dz @@ -597,6 +602,14 @@ class SpaceCharge /// \side side of the TPC void dumpToFile(TFile& outf, const Side side) const; + /// dump sc density, potential, electric fields, global distortions/corrections to tree + /// \param outFileName name of the output file + /// \side side of the TPC + /// \param nZPoints number of vertices of the output in z + /// \param nRPoints number of vertices of the output in r + /// \param nPhiPoints number of vertices of the output in phi + void dumpToTree(const char* outFileName, const Side side, const int nZPoints = 50, const int nRPoints = 150, const int nPhiPoints = 180) const; + /// write electric fields to root file /// \param outf output file where the electrical fields will be written to /// \side side of the TPC @@ -716,7 +729,13 @@ class SpaceCharge /// \param nBinsZNew number of z bins of rebinned histogram /// \param nBinsRNew number of r bins of rebinned histogram /// \param nBinsPhiNew number of phi bins of rebinned histogram - static TH3D rebinDensityHisto(const TH3& hOrig, const unsigned short nBinsZNew, const unsigned short nBinsRNew, const unsigned short nBinsPhiNew); + static TH3DataT rebinDensityHisto(const TH3& hOrig, const unsigned short nBinsZNew, const unsigned short nBinsRNew, const unsigned short nBinsPhiNew); + + /// function to calculate the drift paths of the electron whose starting position is delivered. Electric fields must be set! + /// \param elePos global position of the start position of the electron + /// \param nSamplingPoints number of output points of the electron drift path + /// \return returns the input electron a vector of 3D-points describing the drift path of the electron + std::vector calculateElectronDriftPath(const GlobalPosition3D& elePos, const int nSamplingPoints) const; private: inline static auto& mParamGrid = ParameterSpaceCharge::Instance(); ///< parameters of the grid on which the calculations are performed diff --git a/Detectors/TPC/spacecharge/include/TPCSpaceCharge/Vector.h b/Detectors/TPC/spacecharge/include/TPCSpaceCharge/Vector.h index a98d0d3050218..e09e1d80b5392 100644 --- a/Detectors/TPC/spacecharge/include/TPCSpaceCharge/Vector.h +++ b/Detectors/TPC/spacecharge/include/TPCSpaceCharge/Vector.h @@ -17,7 +17,8 @@ #ifndef ALICEO2_TPC_VECTOR_H_ #define ALICEO2_TPC_VECTOR_H_ -#include +#include +#include namespace o2 { diff --git a/Detectors/TPC/spacecharge/macro/calculateDistortionsCorrections.C b/Detectors/TPC/spacecharge/macro/calculateDistortionsCorrections.C index 7e3c6c3fd028b..4337f355fcf67 100644 --- a/Detectors/TPC/spacecharge/macro/calculateDistortionsCorrections.C +++ b/Detectors/TPC/spacecharge/macro/calculateDistortionsCorrections.C @@ -6,11 +6,15 @@ #include "CommonUtils/TreeStreamRedirector.h" template -void calculateDistortionsAnalytical(const int globalEFieldTypeAna = 1, const int globalDistTypeAna = 1, const int eFieldTypeAna = 1, const int usePoissonSolverAna = 1, const int nSteps = 1, const int simpsonIterations = 3, const int nThreads = -1); +void calculateDistortionsAnalytical(const int sides = 0, const int globalEFieldTypeAna = 1, const int globalDistTypeAna = 1, const int eFieldTypeAna = 1, const int usePoissonSolverAna = 1, const int nSteps = 1, const int simpsonIterations = 3, const int nThreads = -1); template void calculateDistortionsFromHist(const char* path, const char* histoName, const int sides, const int globalEFieldType = 1, const int globalDistType = 1, const int nSteps = 1, const int simpsonIterations = 3, const int nThreads = -1); +int getSideStart(const int sides); +int getSideEnd(const int sides); + +/// \param sides set which sides will be simulated. sides=0: A- and C-Side, sides=1: A-Side only, sides=2: C-Side only /// \param nZVertices number of vertices of the grid in z direction /// \param nRVertices number of vertices of the grid in z direction /// \param nPhiVertices number of vertices of the grid in z direction @@ -21,10 +25,11 @@ void calculateDistortionsFromHist(const char* path, const char* histoName, const /// \param nSteps number of which are used for calculation of distortions/corrections per z-bin /// \param simpsonIterations number of iterations used in the simpson intergration /// \param nThreads number of threads which are used (if the value is -1 all threads should be used) -void calcDistAna(const unsigned short nZVertices = 129, const unsigned short nRVertices = 129, const unsigned short nPhiVertices = 180, const int globalEFieldTypeAna = 1, const int globalDistTypeAna = 1, const int eFieldTypeAna = 1, const int usePoissonSolverAna = 1, const int nSteps = 1, const int simpsonIterations = 3, const int nThreads = -1) +template +void calcDistAna(const int sides = 0, const unsigned short nZVertices = 129, const unsigned short nRVertices = 129, const unsigned short nPhiVertices = 180, const int globalEFieldTypeAna = 1, const int globalDistTypeAna = 1, const int eFieldTypeAna = 1, const int usePoissonSolverAna = 1, const int nSteps = 1, const int simpsonIterations = 3, const int nThreads = -1) { - o2::tpc::SpaceCharge::setGrid(nZVertices, nRVertices, nPhiVertices); - calculateDistortionsAnalytical(globalEFieldTypeAna, globalDistTypeAna, eFieldTypeAna, usePoissonSolverAna, nSteps, simpsonIterations, nThreads); + o2::tpc::SpaceCharge::setGrid(nZVertices, nRVertices, nPhiVertices); + calculateDistortionsAnalytical(sides, globalEFieldTypeAna, globalDistTypeAna, eFieldTypeAna, usePoissonSolverAna, nSteps, simpsonIterations, nThreads); } /// \param path path to the root file containing the 3D density histogram @@ -271,7 +276,7 @@ void writeToTree(o2::tpc::SpaceCharge& spaceCharge3D, o2::utils::TreeStre /// \param nSteps number of which are used for calculation of distortions/corrections per z-bin /// \param simpsonIterations number of iterations used in the simpson intergration template -void calculateDistortionsAnalytical(const int globalEFieldTypeAna, const int globalDistTypeAna, const int eFieldTypeAna, const int usePoissonSolverAna, const int nSteps, const int simpsonIterations, const int nThreads) +void calculateDistortionsAnalytical(const int sides, const int globalEFieldTypeAna, const int globalDistTypeAna, const int eFieldTypeAna, const int usePoissonSolverAna, const int nSteps, const int simpsonIterations, const int nThreads) { const auto integrationStrategy = o2::tpc::SpaceCharge::IntegrationStrategy::SimpsonIterative; o2::tpc::SpaceCharge spaceCharge3D; @@ -286,7 +291,7 @@ void calculateDistortionsAnalytical(const int globalEFieldTypeAna, const int glo // write to root file o2::utils::TreeStreamRedirector pcstream(TString::Format("distortions_ana_nR%hu_nZ%hu_nPhi%hu_SimpsonsIter%i.root", spaceCharge3D.getNRVertices(), spaceCharge3D.getNZVertices(), spaceCharge3D.getNPhiVertices(), simpsonIterations).Data(), "RECREATE"); - for (int iside = 0; iside < 2; ++iside) { + for (int iside = getSideStart(sides); iside < getSideEnd(sides); ++iside) { std::cout << "side: " << iside << std::endl; o2::tpc::Side side = (iside == 0) ? o2::tpc::Side::A : o2::tpc::Side::C; o2::tpc::AnalyticalFields anaFields(side); @@ -327,16 +332,7 @@ void calculateDistortionsFromHist(const char* path, const char* histoName, const fileOutSCDensity.Close(); o2::utils::TreeStreamRedirector pcstream(TString::Format("distortions_real_nR%hu_nZ%hu_nPhi%hu_SimpsonsIter%i.root", spaceCharge3D.getNRVertices(), spaceCharge3D.getNZVertices(), spaceCharge3D.getNPhiVertices(), simpsonIterations).Data(), "RECREATE"); - int iSideStart = 0; - int iSideEnd = 2; - if (sides == 1) { - // a side only - iSideEnd = 1; - } else if (sides == 2) { - // c side only - iSideStart = 1; - } - for (int iside = iSideStart; iside < iSideEnd; ++iside) { + for (int iside = getSideStart(sides); iside < getSideEnd(sides); ++iside) { std::cout << "side: " << iside << std::endl; o2::tpc::Side side = (iside == 0) ? o2::tpc::Side::A : o2::tpc::Side::C; const auto distType = globalDistType == 0 ? SC::GlobalDistType::Standard : SC ::GlobalDistType::Fast; @@ -364,3 +360,23 @@ void calculateDistortionsFromHist(const char* path, const char* histoName, const } fOut.Close(); } + +/// helper function to set the loop over the sides for the tpc +/// \param sides set for which sides the distortions/corrections will be calculated. sides=0: A- and C-Side, sides=1: A-Side only, sides=2: C-Side only +int getSideStart(const int sides) +{ + if (sides == 2) { + return 1; + } + return 0; +} + +/// helper function to set the loop over the sides for the tpc +/// \param sides set for which sides the distortions/corrections will be calculated. sides=0: A- and C-Side, sides=1: A-Side only, sides=2: C-Side only +int getSideEnd(const int sides) +{ + if (sides == 1) { + return 1; + } + return 2; +} diff --git a/Detectors/TPC/spacecharge/src/PoissonSolver.cxx b/Detectors/TPC/spacecharge/src/PoissonSolver.cxx index e20e26aeedf95..b2f1f82d26b23 100644 --- a/Detectors/TPC/spacecharge/src/PoissonSolver.cxx +++ b/Detectors/TPC/spacecharge/src/PoissonSolver.cxx @@ -74,7 +74,7 @@ void PoissonSolver::poissonMultiGrid2D(DataContainer& matricesV, const Da return; } if (!isPowerOfTwo(mParamGrid.NZVertices - 1)) { - LOGP(ERROR, "PoissonMultiGrid2D: PoissonMultiGrid - Error in the number of mParamGrid.NZVertices. Must be 2**N - 1"); + LOGP(ERROR, "PoissonMultiGrid2D: PoissonMultiGrid - Error in the number of mParamGrid.NZVertices. Must be 2**N + 1"); return; } @@ -82,8 +82,8 @@ void PoissonSolver::poissonMultiGrid2D(DataContainer& matricesV, const Da LOGP(info, "{}", fmt::format("PoissonMultiGrid2D: nGridRow={}, nGridCol={}, nLoop={}, nMGCycle={}", nGridRow, nGridCol, nLoop, MGParameters::nMGCycle)); - int iOne = 1; // in/dex - int jOne = 1; // index + unsigned int iOne = 1; // index + unsigned int jOne = 1; // index int tnRRow = mParamGrid.NRVertices; int tnZColumn = mParamGrid.NZVertices; @@ -98,26 +98,27 @@ void PoissonSolver::poissonMultiGrid2D(DataContainer& matricesV, const Da tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1; tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1; // if one just address to matrixV - tvResidue[count - 1].resize(tnRRow, tnZColumn, 1); - tvChargeFMG[count - 1].resize(tnRRow, tnZColumn, 1); - tvArrayV[count - 1].resize(tnRRow, tnZColumn, 1); - tvCharge[count - 1].resize(tnRRow, tnZColumn, 1); + const int index = count - 1; + tvResidue[index].resize(tnRRow, tnZColumn, 1); + tvChargeFMG[index].resize(tnRRow, tnZColumn, 1); + tvArrayV[index].resize(tnRRow, tnZColumn, 1); + tvCharge[index].resize(tnRRow, tnZColumn, 1); if (count == 1) { for (int iphi = iPhi; iphi <= iPhi; ++iphi) { for (int ir = 0; ir < mParamGrid.NRVertices; ++ir) { for (int iz = 0; iz < mParamGrid.NZVertices; ++iz) { - tvChargeFMG[count - 1](ir, iz, iphi) = matricesCharge(iz, ir, iphi); - tvCharge[count - 1](ir, iz, iphi) = matricesCharge(iz, ir, iphi); - tvArrayV[count - 1](ir, iz, iphi) = matricesV(iz, ir, iphi); + tvChargeFMG[index](ir, iz, iphi) = matricesCharge(iz, ir, iphi); + tvCharge[index](ir, iz, iphi) = matricesCharge(iz, ir, iphi); + tvArrayV[index](ir, iz, iphi) = matricesV(iz, ir, iphi); } } } } else { - restrict2D(tvChargeFMG[count - 1], tvChargeFMG[count - 2], tnRRow, tnZColumn, 0); + restrict2D(tvChargeFMG[index], tvChargeFMG[count - 2], tnRRow, tnZColumn, 0); } - iOne = 2 * iOne; - jOne = 2 * jOne; + iOne *= 2; + jOne *= 2; } /// full multi grid @@ -126,14 +127,14 @@ void PoissonSolver::poissonMultiGrid2D(DataContainer& matricesV, const Da LOGP(info, "PoissonMultiGrid2D: Do full cycle"); // FMG // 1) Relax on the coarsest grid - iOne = iOne * 0.5; - jOne = jOne * 0.5; + iOne /= 2; + jOne /= 2; tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1; tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1; const DataT h = gridSpacingR * nLoop; const DataT h2 = h * h; const DataT tempRatio = ratioZ * iOne * iOne / (jOne * jOne); - const DataT tempFourth = 1.0 / (2.0 + 2.0 * tempRatio); + const DataT tempFourth = 1 / (2 + 2 * tempRatio); std::vector coefficient1(tnRRow); std::vector coefficient2(tnRRow); @@ -143,9 +144,8 @@ void PoissonSolver::poissonMultiGrid2D(DataContainer& matricesV, const Da // Do VCycle from nLoop H to h for (int count = nLoop - 2; count >= 0; --count) { - - iOne = iOne * 0.5; - jOne = jOne * 0.5; + iOne /= 2; + jOne /= 2; tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1; tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1; @@ -202,17 +202,13 @@ void PoissonSolver::poissonMultiGrid3D2D(DataContainer& matricesV, const return; } if (!isPowerOfTwo((mParamGrid.NZVertices - 1))) { - LOGP(ERROR, "PoissonMultiGrid3D2D: Poisson3DMultiGrid - Error in the number of mParamGrid.NZVertices. Must be 2**N - 1"); + LOGP(ERROR, "PoissonMultiGrid3D2D: Poisson3DMultiGrid - Error in the number of mParamGrid.NZVertices. Must be 2**N + 1"); return; } if (mParamGrid.NPhiVertices <= 3) { LOGP(ERROR, "PoissonMultiGrid3D2D: Poisson3DMultiGrid - Error in the number of mParamGrid.NPhiVertices. Must be larger than 3"); return; } - if (mParamGrid.NPhiVertices > 1000) { - LOGP(ERROR, "PoissonMultiGrid3D2D: Poisson3D mParamGrid.NPhiVertices > 1000 is not allowed (nor wise)"); - return; - } const DataT gridSpacingR = getSpacingR(); const DataT gridSpacingZ = getSpacingZ(); @@ -248,29 +244,30 @@ void PoissonSolver::poissonMultiGrid3D2D(DataContainer& matricesV, const for (unsigned int count = 1; count <= nLoop; count++) { const unsigned int tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1; const unsigned int tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1; - tvResidue[count - 1].resize(tnRRow, tnZColumn, mParamGrid.NPhiVertices); - tvPrevArrayV[count - 1].resize(tnRRow, tnZColumn, mParamGrid.NPhiVertices); + const unsigned int index = count - 1; + tvResidue[index].resize(tnRRow, tnZColumn, mParamGrid.NPhiVertices); + tvPrevArrayV[index].resize(tnRRow, tnZColumn, mParamGrid.NPhiVertices); // memory for the finest grid is from parameters - tvChargeFMG[count - 1].resize(tnRRow, tnZColumn, mParamGrid.NPhiVertices); - tvArrayV[count - 1].resize(tnRRow, tnZColumn, mParamGrid.NPhiVertices); - tvCharge[count - 1].resize(tnRRow, tnZColumn, mParamGrid.NPhiVertices); + tvChargeFMG[index].resize(tnRRow, tnZColumn, mParamGrid.NPhiVertices); + tvArrayV[index].resize(tnRRow, tnZColumn, mParamGrid.NPhiVertices); + tvCharge[index].resize(tnRRow, tnZColumn, mParamGrid.NPhiVertices); if (count == 1) { for (int iphi = 0; iphi < mParamGrid.NPhiVertices; ++iphi) { for (int ir = 0; ir < mParamGrid.NRVertices; ++ir) { for (int iz = 0; iz < mParamGrid.NZVertices; ++iz) { - tvChargeFMG[count - 1](ir, iz, iphi) = matricesCharge(iz, ir, iphi); - tvArrayV[count - 1](ir, iz, iphi) = matricesV(iz, ir, iphi); + tvChargeFMG[index](ir, iz, iphi) = matricesCharge(iz, ir, iphi); + tvArrayV[index](ir, iz, iphi) = matricesV(iz, ir, iphi); } } } } else { - restrict3D(tvChargeFMG[count - 1], tvChargeFMG[count - 2], tnRRow, tnZColumn, mParamGrid.NPhiVertices, mParamGrid.NPhiVertices); - restrictBoundary3D(tvArrayV[count - 1], tvArrayV[count - 2], tnRRow, tnZColumn, mParamGrid.NPhiVertices, mParamGrid.NPhiVertices); + restrict3D(tvChargeFMG[index], tvChargeFMG[count - 2], tnRRow, tnZColumn, mParamGrid.NPhiVertices, mParamGrid.NPhiVertices); + restrictBoundary3D(tvArrayV[index], tvArrayV[count - 2], tnRRow, tnZColumn, mParamGrid.NPhiVertices, mParamGrid.NPhiVertices); } - iOne = 2 * iOne; // doubling - jOne = 2 * jOne; // doubling + iOne *= 2; // doubling + jOne *= 2; // doubling } std::vector coefficient1(mParamGrid.NRVertices); // coefficient1(mParamGrid.NRVertices) for storing (1 + h_{r}/2r_{i}) from central differences in r direction @@ -282,8 +279,8 @@ void PoissonSolver::poissonMultiGrid3D2D(DataContainer& matricesV, const // Case full multi grid (FMG) if (MGParameters::cycleType == CycleType::FCycle) { // 1) Relax on the coarsest grid - iOne = iOne * 0.5; - jOne = jOne * 0.5; + iOne /= 2; + jOne /= 2; int tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1; int tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1; @@ -301,8 +298,8 @@ void PoissonSolver::poissonMultiGrid3D2D(DataContainer& matricesV, const // 2) Do multiGrid v-cycle from coarsest to finest for (int count = nLoop - 2; count >= 0; --count) { // move to finer grid - iOne = iOne * 0.5; - jOne = jOne * 0.5; + iOne /= 2; + jOne /= 2; tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1; tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1; @@ -355,17 +352,13 @@ void PoissonSolver::poissonMultiGrid3D(DataContainer& matricesV, const Da return; } if (!isPowerOfTwo((mParamGrid.NZVertices - 1))) { - LOGP(ERROR, "PoissonMultiGrid3D: Poisson3DMultiGrid - Error in the number of mParamGrid.NZVertices. Must be 2**N - 1"); + LOGP(ERROR, "PoissonMultiGrid3D: Poisson3DMultiGrid - Error in the number of mParamGrid.NZVertices. Must be 2**N + 1"); return; } if (mParamGrid.NPhiVertices <= 3) { LOGP(ERROR, "PoissonMultiGrid3D: Poisson3DMultiGrid - Error in the number of mParamGrid.NPhiVertices. Must be larger than 3"); return; } - if (mParamGrid.NPhiVertices > 1000) { - LOGP(ERROR, "PoissonMultiGrid3D: Poisson3D mParamGrid.NPhiVertices > 1000 is not allowed (nor wise)"); - return; - } // Solve Poisson's equation in cylindrical coordinates by multi grid technique // Allow for different size grid spacing in R and Z directions @@ -386,16 +379,16 @@ void PoissonSolver::poissonMultiGrid3D(DataContainer& matricesV, const Da int nnPhi = mParamGrid.NPhiVertices; while (nnPhi % 2 == 0) { ++nGridPhi; - nnPhi *= 0.5; + nnPhi /= 2; } LOGP(info, "{}", fmt::format("PoissonMultiGrid3D: nGridRow={}, nGridCol={}, nGridPhi={}", nGridRow, nGridCol, nGridPhi)); const int nLoop = std::max({nGridRow, nGridCol, nGridPhi}); // Calculate the number of nLoop for the binary expansion // Vector for storing multi grid array - int iOne = 1; // index i in gridSize r (original) - int jOne = 1; // index j in gridSize z (original) - int kOne = 1; // index k in gridSize phi + unsigned int iOne = 1; // index i in gridSize r (original) + unsigned int jOne = 1; // index j in gridSize z (original) + unsigned int kOne = 1; // index k in gridSize phi int tnRRow = mParamGrid.NRVertices; int tnZColumn = mParamGrid.NZVertices; int tPhiSlice = mParamGrid.NPhiVertices; @@ -421,27 +414,28 @@ void PoissonSolver::poissonMultiGrid3D(DataContainer& matricesV, const Da tPhiSlice = tPhiSlice < nnPhi ? nnPhi : tPhiSlice; // allocate memory for residue - tvResidue[count - 1].resize(tnRRow, tnZColumn, tPhiSlice); - tvPrevArrayV[count - 1].resize(tnRRow, tnZColumn, tPhiSlice); - tvChargeFMG[count - 1].resize(tnRRow, tnZColumn, tPhiSlice); - tvArrayV[count - 1].resize(tnRRow, tnZColumn, tPhiSlice); - tvCharge[count - 1].resize(tnRRow, tnZColumn, tPhiSlice); + const int index = count - 1; + tvResidue[index].resize(tnRRow, tnZColumn, tPhiSlice); + tvPrevArrayV[index].resize(tnRRow, tnZColumn, tPhiSlice); + tvChargeFMG[index].resize(tnRRow, tnZColumn, tPhiSlice); + tvArrayV[index].resize(tnRRow, tnZColumn, tPhiSlice); + tvCharge[index].resize(tnRRow, tnZColumn, tPhiSlice); // memory for the finest grid is from parameters if (count == 1) { for (int iphi = 0; iphi < mParamGrid.NPhiVertices; ++iphi) { for (int ir = 0; ir < mParamGrid.NRVertices; ++ir) { for (int iz = 0; iz < mParamGrid.NZVertices; ++iz) { - tvChargeFMG[count - 1](ir, iz, iphi) = matricesCharge(iz, ir, iphi); - tvArrayV[count - 1](ir, iz, iphi) = matricesV(iz, ir, iphi); + tvChargeFMG[index](ir, iz, iphi) = matricesCharge(iz, ir, iphi); + tvArrayV[index](ir, iz, iphi) = matricesV(iz, ir, iphi); } } } - tvCharge[count - 1] = tvChargeFMG[count - 1]; + tvCharge[index] = tvChargeFMG[index]; } - iOne = 2 * iOne; // doubling - jOne = 2 * jOne; // doubling - kOne = 2 * kOne; + iOne *= 2; // doubling + jOne *= 2; // doubling + kOne *= 2; } // Case full multi grid (FMG) @@ -466,18 +460,18 @@ void PoissonSolver::poissonMultiGrid3D(DataContainer& matricesV, const Da restrictBoundary3D(tvArrayV[count - 1], tvArrayV[count - 2], tnRRow, tnZColumn, tPhiSlice, otPhiSlice); otPhiSlice = tPhiSlice; - iOne = 2 * iOne; // doubling - jOne = 2 * jOne; // doubling - kOne = 2 * kOne; + iOne *= 2; // doubling + jOne *= 2; // doubling + kOne *= 2; } // Relax on the coarsest grid // FMG // 2) Relax on the coarsest grid // move to the coarsest + 1 - iOne = iOne * 0.5; - jOne = jOne * 0.5; - kOne = kOne * 0.5; + iOne /= 2; + jOne /= 2; + kOne /= 2; tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1; tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1; tPhiSlice = kOne == 1 ? mParamGrid.NPhiVertices : mParamGrid.NPhiVertices / kOne; @@ -504,9 +498,9 @@ void PoissonSolver::poissonMultiGrid3D(DataContainer& matricesV, const Da std::fill(std::begin(coefficient4), std::end(coefficient4), 0); std::fill(std::begin(inverseCoefficient4), std::end(inverseCoefficient4), 0); - iOne = iOne * 0.5; - jOne = jOne * 0.5; - kOne = kOne * 0.5; + iOne /= 2; + jOne /= 2; + kOne /= 2; tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1; tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1; @@ -572,8 +566,8 @@ template void PoissonSolver::wCycle2D(const int gridFrom, const int gridTo, const int gamma, const int nPre, const int nPost, const DataT gridSizeR, const DataT ratio, std::vector& tvArrayV, std::vector& tvCharge, std::vector& tvResidue) { - int iOne = 1 << (gridFrom - 1); - int jOne = 1 << (gridFrom - 1); + unsigned int iOne = 1 << (gridFrom - 1); + unsigned int jOne = 1 << (gridFrom - 1); int tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1; int tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1; @@ -583,16 +577,17 @@ void PoissonSolver::wCycle2D(const int gridFrom, const int gridTo, const // 1) Go to coarsest level for (int count = gridFrom; count <= gridTo - 2; ++count) { + const int index = count - 1; const DataT h = gridSizeR * iOne; const DataT h2 = h * h; - const DataT ih2 = 1.0 / h2; + const DataT ih2 = 1 / h2; const DataT tempRatio = ratio * iOne * iOne / (jOne * jOne); - const DataT tempFourth = 1.0 / (2.0 + 2.0 * tempRatio); - const DataT inverseTempFourth = 1.0 / tempFourth; + const DataT tempFourth = 1 / (2 + 2 * tempRatio); + const DataT inverseTempFourth = 1 / tempFourth; calcCoefficients2D(1, tnRRow - 1, h, coefficient1, coefficient2); - Vector matricesCurrentV = tvArrayV[count - 1]; - Vector matricesCurrentCharge = tvCharge[count - 1]; - Vector residue = tvResidue[count - 1]; + Vector matricesCurrentV = tvArrayV[index]; + Vector matricesCurrentCharge = tvCharge[index]; + Vector residue = tvResidue[index]; // 1) Pre-Smoothing: Gauss-Seidel Relaxation or Jacobi for (int jPre = 1; jPre <= nPre; ++jPre) { @@ -602,8 +597,8 @@ void PoissonSolver::wCycle2D(const int gridFrom, const int gridTo, const // 2) Residue calculation residue2D(residue, matricesCurrentV, matricesCurrentCharge, tnRRow, tnZColumn, ih2, inverseTempFourth, tempRatio, coefficient1, coefficient2); - iOne = 2 * iOne; - jOne = 2 * jOne; + iOne *= 2; + jOne *= 2; tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1; tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1; @@ -621,13 +616,13 @@ void PoissonSolver::wCycle2D(const int gridFrom, const int gridTo, const // Go to finest grid for (int count = gridTo - 2; count >= gridFrom; --count) { - iOne = iOne * 0.5; - jOne = jOne * 0.5; + iOne /= 2; + jOne /= 2; const DataT h = gridSizeR * iOne; const DataT h2 = h * h; const DataT tempRatio = ratio * iOne * iOne / (jOne * jOne); - const DataT tempFourth = 1.0 / (2.0 + 2.0 * tempRatio); + const DataT tempFourth = 1 / (2 + 2 * tempRatio); tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1; tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1; @@ -651,8 +646,8 @@ template void PoissonSolver::vCycle2D(const int gridFrom, const int gridTo, const int nPre, const int nPost, const DataT gridSizeR, const DataT ratio, std::vector& tvArrayV, std::vector& tvCharge, std::vector& tvResidue) { - int iOne = 1 << (gridFrom - 1); - int jOne = 1 << (gridFrom - 1); + unsigned int iOne = 1 << (gridFrom - 1); + unsigned int jOne = 1 << (gridFrom - 1); int tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1; int tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1; @@ -662,28 +657,29 @@ void PoissonSolver::vCycle2D(const int gridFrom, const int gridTo, const // 1) Go to coarsest level for (int count = gridFrom; count <= gridTo - 1; ++count) { + const int index = count - 1; const DataT h = gridSizeR * iOne; const DataT h2 = h * h; - const DataT ih2 = 1.0 / h2; + const DataT ih2 = 1 / h2; const DataT tempRatio = ratio * iOne * iOne / (jOne * jOne); - const DataT tempFourth = 1.0 / (2.0 + 2.0 * tempRatio); - const DataT inverseTempFourth = 1.0 / tempFourth; + const DataT tempFourth = 1 / (2 + 2 * tempRatio); + const DataT inverseTempFourth = 1 / tempFourth; calcCoefficients2D(1, tnRRow - 1, h, coefficient1, coefficient2); for (int jPre = 1; jPre <= nPre; ++jPre) { - relax2D(tvArrayV[count - 1], tvCharge[count - 1], tnRRow, tnZColumn, h2, tempFourth, tempRatio, coefficient1, coefficient2); + relax2D(tvArrayV[index], tvCharge[index], tnRRow, tnZColumn, h2, tempFourth, tempRatio, coefficient1, coefficient2); } // 2) Residue calculation - residue2D(tvResidue[count - 1], tvArrayV[count - 1], tvCharge[count - 1], tnRRow, tnZColumn, ih2, inverseTempFourth, tempRatio, coefficient1, coefficient2); + residue2D(tvResidue[index], tvArrayV[index], tvCharge[index], tnRRow, tnZColumn, ih2, inverseTempFourth, tempRatio, coefficient1, coefficient2); - iOne = 2 * iOne; - jOne = 2 * jOne; + iOne *= 2; + jOne *= 2; tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1; tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1; //3) Restriction - restrict2D(tvCharge[count], tvResidue[count - 1], tnRRow, tnZColumn, 0); + restrict2D(tvCharge[count], tvResidue[index], tnRRow, tnZColumn, 0); //4) Zeroing coarser V std::fill(tvArrayV[count].begin(), tvArrayV[count].end(), 0); // is this necessary??? @@ -693,33 +689,33 @@ void PoissonSolver::vCycle2D(const int gridFrom, const int gridTo, const const DataT h = gridSizeR * iOne; const DataT h2 = h * h; const DataT tempRatio = ratio * iOne * iOne / (jOne * jOne); - const DataT tempFourth = 1.0 / (2.0 + 2.0 * tempRatio); + const DataT tempFourth = 1 / (2 + 2 * tempRatio); calcCoefficients2D(1, tnRRow - 1, h, coefficient1, coefficient2); relax2D(tvArrayV[gridTo - 1], tvCharge[gridTo - 1], tnRRow, tnZColumn, h2, tempFourth, tempRatio, coefficient1, coefficient2); // Go to finest grid for (int count = gridTo - 1; count >= gridFrom; count--) { + const int index = count - 1; + iOne /= 2; + jOne /= 2; - iOne = iOne * 0.5; - jOne = jOne * 0.5; - - const DataT h = gridSizeR * iOne; - const DataT h2 = h * h; - const DataT tempRatio = ratio * iOne * iOne / (jOne * jOne); - const DataT tempFourth = 1.0 / (2.0 + 2.0 * tempRatio); + const DataT hTmp = gridSizeR * iOne; + const DataT h2Tmp = hTmp * hTmp; + const DataT tempRatioTmp = ratio * iOne * iOne / (jOne * jOne); + const DataT tempFourthTmp = 1 / (2 + 2 * tempRatioTmp); tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1; tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1; // 6) Interpolation/Prolongation - addInterp2D(tvArrayV[count - 1], tvArrayV[count], tnRRow, tnZColumn, 0); + addInterp2D(tvArrayV[index], tvArrayV[count], tnRRow, tnZColumn, 0); - calcCoefficients2D(1, tnRRow - 1, h, coefficient1, coefficient2); + calcCoefficients2D(1, tnRRow - 1, hTmp, coefficient1, coefficient2); // 7) Post-Smoothing: Gauss-Seidel Relaxation for (int jPost = 1; jPost <= nPost; ++jPost) { - relax2D(tvArrayV[count - 1], tvCharge[count - 1], tnRRow, tnZColumn, h2, tempFourth, tempRatio, coefficient1, coefficient2); + relax2D(tvArrayV[index], tvCharge[index], tnRRow, tnZColumn, h2Tmp, tempFourthTmp, tempRatioTmp, coefficient1, coefficient2); } // end post smoothing } } @@ -729,62 +725,64 @@ void PoissonSolver::vCycle3D2D(const int symmetry, const int gridFrom, co std::vector& tvArrayV, std::vector& tvCharge, std::vector& tvResidue, std::vector& coefficient1, std::vector& coefficient2, std::vector& coefficient3, std::vector& coefficient4, std::vector& inverseCoefficient4) const { - int iOne = 1 << (gridFrom - 1); - int jOne = 1 << (gridFrom - 1); + unsigned int iOne = 1 << (gridFrom - 1); + unsigned int jOne = 1 << (gridFrom - 1); int tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1; int tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1; for (int count = gridFrom; count <= gridTo - 1; ++count) { + const int index = count - 1; const DataT h = getSpacingR() * iOne; const DataT h2 = h * h; - const DataT ih2 = 1.0 / h2; + const DataT ih2 = 1 / h2; const int iOne2 = iOne * iOne; const DataT tempRatioPhi = ratioPhi * iOne2; // Used tobe divided by ( m_one * m_one ) when m_one was != 1 const DataT tempRatioZ = ratioZ * iOne2 / (jOne * jOne); calcCoefficients(1, tnRRow - 1, h, tempRatioZ, tempRatioPhi, coefficient1, coefficient2, coefficient3, coefficient4); for (unsigned int i = 1; i < tnRRow - 1; ++i) { - inverseCoefficient4[i] = 1.0 / coefficient4[i]; + inverseCoefficient4[i] = 1 / coefficient4[i]; } //Info("VCycle3D2D","Before Pre-smoothing"); // 1) Pre-Smoothing: Gauss-Seidel Relaxation or Jacobi for (int jPre = 1; jPre <= nPre; ++jPre) { - relax3D(tvArrayV[count - 1], tvCharge[count - 1], tnRRow, tnZColumn, mParamGrid.NPhiVertices, symmetry, h2, tempRatioZ, coefficient1, coefficient2, coefficient3, coefficient4); + relax3D(tvArrayV[index], tvCharge[index], tnRRow, tnZColumn, mParamGrid.NPhiVertices, symmetry, h2, tempRatioZ, coefficient1, coefficient2, coefficient3, coefficient4); } // end pre smoothing // 2) Residue calculation - residue3D(tvResidue[count - 1], tvArrayV[count - 1], tvCharge[count - 1], tnRRow, tnZColumn, mParamGrid.NPhiVertices, symmetry, ih2, tempRatioZ, coefficient1, coefficient2, coefficient3, inverseCoefficient4); + residue3D(tvResidue[index], tvArrayV[index], tvCharge[index], tnRRow, tnZColumn, mParamGrid.NPhiVertices, symmetry, ih2, tempRatioZ, coefficient1, coefficient2, coefficient3, inverseCoefficient4); - iOne = 2 * iOne; - jOne = 2 * jOne; + iOne *= 2; + jOne *= 2; tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1; tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1; //3) Restriction - restrict3D(tvCharge[count], tvResidue[count - 1], tnRRow, tnZColumn, mParamGrid.NPhiVertices, mParamGrid.NPhiVertices); + restrict3D(tvCharge[count], tvResidue[index], tnRRow, tnZColumn, mParamGrid.NPhiVertices, mParamGrid.NPhiVertices); //4) Zeroing coarser V std::fill(tvArrayV[count].begin(), tvArrayV[count].end(), 0); } // coarsest grid - const DataT h = getSpacingR() * iOne; - const DataT h2 = h * h; + const DataT hTmp = getSpacingR() * iOne; + const DataT h2Tmp = hTmp * hTmp; - const int iOne2 = iOne * iOne; - const DataT tempRatioPhi = ratioPhi * iOne2; - const DataT tempRatioZ = ratioZ * iOne2 / (jOne * jOne); + const int iOne2Tmp = iOne * iOne; + const DataT tempRatioPhiTmp = ratioPhi * iOne2Tmp; + const DataT tempRatioZTmp = ratioZ * iOne2Tmp / (jOne * jOne); - calcCoefficients(1, tnRRow - 1, h, tempRatioZ, tempRatioPhi, coefficient1, coefficient2, coefficient3, coefficient4); + calcCoefficients(1, tnRRow - 1, hTmp, tempRatioZTmp, tempRatioPhiTmp, coefficient1, coefficient2, coefficient3, coefficient4); // 3) Relax on the coarsest grid - relax3D(tvArrayV[gridTo - 1], tvCharge[gridTo - 1], tnRRow, tnZColumn, mParamGrid.NPhiVertices, symmetry, h2, tempRatioZ, coefficient1, coefficient2, coefficient3, coefficient4); + relax3D(tvArrayV[gridTo - 1], tvCharge[gridTo - 1], tnRRow, tnZColumn, mParamGrid.NPhiVertices, symmetry, h2Tmp, tempRatioZTmp, coefficient1, coefficient2, coefficient3, coefficient4); // back to fine for (int count = gridTo - 1; count >= gridFrom; --count) { - iOne = iOne * 0.5; - jOne = jOne * 0.5; + const int index = count - 1; + iOne /= 2; + jOne /= 2; tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1; tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1; @@ -796,13 +794,13 @@ void PoissonSolver::vCycle3D2D(const int symmetry, const int gridFrom, co const DataT tempRatioZ = ratioZ * iOne2 / (jOne * jOne); // 4) Interpolation/Prolongation - addInterp3D(tvArrayV[count - 1], tvArrayV[count], tnRRow, tnZColumn, mParamGrid.NPhiVertices, mParamGrid.NPhiVertices); + addInterp3D(tvArrayV[index], tvArrayV[count], tnRRow, tnZColumn, mParamGrid.NPhiVertices, mParamGrid.NPhiVertices); calcCoefficients(1, tnRRow - 1, h, tempRatioZ, tempRatioPhi, coefficient1, coefficient2, coefficient3, coefficient4); // 5) Post-Smoothing: Gauss-Seidel Relaxation for (int jPost = 1; jPost <= nPost; ++jPost) { - relax3D(tvArrayV[count - 1], tvCharge[count - 1], tnRRow, tnZColumn, mParamGrid.NPhiVertices, symmetry, h2, tempRatioZ, coefficient1, coefficient2, coefficient3, coefficient4); + relax3D(tvArrayV[index], tvCharge[index], tnRRow, tnZColumn, mParamGrid.NPhiVertices, symmetry, h2, tempRatioZ, coefficient1, coefficient2, coefficient3, coefficient4); } // end post smoothing } } @@ -814,13 +812,13 @@ void PoissonSolver::vCycle3D(const int symmetry, const int gridFrom, cons { const DataT gridSpacingR = getSpacingR(); - int iOne = 1 << (gridFrom - 1); - int jOne = 1 << (gridFrom - 1); - int kOne = 1 << (gridFrom - 1); + unsigned int iOne = 1 << (gridFrom - 1); + unsigned int jOne = 1 << (gridFrom - 1); + unsigned int kOne = 1 << (gridFrom - 1); int nnPhi = mParamGrid.NPhiVertices; while (nnPhi % 2 == 0) { - nnPhi *= 0.5; + nnPhi /= 2; } int tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1; @@ -829,10 +827,11 @@ void PoissonSolver::vCycle3D(const int symmetry, const int gridFrom, cons tPhiSlice = tPhiSlice < nnPhi ? nnPhi : tPhiSlice; for (int count = gridFrom; count <= gridTo - 1; ++count) { + const int index = count - 1; const int otPhiSlice = tPhiSlice; const DataT h = gridSpacingR * iOne; const DataT h2 = h * h; - const DataT ih2 = 1.0 / h2; + const DataT ih2 = 1 / h2; const DataT tempGridSizePhiInv = tPhiSlice * INVTWOPI; // phi now is multiGrid const DataT tempRatioPhi = h2 * tempGridSizePhiInv * tempGridSizePhiInv; // ratio_{phi} = gridSize_{r} / gridSize_{phi} const DataT tempRatioZ = ratioZ * iOne * iOne / (jOne * jOne); @@ -840,49 +839,50 @@ void PoissonSolver::vCycle3D(const int symmetry, const int gridFrom, cons calcCoefficients(1, tnRRow - 1, h, tempRatioZ, tempRatioPhi, coefficient1, coefficient2, coefficient3, coefficient4); for (int i = 1; i < tnRRow - 1; ++i) { - inverseCoefficient4[i] = 1.0 / coefficient4[i]; + inverseCoefficient4[i] = 1 / coefficient4[i]; } // 1) Pre-Smoothing: Gauss-Seidel Relaxation or Jacobi for (int jPre = 1; jPre <= nPre; ++jPre) { - relax3D(tvArrayV[count - 1], tvCharge[count - 1], tnRRow, tnZColumn, tPhiSlice, symmetry, h2, tempRatioZ, coefficient1, coefficient2, coefficient3, coefficient4); + relax3D(tvArrayV[index], tvCharge[index], tnRRow, tnZColumn, tPhiSlice, symmetry, h2, tempRatioZ, coefficient1, coefficient2, coefficient3, coefficient4); } // end pre smoothing // 2) Residue calculation - residue3D(tvResidue[count - 1], tvArrayV[count - 1], tvCharge[count - 1], tnRRow, tnZColumn, tPhiSlice, symmetry, ih2, tempRatioZ, coefficient1, coefficient2, coefficient3, inverseCoefficient4); + residue3D(tvResidue[index], tvArrayV[index], tvCharge[index], tnRRow, tnZColumn, tPhiSlice, symmetry, ih2, tempRatioZ, coefficient1, coefficient2, coefficient3, inverseCoefficient4); - iOne = 2 * iOne; - jOne = 2 * jOne; - kOne = 2 * kOne; + iOne *= 2; + jOne *= 2; + kOne *= 2; tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1; tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1; tPhiSlice = mParamGrid.NPhiVertices / kOne; tPhiSlice = tPhiSlice < nnPhi ? nnPhi : tPhiSlice; //3) Restriction - restrict3D(tvCharge[count], tvResidue[count - 1], tnRRow, tnZColumn, tPhiSlice, otPhiSlice); + restrict3D(tvCharge[count], tvResidue[index], tnRRow, tnZColumn, tPhiSlice, otPhiSlice); //4) Zeroing coarser V std::fill(tvArrayV[count].begin(), tvArrayV[count].end(), 0); } // coarsest grid - const DataT h = gridSpacingR * iOne; - const DataT h2 = h * h; - const DataT tempGridSizePhiInv = tPhiSlice * INVTWOPI; // phi now is multiGrid - const DataT tempRatioPhi = h2 * tempGridSizePhiInv * tempGridSizePhiInv; // ratio_{phi} = gridSize_{r} / gridSize_{phi} - const DataT tempRatioZ = ratioZ * iOne * iOne / (jOne * jOne); + const DataT hTmp = gridSpacingR * iOne; + const DataT h2Tmp = hTmp * hTmp; + const DataT tempGridSizePhiInvTmp = tPhiSlice * INVTWOPI; // phi now is multiGrid + const DataT tempRatioPhiTmp = h2Tmp * tempGridSizePhiInvTmp * tempGridSizePhiInvTmp; // ratio_{phi} = gridSize_{r} / gridSize_{phi} + const DataT tempRatioZTmp = ratioZ * iOne * iOne / (jOne * jOne); - calcCoefficients(1, tnRRow - 1, h, tempRatioZ, tempRatioPhi, coefficient1, coefficient2, coefficient3, coefficient4); + calcCoefficients(1, tnRRow - 1, hTmp, tempRatioZTmp, tempRatioPhiTmp, coefficient1, coefficient2, coefficient3, coefficient4); // 3) Relax on the coarsest grid - relax3D(tvArrayV[gridTo - 1], tvCharge[gridTo - 1], tnRRow, tnZColumn, tPhiSlice, symmetry, h2, tempRatioZ, coefficient1, coefficient2, coefficient3, coefficient4); + relax3D(tvArrayV[gridTo - 1], tvCharge[gridTo - 1], tnRRow, tnZColumn, tPhiSlice, symmetry, h2Tmp, tempRatioZTmp, coefficient1, coefficient2, coefficient3, coefficient4); // back to fine for (int count = gridTo - 1; count >= gridFrom; --count) { + const int index = count - 1; const int otPhiSlice = tPhiSlice; - iOne = iOne * 0.5; - jOne = jOne * 0.5; - kOne = kOne * 0.5; + iOne /= 2; + jOne /= 2; + kOne /= 2; tnRRow = iOne == 1 ? mParamGrid.NRVertices : mParamGrid.NRVertices / iOne + 1; tnZColumn = jOne == 1 ? mParamGrid.NZVertices : mParamGrid.NZVertices / jOne + 1; @@ -896,13 +896,13 @@ void PoissonSolver::vCycle3D(const int symmetry, const int gridFrom, cons const DataT tempRatioZ = ratioZ * iOne * iOne / (jOne * jOne); // 4) Interpolation/Prolongation - addInterp3D(tvArrayV[count - 1], tvArrayV[count], tnRRow, tnZColumn, tPhiSlice, otPhiSlice); + addInterp3D(tvArrayV[index], tvArrayV[count], tnRRow, tnZColumn, tPhiSlice, otPhiSlice); calcCoefficients(1, tnRRow - 1, h, tempRatioZ, tempRatioPhi, coefficient1, coefficient2, coefficient3, coefficient4); // 5) Post-Smoothing: Gauss-Seidel Relaxation for (int jPost = 1; jPost <= nPost; ++jPost) { - relax3D(tvArrayV[count - 1], tvCharge[count - 1], tnRRow, tnZColumn, tPhiSlice, symmetry, h2, tempRatioZ, coefficient1, coefficient2, coefficient3, coefficient4); + relax3D(tvArrayV[index], tvCharge[index], tnRRow, tnZColumn, tPhiSlice, symmetry, h2, tempRatioZ, coefficient1, coefficient2, coefficient3, coefficient4); } } } @@ -985,7 +985,7 @@ void PoissonSolver::interp3D(Vector& matricesCurrentV, const Vector& matr if (newPhiSlice == 2 * oldPhiSlice) { for (int m = 0; m < newPhiSlice; m += 2) { // assuming no symmetry - int mm = m * 0.5; + int mm = m / 2; int mmPlus = mm + 1; int mp1 = m + 1; @@ -999,42 +999,43 @@ void PoissonSolver::interp3D(Vector& matricesCurrentV, const Vector& matr for (int j = 2; j < tnZColumn - 1; j += 2) { for (int i = 2; i < tnRRow - 1; i += 2) { - const int iHalf = i * 0.5; - const int jHalf = j * 0.5; + const int iHalf = i / 2; + const int jHalf = j / 2; matricesCurrentV(i, j, m) = matricesCurrentVC(iHalf, jHalf, mm); // point on corner lines at phi direction - matricesCurrentV(i, j, mp1) = 0.5 * (matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf, mmPlus)); + matricesCurrentV(i, j, mp1) = (matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf, mmPlus)) / 2; } } for (int j = 1; j < tnZColumn - 1; j += 2) { for (int i = 2; i < tnRRow - 1; i += 2) { - const int iHalf = i * 0.5; - const int jHalf = j * 0.5; - matricesCurrentV(i, j, m) = 0.5 * (matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf + 1, mm)); + const int iHalf = i / 2; + const int jHalf = j / 2; + matricesCurrentV(i, j, m) = (matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf + 1, mm)) / 2; // point on corner lines at phi direction - matricesCurrentV(i, j, mp1) = 0.25 * (matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf + 1, mm) + matricesCurrentVC(iHalf, jHalf, mmPlus) + matricesCurrentVC(iHalf, jHalf + 1, mmPlus)); + matricesCurrentV(i, j, mp1) = (matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf + 1, mm) + matricesCurrentVC(iHalf, jHalf, mmPlus) + matricesCurrentVC(iHalf, jHalf + 1, mmPlus)) / 4; } } for (int j = 2; j < tnZColumn - 1; j += 2) { for (int i = 1; i < tnRRow - 1; i += 2) { - const int iHalf = i * 0.5; - const int jHalf = j * 0.5; - matricesCurrentV(i, j, m) = 0.5 * (matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf + 1, jHalf, mm)); + const int iHalf = i / 2; + const int jHalf = j / 2; + matricesCurrentV(i, j, m) = (matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf + 1, jHalf, mm)) / 2; // point on line at phi direction - matricesCurrentV(i, j, mp1) = 0.25 * ((matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf, mmPlus)) + (matricesCurrentVC(iHalf + 1, jHalf, mmPlus) + matricesCurrentVC(iHalf + 1, jHalf, mm))); + matricesCurrentV(i, j, mp1) = ((matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf, mmPlus)) + (matricesCurrentVC(iHalf + 1, jHalf, mmPlus) + matricesCurrentVC(iHalf + 1, jHalf, mm))) / 4; } } for (int j = 1; j < tnZColumn - 1; j += 2) { for (int i = 1; i < tnRRow - 1; i += 2) { - const int iHalf = i * 0.5; - const int jHalf = j * 0.5; - matricesCurrentV(i, j, m) = 0.25 * ((matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf + 1, mm)) + (matricesCurrentVC(iHalf + 1, jHalf, mm) + matricesCurrentVC(iHalf + 1, jHalf + 1, mm))); + const int iHalf = i / 2; + const int jHalf = j / 2; + 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; // point at the center at phi direction - matricesCurrentV(i, j, mp1) = 0.125 * ((matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf + 1, mm) + matricesCurrentVC(iHalf, jHalf, mmPlus) + matricesCurrentVC(iHalf, jHalf + 1, mmPlus)) + - (matricesCurrentVC(iHalf + 1, jHalf, mm) + matricesCurrentVC(iHalf + 1, jHalf + 1, mm) + matricesCurrentVC(iHalf + 1, jHalf, mmPlus) + matricesCurrentVC(iHalf + 1, jHalf + 1, mmPlus))); + matricesCurrentV(i, j, mp1) = ((matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf + 1, mm) + matricesCurrentVC(iHalf, jHalf, mmPlus) + matricesCurrentVC(iHalf, jHalf + 1, mmPlus)) + + (matricesCurrentVC(iHalf + 1, jHalf, mm) + matricesCurrentVC(iHalf + 1, jHalf + 1, mm) + matricesCurrentVC(iHalf + 1, jHalf, mmPlus) + matricesCurrentVC(iHalf + 1, jHalf + 1, mmPlus))) / + 8; } } } @@ -1052,23 +1053,24 @@ void PoissonSolver::interp2D(Vector& matricesCurrentV, const Vector& matr { for (int j = 2; j < tnZColumn - 1; j += 2) { for (int i = 2; i < tnRRow - 1; i += 2) { - matricesCurrentV(i, j, iphi) = matricesCurrentVC(i * 0.5, j * 0.5, iphi); + const int jHalf = j / 2; + matricesCurrentV(i, j, iphi) = matricesCurrentVC(i / 2, jHalf, iphi); } } for (int j = 1; j < tnZColumn - 1; j += 2) { for (int i = 2; i < tnRRow - 1; i += 2) { - const int iHalf = i * 0.5; - const int jHalf = j * 0.5; - matricesCurrentV(i, j, iphi) = 0.5 * (matricesCurrentVC(iHalf, jHalf, iphi) + matricesCurrentVC(iHalf, jHalf + 1, iphi)); + const int iHalf = i / 2; + const int jHalf = j / 2; + matricesCurrentV(i, j, iphi) = (matricesCurrentVC(iHalf, jHalf, iphi) + matricesCurrentVC(iHalf, jHalf + 1, iphi)) / 2; } } for (int j = 2; j < tnZColumn - 1; j += 2) { for (int i = 1; i < tnRRow - 1; i += 2) { - const int iHalf = i * 0.5; - const int jHalf = j * 0.5; - matricesCurrentV(i, j, iphi) = 0.5 * (matricesCurrentVC(iHalf, jHalf, iphi) + matricesCurrentVC(iHalf + 1, jHalf, iphi)); + const int iHalf = i / 2; + const int jHalf = j / 2; + matricesCurrentV(i, j, iphi) = (matricesCurrentVC(iHalf, jHalf, iphi) + matricesCurrentVC(iHalf + 1, jHalf, iphi)) / 2; } } @@ -1076,9 +1078,9 @@ void PoissonSolver::interp2D(Vector& matricesCurrentV, const Vector& matr if (MGParameters::gtType == GridTransferType::Full) { for (int j = 1; j < tnZColumn - 1; j += 2) { for (int i = 1; i < tnRRow - 1; i += 2) { - const int iHalf = i * 0.5; - const int jHalf = j * 0.5; - matricesCurrentV(i, j, iphi) = 0.25 * (matricesCurrentVC(iHalf, jHalf, iphi) + matricesCurrentVC(iHalf, jHalf + 1, iphi) + matricesCurrentVC(iHalf + 1, jHalf, iphi) + matricesCurrentVC(iHalf + 1, jHalf + 1, iphi)); + const int iHalf = i / 2; + const int jHalf = j / 2; + 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; } } } @@ -1091,7 +1093,7 @@ void PoissonSolver::addInterp3D(Vector& matricesCurrentV, const Vector& m if (newPhiSlice == 2 * oldPhiSlice) { for (int m = 0; m < newPhiSlice; m += 2) { // assuming no symmetry - int mm = m * 0.5; + int mm = m / 2; int mmPlus = mm + 1; int mp1 = m + 1; @@ -1105,42 +1107,43 @@ void PoissonSolver::addInterp3D(Vector& matricesCurrentV, const Vector& m for (int j = 2; j < tnZColumn - 1; j += 2) { for (int i = 2; i < tnRRow - 1; i += 2) { - const int iHalf = i * 0.5; - const int jHalf = j * 0.5; + const int iHalf = i / 2; + const int jHalf = j / 2; matricesCurrentV(i, j, m) += matricesCurrentVC(iHalf, jHalf, mm); // point on corner lines at phi direction - matricesCurrentV(i, j, mp1) += 0.5 * (matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf, mmPlus)); + matricesCurrentV(i, j, mp1) += (matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf, mmPlus)) / 2; } } for (int j = 1; j < tnZColumn - 1; j += 2) { for (int i = 2; i < tnRRow - 1; i += 2) { - const int iHalf = i * 0.5; - const int jHalf = j * 0.5; - matricesCurrentV(i, j, m) += 0.5 * (matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf + 1, mm)); + const int iHalf = i / 2; + const int jHalf = j / 2; + matricesCurrentV(i, j, m) += (matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf + 1, mm)) / 2; // point on corner lines at phi direction - matricesCurrentV(i, j, mp1) += 0.25 * (matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf + 1, mm) + matricesCurrentVC(iHalf, jHalf, mmPlus) + matricesCurrentVC(iHalf, jHalf + 1, mmPlus)); + matricesCurrentV(i, j, mp1) += (matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf + 1, mm) + matricesCurrentVC(iHalf, jHalf, mmPlus) + matricesCurrentVC(iHalf, jHalf + 1, mmPlus)) / 4; } } for (int j = 2; j < tnZColumn - 1; j += 2) { for (int i = 1; i < tnRRow - 1; i += 2) { - const int iHalf = i * 0.5; - const int jHalf = j * 0.5; - matricesCurrentV(i, j, m) += 0.5 * (matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf + 1, jHalf, mm)); + const int iHalf = i / 2; + const int jHalf = j / 2; + matricesCurrentV(i, j, m) += (matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf + 1, jHalf, mm)) / 2; // point on line at phi direction - matricesCurrentV(i, j, mp1) += 0.25 * ((matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf, mmPlus)) + (matricesCurrentVC(iHalf + 1, jHalf, mmPlus) + matricesCurrentVC(iHalf + 1, jHalf, mm))); + matricesCurrentV(i, j, mp1) += ((matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf, mmPlus)) + (matricesCurrentVC(iHalf + 1, jHalf, mmPlus) + matricesCurrentVC(iHalf + 1, jHalf, mm))) / 4; } } for (int j = 1; j < tnZColumn - 1; j += 2) { for (int i = 1; i < tnRRow - 1; i += 2) { - const int iHalf = i * 0.5; - const int jHalf = j * 0.5; - matricesCurrentV(i, j, m) += 0.25 * ((matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf + 1, mm)) + (matricesCurrentVC(iHalf + 1, jHalf, mm) + matricesCurrentVC(iHalf + 1, jHalf + 1, mm))); + const int iHalf = i / 2; + const int jHalf = j / 2; + 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; // point at the center at phi direction - matricesCurrentV(i, j, mp1) += 0.125 * ((matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf + 1, mm) + matricesCurrentVC(iHalf, jHalf, mmPlus) + matricesCurrentVC(iHalf, jHalf + 1, mmPlus)) + - (matricesCurrentVC(iHalf + 1, jHalf, mm) + matricesCurrentVC(iHalf + 1, jHalf + 1, mm) + matricesCurrentVC(iHalf + 1, jHalf, mmPlus) + matricesCurrentVC(iHalf + 1, jHalf + 1, mmPlus))); + matricesCurrentV(i, j, mp1) += ((matricesCurrentVC(iHalf, jHalf, mm) + matricesCurrentVC(iHalf, jHalf + 1, mm) + matricesCurrentVC(iHalf, jHalf, mmPlus) + matricesCurrentVC(iHalf, jHalf + 1, mmPlus)) + + (matricesCurrentVC(iHalf + 1, jHalf, mm) + matricesCurrentVC(iHalf + 1, jHalf + 1, mm) + matricesCurrentVC(iHalf + 1, jHalf, mmPlus) + matricesCurrentVC(iHalf + 1, jHalf + 1, mmPlus))) / + 8; } } } @@ -1158,23 +1161,23 @@ void PoissonSolver::addInterp2D(Vector& matricesCurrentV, const Vector& m { for (int j = 2; j < tnZColumn - 1; j += 2) { for (int i = 2; i < tnRRow - 1; i += 2) { - matricesCurrentV(i, j, tnPhi) = matricesCurrentV(i, j, tnPhi) + matricesCurrentVC(i * 0.5, j * 0.5, tnPhi); + matricesCurrentV(i, j, tnPhi) = matricesCurrentV(i, j, tnPhi) + matricesCurrentVC(i / 2, j / 2, tnPhi); } } for (int j = 1; j < tnZColumn - 1; j += 2) { for (int i = 2; i < tnRRow - 1; i += 2) { - const int iHalf = 0.5 * i; - const int jHalf = 0.5 * j; - matricesCurrentV(i, j, tnPhi) = matricesCurrentV(i, j, tnPhi) + 0.5 * (matricesCurrentVC(iHalf, jHalf, tnPhi) + matricesCurrentVC(iHalf, jHalf + 1, tnPhi)); + const int iHalf = i / 2; + const int jHalf = j / 2; + matricesCurrentV(i, j, tnPhi) = matricesCurrentV(i, j, tnPhi) + (matricesCurrentVC(iHalf, jHalf, tnPhi) + matricesCurrentVC(iHalf, jHalf + 1, tnPhi)) / 2; } } for (int j = 2; j < tnZColumn - 1; j += 2) { for (int i = 1; i < tnRRow - 1; i += 2) { - const int iHalf = 0.5 * i; - const int jHalf = 0.5 * j; - matricesCurrentV(i, j, tnPhi) = matricesCurrentV(i, j, tnPhi) + 0.5 * (matricesCurrentVC(iHalf, jHalf, tnPhi) + matricesCurrentVC(iHalf + 1, jHalf, tnPhi)); + const int iHalf = i / 2; + const int jHalf = j / 2; + matricesCurrentV(i, j, tnPhi) = matricesCurrentV(i, j, tnPhi) + (matricesCurrentVC(iHalf, jHalf, tnPhi) + matricesCurrentVC(iHalf + 1, jHalf, tnPhi)) / 2; } } @@ -1182,9 +1185,9 @@ void PoissonSolver::addInterp2D(Vector& matricesCurrentV, const Vector& m if (MGParameters::gtType == GridTransferType::Full) { for (int j = 1; j < tnZColumn - 1; j += 2) { for (int i = 1; i < tnRRow - 1; i += 2) { - const int iHalf = 0.5 * i; - const int jHalf = 0.5 * j; - matricesCurrentV(i, j, tnPhi) = matricesCurrentV(i, j, tnPhi) + 0.25 * (matricesCurrentVC(iHalf, jHalf, tnPhi) + matricesCurrentVC(iHalf, jHalf + 1, tnPhi) + matricesCurrentVC(iHalf + 1, jHalf, tnPhi) + matricesCurrentVC(iHalf + 1, jHalf + 1, tnPhi)); + const int iHalf = i / 2; + const int jHalf = j / 2; + 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; } } } @@ -1393,7 +1396,7 @@ void PoissonSolver::restrict3D(Vector& matricesCurrentCharge, const Vecto const DataT s3 = (residue(iip1, jjp1, mp1) + residue(iip1, jjm1, mp1) + residue(iip1, jjp1, mm1) + residue(iip1, jjm1, mm1)) + (residue(iim1, jjm1, mm1) + residue(iim1, jjp1, mm1) + residue(iim1, jjm1, mp1) + residue(iim1, jjp1, mp1)); - matricesCurrentCharge(i, j, m) = 0.125 * residue(ii, jj, mm) + 0.0625 * s1 + 0.03125 * s2 + 0.015625 * s3; + matricesCurrentCharge(i, j, m) = residue(ii, jj, mm) / 8 + s1 / 16 + s2 / 32 + s3 / 64; } // end cols } // end mParamGrid.NRVertices @@ -1428,10 +1431,10 @@ void PoissonSolver::restrict2D(Vector& matricesCurrentCharge, const Vecto const int jjm1 = jj - 1; if (MGParameters::gtType == GridTransferType::Half) { // half - matricesCurrentCharge(i, j, iphi) = 0.5 * residue(ii, jj, iphi) + 0.125 * (residue(iip1, jj, iphi) + residue(iim1, jj, iphi) + residue(ii, jjp1, iphi) + residue(ii, jjm1, iphi)); + 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; } else if (MGParameters::gtType == GridTransferType::Full) { - matricesCurrentCharge(i, j, iphi) = 0.25 * residue(ii, jj, iphi) + 0.125 * (residue(iip1, jj, iphi) + residue(iim1, jj, iphi) + residue(ii, jjp1, iphi) + residue(ii, jjm1, iphi)) + - 0.0625 * (residue(iip1, jjp1, iphi) + residue(iim1, jjp1, iphi) + residue(iip1, jjm1, iphi) + residue(iim1, jjm1, iphi)); + 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 + + (residue(iip1, jjp1, iphi) + residue(iim1, jjp1, iphi) + residue(iip1, jjm1, iphi) + residue(iim1, jjm1, iphi)) / 16; } } // end cols } // end mParamGrid.NRVertices @@ -1462,7 +1465,7 @@ DataT PoissonSolver::getConvergenceError(const Vector& matricesCurrentV, const auto phiStep = prevArrayV.getNr() * prevArrayV.getNz(); // number of points in one phi slice const auto start = prevArrayV.begin() + m * phiStep; const auto end = start + phiStep; - errorArr[m] = std::inner_product(start, end, start, 0.); // inner product "Sum (matrix[a]*matrix[a])" + errorArr[m] = std::inner_product(start, end, start, DataT(0)); // inner product "Sum (matrix[a]*matrix[a])" } // return largest error return *std::max_element(std::begin(errorArr), std::end(errorArr)); @@ -1472,12 +1475,12 @@ template void PoissonSolver::calcCoefficients(unsigned int from, unsigned int to, const DataT h, const DataT tempRatioZ, const DataT tempRatioPhi, std::vector& coefficient1, std::vector& coefficient2, std::vector& coefficient3, std::vector& coefficient4) const { for (unsigned int i = from; i < to; ++i) { - const DataT radiusInv = 1. / (TPCParameters::IFCRADIUS + i * h); - const DataT hRadiusTmp = h * 0.5 * radiusInv; - coefficient1[i] = 1.0 + hRadiusTmp; - coefficient2[i] = 1.0 - hRadiusTmp; + const DataT radiusInv = 1 / (TPCParameters::IFCRADIUS + i * h); + const DataT hRadiusTmp = h * radiusInv / 2; + coefficient1[i] = 1 + hRadiusTmp; + coefficient2[i] = 1 - hRadiusTmp; coefficient3[i] = tempRatioPhi * radiusInv * radiusInv; - coefficient4[i] = 0.5 / (1.0 + tempRatioZ + coefficient3[i]); + coefficient4[i] = 1 / (2 * (1 + tempRatioZ + coefficient3[i])); } } @@ -1485,10 +1488,11 @@ template void PoissonSolver::calcCoefficients2D(unsigned int from, unsigned int to, const DataT h, std::vector& coefficient1, std::vector& coefficient2) const { for (int i = from; i < to; ++i) { - DataT radiusInvHalf = h * 0.5 / (TPCParameters::IFCRADIUS + i * h); - coefficient1[i] = 1.0 + radiusInvHalf; - coefficient2[i] = 1.0 - radiusInvHalf; + DataT radiusInvHalf = h / (2 * (TPCParameters::IFCRADIUS + i * h)); + coefficient1[i] = 1 + radiusInvHalf; + coefficient2[i] = 1 - radiusInvHalf; } } template class o2::tpc::PoissonSolver; +template class o2::tpc::PoissonSolver; diff --git a/Detectors/TPC/spacecharge/src/SpaceCharge.cxx b/Detectors/TPC/spacecharge/src/SpaceCharge.cxx index d337436ab492c..7be0cd12f4da6 100644 --- a/Detectors/TPC/spacecharge/src/SpaceCharge.cxx +++ b/Detectors/TPC/spacecharge/src/SpaceCharge.cxx @@ -23,6 +23,7 @@ #include "TGeoGlobalMagField.h" #include "TPCBase/ParameterGas.h" #include "Field/MagneticField.h" +#include "CommonUtils/TreeStreamRedirector.h" #include #include "TF1.h" @@ -670,14 +671,15 @@ void SpaceCharge::setLocalDistCorrVectorsFromFile(TFile& inpf, const Side template void SpaceCharge::fillChargeDensityFromFile(TFile& fInp, const char* name) { - const TH3* hisSCDensity3D = (TH3*)fInp.Get(name); + TH3* hisSCDensity3D = (TH3*)fInp.Get(name); fillChargeDensityFromHisto(*hisSCDensity3D); + delete hisSCDensity3D; } template void SpaceCharge::fillChargeDensityFromHisto(const TH3& hisSCDensity3D) { - TH3D hRebin = rebinDensityHisto(hisSCDensity3D, mParamGrid.NZVertices, mParamGrid.NRVertices, mParamGrid.NPhiVertices); + TH3DataT hRebin = rebinDensityHisto(hisSCDensity3D, mParamGrid.NZVertices, mParamGrid.NRVertices, mParamGrid.NPhiVertices); for (int side = Side::A; side < SIDES; ++side) { for (size_t iPhi = 0; iPhi < mParamGrid.NPhiVertices; ++iPhi) { for (size_t iR = 0; iR < mParamGrid.NRVertices; ++iR) { @@ -692,9 +694,9 @@ void SpaceCharge::fillChargeDensityFromHisto(const TH3& hisSCDensity3D) } template -TH3D SpaceCharge::rebinDensityHisto(const TH3& hOrig, const unsigned short nBinsZNew, const unsigned short nBinsRNew, const unsigned short nBinsPhiNew) +typename SpaceCharge::TH3DataT SpaceCharge::rebinDensityHisto(const TH3& hOrig, const unsigned short nBinsZNew, const unsigned short nBinsRNew, const unsigned short nBinsPhiNew) { - TH3D hRebin{}; + TH3DataT hRebin{}; const int nBinsZNewTwo = 2 * nBinsZNew; const auto phiLow = hOrig.GetXaxis()->GetBinLowEdge(1); @@ -1246,7 +1248,7 @@ void SpaceCharge::distortElectron(GlobalPosition3D& point) const } template -DataT SpaceCharge::getChargeCyl(const DataT z, const DataT r, const DataT phi, const Side side) const +DataT SpaceCharge::getDensityCyl(const DataT z, const DataT r, const DataT phi, const Side side) const { return mInterpolatorDensity[side](z, r, phi); } @@ -1407,23 +1409,163 @@ void SpaceCharge::integrateEFieldsRoot(const DataT p1r, const DataT p1phi localIntDeltaEz = getSign(formulaStruct.getSide()) * static_cast(fEz.Integral(p1z, p2z)); } +template +std::vector SpaceCharge::calculateElectronDriftPath(const GlobalPosition3D& elePos, const int nSamplingPoints) const +{ + std::vector elePaths; + elePaths.reserve(nSamplingPoints + 1); + + const DataT z0 = elePos.Z(); + const DataT r0 = elePos.Rho(); + const DataT phi0 = elePos.Phi(); + const Side side = getSide(z0); + const NumericalFields numEFields{getElectricFieldsInterpolator(side)}; + const DataT stepSize = std::abs(getZMax(side) - z0) / nSamplingPoints; + + DataT drDist = 0.0; // global distortion dR + DataT dPhiDist = 0.0; // global distortion dPhi (multiplication with R has to be done at the end) + DataT dzDist = 0.0; // global distortion dZ + int iter = 0; + for (;;) { + const DataT z0Tmp = z0 + dzDist + iter * stepSize; // starting z position + const DataT z1Tmp = regulateZ(z0Tmp + stepSize, side); // electron drifts from z0Tmp to z1Tmp + const DataT radius = regulateR(r0 + drDist, side); // current radial position of the electron + const DataT phi = regulatePhi(phi0 + dPhiDist, side); // current phi position of the electron + elePaths.emplace_back(GlobalPosition3D(radius * std::cos(phi), radius * std::sin(phi), z0Tmp)); + + DataT ddR = 0; // distortion dR for drift from z0Tmp to z1Tmp + DataT ddPhi = 0; // distortion dPhi for drift from z0Tmp to z1Tmp + DataT ddZ = 0; // distortion dZ for drift from z0Tmp to z1Tmp + + // get the distortion from interpolation of local distortions or calculate distortions with the electric field + processGlobalDistCorr(radius, phi, z0Tmp, z1Tmp, ddR, ddPhi, ddZ, numEFields); + + // add local distortions to global distortions + drDist += ddR; + dPhiDist += ddPhi; + dzDist += ddZ; + + // if one uses local distortions the interpolated value for the last bin has to be scaled. + // This has to be done because of the interpolated value is defined for a drift length of one z bin, but in the last bin the distance to the readout can be smaller than one z bin. + const bool checkReached = side == Side::A ? z1Tmp >= getZMax(side) : z1Tmp <= getZMax(side); + + // set loop to exit if the readout is reached and approximate distortion of 'missing' (one never ends exactly on the readout: z1Tmp + ddZ != ZMAX) drift distance. + // approximation is done by the current calculated values of the distortions and scaled linear to the 'missing' distance. + if (checkReached) { + const DataT endPoint = z1Tmp + ddZ; + const DataT deltaZ = getZMax(side) - endPoint; // distance from last point to read out + const DataT diff = endPoint - z0Tmp; + const DataT fac = diff != 0 ? std::abs(deltaZ / diff) : 0; // approximate the distortions for the 'missing' distance deltaZ + drDist += ddR * fac; + dPhiDist += ddPhi * fac; + dzDist += ddZ * fac; + const DataT z1TmpEnd = regulateZ(z0Tmp + stepSize, side); // electron drifts from z0Tmp to z1Tmp + const DataT radiusEnd = regulateR(r0 + drDist, side); // current radial position of the electron + const DataT phiEnd = regulatePhi(phi0 + dPhiDist, side); // current phi position of the electron + elePaths.emplace_back(GlobalPosition3D(radiusEnd * std::cos(phiEnd), radiusEnd * std::sin(phiEnd), z1TmpEnd)); + break; + } + ++iter; + } + + return elePaths; +} + +template +void SpaceCharge::dumpToTree(const char* outFileName, const Side side, const int nZPoints, const int nRPoints, const int nPhiPoints) const +{ + const DataT phiSpacing = GridProp::getGridSpacingR(nPhiPoints); + const DataT rSpacing = GridProp::getGridSpacingR(nRPoints); + const DataT zSpacing = side == Side::A ? GridProp::getGridSpacingZ(nZPoints) : -GridProp::getGridSpacingZ(nZPoints); + + o2::utils::TreeStreamRedirector pcstream(outFileName, "RECREATE"); + pcstream.GetFile()->cd(); + for (int iPhi = 0; iPhi < nPhiPoints; ++iPhi) { + DataT phiPos = iPhi * phiSpacing; + for (int iR = 0; iR < nRPoints; ++iR) { + DataT rPos = getRMin(side) + iR * rSpacing; + for (int iZ = 0; iZ < nZPoints; ++iZ) { + DataT zPos = getZMin(side) + iZ * zSpacing; + DataT density = getDensityCyl(zPos, rPos, phiPos, side); + DataT potential = getPotentialCyl(zPos, rPos, phiPos, side); + + DataT distZ{}; + DataT distR{}; + DataT distRPhi{}; + getDistortionsCyl(zPos, rPos, phiPos, side, distZ, distR, distRPhi); + + // get average distortions + DataT corrZ{}; + DataT corrR{}; + DataT corrRPhi{}; + getCorrectionsCyl(zPos, rPos, phiPos, side, corrZ, corrR, corrRPhi); + + // get average distortions + DataT eZ{}; + DataT eR{}; + DataT ePhi{}; + getElectricFieldsCyl(zPos, rPos, phiPos, side, eZ, eR, ePhi); + + pcstream << "sc" + << "phi=" << phiPos + << "r=" << rPos + << "z=" << zPos + << "scdensity=" << density + << "potential=" << potential + << "eZ=" << eZ + << "eR=" << eR + << "ePhi=" << ePhi + << "distZ=" << distZ + << "distR=" << distR + << "distRPhi=" << distRPhi + << "corrZ=" << corrZ + << "corrR=" << corrR + << "corrRPhi=" << corrRPhi + << "\n"; + } + } + } + pcstream.Close(); +} + using DataTD = double; template class o2::tpc::SpaceCharge; -using NumFields = NumericalFields; -using AnaFields = AnalyticalFields; -using DistCorrInterp = DistCorrInterpolator; -using O2TPCSpaceCharge3DCalc = SpaceCharge; - -template void O2TPCSpaceCharge3DCalc::integrateEFieldsRoot(const DataTD, const DataTD, const DataTD, const DataTD, DataTD&, DataTD&, DataTD&, const NumFields&) const; -template void O2TPCSpaceCharge3DCalc::integrateEFieldsRoot(const DataTD, const DataTD, const DataTD, const DataTD, DataTD&, DataTD&, DataTD&, const AnaFields&) const; -template void O2TPCSpaceCharge3DCalc::calcLocalDistortionsCorrections(const O2TPCSpaceCharge3DCalc::Type, const NumFields&); -template void O2TPCSpaceCharge3DCalc::calcLocalDistortionsCorrections(const O2TPCSpaceCharge3DCalc::Type, const AnaFields&); -template void O2TPCSpaceCharge3DCalc::calcLocalDistortionCorrectionVector(const NumFields&); -template void O2TPCSpaceCharge3DCalc::calcLocalDistortionCorrectionVector(const AnaFields&); -template void O2TPCSpaceCharge3DCalc::calcGlobalCorrections(const NumFields&); -template void O2TPCSpaceCharge3DCalc::calcGlobalCorrections(const AnaFields&); -template void O2TPCSpaceCharge3DCalc::calcGlobalCorrections(const DistCorrInterp&); -template void O2TPCSpaceCharge3DCalc::calcGlobalDistortions(const NumFields&); -template void O2TPCSpaceCharge3DCalc::calcGlobalDistortions(const AnaFields&); -template void O2TPCSpaceCharge3DCalc::calcGlobalDistortions(const DistCorrInterp&); +using NumFieldsD = NumericalFields; +using AnaFieldsD = AnalyticalFields; +using DistCorrInterpD = DistCorrInterpolator; +using O2TPCSpaceCharge3DCalcD = SpaceCharge; + +template void O2TPCSpaceCharge3DCalcD::integrateEFieldsRoot(const DataTD, const DataTD, const DataTD, const DataTD, DataTD&, DataTD&, DataTD&, const NumFieldsD&) const; +template void O2TPCSpaceCharge3DCalcD::integrateEFieldsRoot(const DataTD, const DataTD, const DataTD, const DataTD, DataTD&, DataTD&, DataTD&, const AnaFieldsD&) const; +template void O2TPCSpaceCharge3DCalcD::calcLocalDistortionsCorrections(const O2TPCSpaceCharge3DCalcD::Type, const NumFieldsD&); +template void O2TPCSpaceCharge3DCalcD::calcLocalDistortionsCorrections(const O2TPCSpaceCharge3DCalcD::Type, const AnaFieldsD&); +template void O2TPCSpaceCharge3DCalcD::calcLocalDistortionCorrectionVector(const NumFieldsD&); +template void O2TPCSpaceCharge3DCalcD::calcLocalDistortionCorrectionVector(const AnaFieldsD&); +template void O2TPCSpaceCharge3DCalcD::calcGlobalCorrections(const NumFieldsD&); +template void O2TPCSpaceCharge3DCalcD::calcGlobalCorrections(const AnaFieldsD&); +template void O2TPCSpaceCharge3DCalcD::calcGlobalCorrections(const DistCorrInterpD&); +template void O2TPCSpaceCharge3DCalcD::calcGlobalDistortions(const NumFieldsD&); +template void O2TPCSpaceCharge3DCalcD::calcGlobalDistortions(const AnaFieldsD&); +template void O2TPCSpaceCharge3DCalcD::calcGlobalDistortions(const DistCorrInterpD&); + +using DataTF = float; +template class o2::tpc::SpaceCharge; + +using NumFieldsF = NumericalFields; +using AnaFieldsF = AnalyticalFields; +using DistCorrInterpF = DistCorrInterpolator; +using O2TPCSpaceCharge3DCalcF = SpaceCharge; + +template void O2TPCSpaceCharge3DCalcF::integrateEFieldsRoot(const DataTF, const DataTF, const DataTF, const DataTF, DataTF&, DataTF&, DataTF&, const NumFieldsF&) const; +template void O2TPCSpaceCharge3DCalcF::integrateEFieldsRoot(const DataTF, const DataTF, const DataTF, const DataTF, DataTF&, DataTF&, DataTF&, const AnaFieldsF&) const; +template void O2TPCSpaceCharge3DCalcF::calcLocalDistortionsCorrections(const O2TPCSpaceCharge3DCalcF::Type, const NumFieldsF&); +template void O2TPCSpaceCharge3DCalcF::calcLocalDistortionsCorrections(const O2TPCSpaceCharge3DCalcF::Type, const AnaFieldsF&); +template void O2TPCSpaceCharge3DCalcF::calcLocalDistortionCorrectionVector(const NumFieldsF&); +template void O2TPCSpaceCharge3DCalcF::calcLocalDistortionCorrectionVector(const AnaFieldsF&); +template void O2TPCSpaceCharge3DCalcF::calcGlobalCorrections(const NumFieldsF&); +template void O2TPCSpaceCharge3DCalcF::calcGlobalCorrections(const AnaFieldsF&); +template void O2TPCSpaceCharge3DCalcF::calcGlobalCorrections(const DistCorrInterpF&); +template void O2TPCSpaceCharge3DCalcF::calcGlobalDistortions(const NumFieldsF&); +template void O2TPCSpaceCharge3DCalcF::calcGlobalDistortions(const AnaFieldsF&); +template void O2TPCSpaceCharge3DCalcF::calcGlobalDistortions(const DistCorrInterpF&); diff --git a/Detectors/TPC/spacecharge/src/TPCSpacechargeLinkDef.h b/Detectors/TPC/spacecharge/src/TPCSpacechargeLinkDef.h index 0517a82d7d884..316bc1dd86256 100644 --- a/Detectors/TPC/spacecharge/src/TPCSpacechargeLinkDef.h +++ b/Detectors/TPC/spacecharge/src/TPCSpacechargeLinkDef.h @@ -22,8 +22,8 @@ #pragma link C++ nestedtypedef; #pragma link C++ struct o2::tpc::ParameterSpaceCharge + ; -#pragma link C++ class o2::tpc::MGParameters + ; -#pragma link C++ class o2::tpc::TPCParameters < double> + ; +#pragma link C++ struct o2::tpc::MGParameters + ; +#pragma link C++ struct o2::tpc::TPCParameters < double> + ; #pragma link C++ class o2::tpc::AnalyticalFields < double> + ; #pragma link C++ class o2::tpc::RegularGrid3D < double> + ; #pragma link C++ class o2::tpc::TriCubicInterpolator < double> + ; @@ -33,5 +33,14 @@ #pragma link C++ class o2::tpc::NumericalFields < double> + ; #pragma link C++ class o2::tpc::DistCorrInterpolator < double> + ; #pragma link C++ class o2::tpc::GridProperties < double> + ; - +#pragma link C++ struct o2::tpc::TPCParameters < float> + ; +#pragma link C++ class o2::tpc::AnalyticalFields < float> + ; +#pragma link C++ class o2::tpc::RegularGrid3D < float> + ; +#pragma link C++ class o2::tpc::TriCubicInterpolator < float> + ; +#pragma link C++ class o2::tpc::DataContainer3D < float> + ; +#pragma link C++ class o2::tpc::PoissonSolver < float> + ; +#pragma link C++ class o2::tpc::SpaceCharge < float> + ; +#pragma link C++ class o2::tpc::NumericalFields < float> + ; +#pragma link C++ class o2::tpc::DistCorrInterpolator < float> + ; +#pragma link C++ class o2::tpc::GridProperties < float> + ; #endif diff --git a/Detectors/TPC/spacecharge/src/TriCubic.cxx b/Detectors/TPC/spacecharge/src/TriCubic.cxx index 5c41543c0b36d..f5ddac8aa7a1d 100644 --- a/Detectors/TPC/spacecharge/src/TriCubic.cxx +++ b/Detectors/TPC/spacecharge/src/TriCubic.cxx @@ -1253,3 +1253,4 @@ void TriCubicInterpolator::setValues(const int iz, const int ir, const in } template class o2::tpc::TriCubicInterpolator; +template class o2::tpc::TriCubicInterpolator; From 7a2e55a53ee1d108a2dadef785c6c679cdd817a9 Mon Sep 17 00:00:00 2001 From: Matthias Kleiner Date: Thu, 19 Aug 2021 18:19:50 +0200 Subject: [PATCH 2/4] fixing function name getChargeCyl -> getDensityCyl --- .../TPC/reconstruction/macro/createTPCSpaceChargeCorrection.C | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Detectors/TPC/reconstruction/macro/createTPCSpaceChargeCorrection.C b/Detectors/TPC/reconstruction/macro/createTPCSpaceChargeCorrection.C index c04dff8f35f67..48d90d933c59e 100644 --- a/Detectors/TPC/reconstruction/macro/createTPCSpaceChargeCorrection.C +++ b/Detectors/TPC/reconstruction/macro/createTPCSpaceChargeCorrection.C @@ -222,7 +222,7 @@ void debugInterpolation(utils::TreeStreamRedirector& pcstream, float gxyzf[3] = {gx, gy, gz}; float pointCyl[3] = {gz, r, phi}; double efield[3] = {0.0, 0.0, 0.0}; - double charge = spaceCharge->getChargeCyl(pointCyl[0], pointCyl[1], pointCyl[2], side); + double charge = spaceCharge->getDensityCyl(pointCyl[0], pointCyl[1], pointCyl[2], side); double potential = spaceCharge->getPotentialCyl(pointCyl[0], pointCyl[1], pointCyl[2], side); spaceCharge->getElectricFieldsCyl(pointCyl[0], pointCyl[1], pointCyl[2], side, efield[0], efield[1], efield[2]); spaceCharge->getLocalDistortionsCyl(pointCyl[0], pointCyl[1], pointCyl[2], side, ldD[0], ldD[1], ldD[2]); @@ -379,7 +379,7 @@ void debugGridpoints(utils::TreeStreamRedirector& pcstream, const o2::gpu::TPCFa double efield[3] = {0.0, 0.0, 0.0}; double distLocal[3] = {0.0, 0.0, 0.0}; double dist[3] = {0.0, 0.0, 0.0}; - double charge = spaceCharge->getChargeCyl(pointCyl[0], pointCyl[1], pointCyl[2], side); + double charge = spaceCharge->getDensityCyl(pointCyl[0], pointCyl[1], pointCyl[2], side); double potential = spaceCharge->getPotentialCyl(pointCyl[0], pointCyl[1], pointCyl[2], side); spaceCharge->getElectricFieldsCyl(pointCyl[0], pointCyl[1], pointCyl[2], side, efield[0], efield[1], efield[2]); spaceCharge->getLocalDistortionsCyl(pointCyl[0], pointCyl[1], pointCyl[2], side, distLocal[0], distLocal[1], distLocal[2]); From 43567bc0a54547cb1f92de2ccf4b272c3444bf81 Mon Sep 17 00:00:00 2001 From: Matthias Kleiner Date: Fri, 20 Aug 2021 15:55:19 +0200 Subject: [PATCH 3/4] fixing typo --- Detectors/TPC/spacecharge/src/SpaceCharge.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Detectors/TPC/spacecharge/src/SpaceCharge.cxx b/Detectors/TPC/spacecharge/src/SpaceCharge.cxx index 7be0cd12f4da6..ea493941da748 100644 --- a/Detectors/TPC/spacecharge/src/SpaceCharge.cxx +++ b/Detectors/TPC/spacecharge/src/SpaceCharge.cxx @@ -1474,7 +1474,7 @@ std::vector SpaceCharge::calculateElectronDriftPath(con template void SpaceCharge::dumpToTree(const char* outFileName, const Side side, const int nZPoints, const int nRPoints, const int nPhiPoints) const { - const DataT phiSpacing = GridProp::getGridSpacingR(nPhiPoints); + const DataT phiSpacing = GridProp::getGridSpacingPhi(nPhiPoints); const DataT rSpacing = GridProp::getGridSpacingR(nRPoints); const DataT zSpacing = side == Side::A ? GridProp::getGridSpacingZ(nZPoints) : -GridProp::getGridSpacingZ(nZPoints); From 346ecd4c110a4fb9c3ed9c80b373ba7ef5b06833 Mon Sep 17 00:00:00 2001 From: Matthias Kleiner Date: Tue, 31 Aug 2021 09:23:06 +0200 Subject: [PATCH 4/4] fixing error during compiling of macro --- GPU/TPCFastTransformation/macro/generateTPCCorrectionNTuple.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GPU/TPCFastTransformation/macro/generateTPCCorrectionNTuple.C b/GPU/TPCFastTransformation/macro/generateTPCCorrectionNTuple.C index 2eb7e1f951af4..214235217d652 100644 --- a/GPU/TPCFastTransformation/macro/generateTPCCorrectionNTuple.C +++ b/GPU/TPCFastTransformation/macro/generateTPCCorrectionNTuple.C @@ -25,9 +25,9 @@ #include "TH3.h" #include "DataFormatsTPC/Defs.h" -#include "TPCSpaceCharge/SpaceCharge.h" #include "GPU/TPCFastTransform.h" #include "TPCReconstruction/TPCFastTransformHelperO2.h" +#include "TPCSpaceCharge/SpaceCharge.h" #endif