From 164f005d064699019cd73f9eec2123fa28d17f09 Mon Sep 17 00:00:00 2001
From: schaffer <R.D.Schaffer@cern.ch>
Date: Sun, 16 Aug 2020 15:11:10 +0200
Subject: [PATCH] adding in unit test for field optimizations. This version is
a detailed framework for testing, which will be simplified in further
commits.
---
MagneticField/MagFieldElements/CMakeLists.txt | 7 +
.../MagFieldElements/BFieldCache.h | 20 ++
.../MagFieldElements/BFieldMesh.h | 4 +
.../MagFieldElements/BFieldMesh.icc | 51 ++++-
.../test/BFieldExample_test.cxx | 178 ++++++++++++++++++
5 files changed, 258 insertions(+), 2 deletions(-)
create mode 100644 MagneticField/MagFieldElements/test/BFieldExample_test.cxx
diff --git a/MagneticField/MagFieldElements/CMakeLists.txt b/MagneticField/MagFieldElements/CMakeLists.txt
index d24244fcfb5..87b3674fda8 100644
--- a/MagneticField/MagFieldElements/CMakeLists.txt
+++ b/MagneticField/MagFieldElements/CMakeLists.txt
@@ -13,3 +13,10 @@ atlas_add_library( MagFieldElements
INCLUDE_DIRS ${ROOT_INCLUDE_DIRS}
LINK_LIBRARIES CxxUtils EventPrimitives GaudiKernel
PRIVATE_LINK_LIBRARIES ${ROOT_LIBRARIES} PathResolver )
+
+
+atlas_add_test( BFieldExample_test
+ SOURCES test/BFieldExample_test.cxx
+ INCLUDE_DIRS ${GTEST_INCLUDE_DIRS} #${GMOCK_INCLUDE_DIRS}
+ LINK_LIBRARIES TestTools CxxUtils EventPrimitives ${GTEST_LIBRARIES} #${GMOCK_LIBRARIES}
+ ENVIRONMENT "JOBOPTSEARCHPATH=${CMAKE_CURRENT_SOURCE_DIR}/share" )
diff --git a/MagneticField/MagFieldElements/MagFieldElements/BFieldCache.h b/MagneticField/MagFieldElements/MagFieldElements/BFieldCache.h
index e42b999ef3b..67abc8da026 100644
--- a/MagneticField/MagFieldElements/MagFieldElements/BFieldCache.h
+++ b/MagneticField/MagFieldElements/MagFieldElements/BFieldCache.h
@@ -80,6 +80,26 @@ public:
}
}
+ // set field array, filled externally
+ void setField(double field[][8]) {
+ std::copy( &field[0][0], &field[0][0] + 3*8, &m_field[0][0]);
+ }
+
+
+ // set the field values at each corner (rescale for current scale factor)
+ void printField()
+ {
+ // print out field values
+ std::cout << "Field at corner i, for each component j (Bz, Br, Bphi)" << std::endl;
+
+ for (int i = 0; i < 8; ++i) {
+ for (int j = 0; j < 3; ++j) {
+ std::cout << i << "," << j << ": " << m_field[j][i] << ", ";
+ }
+ std::cout << std::endl;
+ }
+ }
+
// set the multiplicative factor for the field vectors
void setBscale(double bscale) { m_scale = bscale; }
float bscale() { return m_scale; }
diff --git a/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.h b/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.h
index 9c85785a4af..0ba1f8a02f1 100644
--- a/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.h
+++ b/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.h
@@ -98,6 +98,9 @@ public:
double bscale() const { return m_scale; }
int memSize() const;
+ bool m_doNew{false};
+
+
protected:
std::array<double, 3> m_min;
std::array<double, 3> m_max;
@@ -112,6 +115,7 @@ private:
std::array<std::vector<int>,3> m_LUT;
std::array<double,3> m_invUnit; // inverse unit size in the LUT
int m_roff, m_zoff;
+
};
#include "MagFieldElements/BFieldMesh.icc"
#endif
diff --git a/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.icc b/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.icc
index 810c9739fdc..2a59bd862ce 100644
--- a/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.icc
+++ b/MagneticField/MagFieldElements/MagFieldElements/BFieldMesh.icc
@@ -70,10 +70,14 @@ BFieldMesh<T>::getCache(double z,
++iphi;
}
// store the bin edges
- cache.setRange(
- mz[iz], mz[iz + 1], mr[ir], mr[ir + 1], mphi[iphi], mphi[iphi + 1]);
+ cache.setRange(mz[iz], mz[iz + 1], mr[ir], mr[ir + 1], mphi[iphi], mphi[iphi + 1]);
+
+
// store the B field at the 8 corners
const int im0 = iz * m_zoff + ir * m_roff + iphi; // index of the first corner
+
+ if (!m_doNew) {
+
std::array<BFieldVector<T>, 8> field = {
m_field[im0],
m_field[im0 + 1],
@@ -84,7 +88,50 @@ BFieldMesh<T>::getCache(double z,
m_field[im0 + m_zoff + m_roff],
m_field[im0 + m_zoff + m_roff + 1]
};
+
cache.setField(field, scaleFactor);
+
+ }
+ else {
+
+ // cache.printField();
+
+ const double sf = scaleFactor;
+
+ double field1[3][8] = {
+ { sf * m_field[im0][0],
+ sf * m_field[im0 + 1][0],
+ sf * m_field[im0 + m_roff][0],
+ sf * m_field[im0 + m_roff + 1][0],
+ sf * m_field[im0 + m_zoff][0],
+ sf * m_field[im0 + m_zoff + 1][0],
+ sf * m_field[im0 + m_zoff + m_roff][0],
+ sf * m_field[im0 + m_zoff + m_roff + 1][0] },
+ { sf * m_field[im0][1],
+ sf * m_field[im0 + 1][1],
+ sf * m_field[im0 + m_roff][1],
+ sf * m_field[im0 + m_roff + 1][1],
+ sf * m_field[im0 + m_zoff][1],
+ sf * m_field[im0 + m_zoff + 1][1],
+ sf * m_field[im0 + m_zoff + m_roff][1],
+ sf * m_field[im0 + m_zoff + m_roff + 1][1] },
+ { sf * m_field[im0][2],
+ sf * m_field[im0 + 1][2],
+ sf * m_field[im0 + m_roff][2],
+ sf * m_field[im0 + m_roff + 1][2],
+ sf * m_field[im0 + m_zoff][2],
+ sf * m_field[im0 + m_zoff + 1][2],
+ sf * m_field[im0 + m_zoff + m_roff][2],
+ sf * m_field[im0 + m_zoff + m_roff + 1][2] }
+ };
+ cache.setField(field1);
+
+ // cache.printField();
+
+ }
+
+
+
// store the B scale
cache.setBscale(m_scale);
}
diff --git a/MagneticField/MagFieldElements/test/BFieldExample_test.cxx b/MagneticField/MagFieldElements/test/BFieldExample_test.cxx
new file mode 100644
index 00000000000..4ecf1079576
--- /dev/null
+++ b/MagneticField/MagFieldElements/test/BFieldExample_test.cxx
@@ -0,0 +1,178 @@
+/*
+ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include <iostream>
+#include "MagFieldElements/BFieldCache.h"
+#include "MagFieldElements/BFieldZone.h"
+#include <unistd.h>
+
+/**
+ * @brief Specimen of service caching expensive calculations.
+ * It uses internally std::map which is not thread-safely expandable.
+ **/
+class BFieldTest {
+public:
+ int runTest(bool doDebug = false ) {
+
+
+ // Create a field zone
+
+
+ double fieldz[] = { 19487, 19487, 19488, 19488, 19487, 19487, 19531, 19531, 19532, 19532, 19531, 19531, 6399, 6400, 6400, 6400, 6399, -1561, -1561, -1560, -1560, -1560, -1561, -1516, -1516, -1515, -1515, -1516, -1516, 20310, 20310, 20311, 20311, 20311, 20310, 20329, 20329, 20329, 20330, 20329, 20329, 7172, 7173, 7173, 7172, 7172, -814, -814, -813, -812, -813, -814, -795, -795, -794, -793, -794, -795, 20310, 20310, 20311, 20312, 20311, 20310, 20329, 20329, 20330, 20330, 20330, 20329, 7172, 7172, 7173, 7173, 7173, 7172, -813, -813, -812, -812, -813, -813, -794, -794, -793, -793, -793, -794, 19487, 19487, 19488, 19488, 19488, 19487, 19531, 19531, 19532, 19532, 19532, 19531, 6400, 6399, 6400, 6401, 6400, 6400, -1561, -1561, -1560, -1559, -1560, -1561, -1516, -1516, -1515, -1515, -1515, -1516, -1516, -1516 };
+
+ double fieldr[] = { -1357, -1356, -1353, -1354, -1354, -1357, -1366, -1366, -1362, -1363, -1363, -1366, -1378, -1374, -1375, -1375, -1378, -1388, -1388, -1385, -1386, -1386, -1388, -1394, -1394, -1390, -1391, -1391, -1394, -318, -318, -314, -315, -316, -318, -321, -321, -317, -318, -319, -321, -325, -321, -322, -322, -325, -328, -328, -324, -325, -326, -328, -330, -331, -326, -327, -328, -330, 312, 312, 316, 315, 315, 312, 315, 315, 319, 318, 318, 315, 319, 318, 323, 322, 321, 319, 322, 322, 326, 325, 325, 322, 324, 324, 328, 327, 327, 324, 1351, 1351, 1356, 1354, 1354, 1351, 1360, 1360, 1365, 1363, 1363, 1360, 1372, 1372, 1377, 1375, 1375, 1372, 1383, 1383, 1387, 1386, 1386, 1383, 1388, 1388, 1393, 1391, 1391, 1388, 1388, 1388 };
+
+ double fieldphi[] = { -2, 7, 3, 1, 6, -2, -2, 7, 3, 1, 6, -2, -2 , 3, 1, 6, -2, -2, 7, 3, 1, 6, -2, -2, 7, 3, 1, 6, -2, -1, 7, 3, 1, 6, -1, -1, 7, 3, 1, 6, -1, -1, 3, 1, 6, -1, -1, 7, 3, 1, 6, -1, -1, 8, 3, 1, 6, -1, 1, 7, 3, 2, 6, 1, 1, 7, 3, 2, 6, 1, 1, 7, 3, 2, 6, 1, 1, 7, 3, 2, 6, 1, 0, 8, 3, 2, 6, 0, 2, 7, 3, 2, 6, 2, 2, 7, 3, 2, 6, 2, 2, 7, 3, 2, 6, 2, 2, 7, 3, 2, 6, 2, 2, 8, 3, 2, 6, 2, 2, 2 };
+
+
+ int id{5};
+ double zmin{-1400}, zmax{1400}, rmin{1200}, rmax{1300}, phimin{0}, phimax{6.28319}, bscale{1e-07};
+ BFieldZone zone(id, zmin, zmax, rmin, rmax, phimin, phimax, bscale);
+ // add in field
+ int nmeshz{4}, nmeshr{5}, nmeshphi{6};
+ int nfield = nmeshz * nmeshr * nmeshphi;
+
+ zone.reserve(nmeshz, nmeshr, nmeshphi);
+
+ double meshz[] = { -1400, -466.93, 466.14, 1400 };
+
+ for (int j = 0; j < nmeshz; j++) {
+ zone.appendMesh(0, meshz[j]);
+
+ if (doDebug) std::cout << "j, meshz " << j << ", " << meshz[j] << std::endl;
+
+ }
+
+ if (doDebug) std::cout << "did meshz " << std::endl;
+
+ double meshr[] = { 1200, 1225, 1250, 1275, 1300 };
+
+ for (int j = 0; j < nmeshr; j++) {
+ zone.appendMesh(1, meshr[j]);
+ }
+
+ if (doDebug) std::cout << "did meshr " << std::endl;
+
+ double meshphi[] = { 0, 1.25664, 2.51327, 3.76991, 5.02655, 6.28318 };
+
+ for (int j = 0; j < nmeshphi; j++) {
+ zone.appendMesh(2, meshphi[j]);
+ }
+
+ if (doDebug) std::cout << "did meshphi " << std::endl;
+
+ for (int j = 0; j < nfield; j++) {
+ BFieldVector<short> field(fieldz[j], fieldr[j], fieldphi[j]);
+ zone.appendField(field);
+ }
+
+ if (doDebug) std::cout << "did field " << std::endl;
+
+ // build (trivial) look up table for zone
+ zone.buildLUT();
+
+ if (doDebug) std::cout << "did buildLUT " << std::endl;
+
+
+ // fill the cache, pass in current scale factor
+ BFieldCache cache3d;
+ double z{0}, r{1250}, phi{1.6};
+
+ if (doDebug) std::cout << "do getCache old " << std::endl;
+
+ zone.m_doNew = false;
+
+ zone.getCache(z, r, phi, cache3d, 1);
+
+ cache3d.printField();
+
+ if (doDebug) std::cout << "do getCache new " << std::endl;
+
+ zone.m_doNew = true;
+
+ zone.getCache(z, r, phi, cache3d, 1);
+
+ cache3d.printField();
+
+
+ // get field at steps of 10 mm from 1200 to 1300
+
+ int status{0};
+
+ double z0 = z;
+ double r0 = 1200;
+ double phi0 = phi;
+ double xyz[3] = {0, 0, 0};
+ double bxyz_std[3] = {0, 0, 0};
+ double bxyz_new[3] = {0, 0, 0};
+
+ for (unsigned int i = 0; i < 10; ++i) {
+ double r1 = r0 + 5 + i*10.;
+
+ xyz[0] = r1 * cos(phi0);
+ xyz[1] = r1 * sin(phi0);
+ xyz[2] = z0;
+
+ // do interpolation (cache3d has correct scale factor)
+ zone.m_doNew = false;
+ zone.getCache(z, r, phi, cache3d, 1);
+ cache3d.getB(xyz, r1, phi, bxyz_std, 0);
+
+ std::cout << "get field std: i, bxyz_std " << i << " "
+ << bxyz_std[0] << ", " << bxyz_std[1] << ", "
+ << bxyz_std[2] << std::endl;
+
+ zone.m_doNew = true;
+ zone.getCache(z, r, phi, cache3d, 1);
+ cache3d.getB(xyz, r1, phi, bxyz_new, 0);
+
+ std::cout << "get field new: i, bxyz_new " << i << " "
+ << bxyz_new[0] << ", " << bxyz_new[1] << ", "
+ << bxyz_new[2] << std::endl;
+
+ if (bxyz_std[0] != bxyz_new[0]) {
+ std::cout << "failed bz comparison" << std::endl;
+ status = 1;
+ }
+ if (bxyz_std[1] != bxyz_new[1]) {
+ std::cout << "failed br comparison" << std::endl;
+ status = 1;
+ }
+ if (bxyz_std[2] != bxyz_new[2]) {
+ std::cout << "failed bphi comparison" << std::endl;
+ status = 1;
+ }
+
+ }
+
+ std::cout << "runTest: status " << status << std::endl;
+
+ return status;
+
+ }
+
+private:
+};
+
+
+
+
+int main() {
+
+ std::cout << "start BFieldExample test" << std::endl;
+
+
+ BFieldTest bftest;
+ int status = bftest.runTest(true);
+
+ if (status == 0) {
+ std::cout << "Test passed OK" << std::endl;
+ return 0;
+ }
+ else {
+ std::cout << "Test failed" << std::endl;
+ return 1;
+ }
+
+}
--
GitLab