OPALX (Object Oriented Parallel Accelerator Library for Exascal) master (dc2a29eed580)
OPALX
Loading...
Searching...
No Matches
TestMultipoleTCurvedConstRadius.cpp
Go to the documentation of this file.
1//
2// Copyright (c) 2026, Paul Scherrer Institute, Villigen PSI, Switzerland
3// All rights reserved
4//
5// This file is part of OPAL.
6//
7// OPAL is free software: you can redistribute it and/or modify
8// it under the terms of the GNU General Public License as published by
9// the Free Software Foundation, either version 3 of the License, or
10// (at your option) any later version.
11//
12// You should have received a copy of the GNU General License
13// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
14//
15
16#include <vector>
17#include "AbsBeamline/MultipoleT.h"
18#include "AbstractObjects/OpalData.h"
19#include "Structure/Beam.h"
20#include "Structure/DataSink.h"
21#include "gtest/gtest.h"
22
23class TestMultipoleTCurvedConstRadius : public testing::Test, public MultipoleT {
24public:
25 TestMultipoleTCurvedConstRadius() : MultipoleT("Magnet") {}
26
27 static void SetUpTestSuite() {
28 int argc = 0;
29 char** argv = nullptr;
30 ippl::initialize(argc, argv);
31 // DataSink requires a basename to create *.stat / *.lbal writers.
32 OpalData::getInstance()->storeInputFn("unit_test.opal");
33 // Many OPAL writers assume `gmsg` is initialized (see SDDSWriter/StatWriter).
34 // Unit tests normally don't set this up via Main().
35 gmsg = new Inform(nullptr, -1);
36 // DataSink::DataSink() constructs HDF5 writers when enabled, but the unit
37 // test doesn't have an H5PartWrapper. Disable HDF5 for this smoke test.
38 Options::enableHDF5 = false;
39 }
40 static void TearDownTestSuite() { ippl::finalize(); }
41
42 // Test helper functions
43 static Vector_t<double, 3> curvilinearToGlobal(
44 const Vector_t<double, 3>& local, const Vector_t<double, 3>& elementEntry,
45 const double elementLength, const double bendAngle) {
46 const auto radius = elementLength / bendAngle;
47 const double s = local[2] + elementLength / 2.0;
48 const auto cosSbyR = std::cos(s / radius);
49 const auto sinSbyR = std::sin(s / radius);
50 const double x = radius - (local[0] + radius) * cosSbyR + elementEntry[0];
51 const double y = local[1] + elementEntry[1];
52 const double z = (local[0] + radius) * sinSbyR + elementEntry[2];
53 return {x, y, z};
54 }
55
56 void grabTransverseDataLine(
57 std::vector<double>& line, const double s, const double width,
58 const Vector_t<double, 3>& elementEntry, const double elementLength,
59 const double bendAngle) {
60 // Make the bunch
61 const auto bunch = makeBunch(line.size());
62 const auto pc = bunch->getParticleContainer();
63 // Create the local views and data
64 std::vector<Vector_t<double, 3>> localR(line.size());
65 const auto hostR = Kokkos::create_mirror_view(pc->R.getView());
66 const auto hostB = Kokkos::create_mirror_view(pc->B.getView());
67 // Set the particle positions
68 const double stepSize = width / static_cast<double>(line.size() - 1);
69 for (size_t i = 0; i < line.size(); ++i) {
70 localR[i] = {static_cast<double>(i) * stepSize - width / 2, 0, s};
71 hostR(i) = curvilinearToGlobal(localR[i], elementEntry, elementLength, bendAngle);
72 }
73 Kokkos::deep_copy(pc->R.getView(), hostR);
74 pc->setQ(pc->getChargePerParticle());
75 ippl::Comm->barrier();
76 Kokkos::fence();
77 // Register the bunch with the element
78 bunch->setT(0.0);
79 double startField, endField;
80 initialise(bunch.get(), startField, endField);
81 // Get the fields
82 apply(pc);
83 // Return the fields
84 Kokkos::deep_copy(hostB, pc->B.getView());
85 Kokkos::fence();
86 for (size_t i = 0; i < line.size(); ++i) {
87 line[i] = std::hypot(hostB(i)[0], hostB(i)[1], hostB(i)[2]);
88 // std::cout << i << ": Local=" << local[i] << ", Global=" << R[i]
89 // << ", mag(B)=" << line[i] << std::endl;
90 }
91 // std::cout << std::endl;
92 }
93
94 void grabVerticalDataLine(
95 std::vector<double>& line, const double s, const double height,
96 const Vector_t<double, 3>& elementEntry, const double elementLength,
97 const double bendAngle) {
98 // Make the bunch
99 const auto bunch = makeBunch(line.size());
100 const auto pc = bunch->getParticleContainer();
101 // Create the local views and data
102 std::vector<Vector_t<double, 3>> localR(line.size());
103 const auto hostR = Kokkos::create_mirror_view(pc->R.getView());
104 const auto hostB = Kokkos::create_mirror_view(pc->B.getView());
105 // Set the particle positions
106 const double stepSize = height / static_cast<double>(line.size() - 1);
107 for (size_t i = 0; i < line.size(); ++i) {
108 localR[i] = {0, static_cast<double>(i) * stepSize - height / 2, s};
109 hostR(i) = curvilinearToGlobal(localR[i], elementEntry, elementLength, bendAngle);
110 }
111 Kokkos::deep_copy(pc->R.getView(), hostR);
112 pc->setQ(pc->getChargePerParticle());
113 ippl::Comm->barrier();
114 Kokkos::fence();
115 // Register the bunch with the element
116 bunch->setT(0.0);
117 double startField, endField;
118 initialise(bunch.get(), startField, endField);
119 // Get the fields
120 apply(pc);
121 // Return the fields
122 Kokkos::deep_copy(hostB, pc->B.getView());
123 Kokkos::fence();
124 for (size_t i = 0; i < line.size(); ++i) {
125 line[i] = std::hypot(hostB(i)[0], hostB(i)[1], hostB(i)[2]);
126 // std::cout << i << ": Local=" << local[i] << ", Global=" << R[i]
127 // << ", mag(B)=" << line[i] << std::endl;
128 }
129 // std::cout << std::endl;
130 }
131
132 void grabLongitudinalDivCurlLine(
133 std::vector<double>& fieldLine, std::vector<double>& divLine,
134 std::vector<Vector_t<double, 3>>& curlLine, const double x, const double length,
135 const Vector_t<double, 3>& elementEntry, const double elementLength,
136 const double bendAngle, const double dr) {
137 const double stepSize = length / static_cast<double>(divLine.size() - 1);
138 const double startS = 0 - (elementLength - length) / 2;
139 for (size_t i = 0; i < divLine.size(); ++i) {
140 // Get the sourrounding 6 B fields
141 Vector_t<double, 3> local{x, 0, static_cast<double>(i) * stepSize + startS};
142 Vector_t<double, 3> R =
143 curvilinearToGlobal(local, elementEntry, elementLength, bendAngle);
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 class TestableFieldSolverCmd : public FieldSolverCmd {
171 public:
172 void setType(const std::string& t) {
173 Attributes::setPredefinedString(this->itsAttr[FIELDSOLVER::TYPE], t);
174 }
175
176 void setBCX(const std::string& bc) {
177 Attributes::setPredefinedString(this->itsAttr[FIELDSOLVER::BCFFTX], bc);
178 }
179 void setBCY(const std::string& bc) {
180 Attributes::setPredefinedString(this->itsAttr[FIELDSOLVER::BCFFTY], bc);
181 }
182 void setBCZ(const std::string& bc) {
183 Attributes::setPredefinedString(this->itsAttr[FIELDSOLVER::BCFFTZ], bc);
184 }
185 };
186
187 std::shared_ptr<FieldSolverCmd> fsCmdBase_m;
188 std::shared_ptr<DataSink> dataSink_m;
189
190 std::shared_ptr<PartBunch_t> makeBunch(const size_t numParticles) {
191 dataSink_m = std::make_shared<DataSink>();
192 const auto fsCmd = std::make_shared<TestableFieldSolverCmd>();
193 fsCmdBase_m = fsCmd;
194 fsCmd->setType("NONE");
195 fsCmd->setNX(8);
196 fsCmd->setNY(8);
197 fsCmd->setNZ(8);
198 fsCmd->setBCX("PERIODIC");
199 fsCmd->setBCY("PERIODIC");
200 fsCmd->setBCZ("PERIODIC");
201 auto beam = std::make_shared<Beam>();
202 Beam* opBeam = Beam::find("UNNAMED_BEAM");
203 EXPECT_NE(opBeam, nullptr);
204 auto bunch = std::make_shared<PartBunch_t>(
205 /*qi=*/std::vector{1.0}, /*mi=*/std::vector{1.0},
206 /*beams=*/std::vector<Beam*>{opBeam},
207 /*totalParticlesPerBeam=*/std::vector<size_t>{numParticles},
208 /*lbt=*/1.0, /*integration_method=*/"LF2", fsCmdBase_m.get(), dataSink_m.get());
209 bunch->getParticleContainer()->createParticles(numParticles);
210 return bunch;
211 }
212};
213
214// Test transversely across the magnet to check for a dipole (constant) field.
215TEST_F(TestMultipoleTCurvedConstRadius, Dipole) {
216 std::vector<double> line(101);
217 // Set up the magnet
218 constexpr double length = 4.4;
219 constexpr double bendAngle = M_PI / 8.0;
220 constexpr double dipoleField = 1.0;
221 setBendAngle(bendAngle, false);
222 setElementLength(length);
223 setAperture(3.5, 3.5);
224 setFringeField(length / 2, 0.3, 0.3);
225 setRotation(0.0);
226 setEntranceAngle(0.0);
227 setMaxOrder(5, 10);
228 // Check dipole has the correct constant transverse field magnitude at the center
229 setTransProfile({dipoleField});
230 grabTransverseDataLine(line, 0, 3.0, {0, 0, 0}, length, bendAngle);
231 for (const double val : line) {
232 EXPECT_NEAR(val, dipoleField, 1e-2);
233 }
234 // Check dipole has the correct constant transverse field magnitude at the edge
235 setTransProfile({dipoleField});
236 grabTransverseDataLine(line, -length / 2, 3.0, {0, 0, 0}, length, bendAngle);
237 for (const double val : line) {
238 EXPECT_NEAR(val, dipoleField / 2, 1e-2);
239 }
240}
241
242// Test transversely across the magnet to check for a quadrupole (linear) field.
243TEST_F(TestMultipoleTCurvedConstRadius, Quadrupole) {
244 constexpr unsigned int samplesPerSide = 50;
245 std::vector<double> line(2 * samplesPerSide + 1);
246 // Set up the magnet
247 constexpr double length = 4.4;
248 constexpr double bendAngle = M_PI / 8.0;
249 constexpr double quadrupoleField = 1.0;
250 setBendAngle(bendAngle, false);
251 setElementLength(length);
252 setAperture(3.5, 3.5);
253 setFringeField(length / 2, 0.3, 0.3);
254 setRotation(0.0);
255 setEntranceAngle(0.0);
256 setMaxOrder(5, 10);
257 // Check quadrupole has linear field magnitude at the center
258 setTransProfile({0, quadrupoleField});
259 grabTransverseDataLine(line, 0, 3.0, {0, 0, 0}, length, bendAngle);
260 double expected = 0;
261 double delta = line[samplesPerSide + 1] - line[samplesPerSide];
262 EXPECT_NEAR(delta, 0.03, 1e-5);
263 EXPECT_NEAR(line[samplesPerSide], expected, 1e-5);
264 for (size_t i = 0; i < samplesPerSide; ++i) {
265 expected += delta;
266 EXPECT_NEAR(line[samplesPerSide - i - 1], expected, 1e-2);
267 EXPECT_NEAR(line[samplesPerSide + i + 1], expected, 1e-2);
268 }
269 // Check quadrupole has linear field magnitude at the edge with half the slope
270 setTransProfile({0, quadrupoleField});
271 grabTransverseDataLine(line, -length / 2, 3.0, {0, 0, 0}, length, bendAngle);
272 expected = 0;
273 delta = line[samplesPerSide + 1] - line[samplesPerSide];
274 EXPECT_NEAR(delta, 0.015, 1e-5);
275 EXPECT_NEAR(line[samplesPerSide], expected, 1e-5);
276 for (size_t i = 0; i < samplesPerSide; ++i) {
277 expected += delta;
278 EXPECT_NEAR(line[samplesPerSide - i - 1], expected, 1e-2);
279 EXPECT_NEAR(line[samplesPerSide + i + 1], expected, 1e-2);
280 }
281}
282
283// Test transversely across the magnet to check for a sextupole (quadratic) field.
284TEST_F(TestMultipoleTCurvedConstRadius, Sextupole) {
285 constexpr unsigned int samplesPerSide = 50;
286 std::vector<double> line(2 * samplesPerSide + 1);
287 // Set up the magnet
288 constexpr double length = 4.4;
289 constexpr double bendAngle = M_PI / 8.0;
290 constexpr double sextupoleField = 1.0;
291 setBendAngle(bendAngle, false);
292 setElementLength(length);
293 setAperture(3.5, 3.5);
294 setFringeField(length / 2, 0.3, 0.3);
295 setRotation(0.0);
296 setEntranceAngle(0.0);
297 setMaxOrder(5, 10);
298 // Check sextupole has quadratic field magnitude at the center
299 setTransProfile({0, 0, sextupoleField});
300 grabTransverseDataLine(line, 0, 3.0, {0, 0, 0}, length, bendAngle);
301 double expected = 0;
302 double delta = std::sqrt(line[samplesPerSide + 1] - line[samplesPerSide]);
303 EXPECT_NEAR(delta, 0.03, 1e-5);
304 EXPECT_NEAR(line[samplesPerSide], expected, 1e-5);
305 for (size_t i = 0; i < samplesPerSide; ++i) {
306 expected = std::pow(static_cast<double>(i + 1) * delta, 2);
307 EXPECT_NEAR(line[samplesPerSide - i - 1], expected, 1e-2);
308 EXPECT_NEAR(line[samplesPerSide + i + 1], expected, 1e-2);
309 }
310 // Check sextupole has quadratic field magnitude at the edge with half the slope
311 setTransProfile({0, 0, sextupoleField});
312 grabTransverseDataLine(line, -length / 2, 3.0, {0, 0, 0}, length, bendAngle);
313 expected = 0;
314 delta = std::sqrt(line[samplesPerSide + 1] - line[samplesPerSide]);
315 EXPECT_NEAR(delta, 0.02121, 1e-5);
316 EXPECT_NEAR(line[samplesPerSide], expected, 1e-5);
317 for (size_t i = 0; i < samplesPerSide; ++i) {
318 expected = std::pow(static_cast<double>(i + 1) * delta, 2);
319 EXPECT_NEAR(line[samplesPerSide - i - 1], expected, 1e-2);
320 EXPECT_NEAR(line[samplesPerSide + i + 1], expected, 1e-2);
321 }
322}
323
324// Test transversely across the magnet to check for an octupole (cubic) field.
325TEST_F(TestMultipoleTCurvedConstRadius, Octupole) {
326 constexpr unsigned int samplesPerSide = 50;
327 std::vector<double> line(2 * samplesPerSide + 1);
328 // Set up the magnet
329 constexpr double length = 4.4;
330 constexpr double bendAngle = M_PI / 8.0;
331 constexpr double sextupoleField = 1.0;
332 setBendAngle(bendAngle, false);
333 setElementLength(length);
334 setAperture(3.5, 3.5);
335 setFringeField(length / 2, 0.3, 0.3);
336 setRotation(0.0);
337 setEntranceAngle(0.0);
338 setMaxOrder(5, 10);
339 // Check octupole has cubic field magnitude at the center
340 setTransProfile({0, 0, 0, sextupoleField});
341 grabTransverseDataLine(line, 0, 3.0, {0, 0, 0}, length, bendAngle);
342 double expected = 0;
343 double delta = std::cbrt(line[samplesPerSide + 1] - line[samplesPerSide]);
344 EXPECT_NEAR(delta, 0.03, 1e-5);
345 EXPECT_NEAR(line[samplesPerSide], expected, 1e-5);
346 for (size_t i = 0; i < samplesPerSide; ++i) {
347 expected = std::pow(static_cast<double>(i + 1) * delta, 3);
348 EXPECT_NEAR(line[samplesPerSide - i - 1], expected, 1e-2);
349 EXPECT_NEAR(line[samplesPerSide + i + 1], expected, 1e-2);
350 }
351 // Check octupole has cubic field magnitude at the edge with half the slope
352 setTransProfile({0, 0, 0, sextupoleField});
353 grabTransverseDataLine(line, -length / 2, 3.0, {0, 0, 0}, length, bendAngle);
354 expected = 0;
355 delta = std::cbrt(line[samplesPerSide + 1] - line[samplesPerSide]);
356 EXPECT_NEAR(delta, 0.023811, 1e-5);
357 EXPECT_NEAR(line[samplesPerSide], expected, 1e-5);
358 for (size_t i = 0; i < samplesPerSide; ++i) {
359 expected = std::pow(static_cast<double>(i + 1) * delta, 3);
360 EXPECT_NEAR(line[samplesPerSide - i - 1], expected, 1e-2);
361 EXPECT_NEAR(line[samplesPerSide + i + 1], expected, 1e-2);
362 }
363}
364
365// Test transversely across the magnet to check for a decapole (quartic) field.
366TEST_F(TestMultipoleTCurvedConstRadius, Decapole) {
367 constexpr unsigned int samplesPerSide = 50;
368 std::vector<double> line(2 * samplesPerSide + 1);
369 // Set up the magnet
370 constexpr double length = 4.4;
371 constexpr double bendAngle = M_PI / 8.0;
372 constexpr double decapoleField = 1.0;
373 setBendAngle(bendAngle, false);
374 setElementLength(length);
375 setAperture(3.5, 3.5);
376 setFringeField(length / 2, 0.3, 0.3);
377 setRotation(0.0);
378 setEntranceAngle(0.0);
379 setMaxOrder(5, 10);
380 // Check decapole has quartic field magnitude at the center
381 setTransProfile({0, 0, 0, 0, decapoleField});
382 grabTransverseDataLine(line, 0, 3.0, {0, 0, 0}, length, bendAngle);
383 double expected = 0;
384 double delta = std::pow(line[samplesPerSide + 1] - line[samplesPerSide], 1.0 / 4.0);
385 EXPECT_NEAR(delta, 0.03, 1e-5);
386 EXPECT_NEAR(line[samplesPerSide], expected, 1e-5);
387 for (size_t i = 0; i < samplesPerSide; ++i) {
388 expected = std::pow(static_cast<double>(i + 1) * delta, 4);
389 EXPECT_NEAR(line[samplesPerSide - i - 1], expected, 1e-2);
390 EXPECT_NEAR(line[samplesPerSide + i + 1], expected, 1e-2);
391 }
392 // Check decapole has quartic field magnitude at the edge with half the slope
393 setTransProfile({0, 0, 0, 0, decapoleField});
394 grabTransverseDataLine(line, -length / 2, 3.0, {0, 0, 0}, length, bendAngle);
395 expected = 0;
396 delta = std::pow(line[samplesPerSide + 1] - line[samplesPerSide], 1.0 / 4.0);
397 EXPECT_NEAR(delta, 0.025227, 1e-5);
398 EXPECT_NEAR(line[samplesPerSide], expected, 1e-5);
399 for (size_t i = 0; i < samplesPerSide; ++i) {
400 expected = std::pow(static_cast<double>(i + 1) * delta, 4);
401 EXPECT_NEAR(line[samplesPerSide - i - 1], expected, 1e-2);
402 EXPECT_NEAR(line[samplesPerSide + i + 1], expected, 1e-2);
403 }
404}
405
406// Check that the field along the center line is Maxwellian (zero div and curl)
407TEST_F(TestMultipoleTCurvedConstRadius, DivCurl) {
408 constexpr unsigned int samplesPerSide = 20;
409 std::vector<double> divLine(2 * samplesPerSide + 1);
410 std::vector<double> fieldLine(2 * samplesPerSide + 1);
411 std::vector<Vector_t<double, 3>> curlLine(2 * samplesPerSide + 1);
412 constexpr double dr = 0.001;
413 // Set up the magnet
414 constexpr double length = 4.4;
415 constexpr double bendAngle = M_PI / 8.0;
416 setBendAngle(bendAngle, false);
417 setElementLength(length);
418 setAperture(3.5, 3.5);
419 setFringeField(length / 2, 0.3, 0.3);
420 setRotation(0.0);
421 setEntranceAngle(0.0);
422 setMaxOrder(5, 10);
423 setTransProfile({1, 1, 1, 1, 1});
424 // Get the data
425 grabLongitudinalDivCurlLine(
426 fieldLine, divLine, curlLine, 0, length * 1.5, {}, length, bendAngle, 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(TestMultipoleTCurvedConstRadius, Geometry) {
439 TestMultipoleTCurvedConstRadius* magnet = this;
440 auto* geom = &magnet->getGeometry();
441 EXPECT_NE(geom, nullptr);
442 auto* constGeom = &const_cast<const TestMultipoleTCurvedConstRadius*>(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(TestMultipoleTCurvedConstRadius, BoundingBox) {
449 std::vector<double> line(3);
450 // Set up the magnet
451 constexpr double length = 4;
452 constexpr double bendAngle = M_PI / 8.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, bendAngle);
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, bendAngle);
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, bendAngle);
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, bendAngle);
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, bendAngle);
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(TestMultipoleTCurvedConstRadius, HorzAperture) {
492 std::vector<double> line(9);
493 // Set up the magnet
494 constexpr double length = 4;
495 constexpr double bendAngle = M_PI / 8.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, bendAngle);
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(TestMultipoleTCurvedConstRadius, 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, bendAngle);
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}
gmsg
Inform * gmsg
Definition changes.cpp:7
Vector_t
ippl::Vector< T, Dim > Vector_t
Definition DistributionMoments.h:31
TEST_F
TEST_F(TestMultipoleTCurvedConstRadius, Dipole)
Definition TestMultipoleTCurvedConstRadius.cpp:215
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
TestMultipoleTCurvedConstRadius::curvilinearToGlobal
static Vector_t< double, 3 > curvilinearToGlobal(const Vector_t< double, 3 > &local, const Vector_t< double, 3 > &elementEntry, const double elementLength, const double bendAngle)
Definition TestMultipoleTCurvedConstRadius.cpp:43
TestMultipoleTCurvedConstRadius::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 bendAngle, const double dr)
Definition TestMultipoleTCurvedConstRadius.cpp:132
TestMultipoleTCurvedConstRadius::grabTransverseDataLine
void grabTransverseDataLine(std::vector< double > &line, const double s, const double width, const Vector_t< double, 3 > &elementEntry, const double elementLength, const double bendAngle)
Definition TestMultipoleTCurvedConstRadius.cpp:56
TestMultipoleTCurvedConstRadius::grabVerticalDataLine
void grabVerticalDataLine(std::vector< double > &line, const double s, const double height, const Vector_t< double, 3 > &elementEntry, const double elementLength, const double bendAngle)
Definition TestMultipoleTCurvedConstRadius.cpp:94
TestMultipoleTCurvedConstRadius::fsCmdBase_m
std::shared_ptr< FieldSolverCmd > fsCmdBase_m
Definition TestMultipoleTCurvedConstRadius.cpp:187
TestMultipoleTCurvedConstRadius::makeBunch
std::shared_ptr< PartBunch_t > makeBunch(const size_t numParticles)
Definition TestMultipoleTCurvedConstRadius.cpp:190
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