Skip to content

Commit 36c495e

Browse files
cmsbuildcmsbuild
authored
Merge pull request #44360 from lecriste/para-magnetic-field-alpaka_1410
Parametrized magnetic field in alpaka
2 parents bc22d28 + 0828b85 commit 36c495e

File tree

3 files changed

+187
-0
lines changed

3 files changed

+187
-0
lines changed

‎MagneticField/ParametrizedEngine/interface/ParabolicParametrizedMagneticField.h‎

Lines changed: 47 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,47 @@
1+
/**
2+
Description: Utility function to calculate the Magnetic Field on the GPU. The Vec3 argument of the functions must support access to its components via (), note that e.g. Eigen::Matrix provides such an interface.
3+
*/
4+
5+
#ifndef MagneticField_ParametrizedEngine_interface_ParabolicParametrizedMagneticField_h
6+
#define MagneticField_ParametrizedEngine_interface_ParabolicParametrizedMagneticField_h
7+
8+
namespace magneticFieldParabolicPortable {
9+
10+
struct Parameters {
11+
// These parameters are the best fit of 3.8T to the OAEParametrizedMagneticField parametrization.
12+
// See MagneticField/ParametrizedEngine/src/ParabolicParametrizedMagneticField.cc
13+
static constexpr float c1 = 3.8114;
14+
static constexpr float b0 = -3.94991e-06;
15+
static constexpr float b1 = 7.53701e-06;
16+
static constexpr float a = 2.43878e-11;
17+
static constexpr float max_radius2 = 13225.f; // tracker radius
18+
static constexpr float max_z = 280.f; // tracker z
19+
};
20+
21+
template <typename Vec3>
22+
constexpr float Kr(Vec3 const& vec) {
23+
return Parameters::a * (vec(0) * vec(0) + vec(1) * vec(1)) + 1.f;
24+
}
25+
26+
template <typename Vec3>
27+
constexpr float B0Z(Vec3 const& vec) {
28+
return Parameters::b0 * vec(2) * vec(2) + Parameters::b1 * vec(2) + Parameters::c1;
29+
}
30+
31+
template <typename Vec3>
32+
constexpr bool isValid(Vec3 const& vec) {
33+
return ((vec(0) * vec(0) + vec(1) * vec(1)) < Parameters::max_radius2 && fabs(vec(2)) < Parameters::max_z);
34+
}
35+
36+
template <typename Vec3>
37+
constexpr float magneticFieldAtPoint(Vec3 const& vec) {
38+
if (isValid(vec)) {
39+
return B0Z(vec) * Kr(vec);
40+
} else {
41+
return 0;
42+
}
43+
}
44+
45+
} // namespace magneticFieldParabolicPortable
46+
47+
#endif

‎MagneticField/ParametrizedEngine/test/BuildFile.xml‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,9 @@
1+
<bin file="alpaka/testParabolicParametrizedMagneticField.dev.cc">
2+
<use name="alpaka"/>
3+
<use name="eigen"/>
4+
<use name="DataFormats/GeometryVector"/>
5+
<use name="FWCore/Utilities"/>
6+
<use name="HeterogeneousCore/AlpakaInterface"/>
7+
<use name="MagneticField/ParametrizedEngine" source_only="1"/>
8+
<flags ALPAKA_BACKENDS="1"/>
9+
</bin>

‎MagneticField/ParametrizedEngine/test/alpaka/testParabolicParametrizedMagneticField.dev.cc‎

