TPC SC: faster and more robust space-charge distortion calculation

Major improvements:
- Poisson solver: anisotropy aware multigrid with phi line relaxation. Fine
  phi grids converge now (1800 phi bins did not converge before), relative
  stopping criterion independent of units and grid size, parallelised
- tricubic interpolation without Vc, the three components of the fields and
  of the distortions/corrections are interpolated together
- fix of the global distortions (standard method) at the readout and several
  other bug fixes in SpaceCharge, no crash or endless loop for invalid
  interpolation positions
- faster filling of the density from histograms and faster reading of dumped
  files (ZSTD)
- SpaceCharge objects read from a file are directly usable (the grid is set up
  by a read rule, initAfterReadingFromFile() is not needed anymore)
- the slow space-charge correction is removed from TPCFastTransform
- closure test of the global distortions and corrections

Time per TPC side for calculateDistortionsCorrections() (default: fast
global distortions), realistic space-charge density, 32 threads (AMD EPYC
7702P):

  grid (z x r x phi)    old              new       speed-up
  129 x 129 x  180      4.4 s            0.49 s     9x
  257 x 257 x  180     50 s              1.6 s     32x
  129 x 129 x 1800     68 s (*)          4.7 s     15x

  (*) Poisson solver not converged after 200 V-cycles

  Poisson solver only   old              new
  129 x 129 x  180      3.4 s            0.11 s
  257 x 257 x  180     47 s              0.33 s
  129 x 129 x 1800     62 s (*)          1.55 s

  other                                  old       new
  density from TH3 (3600x257x257 -> 129x129x180)    11.7 s    0.3 s
  read all dumped containers of one side (1800 phi) 2.7 s     1.0 s

The global distortions and corrections agree with the old code to a few um
rms (max. ~30 um for r < 248 cm). Larger differences are only at the outer
radial edge (extrapolation) and from the fixes listed below.

Poisson solver
- phi is only coarsened if its coupling (h_rz / (r_min * h_phi))^2 is not
  weak compared to r/z (MGParameters::phiCoarseningThreshold, default 0.5).
  On levels where the phi coupling dominates, each phi ring is solved exactly
  (periodic tridiagonal system, Thomas algorithm + Sherman-Morrison) instead
  of the point-wise relaxation, which does not smooth the error there
  (MGParameters::phiLineRelaxThreshold, default 4, 0 disables). The levels
  are described by a list of MGLevel. All levels converge in 3-15 V-cycles.
- relative stopping criterion (MGParameters::relativeConvergence, default
  true): relative change of the potential between two V-cycles < 1e-8
  (1e-6 for float, which stays at ~1e-7 because of the rounding). The
  previous absolute criterion (sum of dV^2 < 1e-6) depended on the units of
  the potential and on the grid size. The error w.r.t. the converged solution
  is now ~3e-9 of max|V|. relativeConvergence = false with
  setConvergenceError(1e-6) reproduces the previous criterion.
- parallelised: red-black relaxation over phi, line relaxation (factorisation
  and blocks of lines), restriction, interpolation, convergence error,
  allocation of the level buffers and the copies between the V-cycles
- removed unused options (available in the git history): 2D solver, 3D solver
  with constant number of phi slices, W cycle (not functional), Jacobi
  relaxation, half weighting, the symmetry argument (diverged; also removed
  from SpaceCharge::poissonSolver() and calcEField()) and the MGParameters
  isFull3D, cycleType, gtType, relaxType, maxLoop, gamma with their enums
- merged duplicated code (interpolate3D<Add>(), restrictBoundary2D())
- fixes: minimum radius of the grid instead of the hard coded IFC radius for
  the coefficients, Vector3D constructor (list initialisation gave a single
  element) and size()

Tricubic interpolation
- separable Catmull-Rom interpolation (z, then r, then phi) with the values
  taken directly from the grid. The cell (position, index, weights) is
  calculated once for the three components (new
  TriCubicInterpolator::interpolate(z, r, phi, other1, other2, ...), used by
  NumericalFields::evalFields() and DistCorrInterpolator::evalDistCorr()) and
  the three sums are done in one pass. The results change at the rounding
  level (< 2.3e-14 of max|V| for double).
- the cells at the boundary use a generic extrapolation of the stencil points
  outside of the grid instead of 26 hand-written cases. This fixes the
  extrapolation at the (first z, last r) and (last z, first r) edges, which
  was done along the wrong diagonal.
- no Vc: Vector.h is removed (std::array in RegularGrid3D), the public
  headers and the dictionary do not include Vc anymore and TPCSpaceCharge
  does not link Vc
