OPALX (Object Oriented Parallel Accelerator Library for Exascal) master (dc2a29eed580)
OPALX
Loading...
Searching...
No Matches
TestMultipoleTStraight.cpp
Go to the documentation of this file.
1//
2// Cubic Spline Interpolation to replace GSL spline
3//
4// Copyright (c) 2023, Paul Scherrer Institute, Villigen PSI, Switzerland
5// All rights reserved
6//
7// This file is part of OPAL.
8//
9// OPAL is free software: you can redistribute it and/or modify
10// it under the terms of the GNU General Public License as published by
11// the Free Software Foundation, either version 3 of the License, or
12// (at your option) any later version.
13//
14// You should have received a copy of the GNU General License
15// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
16//
17
18#include <vector>
19#include "AbsBeamline/MultipoleT.h"
20#include "AbstractObjects/OpalData.h"
21#include "Structure/Beam.h"
22#include "Structure/DataSink.h"
23#include "gtest/gtest.h"
24
25class TestMultipoleTStraight : public testing::Test, public MultipoleT {
26public:
27 TestMultipoleTStraight() : MultipoleT("Magnet") {}
28
29 static void SetUpTestSuite() {
30 int argc = 0;
31 char** argv = nullptr;
32 ippl::initialize(argc, argv);
33 // DataSink requires a basename to create *.stat / *.lbal writers.
34 OpalData::getInstance()->storeInputFn("unit_test.opal");
35 // Many OPAL writers assume `gmsg` is initialized (see SDDSWriter/StatWriter).
36 // Unit tests normally don't set this up via Main().
37 gmsg = new Inform(nullptr, -1);
38 // DataSink::DataSink() constructs HDF5 writers when enabled, but the unit
39 // test doesn't have an H5PartWrapper. Disable HDF5 for this smoke test.
40 Options::enableHDF5 = false;
41 }
42 static void TearDownTestSuite() {
43 delete gmsg;
44 gmsg = nullptr;
45 ippl::finalize();
46 }
47
48 // Test helper functions
49 static Vector_t<double, 3> curvilinearToGlobal(
50 const Vector_t<double, 3>& local, const Vector_t<double, 3>& elementEntry,
51 const double elementLength) {
52 const double x = local[0] + elementEntry[0];
53 const double y = local[1] + elementEntry[1];
54 const double z = local[2] + elementLength / 2 + elementEntry[2];
55 return {x, y, z};
56 }
57
58 void grabTransverseDataLine(
59 std::vector<double>& line, const double s, const double width,
60 const Vector_t<double, 3>& elementEntry, const double elementLength) {
61 // Make the bunch
62 const auto bunch = makeBunch(line.size());
63 const auto pc = bunch->getParticleContainer();
64 // Create the local views and data
65 std::vector<Vector_t<double, 3>> localR(line.size());
66 const auto hostR = Kokkos::create_mirror_view(pc->R.getView());
67 const auto hostB = Kokkos::create_mirror_view(pc->B.getView());
68 // Set the particle positions
69 const double stepSize = width / static_cast<double>(line.size() - 1);
70 for (size_t i = 0; i < line.size(); ++i) {
71 localR[i] = {static_cast<double>(i) * stepSize - width / 2, 0, s};
72 hostR(i) = curvilinearToGlobal(localR[i], elementEntry, elementLength);
73 }
74 Kokkos::deep_copy(pc->R.getView(), hostR);
75 Kokkos::deep_copy(pc->B.getView(), Vector_t<double, 3>{});
76 pc->setQ(pc->getChargePerParticle());
77 ippl::Comm->barrier();
78 Kokkos::fence();
79 // Register the bunch with the element
80 bunch->setT(0.0);
81 double startField, endField;
82 initialise(bunch.get(), startField, endField);
83 // Get the fields
84 apply(pc);
85 // Return the fields
86 Kokkos::deep_copy(hostB, pc->B.getView());
87 Kokkos::fence();
88 for (size_t i = 0; i < line.size(); ++i) {
89 line[i] = std::hypot(hostB(i)[0], hostB(i)[1], hostB(i)[2]);
90 // std::cout << i << ": Local=" << localR[i] << ", Global=" << hostR[i]
91 // << ", mag(B)=" << line[i] << std::endl;
92 }
93 // std::cout << std::endl;
94 }
95
96 void grabVerticalDataLine(
97 std::vector<double>& line, const double s, const double height,
98 const Vector_t<double, 3>& elementEntry, const double elementLength) {
99 // Make the bunch
100 const auto bunch = makeBunch(line.size());
101 const auto pc = bunch->getParticleContainer();
102 // Create the local views and data
103 std::vector<Vector_t<double, 3>> localR(line.size());
104 const auto hostR = Kokkos::create_mirror_view(pc->R.getView());
105 const auto hostB = Kokkos::create_mirror_view(pc->B.getView());
106 // Set the particle positions
107 const double stepSize = height / static_cast<double>(line.size() - 1);
108 for (size_t i = 0; i < line.size(); ++i) {
109 localR[i] = {0, static_cast<double>(i) * stepSize - height / 2, s};
110 hostR(i) = curvilinearToGlobal(localR[i], elementEntry, elementLength);
111 }
112 Kokkos::deep_copy(pc->R.getView(), hostR);
113 Kokkos::deep_copy(pc->B.getView(), Vector_t<double, 3>{});
114 pc->setQ(pc->getChargePerParticle());
115 ippl::Comm->barrier();
116 Kokkos::fence();
117 // Register the bunch with the element
118 bunch->setT(0.0);
119 double startField, endField;
120 initialise(bunch.get(), startField, endField);
121 // Get the fields
122 apply(pc);
123 // Return the fields
124 Kokkos::deep_copy(hostB, pc->B.getView());
125 Kokkos::fence();
126 for (size_t i = 0; i < line.size(); ++i) {
127 line[i] = std::hypot(hostB(i)[0], hostB(i)[1], hostB(i)[2]);
128 // std::cout << i << ": Local=" << localR[i] << ", Global=" << hostR[i]
129 // << ", mag(B)=" << line[i] << std::endl;
130 }
131 // std::cout << std::endl;
132 }
133
134 void grabLongitudinalDivCurlLine(
135 std::vector<double>& fieldLine, std::vector<double>& divLine,
136 std::vector<Vector_t<double, 3>>& curlLine, const double x, const double length,
137 const Vector_t<double, 3>& elementEntry, const double elementLength, const double dr) {
138 const double stepSize = length / static_cast<double>(divLine.size() - 1);
139 const double startS = 0 - (elementLength - length) / 2;
140 for (size_t i = 0; i < divLine.size(); ++i) {
141 // Get the sourrounding 6 B fields
142 Vector_t<double, 3> local{x, 0, static_cast<double>(i) * stepSize + startS};
143 Vector_t<double, 3> R = curvilinearToGlobal(local, elementEntry, elementLength);
144 Vector_t<double, 3> B;
145 Vector_t<double, 3> Bxp;
146 Vector_t<double, 3> Bxm;
147 Vector_t<double, 3> Byp;
148 Vector_t<double, 3> Bym;
149 Vector_t<double, 3> Bzp;
150 Vector_t<double, 3> Bzm;
151 Vector_t<double, 3> E;
152 apply(R, {}, 0.0, E, B);
153 apply(R + Vector_t<double, 3>{dr, 0, 0}, {}, 0.0, E, Bxp);
154 apply(R - Vector_t<double, 3>{dr, 0, 0}, {}, 0.0, E, Bxm);
155 apply(R + Vector_t<double, 3>{0, dr, 0}, {}, 0.0, E, Byp);
156 apply(R - Vector_t<double, 3>{0, dr, 0}, {}, 0.0, E, Bym);
157 apply(R + Vector_t<double, 3>{0, 0, dr}, {}, 0.0, E, Bzp);
158 apply(R - Vector_t<double, 3>{0, 0, dr}, {}, 0.0, E, Bzm);
159 // Calculate the divergence and curl using central differences
160 divLine[i] = (Bxp[0] - Bxm[0] + Byp[1] - Bym[1] + Bzp[2] - Bzm[2]) / 2 / dr;
161 curlLine[i][0] = (Byp[2] - Bym[2] - Bzp[1] + Bzm[1]) / 2 / dr;
162 curlLine[i][1] = (Bzp[0] - Bzm[0] - Bxp[2] + Bxm[2]) / 2 / dr;
163 curlLine[i][2] = (Bxp[1] - Bxm[1] - Byp[0] + Bym[0]) / 2 / dr;
164 fieldLine[i] = std::hypot(B[0], B[1], B[2]);
165 // std::cout << i << ": Local=" << local << ", Global=" << R << ", div(B)=" <<
166 // divLine[i]
167 // << ", curl(B)=" << curlLine[i] << ", B=" << fieldLine[i] << std::endl;
168 }
169 }
170
171 class TestableFieldSolverCmd : public FieldSolverCmd {
172 public:
173 void setType(const std::string& t) {
174 Attributes::setPredefinedString(this->itsAttr[FIELDSOLVER::TYPE], t);
175 }
176
177 void setBCX(const std::string& bc) {
178 Attributes::setPredefinedString(this->itsAttr[FIELDSOLVER::BCFFTX], bc);
179 }
180 void setBCY(const std::string& bc) {
181 Attributes::setPredefinedString(this->itsAttr[FIELDSOLVER::BCFFTY], bc);
182 }
183 void setBCZ(const std::string& bc) {
184 Attributes::setPredefinedString(this->itsAttr[FIELDSOLVER::BCFFTZ], bc);
185 }
186 };
187
188 std::shared_ptr<FieldSolverCmd> fsCmdBase_m;
189 std::shared_ptr<DataSink> dataSink_m;
190
191 std::shared_ptr<PartBunch_t> makeBunch(const size_t numParticles) {
192 dataSink_m = std::make_shared<DataSink>();
193 const auto fsCmd = std::make_shared<TestableFieldSolverCmd>();
194 fsCmdBase_m = fsCmd;
195 fsCmd->setType("NONE");
196 fsCmd->setNX(8);
197 fsCmd->setNY(8);
198 fsCmd->setNZ(8);
199 fsCmd->setBCX("PERIODIC");
200 fsCmd->setBCY("PERIODIC");
201 fsCmd->setBCZ("PERIODIC");
202 auto beam = std::make_shared<Beam>();
203 Beam* opBeam = Beam::find("UNNAMED_BEAM");
204 EXPECT_NE(opBeam, nullptr);
205 auto bunch = std::make_shared<PartBunch_t>(
206 /*qi=*/std::vector{1.0}, /*mi=*/std::vector{1.0},
207 /*beams=*/std::vector<Beam*>{opBeam},
208 /*totalParticlesPerBeam=*/std::vector<size_t>{numParticles},
209 /*lbt=*/1.0, /*integration_method=*/"LF2", fsCmdBase_m.get(), dataSink_m.get());
210 bunch->getParticleContainer()->createParticles(numParticles);
211 return bunch;
212 }
213};
214
215// Test transversely across the magnet to check for a dipole (constant) field.
216TEST_F(TestMultipoleTStraight, Dipole) {
217 std::vector<double> line(101);
218 // Set up the magnet
219 constexpr double length = 4.4;
220 constexpr double bendAngle = 0.0;
221 constexpr double dipoleField = 1.0;
222 setBendAngle(bendAngle, false);
223 setElementLength(length);
224 setAperture(3.5, 3.5);
225 setFringeField(length / 2, 0.3, 0.3);
226 setRotation(0.0);
227 setEntranceAngle(0.0);
228 setMaxOrder(5, 10);
229 // Check dipole has the correct constant transverse field magnitude at the center
230 setTransProfile({dipoleField});
231 grabTransverseDataLine(line, 0, 3.0, {0, 0, 0}, length);
232 for (const double val : line) {
233 EXPECT_NEAR(val, dipoleField, 1e-2);
234 }
235 // Check dipole has the correct constant transverse field magnitude at the edge
236 setTransProfile({dipoleField});
237 grabTransverseDataLine(line, -length / 2, 3.0, {0, 0, 0}, length);
238 for (const double val : line) {
239 EXPECT_NEAR(val, dipoleField / 2, 1e-2);
240 }
241}
242
243// Test transversely across the magnet to check for a quadrupole (linear) field.
244TEST_F(TestMultipoleTStraight, Quadrupole) {
245 constexpr unsigned int samplesPerSide = 50;
246 std::vector<double> line(2 * samplesPerSide + 1);
247 // Set up the magnet
248 constexpr double length = 4.4;
249 constexpr double bendAngle = 0.0;
250 constexpr double quadrupoleField = 1.0;
251 setBendAngle(bendAngle, false);
252 setElementLength(length);
253 setAperture(3.5, 3.5);
254 setFringeField(length / 2, 0.3, 0.3);
255 setRotation(0.0);
256 setEntranceAngle(0.0);
257 setMaxOrder(5, 10);
258 // Check quadrupole has linear field magnitude at the center
259 setTransProfile({0, quadrupoleField});
260 grabTransverseDataLine(line, 0, 3.0, {0, 0, 0}, length);
261 double expected = 0;
262 double delta = line[samplesPerSide + 1] - line[samplesPerSide];
263 EXPECT_NEAR(delta, 0.03, 1e-5);
264 EXPECT_NEAR(line[samplesPerSide], expected, 1e-5);
265 for (size_t i = 0; i < samplesPerSide; ++i) {
266 expected += delta;
267 EXPECT_NEAR(line[samplesPerSide - i - 1], expected, 1e-2);
268 EXPECT_NEAR(line[samplesPerSide + i + 1], expected, 1e-2);
269 }
270 // Check quadrupole has linear field magnitude at the edge with half the slope
271 setTransProfile({0, quadrupoleField});
272 grabTransverseDataLine(line, -length / 2, 3.0, {0, 0, 0}, length);
273 expected = 0;
274 delta = line[samplesPerSide + 1] - line[samplesPerSide];
275 EXPECT_NEAR(delta, 0.015, 1e-5);
276 EXPECT_NEAR(line[samplesPerSide], expected, 1e-5);
277 for (size_t i = 0; i < samplesPerSide; ++i) {
278 expected += delta;
279 EXPECT_NEAR(line[samplesPerSide - i - 1], expected, 1e-2);
280 EXPECT_NEAR(line[samplesPerSide + i + 1], expected, 1e-2);
281 }
282}
283
284// Test transversely across the magnet to check for a sextupole (quadratic) field.
285TEST_F(TestMultipoleTStraight, Sextupole) {
286 constexpr unsigned int samplesPerSide = 50;
287 std::vector<double> line(2 * samplesPerSide + 1);
288 // Set up the magnet
289 constexpr double length = 4.4;
290 constexpr double bendAngle = 0.0;
291 constexpr double sextupoleField = 1.0;
292 setBendAngle(bendAngle, false);
293 setElementLength(length);
294 setAperture(3.5, 3.5);
295 setFringeField(length / 2, 0.3, 0.3);
296 setRotation(0.0);
297 setEntranceAngle(0.0);
298 setMaxOrder(5, 10);
299 // Check sextupole has quadratic field magnitude at the center
300 setTransProfile({0, 0, sextupoleField});
301 grabTransverseDataLine(line, 0, 3.0, {0, 0, 0}, length);
302 double expected = 0;
303 double delta = std::sqrt(line[samplesPerSide + 1] - line[samplesPerSide]);
304 EXPECT_NEAR(delta, 0.03, 1e-5);
305 EXPECT_NEAR(line[samplesPerSide], expected, 1e-5);
306 for (size_t i = 0; i < samplesPerSide; ++i) {
307 expected = std::pow(static_cast<double>(i + 1) * delta, 2);
308 EXPECT_NEAR(line[samplesPerSide - i - 1], expected, 1e-2);
309 EXPECT_NEAR(line[samplesPerSide + i + 1], expected, 1e-2);
310 }
311 // Check sextupole has quadratic field magnitude at the edge with half the slope
312 setTransProfile({0, 0, sextupoleField});
313 grabTransverseDataLine(line, -length / 2, 3.0, {0, 0, 0}, length);
314 expected = 0;
315 delta = std::sqrt(line[samplesPerSide + 1] - line[samplesPerSide]);
316 EXPECT_NEAR(delta, 0.02121, 1e-5);
317 EXPECT_NEAR(line[samplesPerSide], expected, 1e-5);
318 for (size_t i = 0; i < samplesPerSide; ++i) {
319 expected = std::pow(static_cast<double>(i + 1) * delta, 2);
320 EXPECT_NEAR(line[samplesPerSide - i - 1], expected, 1e-2);
321 EXPECT_NEAR(line[samplesPerSide + i + 1], expected, 1e-2);
322 }
323}
324
325// Test transversely across the magnet to check for an octupole (cubic) field.
326TEST_F(TestMultipoleTStraight, Octupole) {
327 constexpr unsigned int samplesPerSide = 50;
328 std::vector<double> line(2 * samplesPerSide + 1);
329 // Set up the magnet
330 constexpr double length = 4.4;
331 constexpr double bendAngle = 0.0;
332 constexpr double octupoleField = 1.0;
333 setBendAngle(bendAngle, false);
334 setElementLength(length);
335 setAperture(3.5, 3.5);
336 setFringeField(length / 2, 0.3, 0.3);
337 setRotation(0.0);
338 setEntranceAngle(0.0);
339 setMaxOrder(5, 10);
340 // Check octupole has cubic field magnitude at the center
341 setTransProfile({0, 0, 0, octupoleField});
342 grabTransverseDataLine(line, 0, 3.0, {0, 0, 0}, length);
343 double expected = 0;
344 double delta = std::cbrt(line[samplesPerSide + 1] - line[samplesPerSide]);
345 EXPECT_NEAR(delta, 0.03, 1e-5);
346 EXPECT_NEAR(line[samplesPerSide], expected, 1e-5);
347 for (size_t i = 0; i < samplesPerSide; ++i) {
348 expected = std::pow(static_cast<double>(i + 1) * delta, 3);
349 EXPECT_NEAR(line[samplesPerSide - i - 1], expected, 1e-2);
350 EXPECT_NEAR(line[samplesPerSide + i + 1], expected, 1e-2);
351 }
352 // Check octupole has cubic field magnitude at the edge with half the slope
353 setTransProfile({0, 0, 0, octupoleField});
354 grabTransverseDataLine(line, -length / 2, 3.0, {0, 0, 0}, length);
355 expected = 0;
356 delta = std::cbrt(line[samplesPerSide + 1] - line[samplesPerSide]);
357 EXPECT_NEAR(delta, 0.023811, 1e-5);
358 EXPECT_NEAR(line[samplesPerSide], expected, 1e-5);
359 for (size_t i = 0; i < samplesPerSide; ++i) {
360 expected = std::pow(static_cast<double>(i + 1) * delta, 3);
361 EXPECT_NEAR(line[samplesPerSide - i - 1], expected, 1e-2);
362 EXPECT_NEAR(line[samplesPerSide + i + 1], expected, 1e-2);
363 }
364}
365
366// Test transversely across the magnet to check for a decapole (quartic) field.
367TEST_F(TestMultipoleTStraight, Decapole) {
368 constexpr unsigned int samplesPerSide = 50;
369 std::vector<double> line(2 * samplesPerSide + 1);
370 // Set up the magnet
371 constexpr double length = 4.4;
372 constexpr double bendAngle = 0.0;
373 constexpr double decapoleField = 1.0;
374 setBendAngle(bendAngle, false);
375 setElementLength(length);
376 setAperture(3.5, 3.5);
377 setFringeField(length / 2, 0.3, 0.3);
378 setRotation(0.0);
379 setEntranceAngle(0.0);
380 setMaxOrder(5, 10);
381 // Check decapole has quartic field magnitude at the center
382 setTransProfile({0, 0, 0, 0, decapoleField});
383 grabTransverseDataLine(line, 0, 3.0, {0, 0, 0}, length);
384 double expected = 0;
385 double delta = std::pow(line[samplesPerSide + 1] - line[samplesPerSide], 1.0 / 4.0);
386 EXPECT_NEAR(delta, 0.03, 1e-5);
387 EXPECT_NEAR(line[samplesPerSide], expected, 1e-5);
388 for (size_t i = 0; i < samplesPerSide; ++i) {
389 expected = std::pow(static_cast<double>(i + 1) * delta, 4);
390 EXPECT_NEAR(line[samplesPerSide - i - 1], expected, 1e-2);
391 EXPECT_NEAR(line[samplesPerSide + i + 1], expected, 1e-2);
392 }
393 // Check decapole has quartic field magnitude at the edge with half the slope
394 setTransProfile({0, 0, 0, 0, decapoleField});
395 grabTransverseDataLine(line, -length / 2, 3.0, {0, 0, 0}, length);
396 expected = 0;
397 delta = std::pow(line[samplesPerSide + 1] - line[samplesPerSide], 1.0 / 4.0);
398 EXPECT_NEAR(delta, 0.025227, 1e-5);
399 EXPECT_NEAR(line[samplesPerSide], expected, 1e-5);
400 for (size_t i = 0; i < samplesPerSide; ++i) {
401 expected = std::pow(static_cast<double>(i + 1) * delta, 4);
402 EXPECT_NEAR(line[samplesPerSide - i - 1], expected, 1e-2);
403 EXPECT_NEAR(line[samplesPerSide + i + 1], expected, 1e-2);
404 }
405}
406
407// Check that the field along the center line is Maxwellian (zero div and curl)
408TEST_F(TestMultipoleTStraight, DivCurl) {
409 constexpr unsigned int samplesPerSide = 20;
410 std::vector<double> divLine(2 * samplesPerSide + 1);
411 std::vector<double> fieldLine(2 * samplesPerSide + 1);
412 std::vector<Vector_t<double, 3>> curlLine(2 * samplesPerSide + 1);
413 constexpr double dr = 0.001;
414 // Set up the magnet
415 constexpr double length = 4.4;
416 constexpr double bendAngle = 0.0;
417 setBendAngle(bendAngle, false);
418 setElementLength(length);
419 setAperture(3.5, 3.5);
420 setFringeField(length / 2, 0.3, 0.3);
421 setRotation(0.0);
422 setEntranceAngle(0.0);
423 setMaxOrder(5, 10);
424 setTransProfile({1, 1, 1, 1, 1});
425 // Get the data
426 grabLongitudinalDivCurlLine(fieldLine, divLine, curlLine, 0, length * 1.5, {}, length, dr);
427 // Analyse the results
428 for (size_t i = 0; i < divLine.size(); ++i) {
429 const double divError = std::abs(divLine[i]) * dr / fieldLine[i];
430 const double curlError =
431 std::hypot(curlLine[i][0], curlLine[i][1], curlLine[i][2]) * dr / fieldLine[i];
432 EXPECT_NEAR(divError, 0, 1e-5);
433 EXPECT_NEAR(curlError, 0, 1e-2);
434 }
435}
436
437// Check the const and non-const geometry access APIs return the same object
438TEST_F(TestMultipoleTStraight, Geometry) {
439 TestMultipoleTStraight* magnet = this;
440 auto* geom = &magnet->getGeometry();
441 EXPECT_NE(geom, nullptr);
442 auto* constGeom = &const_cast<const TestMultipoleTStraight*>(magnet)->getGeometry();
443 EXPECT_NE(constGeom, nullptr);
444 EXPECT_EQ(geom, constGeom);
445}
446
447// Check that the field outside the bounding box is not calculated
448TEST_F(TestMultipoleTStraight, BoundingBox) {
449 std::vector<double> line(3);
450 // Set up the magnet
451 constexpr double length = 4;
452 constexpr double bendAngle = 0.0;
453 constexpr double dipoleField = 1.0;
454 setBendAngle(bendAngle, false);
455 setElementLength(length);
456 setAperture(3.5, 3.5);
457 setFringeField(length / 2, 3, 3);
458 setRotation(0.0);
459 setEntranceAngle(0.0);
460 setMaxOrder(5, 10);
461 setBoundingBoxLength(6);
462 setTransProfile({dipoleField});
463 // Check field vanishes outside the bounding box
464 grabTransverseDataLine(line, -3.1, 3.0, {0, 0, 0}, length);
465 for (const double val : line) {
466 EXPECT_EQ(val, 0.0);
467 }
468 // Check field is present inside the bounding box
469 grabTransverseDataLine(line, -2.9, 3.0, {0, 0, 0}, length);
470 for (const double val : line) {
471 EXPECT_NE(val, 0.0);
472 }
473 // Check field is present inside the bounding box
474 grabTransverseDataLine(line, 0.0, 3.0, {0, 0, 0}, length);
475 for (const double val : line) {
476 EXPECT_NE(val, 0.0);
477 }
478 // Check field is present inside the bounding box
479 grabTransverseDataLine(line, 2.9, 3.0, {0, 0, 0}, length);
480 for (const double val : line) {
481 EXPECT_NE(val, 0.0);
482 }
483 // Check field vanishes outside the bounding box
484 grabTransverseDataLine(line, 3.1, 3.0, {0, 0, 0}, length);
485 for (const double val : line) {
486 EXPECT_EQ(val, 0.0);
487 }
488}
489
490// Check that the field outside the aperture is not calculated
491TEST_F(TestMultipoleTStraight, HorzAperture) {
492 std::vector<double> line(9);
493 // Set up the magnet
494 constexpr double length = 4;
495 constexpr double bendAngle = 0.0;
496 constexpr double dipoleField = 1.0;
497 setBendAngle(bendAngle, false);
498 setElementLength(length);
499 setAperture(3.5, 3.5);
500 setFringeField(length / 2, 3, 3);
501 setRotation(0.0);
502 setEntranceAngle(0.0);
503 setMaxOrder(5, 10);
504 setTransProfile({dipoleField});
505 // Check field vanishes outside the aperture
506 grabTransverseDataLine(line, 0, 4.0, {0, 0, 0}, length);
507 EXPECT_EQ(line[0], 0.0); // -4.0
508 EXPECT_NE(line[1], 0.0); // -3.0
509 EXPECT_NE(line[2], 0.0); // -2.0
510 EXPECT_NE(line[3], 0.0); // -1.0
511 EXPECT_NE(line[4], 0.0); // 0.0
512 EXPECT_NE(line[5], 0.0); // 1.0
513 EXPECT_NE(line[6], 0.0); // 2.0
514 EXPECT_NE(line[7], 0.0); // 3.0
515 EXPECT_EQ(line[8], 0.0); // 3.0
516}
517
518// Check that the field outside the aperture is not calculated
519TEST_F(TestMultipoleTStraight, VertAperture) {
520 std::vector<double> line(9);
521 // Set up the magnet
522 constexpr double length = 4;
523 constexpr double bendAngle = M_PI / 8.0;
524 constexpr double dipoleField = 1.0;
525 setBendAngle(bendAngle, false);
526 setElementLength(length);
527 setAperture(3.5, 3.5);
528 setFringeField(length / 2, 3, 3);
529 setRotation(0.0);
530 setEntranceAngle(0.0);
531 setMaxOrder(5, 10);
532 setTransProfile({dipoleField});
533 // Check field vanishes outside the aperture
534 grabVerticalDataLine(line, 0, 4.0, {0, 0, 0}, length);
535 EXPECT_EQ(line[0], 0.0); // -4.0
536 EXPECT_NE(line[1], 0.0); // -3.0
537 EXPECT_NE(line[2], 0.0); // -2.0
538 EXPECT_NE(line[3], 0.0); // -1.0
539 EXPECT_NE(line[4], 0.0); // 0.0
540 EXPECT_NE(line[5], 0.0); // 1.0
541 EXPECT_NE(line[6], 0.0); // 2.0
542 EXPECT_NE(line[7], 0.0); // 3.0
543 EXPECT_EQ(line[8], 0.0); // 3.0
544}
545
546// Test the single particle version of the apply function
547TEST_F(TestMultipoleTStraight, FieldAtSingleParticlePosition) {
548 // Set up the magnet
549 constexpr double length = 4.4;
550 constexpr double bendAngle = 0.0;
551 constexpr double dipoleField = 1.0;
552 setBendAngle(bendAngle, false);
553 setElementLength(length);
554 setAperture(3.5, 3.5);
555 setFringeField(length / 2, 0.3, 0.3);
556 setRotation(0.0);
557 setEntranceAngle(0.0);
558 setMaxOrder(5, 10);
559 setTransProfile({dipoleField});
560 // Make the bunch
561 const auto bunch = makeBunch(1);
562 const auto pc = bunch->getParticleContainer();
563 // Create the local views and data
564 std::vector<Vector_t<double, 3>> localR(1);
565 const auto hostR = Kokkos::create_mirror_view(pc->R.getView());
566 const auto hostB = Kokkos::create_mirror_view(pc->B.getView());
567 // Set the particle position
568 localR[0] = {0, 0, 0};
569 hostR(0) = curvilinearToGlobal(localR[0], {0, 0, 0}, length);
570 Kokkos::deep_copy(pc->R.getView(), hostR);
571 pc->setQ(pc->getChargePerParticle());
572 ippl::Comm->barrier();
573 Kokkos::fence();
574 // Register the bunch with the element
575 bunch->setT(0.0);
576 double startField, endField;
577 initialise(bunch.get(), startField, endField);
578 // Get the fields
579 Vector_t<double, 3> E{}, B{};
580 apply(0, 0.0, E, B);
581 // Check the field
582 const auto val = std::hypot(B[0], B[1], B[2]);
583 EXPECT_NEAR(val, dipoleField, 1e-2);
584 EXPECT_EQ(E[0], 0.0);
585 EXPECT_EQ(E[1], 0.0);
586 EXPECT_EQ(E[2], 0.0);
587}
gmsg
Inform * gmsg
Definition changes.cpp:7
Vector_t
ippl::Vector< T, Dim > Vector_t
Definition DistributionMoments.h:31
TEST_F
TEST_F(TestMultipoleTStraight, Dipole)
Definition TestMultipoleTStraight.cpp:216
Beam
Definition Beam.h:32
Beam::find
static Beam * find(const std::string &name)
Find named BEAM.
Definition Beam.cpp:290
MultipoleT::initialise
void initialise(PartBunch_t *bunch, double &startField, double &endField) override
Definition MultipoleT.cpp:167
MultipoleT::getGeometry
BGeometryBase & getGeometry() override
Definition MultipoleT.cpp:187
MultipoleT::apply
bool apply(const std::shared_ptr< ParticleContainer_t > &pc) override
Definition MultipoleT.cpp:49
Object::itsAttr
std::vector< Attribute > itsAttr
The object attributes.
Definition Object.h:210
OpalData::storeInputFn
void storeInputFn(const std::string &fn)
store opals input filename
Definition OpalData.cpp:561
OpalData::getInstance
static OpalData * getInstance()
Definition OpalData.cpp:193
TestMultipoleTStraight::dataSink_m
std::shared_ptr< DataSink > dataSink_m
Definition TestMultipoleTStraight.cpp:189
TestMultipoleTStraight::makeBunch
std::shared_ptr< PartBunch_t > makeBunch(const size_t numParticles)
Definition TestMultipoleTStraight.cpp:191
TestMultipoleTStraight::grabTransverseDataLine
void grabTransverseDataLine(std::vector< double > &line, const double s, const double width, const Vector_t< double, 3 > &elementEntry, const double elementLength)
Definition TestMultipoleTStraight.cpp:58
TestMultipoleTStraight::curvilinearToGlobal
static Vector_t< double, 3 > curvilinearToGlobal(const Vector_t< double, 3 > &local, const Vector_t< double, 3 > &elementEntry, const double elementLength)
Definition TestMultipoleTStraight.cpp:49
TestMultipoleTStraight::fsCmdBase_m
std::shared_ptr< FieldSolverCmd > fsCmdBase_m
Definition TestMultipoleTStraight.cpp:188
TestMultipoleTStraight::grabVerticalDataLine
void grabVerticalDataLine(std::vector< double > &line, const double s, const double height, const Vector_t< double, 3 > &elementEntry, const double elementLength)
Definition TestMultipoleTStraight.cpp:96
TestMultipoleTStraight::grabLongitudinalDivCurlLine
void grabLongitudinalDivCurlLine(std::vector< double > &fieldLine, std::vector< double > &divLine, std::vector< Vector_t< double, 3 > > &curlLine, const double x, const double length, const Vector_t< double, 3 > &elementEntry, const double elementLength, const double dr)
Definition TestMultipoleTStraight.cpp:134
Attributes::setPredefinedString
void setPredefinedString(Attribute &attr, const std::string &val)
Set predefined string value.
Definition Attributes.cpp:413
Options::enableHDF5
bool enableHDF5
If true HDF5 files are written.
Definition Options.cpp:83