Skip to content

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

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
matthias-kleiner wants to merge 1 commit into
base: dev
Choose a base branch
from
+1,546 −2,982
Open

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

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 9 additions & 3 deletions Detectors/TPC/spacecharge/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,6 @@ o2_add_library(TPCSpaceCharge
src/DataContainer3D.cxx
PUBLIC_LINK_LIBRARIES O2::TPCBaseRecSim
O2::Field
Vc::Vc
ROOT::Core
O2::DataFormatsParameters)

Expand All @@ -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)
Expand All @@ -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
Expand Down Expand Up @@ -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)
Expand Down
2 changes: 2 additions & 0 deletions Detectors/TPC/spacecharge/include/TPCSpaceCharge/DataContainer3D.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,8 @@

#include "Rtypes.h"
#include <vector>
#include <string_view>
#include <utility>

class TFile;
class TTree;
Expand Down
335 changes: 88 additions & 247 deletions Detectors/TPC/spacecharge/include/TPCSpaceCharge/PoissonSolver.h
View file
Open in desktop

Large diffs are not rendered by default.

39 changes: 8 additions & 31 deletions Detectors/TPC/spacecharge/include/TPCSpaceCharge/PoissonSolverHelpers.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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>
Expand Down
51 changes: 31 additions & 20 deletions Detectors/TPC/spacecharge/include/TPCSpaceCharge/RegularGrid3D.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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]; }
Expand Down Expand Up @@ -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
Expand All @@ -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();

Expand Down Expand Up @@ -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];
}
Expand All @@ -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];
}
Expand Down
Loading
Loading