Lines changed: 131 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,131 @@
1+
#include <iostream>
2+
#include <fstream>
3+
#include <Eigen/Core>
4+
#include <alpaka/alpaka.hpp>
5+
6+
#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
7+
#include "DataFormats/GeometryVector/interface/GlobalVector.h"
8+
#include "FWCore/Utilities/interface/FileInPath.h"
9+
#include "FWCore/Utilities/interface/Exception.h"
10+
#include "FWCore/Utilities/interface/stringize.h"
11+
#include "HeterogeneousCore/AlpakaInterface/interface/config.h"
12+
#include "HeterogeneousCore/AlpakaInterface/interface/memory.h"
13+
#include "HeterogeneousCore/AlpakaInterface/interface/workdivision.h"
14+
#include "MagneticField/ParametrizedEngine/interface/ParabolicParametrizedMagneticField.h"
15+
16+
using namespace edm;
17+
using namespace std;
18+
using namespace ALPAKA_ACCELERATOR_NAMESPACE;
19+
using namespace magneticFieldParabolicPortable;
20+
using Vector3f = Eigen::Matrix<float, 3, 1>;
21+
22+
struct MagneticFieldKernel {
23+
template <typename TAcc, typename T>
24+
ALPAKA_FN_ACC void operator()(TAcc const& acc, T const* __restrict__ in, T* __restrict__ out, size_t size) const {
25+
for (auto index : cms::alpakatools::uniform_elements(acc, size)) {
26+
out[index](0) = 0;
27+
out[index](1) = 0;
28+
out[index](2) = magneticFieldAtPoint(in[index]);
29+
}
30+
}
31+
};
32+
33+
int main() {
34+
// get the list of devices on the current platform
35+
auto const& devices = cms::alpakatools::devices<Platform>();
36+
if (devices.empty()) {
37+
std::cerr << "No devices available for the " EDM_STRINGIZE(ALPAKA_ACCELERATOR_NAMESPACE) " backend, "
38+
"the test will be skipped.\n";
39+
exit(EXIT_FAILURE);
40+
}
41+
42+
ifstream file;
43+
edm::FileInPath mydata("MagneticField/Engine/data/Regression/referenceField_160812_RII_3_8T.bin");
44+
file.open(mydata.fullPath().c_str(), ios::binary);
45+
46+
int count = 0;
47+
float px, py, pz;
48+
float bx, by, bz;
49+
vector<Vector3f> points;
50+
vector<GlobalVector> referenceB_vec;
51+
52+
int numberOfPoints = 100;
53+
points.reserve(numberOfPoints);
54+
referenceB_vec.reserve(numberOfPoints);
55+
do {
56+
if (!(file.read((char*)&px, sizeof(float)) && file.read((char*)&py, sizeof(float)) &&
57+
file.read((char*)&pz, sizeof(float)) && file.read((char*)&bx, sizeof(float)) &&
58+
file.read((char*)&by, sizeof(float)) && file.read((char*)&bz, sizeof(float))))
59+
break;
60+
61+
const auto point = Vector3f(px, py, pz);
62+
if (!isValid(point))
63+
continue;
64+
65+
points.push_back(Vector3f(px, py, pz));
66+
referenceB_vec.push_back(GlobalVector(bx, by, bz));
67+
count++;
68+
} while (count < numberOfPoints);
69+
70+
const size_t size = points.size();
71+
// allocate the input and output host buffer in pinned memory accessible by the Platform devices
72+
auto points_host = cms::alpakatools::make_host_buffer<Vector3f[], Platform>(size);
73+
auto field_host = cms::alpakatools::make_host_buffer<Vector3f[], Platform>(size);
74+
// fill the input buffers, and the output buffer with zeros
75+
for (size_t i = 0; i < size; ++i) {
76+
points_host[i] = points[i];
77+
field_host[i] = Vector3f::Zero();
78+
}
79+
80+
float resolution = 0.2;
81+
float maxdelta = 0.;
82+
int fail = 0;
83+
84+
// run the test on each device
85+
for (auto const& device : devices) {
86+
auto queue = Queue(device);
87+
// allocate input and output buffers on the device
88+
auto points_dev = cms::alpakatools::make_device_buffer<Vector3f[]>(queue, size);
89+
auto field_dev = cms::alpakatools::make_device_buffer<Vector3f[]>(queue, size);
90+
91+
// copy the input data to the device; the size is known from the buffer objects
92+
alpaka::memcpy(queue, points_dev, points_host);
93+
// fill the output buffer with zeros; the size is known from the buffer objects
94+
alpaka::memset(queue, field_dev, 0.);
95+
96+
auto workDiv = cms::alpakatools::make_workdiv<Acc1D>(1, size);
97+
alpaka::exec<Acc1D>(queue, workDiv, MagneticFieldKernel{}, points_dev.data(), field_dev.data(), size);
98+
99+
// copy the results from the device to the host
100+
alpaka::memcpy(queue, field_host, field_dev);
101+
102+
// wait for the kernel and the potential copy to complete
103+
alpaka::wait(queue);
104+
105+
// check the results
106+
for (uint i = 0; i < points.size(); i++) {
107+
const auto& point = points[i];
108+
const auto& referenceB = referenceB_vec[i];
109+
GlobalVector parametricB(field_host[i](0), field_host[i](1), field_host[i](2));
110+
float delta = (referenceB - parametricB).mag();
111+
if (delta > resolution) {
112+
++fail;
113+
if (delta > maxdelta)
114+
maxdelta = delta;
115+
if (fail < 10) {
116+
const GlobalPoint gp(point(0), point(1), point(2));
117+
cout << " Discrepancy at point # " << i + 1 << ": " << gp << ", R " << gp.perp() << ", Phi " << gp.phi()
118+
<< ", delta: " << referenceB - parametricB << " " << delta << endl;
119+
cout << " Reference: " << referenceB << ", Approximation: " << parametricB << endl;
120+
} else if (fail == 10) {
121+
cout << "..." << endl;
122+
}
123+
}
124+
}
125+
126+
if (fail != 0) {
127+
cout << "MF regression found: " << fail << " failures; max delta = " << maxdelta << endl;
128+
exit(EXIT_FAILURE);
129+
}
130+
}
131+
}

0 commit comments

Comments
 (0)