Skip to content

Make the magnetic field evaluation reentrant #16057

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

Merged
sawenzel merged 1 commit into from
Nov 7, 2026
Merged

Make the magnetic field evaluation reentrant #16057

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
4 changes: 2 additions & 2 deletions Common/Field/src/MagneticWrapperChebyshev.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -690,7 +690,7 @@ Int_t MagneticWrapperChebyshev::findTPCRatSegment(const Double_t* rpz) const

void MagneticWrapperChebyshev::getTPCIntegral(const Double_t* xyz, Double_t* b) const
{
static Double_t rphiz[3];
Double_t rphiz[3];

// TPCInt region
// convert coordinates to cyl system
Expand All @@ -712,7 +712,7 @@ void MagneticWrapperChebyshev::getTPCIntegral(const Double_t* xyz, Double_t* b)

void MagneticWrapperChebyshev::getTPCRatIntegral(const Double_t* xyz, Double_t* b) const
{
static Double_t rphiz[3];
Double_t rphiz[3];

// TPCRatIntegral region
// convert coordinates to cylindrical system
Expand Down
41 changes: 41 additions & 0 deletions Common/Field/test/testMagneticField.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,10 @@
#include <TFile.h>
#include <TStopwatch.h>
#include <TRandom.h>
#include <algorithm>
#include <array>
#include <thread>
#include <vector>

using namespace o2::field;

Expand Down Expand Up @@ -144,3 +148,40 @@ BOOST_AUTO_TEST_CASE(MagneticField_reinitialization_test)
}
}
}

BOOST_AUTO_TEST_CASE(MagneticField_threads)
{
// evaluating one field from several threads at once gives bit-identical values to a sequential pass
MagneticField fld("Maps", "Maps", 1., 1., o2::field::MagFieldParam::k5kG);
constexpr int npoints = 20000;
constexpr int nthreads = 8;
std::vector<std::array<double, 3>> xyz(npoints);
for (auto& p : xyz) {
p = {gRandom->Uniform(-450., 450.), gRandom->Uniform(-450., 450.), gRandom->Uniform(-1500., 1500.)};
}
// field, Bz, TPC integral and TPC ratio integral for every point
constexpr int nval = 10;
auto evaluate = [&fld, &xyz](std::vector<double>& out) {
out.resize(npoints * nval);
for (int i = 0; i < npoints; ++i) {
fld.Field(xyz[i].data(), &out[i * nval]);
out[i * nval + 3] = fld.getBz(xyz[i].data());
fld.getTPCIntegral(xyz[i].data(), &out[i * nval + 4]);
fld.getTPCRatIntegral(xyz[i].data(), &out[i * nval + 7]);
}
};
std::vector<double> reference;
evaluate(reference);
std::vector<std::vector<double>> results(nthreads);
std::vector<std::thread> threads;
for (int t = 0; t < nthreads; ++t) {
threads.emplace_back(evaluate, std::ref(results[t]));
}
for (auto& th : threads) {
th.join();
}
for (int t = 0; t < nthreads; ++t) {
auto diff = std::mismatch(reference.begin(), reference.end(), results[t].begin());
BOOST_CHECK_MESSAGE(diff.first == reference.end(), "thread " << t << " differs at value " << (diff.first - reference.begin()));
}
}
50 changes: 30 additions & 20 deletions Common/MathUtils/include/MathUtils/Chebyshev3D.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -292,66 +292,72 @@ inline Bool_t Chebyshev3D::isInside(const Double_t* par) const
/// Evaluates Chebyshev parameterization for 3d->DimOut function
inline void Chebyshev3D::Eval(const Float_t* par, Float_t* res)
{
Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant
for (int i = 3; i--;) {
mTemporaryCoefficient[i] = mapToInternal(par[i], i);
coef[i] = mapToInternal(par[i], i);
}
for (int i = mOutputArrayDimension; i--;) {
res[i] = getChebyshevCalc(i)->Eval(mTemporaryCoefficient);
res[i] = getChebyshevCalc(i)->Eval(coef);
}
}