- fixes: without extrapolation, queries at or beyond the last vertex in r or z
  returned the value one grid spacing inwards; the move constructor did not
  take over the extrapolation setting
- non-finite query positions return NaN instead of reading outside of the
  grid, and positions far outside of the grid in phi are wrapped with fmod
  instead of loops which did not terminate for huge values

SpaceCharge
- calcGlobalDistortions(): in the last step before the readout the local
  distortions were scaled to the remaining drift length, but the correction
  for the missing distance assumed a full z bin, and std::abs turned the
  correction for electrons which passed the readout into an addition. The z
  bias of the standard method w.r.t. a direct integration of the fields
  changes from -83 um to -12 um, the closure with the global corrections in z
  from 110 um to 2.7 um rms. Standard and fast method agree in z to ~2 um.
- correctElectron(): the z coordinate of the corrected point was set to
  y + corrY
- calcGlobalDistCorrIterativeLinearCartesian(): the scaled contribution was
  lost for corrections with a scaling object; float precision of the angles
- calcGlobalDistCorrIterative*(): reference to a member of scSCale also for
  scSCale == nullptr (undefined behaviour)
- calcGlobalCorrections() with type 1: dR was copied to dRPhi and dZ for the
  vertices out of the volume
- scaleChargeDensityStack(): phi of the vertex taken from the radial index
- fillChargeDensityFromHisto(file, ...): nullptr dereference if a histogram
  is missing
- getStepWidth(): integer division
- DataContainer3D: the element-wise operators did not check the grid of the
  other container (out-of-bounds read)
- reading from a file: the grid properties are transient and were only set up
  for the default number of vertices. Without calling
  initAfterReadingFromFile() an object with a different grid gave wrong
  values or read outside of the data. A read rule in the LinkDef now calls it
  for every SpaceCharge object read from a file (calling it explicitly, e.g.
  in the Digitizer, is still possible).
- warning if the global distortions/corrections contain non-finite values
  (e.g. for unphysically large fields)
- the functions processing both sides used only two threads (nested OpenMP
  regions are not active); the sides are processed one after the other with
  all threads, and the number of threads is passed to the Poisson solver
- the density is rebinned from the histograms of both sides directly, without
  a merged copy (1.9 GB for 3600x257x257 bins)
- iterative Simpson integration: the fields at the start point are evaluated
  once per step instead of once per iteration
- calculateDistortionsCorrections() and initAllBuffers() allocate the
  containers concurrently: the first access to the memory of a large
  container by the single thread zeroing it took ~1.2 s per side at 1800 phi
- the dumped containers use ZSTD instead of ZLIB compression: same file size
  (the data are hardly compressible), ~2.7 times faster reading. Old files are
  read as before.

GPU TPCFastTransform
- TPCSlowSpaceChargeCorrection held a SpaceCharge<float> object with the
  global corrections (setSlowTPCSCCorrection()), which replaced the spline
  correction in TransformLocal() on the CPU. It is not used anymore (the
  reconstruction uses the POD correction object) and is removed together with
  setSlowTPCSCCorrection(), getCorrectionSlow() and the dependency of
  TPCFastTransformation on O2::TPCSpaceCharge (O2::DataFormatsTPC is linked
  directly instead).
- the persistent member mCorrectionSlow is removed: ClassDefNV of
  TPCFastTransform is increased to 6 (schema evolution). An empty
  dictionary-only stub of TPCSlowSpaceChargeCorrection is kept, so that
  reading existing objects (ClassDef <= 5, e.g. the CCDB maps) does not print
  a missing-dictionary warning. It can be removed once old maps are no longer
  read.

Tests
- Poisson solver: phi coarsening, phi line relaxation (720 phi bins)
- tricubic interpolation: interpolation at all vertices, combined vs.
  separate interpolation of three components, phi spacing of the query points
- new SpaceCharge test: closure of the global distortions (standard and fast
  method) with the global corrections from analytical fields on both sides,
  large fields and queries at invalid positions, reading of a stored object
  with a non-default grid (float and double)
- createResidualDistortionObject.C is tested compiled only (COMPILE_ONLY, as the
  other space-charge macros and as documented in the macro): loading it in
  the interpreter crashed while auto-loading the CommonUtils dictionary
  (DebugStreamer's map of unique_ptr<TreeStreamRedirector> with an incomplete
  type), independent of this change

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

Author: matthias-kleiner

Detectors/TPC/spacecharge/CMakeLists.txt

