OPAL (Object Oriented Parallel Accelerator Library) 2022.1
OPAL
FM2DMagnetoStatic.cpp
Go to the documentation of this file.
1#include "Fields/FM2DMagnetoStatic.h"
2#include "Fields/Fieldmap.hpp"
3#include "Physics/Units.h"
4#include "Utilities/GeneralClassicException.h"
5#include "Utilities/Util.h"
6
7#include <fstream>
8#include <ios>
9#include <cmath>
10
11FM2DMagnetoStatic::FM2DMagnetoStatic(std::string aFilename):
12 Fieldmap(aFilename),
13 FieldstrengthBz_m(nullptr),
14 FieldstrengthBr_m(nullptr) {
15 std::ifstream file;
16 std::string tmpString;
17 double tmpDouble;
18
19 Type = T2DMagnetoStatic;
20
21 // open field map, parse it and disable element on error
22 file.open(Filename_m.c_str());
23 if(file.good()) {
24 bool parsing_passed = true;
25 try {
26 parsing_passed = interpretLine<std::string, std::string>(file,
27 tmpString,
28 tmpString);
29 } catch (GeneralClassicException &e) {
30 parsing_passed = interpretLine<std::string, std::string, std::string>(file,
31 tmpString,
32 tmpString,
33 tmpString);
34
35 tmpString = Util::toUpper(tmpString);
36 if (tmpString != "TRUE" &&
37 tmpString != "FALSE")
38 throw GeneralClassicException("FM2DMagnetoStatic::FM2DMagnetoStatic",
39 "The third string on the first line of 2D field "
40 "maps has to be either TRUE or FALSE");
41
42 normalize_m = (tmpString == "TRUE");
43 }
44
45 if(tmpString == "ZX") {
46 swap_m = true;
47 parsing_passed = parsing_passed &&
48 interpretLine<double, double, int>(file, rbegin_m, rend_m, num_gridpr_m);
49 parsing_passed = parsing_passed &&
50 interpretLine<double, double, int>(file, zbegin_m, zend_m, num_gridpz_m);
51 } else if(tmpString == "XZ") {
52 swap_m = false;
53 parsing_passed = parsing_passed &&
54 interpretLine<double, double, int>(file, zbegin_m, zend_m, num_gridpz_m);
55 parsing_passed = parsing_passed &&
56 interpretLine<double, double, int>(file, rbegin_m, rend_m, num_gridpr_m);
57 } else {
58 std::cerr << "unknown orientation of 2D magnetostatic fieldmap" << std::endl;
59 parsing_passed = false;
60 }
61
62 for(long i = 0; (i < (num_gridpz_m + 1) * (num_gridpr_m + 1)) && parsing_passed; ++ i) {
63 parsing_passed = parsing_passed && interpretLine<double, double>(file, tmpDouble, tmpDouble);
64 }
65
66 parsing_passed = parsing_passed &&
67 interpreteEOF(file);
68
69 file.close();
70 lines_read_m = 0;
71
72 if(!parsing_passed) {
73 disableFieldmapWarning();
74 zend_m = zbegin_m - 1e-3;
75 throw GeneralClassicException("FM2DMagnetoStatic::FM2DMagnetoStatic",
76 "An error occured when reading the fieldmap '" + Filename_m + "'");
77 } else {
78 // conversion from cm to m
79 rbegin_m *= Units::cm2m;
80 rend_m *= Units::cm2m;
81 zbegin_m *= Units::cm2m;
82 zend_m *= Units::cm2m;
83
84 hr_m = (rend_m - rbegin_m) / num_gridpr_m;
85 hz_m = (zend_m - zbegin_m) / num_gridpz_m;
86
87 // num spacings -> num grid points
88 num_gridpr_m++;
89 num_gridpz_m++;
90 }
91 } else {
92 noFieldmapWarning();
93 zbegin_m = 0.0;
94 zend_m = -1e-3;
95 }
96}
97
98FM2DMagnetoStatic::~FM2DMagnetoStatic() {
99 freeMap();
100}
101
102void FM2DMagnetoStatic::readMap() {
103 if(FieldstrengthBz_m == nullptr) {
104 // declare variables and allocate memory
105 std::ifstream in;
106 std::string tmpString;
107
108 FieldstrengthBz_m = new double[num_gridpz_m * num_gridpr_m];
109 FieldstrengthBr_m = new double[num_gridpz_m * num_gridpr_m];
110
111 // read in and parse field map
112 in.open(Filename_m.c_str());
113 getLine(in, tmpString);
114 getLine(in, tmpString);
115 getLine(in, tmpString);
116
117 if(swap_m) {
118 for(int i = 0; i < num_gridpz_m; i++) {
119 for(int j = 0; j < num_gridpr_m; j++) {
120 interpretLine<double, double>(in,
121 FieldstrengthBr_m[i + j * num_gridpz_m],
122 FieldstrengthBz_m[i + j * num_gridpz_m]);
123 }
124 }
125 } else {
126 for(int j = 0; j < num_gridpr_m; j++) {
127 for(int i = 0; i < num_gridpz_m; i++) {
128 interpretLine<double, double>(in,
129 FieldstrengthBz_m[i + j * num_gridpz_m],
130 FieldstrengthBr_m[i + j * num_gridpz_m]);
131 }
132 }
133 }
134 in.close();
135
136 if (normalize_m) {
137 double Bzmax = 0.0;
138 // find maximum field
139 for(int i = 0; i < num_gridpz_m; ++ i) {
140 if(std::abs(FieldstrengthBz_m[i]) > Bzmax) {
141 Bzmax = std::abs(FieldstrengthBz_m[i]);
142 }
143 }
144
145 // normalize field
146 for(int i = 0; i < num_gridpr_m * num_gridpz_m; ++ i) {
147 FieldstrengthBz_m[i] /= Bzmax;
148 FieldstrengthBr_m[i] /= Bzmax;
149 }
150 }
151
152 INFOMSG(level3 << typeset_msg("read in fieldmap '" + Filename_m + "'", "info") << endl);
153 }
154}
155
156void FM2DMagnetoStatic::freeMap() {
157 if(FieldstrengthBz_m != nullptr) {
158 delete[] FieldstrengthBz_m;
159 delete[] FieldstrengthBr_m;
160
161 FieldstrengthBz_m = nullptr;
162 FieldstrengthBr_m = nullptr;
163
164 INFOMSG(level3 << typeset_msg("freed fieldmap '" + Filename_m + "'", "info") << endl);
165 }
166}
167
168bool FM2DMagnetoStatic::getFieldstrength(const Vector_t &R, Vector_t &/*E*/, Vector_t &B) const {
169 // do bi-linear interpolation
170 const double RR = std::sqrt(R(0) * R(0) + R(1) * R(1));
171
172 const int indexr = std::abs((int)std::floor(RR / hr_m));
173 const double leverr = (RR / hr_m) - indexr;
174
175 const int indexz = std::abs((int)std::floor((R(2) - zbegin_m) / hz_m));
176 const double leverz = (R(2) - zbegin_m) / hz_m - indexz;
177
178 if((indexz < 0) || (indexz + 2 > num_gridpz_m))
179 return false;
180
181 if(indexr + 2 > num_gridpr_m)
182 return true;
183
184 const int index1 = indexz + indexr * num_gridpz_m;
185 const int index2 = index1 + num_gridpz_m;
186 const double BfieldR = (1.0 - leverz) * (1.0 - leverr) * FieldstrengthBr_m[index1]
187 + leverz * (1.0 - leverr) * FieldstrengthBr_m[index1 + 1]
188 + (1.0 - leverz) * leverr * FieldstrengthBr_m[index2]
189 + leverz * leverr * FieldstrengthBr_m[index2 + 1];
190
191 const double BfieldZ = (1.0 - leverz) * (1.0 - leverr) * FieldstrengthBz_m[index1]
192 + leverz * (1.0 - leverr) * FieldstrengthBz_m[index1 + 1]
193 + (1.0 - leverz) * leverr * FieldstrengthBz_m[index2]
194 + leverz * leverr * FieldstrengthBz_m[index2 + 1];
195
196 if(RR != 0) {
197 B(0) += BfieldR * R(0) / RR;
198 B(1) += BfieldR * R(1) / RR;
199 }
200 B(2) += BfieldZ;
201 return false;
202}
203
204bool FM2DMagnetoStatic::getFieldDerivative(const Vector_t &R, Vector_t &/*E*/, Vector_t &B, const DiffDirection &dir) const {
205
206 double BfieldR, BfieldZ;
207
208 const double RR = std::sqrt(R(0) * R(0) + R(1) * R(1));
209
210 const int indexr = (int)std::floor(RR / hr_m);
211 const double leverr = (RR / hr_m) - indexr;
212
213 const int indexz = (int)std::floor((R(2)) / hz_m);
214 const double leverz = (R(2) / hz_m) - indexz;
215
216 if((indexz < 0) || (indexz + 2 > num_gridpz_m))
217 return false;
218
219 if(indexr + 2 > num_gridpr_m)
220 return true;
221
222 const int index1 = indexz + indexr * num_gridpz_m;
223 const int index2 = index1 + num_gridpz_m;
224 if(dir == DZ) {
225 if(indexz > 0) {
226 if(indexz < num_gridpz_m - 1) {
227 BfieldR = (1.0 - leverr) * ((FieldstrengthBr_m[index1 - 1] - 2. * FieldstrengthBr_m[index1] + FieldstrengthBr_m[index1 + 1]) * (R(2) - indexz * hz_m) / (hz_m * hz_m)
228 + (FieldstrengthBr_m[index1 + 1] - FieldstrengthBr_m[index1 - 1]) / (2.*hz_m))
229 + leverr * ((FieldstrengthBr_m[index2 - 1] - 2. * FieldstrengthBr_m[index2] + FieldstrengthBr_m[index2 + 1]) * (R(2) - indexz * hz_m) / (hz_m * hz_m)
230 + (FieldstrengthBr_m[index2 + 1] - FieldstrengthBr_m[index2 - 1]) / (2.*hz_m));
231
232 BfieldZ = (1.0 - leverr) * ((FieldstrengthBz_m[index1 - 1] - 2. * FieldstrengthBz_m[index1] + FieldstrengthBz_m[index1 + 1]) * (R(2) - indexz * hz_m) / (hz_m * hz_m)
233 + (FieldstrengthBz_m[index1 + 1] - FieldstrengthBz_m[index1 - 1]) / (2.*hz_m))
234 + leverr * ((FieldstrengthBz_m[index2 - 1] - 2. * FieldstrengthBz_m[index2] + FieldstrengthBz_m[index2 + 1]) * (R(2) - indexz * hz_m) / (hz_m * hz_m)
235 + (FieldstrengthBz_m[index2 + 1] - FieldstrengthBz_m[index2 - 1]) / (2.*hz_m));
236 } else {
237 BfieldR = (1.0 - leverr) * ((FieldstrengthBr_m[index1] - 2. * FieldstrengthBr_m[index1 - 1] + FieldstrengthBr_m[index1 - 2]) * (R(2) - indexz * hz_m) / (hz_m * hz_m)
238 - (3. * FieldstrengthBr_m[index1] - 4. * FieldstrengthBr_m[index1 - 1] + FieldstrengthBr_m[index1 - 2]) / (2.*hz_m))
239 + leverr * ((FieldstrengthBr_m[index2] - 2. * FieldstrengthBr_m[index2 - 1] + FieldstrengthBr_m[index2 - 2]) * (R(2) - indexz * hz_m) / (hz_m * hz_m)
240 - (3. * FieldstrengthBr_m[index2] - 4. * FieldstrengthBr_m[index1 - 1] + FieldstrengthBr_m[index2 - 2]) / (2.*hz_m));
241
242 BfieldZ = (1.0 - leverr) * ((FieldstrengthBz_m[index1] - 2. * FieldstrengthBz_m[index1 - 1] + FieldstrengthBz_m[index1 - 2]) * (R(2) - indexz * hz_m) / (hz_m * hz_m)
243 - (3. * FieldstrengthBz_m[index1] - 4. * FieldstrengthBr_m[index1 - 1] + FieldstrengthBz_m[index1 - 2]) / (2.*hz_m))
244 + leverr * ((FieldstrengthBz_m[index2] - 2. * FieldstrengthBz_m[index2 - 1] + FieldstrengthBz_m[index2 - 2]) * (R(2) - indexz * hz_m) / (hz_m * hz_m)
245 - (3. * FieldstrengthBz_m[index2] - 4. * FieldstrengthBr_m[index1 - 1] + FieldstrengthBz_m[index2 - 2]) / (2.*hz_m));
246 }
247 } else {
248 BfieldR = (1.0 - leverr) * ((FieldstrengthBr_m[index1] - 2. * FieldstrengthBr_m[index1 + 1] + FieldstrengthBr_m[index1 + 2]) * (R(2) - indexz * hz_m) / (hz_m * hz_m)
249 - (3. * FieldstrengthBr_m[index1] - 4. * FieldstrengthBr_m[index1 + 1] + FieldstrengthBr_m[index1 + 2]) / (2.*hz_m))
250 + leverr * ((FieldstrengthBr_m[index2] - 2. * FieldstrengthBr_m[index2 + 1] + FieldstrengthBr_m[index2 + 2]) * (R(2) - indexz * hz_m) / (hz_m * hz_m)
251 - (3. * FieldstrengthBr_m[index2] - 4. * FieldstrengthBr_m[index1 + 1] + FieldstrengthBr_m[index2 + 2]) / (2.*hz_m));
252
253 BfieldZ = (1.0 - leverr) * ((FieldstrengthBz_m[index1] - 2. * FieldstrengthBz_m[index1 + 1] + FieldstrengthBz_m[index1 + 2]) * (R(2) - indexz * hz_m) / (hz_m * hz_m)
254 - (3. * FieldstrengthBz_m[index1] - 4. * FieldstrengthBr_m[index1 + 1] + FieldstrengthBz_m[index1 + 2]) / (2.*hz_m))
255 + leverr * ((FieldstrengthBz_m[index2] - 2. * FieldstrengthBz_m[index2 + 1] + FieldstrengthBz_m[index2 + 2]) * (R(2) - indexz * hz_m) / (hz_m * hz_m)
256 - (3. * FieldstrengthBz_m[index2] - 4. * FieldstrengthBr_m[index1 + 1] + FieldstrengthBz_m[index2 + 2]) / (2.*hz_m));
257 }
258 } else {
259 if(indexr > 0) {
260 const int index_1 = index1 - num_gridpz_m;
261 if(indexr < num_gridpr_m - 1) {
262 BfieldR = (1.0 - leverz) * ((FieldstrengthBr_m[index_1 + 1] - 2. * FieldstrengthBr_m[index1 + 1] + FieldstrengthBr_m[index2 + 1]) * (RR - indexr * hr_m) / (hr_m * hr_m)
263 + (FieldstrengthBr_m[index2 + 1] - FieldstrengthBr_m[index_1 + 1]) / (2.*hr_m))
264 + leverz * ((FieldstrengthBr_m[index_1] - 2. * FieldstrengthBr_m[index1] + FieldstrengthBr_m[index2]) * (RR - indexr * hr_m) / (hr_m * hr_m)
265 + (FieldstrengthBr_m[index2] - FieldstrengthBr_m[index_1]) / (2.*hr_m));
266
267 BfieldZ = (1.0 - leverz) * ((FieldstrengthBz_m[index_1 + 1] - 2. * FieldstrengthBz_m[index1 + 1] + FieldstrengthBz_m[index2 + 1]) * (RR - indexr * hr_m) / (hr_m * hr_m)
268 + (FieldstrengthBz_m[index2 + 1] - FieldstrengthBz_m[index_1 + 1]) / (2.*hr_m))
269 + leverz * ((FieldstrengthBz_m[index_1] - 2. * FieldstrengthBz_m[index1] + FieldstrengthBz_m[index2]) * (RR - indexr * hr_m) / (hr_m * hr_m)
270 + (FieldstrengthBz_m[index2] - FieldstrengthBz_m[index_1]) / (2.*hr_m));
271 } else {
272 const int index_2 = index_1 - num_gridpz_m;
273 BfieldR = (1.0 - leverz) * ((FieldstrengthBr_m[index1 + 1] - 2. * FieldstrengthBr_m[index_1 + 1] + FieldstrengthBr_m[index_2 + 1]) * (RR - indexr * hr_m) / (hr_m * hr_m)
274 - (3. * FieldstrengthBr_m[index1 + 1] - 4. * FieldstrengthBr_m[index_1 + 1] + FieldstrengthBr_m[index_2 + 1]) / (2.*hr_m))
275 + leverz * ((FieldstrengthBr_m[index1] - 2. * FieldstrengthBr_m[index_1] + FieldstrengthBr_m[index_2]) * (RR - indexr * hr_m) / (hr_m * hr_m)
276 - (3. * FieldstrengthBr_m[index1] - 4. * FieldstrengthBr_m[index_1] + FieldstrengthBr_m[index_2]) / (2.*hr_m));
277
278 BfieldZ = (1.0 - leverz) * ((FieldstrengthBz_m[index1 + 1] - 2. * FieldstrengthBz_m[index_1 + 1] + FieldstrengthBz_m[index_2 + 1]) * (RR - indexr * hr_m) / (hr_m * hr_m)
279 - (3. * FieldstrengthBz_m[index1 + 1] - 4. * FieldstrengthBr_m[index_1 + 1] + FieldstrengthBz_m[index_2 + 1]) / (2.*hr_m))
280 + leverz * ((FieldstrengthBz_m[index1] - 2. * FieldstrengthBz_m[index_1] + FieldstrengthBz_m[index_2]) * (RR - indexr * hr_m) / (hr_m * hr_m)
281 - (3. * FieldstrengthBz_m[index1] - 4. * FieldstrengthBr_m[index_1] + FieldstrengthBz_m[index_2]) / (2.*hr_m));
282 }
283 } else {
284 const int index3 = index2 + num_gridpz_m;
285 BfieldR = (1.0 - leverz) * ((FieldstrengthBr_m[index1 + 1] - 2. * FieldstrengthBr_m[index2 + 1] + FieldstrengthBr_m[index3 + 1]) * (RR - indexr * hr_m) / (hr_m * hr_m)
286 - (3. * FieldstrengthBr_m[index1 + 1] - 4. * FieldstrengthBr_m[index2 + 1] + FieldstrengthBr_m[index3 + 1]) / (2.*hr_m))
287 + leverz * ((FieldstrengthBr_m[index1] - 2. * FieldstrengthBr_m[index2] + FieldstrengthBr_m[index3]) * (RR - indexr * hr_m) / (hr_m * hr_m)
288 - (3. * FieldstrengthBr_m[index1] - 4. * FieldstrengthBr_m[index2] + FieldstrengthBr_m[index3]) / (2.*hr_m));
289
290 BfieldZ = (1.0 - leverz) * ((FieldstrengthBz_m[index1 + 1] - 2. * FieldstrengthBz_m[index2 + 1] + FieldstrengthBz_m[index3 + 1]) * (RR - indexr * hr_m) / (hr_m * hr_m)
291 - (3. * FieldstrengthBz_m[index1 + 1] - 4. * FieldstrengthBr_m[index2 + 1] + FieldstrengthBz_m[index3 + 1]) / (2.*hr_m))
292 + leverz * ((FieldstrengthBz_m[index1] - 2. * FieldstrengthBz_m[index2] + FieldstrengthBz_m[index3]) * (RR - indexr * hr_m) / (hr_m * hr_m)
293 - (3. * FieldstrengthBz_m[index1] - 4. * FieldstrengthBr_m[index2] + FieldstrengthBz_m[index3]) / (2.*hr_m));
294 }
295
296 }
297
298 if(RR > 1.e-12) {
299 B(0) += BfieldR * R(0) / RR;
300 B(1) += BfieldR * R(1) / RR;
301 }
302 B(2) += BfieldZ;
303 return false;
304}
305
306void FM2DMagnetoStatic::getFieldDimensions(double &zBegin, double &zEnd) const {
307 zBegin = zbegin_m;
308 zEnd = zend_m;
309}
310
311void FM2DMagnetoStatic::getFieldDimensions(double &/*xIni*/, double &/*xFinal*/, double &/*yIni*/, double &/*yFinal*/, double &/*zIni*/, double &/*zFinal*/) const {}
312
313void FM2DMagnetoStatic::swap() {
314 if(swap_m) swap_m = false;
315 else swap_m = true;
316}
317
318void FM2DMagnetoStatic::getInfo(Inform *msg) {
319 (*msg) << Filename_m << " (2D magnetostatic); zini= " << zbegin_m << " m; zfinal= " << zend_m << " m;" << endl;
320}
321
322double FM2DMagnetoStatic::getFrequency() const {
323 return 0.0;
324}
325
326void FM2DMagnetoStatic::setFrequency(double /*freq*/)
327{ ;}
T2DMagnetoStatic
@ T2DMagnetoStatic
Definition: Fieldmap.h:28
DiffDirection
DiffDirection
Definition: Fieldmap.h:54
DZ
@ DZ
Definition: Fieldmap.h:57
sqrt
Tps< T > sqrt(const Tps< T > &x)
Square root.
Definition: TpsMath.h:91
floor
PETE_TUTree< FnFloor, typename T::PETE_Expr_t > floor(const PETE_Expr< T > &l)
Definition: PETE.h:733
abs
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
Definition: IpplExpressions.h:62
endl
Inform & endl(Inform &inf)
Definition: Inform.cpp:42
level3
Inform & level3(Inform &inf)
Definition: Inform.cpp:47
INFOMSG
#define INFOMSG(msg)
Definition: IpplInfo.h:348
Physics::e
constexpr double e
The value of.
Definition: Physics.h:39
Units::cm2m
constexpr double cm2m
Definition: Units.h:35
Util::toUpper
std::string toUpper(const std::string &str)
Definition: Util.cpp:146
Fieldmap::Type
MapType Type
Definition: Fieldmap.h:114
Fieldmap::interpreteEOF
bool interpreteEOF(std::ifstream &in)
Definition: Fieldmap.cpp:558
Fieldmap::disableFieldmapWarning
void disableFieldmapWarning()
Definition: Fieldmap.cpp:613
Fieldmap::normalize_m
bool normalize_m
Definition: Fieldmap.h:120
Fieldmap::lines_read_m
int lines_read_m
Definition: Fieldmap.h:118
Fieldmap::typeset_msg
static std::string typeset_msg(const std::string &msg, const std::string &title)
Definition: Fieldmap.cpp:652
Fieldmap::Filename_m
std::string Filename_m
Definition: Fieldmap.h:117
Fieldmap::getLine
void getLine(std::ifstream &in, std::string &buffer)
Definition: Fieldmap.h:121
Fieldmap::noFieldmapWarning
void noFieldmapWarning()
Definition: Fieldmap.cpp:621
FM2DMagnetoStatic::readMap
virtual void readMap()
Definition: FM2DMagnetoStatic.cpp:102
FM2DMagnetoStatic::getFieldstrength
virtual bool getFieldstrength(const Vector_t &R, Vector_t &E, Vector_t &B) const
Definition: FM2DMagnetoStatic.cpp:168
FM2DMagnetoStatic::getFrequency
virtual double getFrequency() const
Definition: FM2DMagnetoStatic.cpp:322
FM2DMagnetoStatic::getInfo
virtual void getInfo(Inform *msg)
Definition: FM2DMagnetoStatic.cpp:318
FM2DMagnetoStatic::getFieldDimensions
virtual void getFieldDimensions(double &zBegin, double &zEnd) const
Definition: FM2DMagnetoStatic.cpp:306
FM2DMagnetoStatic::setFrequency
virtual void setFrequency(double freq)
Definition: FM2DMagnetoStatic.cpp:326
FM2DMagnetoStatic::getFieldDerivative
virtual bool getFieldDerivative(const Vector_t &R, Vector_t &E, Vector_t &B, const DiffDirection &dir) const
Definition: FM2DMagnetoStatic.cpp:204
FM2DMagnetoStatic::FM2DMagnetoStatic
FM2DMagnetoStatic(std::string aFilename)
Definition: FM2DMagnetoStatic.cpp:11
FM2DMagnetoStatic::freeMap
virtual void freeMap()
Definition: FM2DMagnetoStatic.cpp:156
Inform
Definition: Inform.h:42
Generated on Thu Oct 20 2022 17:40:10 for OPAL (Object Oriented Parallel Accelerator Library) by doxygen 1.9.3