/// Evaluates Chebyshev parameterization for 3d->DimOut function
inline void Chebyshev3D::Eval(const Double_t* par, Double_t* res)
{
Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant
for (int i = 3; i--;) {
mTemporaryCoefficient[i] = mapToInternal(par[i], i);
coef[i] = mapToInternal(par[i], i);
}
for (int i = mOutputArrayDimension; i--;) {
res[i] = getChebyshevCalc(i)->Eval(mTemporaryCoefficient);
res[i] = getChebyshevCalc(i)->Eval(coef);
}
}

/// Evaluates Chebyshev parameterization for idim-th output dimension of 3d->DimOut function
inline Double_t Chebyshev3D::Eval(const Double_t* par, int idim)
{
Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant
for (int i = 3; i--;) {
mTemporaryCoefficient[i] = mapToInternal(par[i], i);
coef[i] = mapToInternal(par[i], i);
}
return getChebyshevCalc(idim)->Eval(mTemporaryCoefficient);
return getChebyshevCalc(idim)->Eval(coef);
}

/// Evaluates Chebyshev parameterization for idim-th output dimension of 3d->DimOut function
inline Float_t Chebyshev3D::Eval(const Float_t* par, int idim)
{
Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant
for (int i = 3; i--;) {
mTemporaryCoefficient[i] = mapToInternal(par[i], i);
coef[i] = mapToInternal(par[i], i);
}
return getChebyshevCalc(idim)->Eval(mTemporaryCoefficient);
return getChebyshevCalc(idim)->Eval(coef);
}

/// Returns the gradient matrix
inline void Chebyshev3D::evaluateDerivative3D(const Float_t* par, Float_t dbdr[3][3])
{
Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant
for (int i = 3; i--;) {
mTemporaryCoefficient[i] = mapToInternal(par[i], i);
coef[i] = mapToInternal(par[i], i);
}
for (int ib = 3; ib--;) {
for (int id = 3; id--;) {
dbdr[ib][id] = getChebyshevCalc(ib)->evaluateDerivative(id, mTemporaryCoefficient) * mBoundaryMappingScale[id];
dbdr[ib][id] = getChebyshevCalc(ib)->evaluateDerivative(id, coef) * mBoundaryMappingScale[id];
}
}
}