@@ -17,7 +17,6 @@ o2_add_library(TPCSpaceCharge
                        src/DataContainer3D.cxx
                PUBLIC_LINK_LIBRARIES O2::TPCBaseRecSim
                                      O2::Field
-                                     Vc::Vc
                                      ROOT::Core
                                      O2::DataFormatsParameters)
 
@@ -35,7 +34,6 @@ o2_target_root_dictionary(TPCSpaceCharge
                                   include/TPCSpaceCharge/PoissonSolverHelpers.h
                                   include/TPCSpaceCharge/SpaceChargeHelpers.h
                                   include/TPCSpaceCharge/TriCubic.h
-                                  include/TPCSpaceCharge/Vector.h
                                   include/TPCSpaceCharge/Vector3D.h
                                   include/TPCSpaceCharge/SpaceChargeParameter.h
                           LINKDEF src/TPCSpacechargeLinkDef.h)
@@ -49,7 +47,7 @@ o2_add_test_root_macro(macro/createResidualDistortionObject.C
                       PUBLIC_LINK_LIBRARIES O2::TPCSpaceCharge
                                             O2::TPCSpaceChargeIO
                                             O2::CommonUtils
-                      LABELS tpc)
+                      LABELS tpc COMPILE_ONLY)
 
 
 install(FILES macro/createSCHistosFromHits.C
@@ -81,6 +79,14 @@ o2_add_test(TriCubic
             LABELS tpc
             CONFIGURATIONS RelWithDebInfo Release MinSizeRel)
 
+o2_add_test(SpaceCharge
+            COMPONENT_NAME spacecharge
+            PUBLIC_LINK_LIBRARIES O2::TPCSpaceCharge
+            SOURCES test/testO2TPCSpaceCharge.cxx
+            ENVIRONMENT O2_ROOT=${CMAKE_BINARY_DIR}/stage
+            LABELS tpc
+            CONFIGURATIONS RelWithDebInfo Release MinSizeRel)
+
 if (OpenMP_CXX_FOUND)
     target_compile_definitions(${targetName} PRIVATE WITH_OPENMP)
     target_link_libraries(${targetName} PRIVATE OpenMP::OpenMP_CXX)

Detectors/TPC/spacecharge/include/TPCSpaceCharge/DataContainer3D.h

@@ -20,6 +20,8 @@
 
 #include "Rtypes.h"
 #include <vector>
+#include <string_view>
+#include <utility>
 
 class TFile;
 class TTree;

Detectors/TPC/spacecharge/include/TPCSpaceCharge/PoissonSolver.h

@@ -26,37 +26,14 @@ namespace o2
 namespace tpc
 {
 
-///< Enumeration of Cycles Type
-enum class CycleType {
-  VCycle = 0, ///< V Cycle
-  WCycle = 1, ///< W Cycle (TODO)
-  FCycle = 2  ///< Full Cycle
-};
-
-///< Fine -> Coarse Grid transfer operator types
-enum class GridTransferType {
-  Half = 0, ///< Half weighting
-  Full = 1, ///< Full weighting
-};
-
-///< Smoothing (Relax) operator types
-enum class RelaxType {
-  Jacobi = 0,         ///< Jacobi (5 Stencil 2D, 7 Stencil 3D_
-  WeightedJacobi = 1, ///< (TODO)
-  GaussSeidel = 2     ///< Gauss Seidel 2D (2 Color, 5 Stencil), 3D (7 Stencil)
-};
-
-struct MGParameters {                                             ///< Parameters choice for MultiGrid algorithm
-  inline static bool isFull3D = true;                             ///<  TRUE: full coarsening, FALSE: semi coarsening
-  inline static CycleType cycleType = CycleType::FCycle;          ///< cycleType follow  CycleType
-  inline static GridTransferType gtType = GridTransferType::Full; ///< gtType grid transfer type follow GridTransferType
-  inline static RelaxType relaxType = RelaxType::GaussSeidel;     ///< relaxType follow RelaxType
-  inline static int nPre = 2;                                     ///< number of iteration for pre smoothing
-  inline static int nPost = 2;                                    ///< number of iteration for post smoothing
-  inline static int nMGCycle = 200;                               ///< number of multi grid cycle (V type)
-  inline static int maxLoop = 7;                                  ///< the number of tree-deep of multi grid
-  inline static int gamma = 1;                                    ///< number of iteration at coarsest level !TODO SET TO REASONABLE VALUE!
-  inline static int normalizeGridToNSector = SECTORSPERSIDE;      ///< the grid in phi direction is squashed from 2 Pi to (2 Pi / SECTORSPERSIDE). This can used to get the potential for phi symmetric sc density or boundary potentials
+struct MGParameters {                                        ///< Parameters of the multi grid algorithm
+  inline static int nPre = 2;                                ///< number of iteration for pre smoothing
+  inline static int nPost = 2;                               ///< number of iteration for post smoothing
+  inline static int nMGCycle = 200;                          ///< maximum number of V cycles on each level of the full multi grid
+  inline static int normalizeGridToNSector = SECTORSPERSIDE; ///< the grid in phi direction is squashed from 2 Pi to (2 Pi / SECTORSPERSIDE). This can used to get the potential for phi symmetric sc density or boundary potentials
+  inline static float phiCoarseningThreshold = 0.5;          ///< full 3D: the phi direction is only coarsened if the relative coupling in phi (h_{r,z} / (r_{min} * h_{phi}))^2 is larger than this value (weakly coupled directions should not be coarsened). 0: always coarsen phi if possible
+  inline static float phiLineRelaxThreshold = 4;             ///< full 3D: line relaxation in phi instead of point relaxation on the multi grid levels where the relative coupling in phi (h_{r,z} / (r_{min} * h_{phi}))^2 is larger than this value. 0: always point relaxation
+  inline static bool relativeConvergence = true;             ///< stopping criterion of the multi grid V cycles: true: relative change of the potential between two V cycles sqrt(sum dV^2 / sum V^2) (independent of the units of the potential and the grid size), false: maximum over the phi slices of the sum of dV^2 (depends on the units and on the grid size)
 };
 
 template <typename DataT = double>

Detectors/TPC/spacecharge/include/TPCSpaceCharge/PoissonSolverHelpers.h

@@ -17,9 +17,10 @@
 #ifndef ALICEO2_TPC_REGULARGRID3D_H_
 #define ALICEO2_TPC_REGULARGRID3D_H_
 
-#include "TPCSpaceCharge/Vector.h"
 #include "Rtypes.h" // for ClassDefNV
 #include "TPCSpaceCharge/SpaceChargeParameter.h"
+#include <array>
+#include <cmath>
 #include <vector>
 
 namespace o2
@@ -90,17 +91,17 @@ struct RegularGrid3D {
   static constexpr unsigned int getFR() { return FR; }     /// \return returns the index for dimension y (1)
   static constexpr unsigned int getFPhi() { return FPHI; } /// \return returns the index for dimension z (2)
 
-  const Vector<DataT, 3>& getGridMin() const { return mMin; } /// \return returns the minimum coordinates of the grid in all dimensions
-  DataT getGridMinZ() const { return mMin[FZ]; }              /// \return returns the minimum coordinate of the grid in x dimension
-  DataT getGridMinR() const { return mMin[FR]; }              /// \return returns the minimum coordinate of the grid in y dimension
-  DataT getGridMinPhi() const { return mMin[FPHI]; }          /// \return returns the minimum coordinate of the grid in z dimension
+  const std::array<DataT, 3>& getGridMin() const { return mMin; } /// \return returns the minimum coordinates of the grid in all dimensions
+  DataT getGridMinZ() const { return mMin[FZ]; }                  /// \return returns the minimum coordinate of the grid in x dimension
+  DataT getGridMinR() const { return mMin[FR]; }                  /// \return returns the minimum coordinate of the grid in y dimension
+  DataT getGridMinPhi() const { return mMin[FPHI]; }              /// \return returns the minimum coordinate of the grid in z dimension
 
   DataT getGridMaxZ() const { return mMax[FZ]; }
   DataT getGridMaxR() const { return mMax[FR]; }
   DataT getGridMaxPhi() const { return mMax[FPHI]; }
 
   /// \return returns the inversed spacing of the grid for all dimensions
-  const Vector<DataT, 3>& getInvSpacing() const { return mInvSpacing; }
+  const std::array<DataT, 3>& getInvSpacing() const { return mInvSpacing; }
   DataT getInvSpacingZ() const { return mInvSpacing[FZ]; }
   DataT getInvSpacingR() const { return mInvSpacing[FR]; }
   DataT getInvSpacingPhi() const { return mInvSpacing[FPHI]; }
@@ -141,7 +142,7 @@ struct RegularGrid3D {
   /// \return returns the z positon for given vertex
   DataT getPhiVertex(const size_t vertexZ) const { return mPhiVertices[vertexZ]; }
 
-  const Vector<DataT, 3>& getMaxIndices() const { return sMaxIndex; } /// get max indices for all dimensions
+  const std::array<DataT, 3>& getMaxIndices() const { return sMaxIndex; } /// get max indices for all dimensions
 
   DataT getMaxIndexZ() const { return sMaxIndex[0]; }   /// get max index in x direction
   DataT getMaxIndexR() const { return sMaxIndex[1]; }   /// get max index in y direction
@@ -159,19 +160,19 @@ struct RegularGrid3D {
 
  private:
   ParamSpaceCharge mParamGrid{};
-  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
-  Vector<DataT, FDIM> mMin{};                                                                                                                                                      ///<! min vertices positions of the grid
-  Vector<DataT, FDIM> mMax{};                                                                                                                                                      ///<! max vertices positions of the grid
-  Vector<DataT, FDIM> mSpacing{};                                                                                                                                                  ///<! spacing of the grid
-  Vector<DataT, FDIM> mInvSpacing{};                                                                                                                                               ///<! inverse spacing of grid
-  Vector<DataT, FDIM> sMaxIndex{{static_cast<DataT>(mParamGrid.NZVertices - 1.), static_cast<DataT>(mParamGrid.NRVertices - 1), static_cast<DataT>(mParamGrid.NPhiVertices - 1)}}; ///<! max index which is on the grid in all dimensions
-  Vector<int, FDIM> sNdim{{static_cast<int>(mParamGrid.NZVertices), static_cast<int>(mParamGrid.NRVertices), static_cast<int>(mParamGrid.NPhiVertices)}};                          ///<! number of vertices for each dimension
-  std::vector<DataT> mZVertices{};                                                                                                                                                 ///< positions of vertices in x direction
-  std::vector<DataT> mRVertices{};                                                                                                                                                 ///< positions of vertices in y direction
-  std::vector<DataT> 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
+  std::array<DataT, FDIM> mMin{};                                                                                                                                                      ///<! min vertices positions of the grid
+  std::array<DataT, FDIM> mMax{};                                                                                                                                                      ///<! max vertices positions of the grid
+  std::array<DataT, FDIM> mSpacing{};                                                                                                                                                  ///<! spacing of the grid
+  std::array<DataT, FDIM> mInvSpacing{};                                                                                                                                               ///<! inverse spacing of grid
+  std::array<DataT, FDIM> sMaxIndex{{static_cast<DataT>(mParamGrid.NZVertices - 1.), static_cast<DataT>(mParamGrid.NRVertices - 1), static_cast<DataT>(mParamGrid.NPhiVertices - 1)}}; ///<! max index which is on the grid in all dimensions
+  std::array<int, FDIM> sNdim{{static_cast<int>(mParamGrid.NZVertices), static_cast<int>(mParamGrid.NRVertices), static_cast<int>(mParamGrid.NPhiVertices)}};                          ///<! number of vertices for each dimension
+  std::vector<DataT> mZVertices{};                                                                                                                                                     ///< positions of vertices in x direction
+  std::vector<DataT> mRVertices{};                                                                                                                                                     ///< positions of vertices in y direction
+  std::vector<DataT> mPhiVertices{};                                                                                                                                                   ///< positions of vertices in z direction
 
   void initLists();
 
@@ -217,6 +218,11 @@ DataT RegularGrid3D<DataT>::clampToGridRel(const DataT pos, const unsigned int d
 template <typename DataT>
 DataT RegularGrid3D<DataT>::clampToGridCircular(DataT pos, const unsigned int dim) const
 {
+  const DataT period = mMax[dim] - mMin[dim] + mSpacing[dim];
+  if (!(std::abs(pos - mMin[dim]) < 4 * period)) {
+    // far outside of the grid or not finite: wrap directly with the exact remainder instead of looping (a non-finite position becomes NaN)
+    pos = mMin[dim] + std::fmod(pos - mMin[dim], period);
+  }
   while (pos < mMin[dim]) {
     pos += mMax[dim] - mMin[dim] + mSpacing[dim];
   }
@@ -229,6 +235,11 @@ DataT RegularGrid3D<DataT>::clampToGridCircular(DataT pos, const unsigned int di
 template <typename DataT>
 DataT RegularGrid3D<DataT>::clampToGridCircularRel(DataT pos, const unsigned int dim) const
 {
+  const DataT period = sNdim[dim];
+  if (!(std::abs(pos) < 4 * period)) {
+    // far outside of the grid or not finite: wrap directly with the exact remainder instead of looping (a non-finite position becomes NaN)
+    pos = std::fmod(pos, period);
+  }
   while (pos < 0) {
     pos += sNdim[dim];
   }