test/alpaka/testParabolicParametrizedMagneticField.dev.cc

0001 #include <iostream>
0002 #include <fstream>
0003 #include <Eigen/Core>
0004 #include <alpaka/alpaka.hpp>
0005
0006 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0007 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0008 #include "FWCore/Utilities/interface/FileInPath.h"
0009 #include "FWCore/Utilities/interface/Exception.h"
0010 #include "FWCore/Utilities/interface/stringize.h"
0011 #include "HeterogeneousCore/AlpakaInterface/interface/config.h"
0012 #include "HeterogeneousCore/AlpakaInterface/interface/memory.h"
0013 #include "HeterogeneousCore/AlpakaInterface/interface/workdivision.h"
0014 #include "MagneticField/ParametrizedEngine/interface/ParabolicParametrizedMagneticField.h"
0015
0016 using namespace edm;
0017 using namespace std;
0018 using namespace ALPAKA_ACCELERATOR_NAMESPACE;
0019 using namespace magneticFieldParabolicPortable;
0020 using Vector3f = Eigen::Matrix<float, 3, 1>;
0021
0022 struct MagneticFieldKernel {
0023   template <typename TAcc, typename T>
0024   ALPAKA_FN_ACC void operator()(TAcc const& acc, T const* __restrict__ in, T* __restrict__ out, size_t size) const {
0025     for (auto index : cms::alpakatools::uniform_elements(acc, size)) {
0026       out[index](0) = 0;
0027       out[index](1) = 0;
0028       out[index](2) = magneticFieldAtPoint(in[index]);
0029     }
0030   }
0031 };
0032
0033 int main() {
0034   // get the list of devices on the current platform
0035   auto const& devices = cms::alpakatools::devices<Platform>();
0036   if (devices.empty()) {
0037     std::cerr << "No devices available for the " EDM_STRINGIZE(ALPAKA_ACCELERATOR_NAMESPACE) " backend, "
0038       "the test will be skipped.\n";
0039     exit(EXIT_FAILURE);
0040   }
0041
0042   ifstream file;
0043   edm::FileInPath mydata("MagneticField/Engine/data/Regression/referenceField_160812_RII_3_8T.bin");
0044   file.open(mydata.fullPath().c_str(), ios::binary);
0045
0046   int count = 0;
0047   float px, py, pz;
0048   float bx, by, bz;
0049   vector<Vector3f> points;
0050   vector<GlobalVector> referenceB_vec;
0051
0052   int numberOfPoints = 100;
0053   points.reserve(numberOfPoints);
0054   referenceB_vec.reserve(numberOfPoints);
0055   do {
0056     if (!(file.read((char*)&px, sizeof(float)) && file.read((char*)&py, sizeof(float)) &&
0057           file.read((char*)&pz, sizeof(float)) && file.read((char*)&bx, sizeof(float)) &&
0058           file.read((char*)&by, sizeof(float)) && file.read((char*)&bz, sizeof(float))))
0059       break;
0060
0061     const auto point = Vector3f(px, py, pz);
0062     if (!isValid(point))
0063       continue;
0064
0065     points.push_back(Vector3f(px, py, pz));
0066     referenceB_vec.push_back(GlobalVector(bx, by, bz));
0067     count++;
0068   } while (count < numberOfPoints);
0069
0070   const size_t size = points.size();
0071   // allocate the input and output host buffer in pinned memory accessible by the Platform devices
0072   auto points_host = cms::alpakatools::make_host_buffer<Vector3f[], Platform>(size);
0073   auto field_host = cms::alpakatools::make_host_buffer<Vector3f[], Platform>(size);
0074   // fill the input buffers, and the output buffer with zeros
0075   for (size_t i = 0; i < size; ++i) {
0076     points_host[i] = points[i];
0077     field_host[i] = Vector3f::Zero();
0078   }
0079
0080   float resolution = 0.2;
0081   float maxdelta = 0.;
0082   int fail = 0;
0083
0084   // run the test on each device
0085   for (auto const& device : devices) {
0086     auto queue = Queue(device);
0087     // allocate input and output buffers on the device
0088     auto points_dev = cms::alpakatools::make_device_buffer<Vector3f[]>(queue, size);
0089     auto field_dev = cms::alpakatools::make_device_buffer<Vector3f[]>(queue, size);
0090
0091     // copy the input data to the device; the size is known from the buffer objects
0092     alpaka::memcpy(queue, points_dev, points_host);
0093     // fill the output buffer with zeros; the size is known from the buffer objects
0094     alpaka::memset(queue, field_dev, 0.);
0095
0096     auto workDiv = cms::alpakatools::make_workdiv<Acc1D>(1, size);
0097     alpaka::exec<Acc1D>(queue, workDiv, MagneticFieldKernel{}, points_dev.data(), field_dev.data(), size);
0098
0099     // copy the results from the device to the host
0100     alpaka::memcpy(queue, field_host, field_dev);
0101
0102     // wait for the kernel and the potential copy to complete
0103     alpaka::wait(queue);
0104
0105     // check the results
0106     for (uint i = 0; i < points.size(); i++) {
0107       const auto& point = points[i];
0108       const auto& referenceB = referenceB_vec[i];
0109       GlobalVector parametricB(field_host[i](0), field_host[i](1), field_host[i](2));
0110       float delta = (referenceB - parametricB).mag();
0111       if (delta > resolution) {
0112         ++fail;
0113         if (delta > maxdelta)
0114           maxdelta = delta;
0115         if (fail < 10) {
0116           const GlobalPoint gp(point(0), point(1), point(2));
0117           cout << " Discrepancy at point  # " << i + 1 << ": " << gp << ", R " << gp.perp() << ", Phi " << gp.phi()
0118                << ", delta: " << referenceB - parametricB << " " << delta << endl;
0119           cout << " Reference: " << referenceB << ", Approximation: " << parametricB << endl;
0120         } else if (fail == 10) {
0121           cout << "..." << endl;
0122         }
0123       }
0124     }
0125
0126     if (fail != 0) {
0127       cout << "MF regression found: " << fail << " failures; max delta = " << maxdelta << endl;
0128       exit(EXIT_FAILURE);
0129     }
0130   }
0131 }