/// Returns the gradient matrix
inline void Chebyshev3D::evaluateDerivative3D2(const Float_t* par, Float_t dbdrdr[3][3][3])
{
Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant
for (int i = 3; i--;) {
mTemporaryCoefficient[i] = mapToInternal(par[i], i);
coef[i] = mapToInternal(par[i], i);
}
for (int ib = 3; ib--;) {
for (int id = 3; id--;) {
for (int id1 = 3; id1--;) {
dbdrdr[ib][id][id1] = getChebyshevCalc(ib)->evaluateDerivative2(id, id1, mTemporaryCoefficient) *
dbdrdr[ib][id][id1] = getChebyshevCalc(ib)->evaluateDerivative2(id, id1, coef) *
mBoundaryMappingScale[id] * mBoundaryMappingScale[id1];
}
}
Expand All @@ -361,22 +367,24 @@ inline void Chebyshev3D::evaluateDerivative3D2(const Float_t* par, Float_t dbdrd
// Evaluates Chebyshev parameterization derivative for 3d->DimOut function
inline void Chebyshev3D::evaluateDerivative(int dimd, const Float_t* par, Float_t* res)
{
Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant
for (int i = 3; i--;) {
mTemporaryCoefficient[i] = mapToInternal(par[i], i);
coef[i] = mapToInternal(par[i], i);
}
for (int i = mOutputArrayDimension; i--;) {
res[i] = getChebyshevCalc(i)->evaluateDerivative(dimd, mTemporaryCoefficient) * mBoundaryMappingScale[dimd];
res[i] = getChebyshevCalc(i)->evaluateDerivative(dimd, coef) * mBoundaryMappingScale[dimd];
};
}

// Evaluates Chebyshev parameterization 2nd derivative over dimd1 and dimd2 dimensions for 3d->DimOut function
inline void Chebyshev3D::evaluateDerivative2(int dimd1, int dimd2, const Float_t* par, Float_t* res)
{
Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant
for (int i = 3; i--;) {
mTemporaryCoefficient[i] = mapToInternal(par[i], i);
coef[i] = mapToInternal(par[i], i);
}
for (int i = mOutputArrayDimension; i--;) {
res[i] = getChebyshevCalc(i)->evaluateDerivative2(dimd1, dimd2, mTemporaryCoefficient) *
res[i] = getChebyshevCalc(i)->evaluateDerivative2(dimd1, dimd2, coef) *
mBoundaryMappingScale[dimd1] * mBoundaryMappingScale[dimd2];
}
}
Expand All @@ -385,20 +393,22 @@ inline void Chebyshev3D::evaluateDerivative2(int dimd1, int dimd2, const Float_t
/// function
inline Float_t Chebyshev3D::evaluateDerivative(int dimd, const Float_t* par, int idim)
{
Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant
for (int i = 3; i--;) {
mTemporaryCoefficient[i] = mapToInternal(par[i], i);
coef[i] = mapToInternal(par[i], i);
}
return getChebyshevCalc(idim)->evaluateDerivative(dimd, mTemporaryCoefficient) * mBoundaryMappingScale[dimd];
return getChebyshevCalc(idim)->evaluateDerivative(dimd, coef) * mBoundaryMappingScale[dimd];
}

/// Evaluates Chebyshev parameterization 2ns derivative over dimd1 and dimd2 dimensions for idim-th output dimension of
/// 3d->DimOut function
inline Float_t Chebyshev3D::evaluateDerivative2(int dimd1, int dimd2, const Float_t* par, int idim)
{
Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant
for (int i = 3; i--;) {
mTemporaryCoefficient[i] = mapToInternal(par[i], i);
coef[i] = mapToInternal(par[i], i);
}
return getChebyshevCalc(idim)->evaluateDerivative2(dimd1, dimd2, mTemporaryCoefficient) *
return getChebyshevCalc(idim)->evaluateDerivative2(dimd1, dimd2, coef) *
mBoundaryMappingScale[dimd1] * mBoundaryMappingScale[dimd2];
}

Expand Down
55 changes: 47 additions & 8 deletions Common/MathUtils/include/MathUtils/Chebyshev3DCalc.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,8 @@
#include <TNamed.h> // for TNamed
#include <cstdio> // for FILE, stdout
#include <cmath> // for std::fma
#include <algorithm>
#include <vector>
#include "Rtypes.h" // for Float_t, UShort_t, Int_t, Double_t, etc

class TString;
Expand Down Expand Up @@ -171,6 +173,37 @@ class Chebyshev3DCalc : public TNamed
Double_t Eval(const Double_t* par) const;

private:
/// Scratch arrays for one evaluation: on the stack for the usual sizes, thread-local beyond them,
/// so that evaluation is reentrant. The ALICE field maps use at most 14 rows and columns.
class EvalScratch
{
public:
EvalScratch(int nrows, int ncols)
{
if (nrows <= kStackSize && ncols <= kStackSize) {
mRows = mStackRows;
mCols = mStackCols;
} else {
thread_local std::vector<Float_t> rows, cols;
rows.resize(std::max<size_t>(nrows, rows.size()));
cols.resize(std::max<size_t>(ncols, cols.size()));
mRows = rows.data();
mCols = cols.data();
}
}
EvalScratch(const EvalScratch&) = delete;
EvalScratch& operator=(const EvalScratch&) = delete;
Float_t* rows() { return mRows; }
Float_t* cols() { return mCols; }

private:
static constexpr int kStackSize = 64;
Float_t mStackRows[kStackSize];
Float_t mStackCols[kStackSize];
Float_t* mRows;
Float_t* mCols;
};

Int_t mNumberOfCoefficients; ///< total number of coeeficients
Int_t mNumberOfRows; ///< number of significant rows in the 3D coeffs matrix
Int_t mNumberOfColumns; ///< max number of significant cols in the 3D coeffs matrix
Expand All @@ -188,8 +221,8 @@ class Chebyshev3DCalc : public TNamed
// coeffs for col/row
Float_t* mCoefficients; //[mNumberOfCoefficients] array of Chebyshev coefficients

Float_t* mTemporaryCoefficients2D; //[mNumberOfColumns] temp. coeffs for 2d summation
Float_t* mTemporaryCoefficients1D; //[mNumberOfRows] temp. coeffs for 1d summation
Float_t* mTemporaryCoefficients2D; //[mNumberOfColumns] no longer used for evaluation
Float_t* mTemporaryCoefficients1D; //[mNumberOfRows] no longer used for evaluation

ClassDefOverride(o2::math_utils::Chebyshev3DCalc,
2) // Class for interpolation of 3D->1 function by Chebyshev parametrization
Expand Down Expand Up @@ -223,32 +256,38 @@ inline Float_t Chebyshev3DCalc::chebyshevEvaluation1D(Float_t x, const Float_t*
/// VERY IMPORTANT: par must contain the function arguments ALREADY MAPPED to [-1:1] interval
inline Float_t Chebyshev3DCalc::Eval(const Float_t* par) const
{
EvalScratch scratch(mNumberOfRows, mNumberOfColumns);
Float_t* coefs2D = scratch.cols();
Float_t* coefs1D = scratch.rows();
for (int id0 = mNumberOfRows; id0--;) {
int nCLoc = mNumberOfColumnsAtRow[id0]; // number of significant coefs on this row
int col0 = mColumnAtRowBeginning[id0]; // beginning of local column in the 2D boundary matrix
for (int id1 = nCLoc; id1--;) {
int id = id1 + col0;
mTemporaryCoefficients2D[id1] = chebyshevEvaluation1D(par[2], mCoefficients + mCoefficientBound2D1[id], mCoefficientBound2D0[id]);
coefs2D[id1] = chebyshevEvaluation1D(par[2], mCoefficients + mCoefficientBound2D1[id], mCoefficientBound2D0[id]);
}
mTemporaryCoefficients1D[id0] = chebyshevEvaluation1D(par[1], mTemporaryCoefficients2D, nCLoc);
coefs1D[id0] = chebyshevEvaluation1D(par[1], coefs2D, nCLoc);
}
return chebyshevEvaluation1D(par[0], mTemporaryCoefficients1D, mNumberOfRows);
return chebyshevEvaluation1D(par[0], coefs1D, mNumberOfRows);
}

/// Evaluates Chebyshev parameterization for 3D function.
/// VERY IMPORTANT: par must contain the function arguments ALREADY MAPPED to [-1:1] interval
inline Double_t Chebyshev3DCalc::Eval(const Double_t* par) const
{
EvalScratch scratch(mNumberOfRows, mNumberOfColumns);
Float_t* coefs2D = scratch.cols();
Float_t* coefs1D = scratch.rows();
for (int id0 = mNumberOfRows; id0--;) {
int nCLoc = mNumberOfColumnsAtRow[id0]; // number of significant coefs on this row
int col0 = mColumnAtRowBeginning[id0]; // beginning of local column in the 2D boundary matrix
for (int id1 = nCLoc; id1--;) {
int id = id1 + col0;
mTemporaryCoefficients2D[id1] = chebyshevEvaluation1D(par[2], mCoefficients + mCoefficientBound2D1[id], mCoefficientBound2D0[id]);
coefs2D[id1] = chebyshevEvaluation1D(par[2], mCoefficients + mCoefficientBound2D1[id], mCoefficientBound2D0[id]);
}
mTemporaryCoefficients1D[id0] = chebyshevEvaluation1D(par[1], mTemporaryCoefficients2D, nCLoc);
coefs1D[id0] = chebyshevEvaluation1D(par[1], coefs2D, nCLoc);
}
return chebyshevEvaluation1D(par[0], mTemporaryCoefficients1D, mNumberOfRows);
return chebyshevEvaluation1D(par[0], coefs1D, mNumberOfRows);
}
} // namespace math_utils
} // namespace o2
Expand Down
Loading
Loading