Vlasiator ebf0dd394 on dev (v5.4.0 + 1054 commits)
Loading...
Searching...
No Matches
test1.cpp
Go to the documentation of this file.
1/*
2Background magnetic field test of GUMICS.
3
4Copyright 1997, 1998, 1999, 2000, 2001, 2010,
52011, 2012 Finnish Meteorological Institute
6*/
7
8#include "boost/program_options.hpp"
9#include "cmath"
10#include "cstdlib"
11#include "limits"
12
13#include "../B0.hpp"
14
15#define R_E 6.3712e6
16
17using namespace std;
18
19void get_dipole_B(
20 const double x, const double y, const double z,
21 double& Bx, double& By, double& Bz
22) {
23 const double k0 = 8e15;
24 Bx = -3 * k0 * x * z / pow(x*x + y*y + z*z, 5.0/2.0);
25 By = -3 * k0 * y * z / pow(x*x + y*y + z*z, 5.0/2.0);
26 Bz = -3 * k0 * (z*z - (x*x + y*y + z*z) / 3.0) / pow(x*x + y*y + z*z, 5.0/2.0);
27}
28
29
30int main(int argc, char* argv[])
31{
32 boost::program_options::options_description options(
33 "Usage: test1 [options (options given on the command line "
34 "override options given everywhere else)], where options are:"
35 );
36 options.add_options()("help", "print this help message")("verbose", "print run time information");
37
38 TB0 background_B(&options);
39
40 /*
41 Option parsing
42 */
43 boost::program_options::variables_map option_variables;
44
45 // read options from command line
46 boost::program_options::store(boost::program_options::parse_command_line(argc, argv, options), option_variables);
47 boost::program_options::notify(option_variables);
48 // read options from environment variables
49 boost::program_options::store(boost::program_options::parse_environment(options, "GUMICS_"), option_variables);
50 boost::program_options::notify(option_variables);
51
52 // print a help message if asked
53 if (option_variables.count("help") > 0) {
54 return EXIT_SUCCESS;
55 }
56
57 bool verbose = false;
58 if (option_variables.count("verbose") > 0) {
59 verbose = true;
60 }
61
62 background_B.initialize();
63
64 if (verbose) {
65 cout << "Testing background magnetic field calculation..." << endl;
66 }
67
68 // calculate volume average for smaller and smaller cube,
69 // result should converge to the value at the center
70 double
71 start[3] = {0, 0, 0},
72 end[3] = {0, 0, 0},
73 center[3] = {0, 0, 0};
74
75 // center in positive x direction
76 center[0] = 5 * R_E;
77 center[1] = 0 * R_E;
78 center[2] = 0 * R_E;
79
80 if (verbose) {
81 cout << "Cube center at "
82 << center[0] / R_E << ", " << center[1] / R_E << ", " << center[2] / R_E << " (R_E):"
83 << endl;
84 }
85
86 double
87 previous_difference = std::numeric_limits<double>::max(),
88 difference = std::numeric_limits<double>::max();
89
90 for (double size = 1 * R_E; size >= 0.01 * R_E; size /= 2) {
91 start[0] = center[0] - size / 2;
92 start[1] = center[1] - size / 2;
93 start[2] = center[2] - size / 2;
94 end[0] = center[0] + size / 2;
95 end[1] = center[1] + size / 2;
96 end[2] = center[2] + size / 2;
97
98 double
99 result_x = 0, result_y = 0, result_z = 0,
100 reference_x = 0, reference_y = 0, reference_z = 0;
101
102 get_dipole_B(center[0], center[1], center[2], reference_x, reference_y, reference_z);
103
104 background_B.BackgroundVolumeAverageFast(start, end, result_x, result_y, result_z);
105
106 const double
107 dx = result_x - reference_x,
108 dy = result_y - reference_y,
109 dz = result_z - reference_z,
110
111 difference = sqrt(dx*dx + dy*dy + dz*dz);
112
113 if (verbose) {
114 cout << "Difference for size " << size / R_E << " (R_E): " << difference << endl;
115 }
116
117 if (difference >= previous_difference) {
118 cerr << __FILE__ << ":" << __LINE__
119 << " New difference (" << difference
120 << ") is larger than previous (" << previous_difference
121 << ")"
122 << endl;
123 abort();
124 }
125
126 previous_difference = difference;
127 }
128
129
130 // center in -y direction
131 center[0] = 0 * R_E;
132 center[1] = -3 * R_E;
133 center[2] = 0 * R_E;
134
135 if (verbose) {
136 cout << "Cube center at "
137 << center[0] / R_E << ", " << center[1] / R_E << ", " << center[2] / R_E << " (R_E):"
138 << endl;
139 }
140
141 previous_difference = std::numeric_limits<double>::max(),
142 difference = std::numeric_limits<double>::max();
143
144 for (double size = 1 * R_E; size >= 0.01 * R_E; size /= 2) {
145 start[0] = center[0] - size / 2;
146 start[1] = center[1] - size / 2;
147 start[2] = center[2] - size / 2;
148 end[0] = center[0] + size / 2;
149 end[1] = center[1] + size / 2;
150 end[2] = center[2] + size / 2;
151
152 double
153 result_x = 0, result_y = 0, result_z = 0,
154 reference_x = 0, reference_y = 0, reference_z = 0;
155
156 get_dipole_B(center[0], center[1], center[2], reference_x, reference_y, reference_z);
157
158 background_B.BackgroundVolumeAverageFast(start, end, result_x, result_y, result_z);
159
160 const double
161 dx = result_x - reference_x,
162 dy = result_y - reference_y,
163 dz = result_z - reference_z,
164
165 difference = sqrt(dx*dx + dy*dy + dz*dz);
166
167 if (verbose) {
168 cout << "Difference for size " << size / R_E << " (R_E): " << difference << endl;
169 }
170
171 if (difference >= previous_difference) {
172 cerr << __FILE__ << ":" << __LINE__
173 << " New difference (" << difference
174 << ") is larger than previous (" << previous_difference
175 << ")"
176 << endl;
177 abort();
178 }
179
180 previous_difference = difference;
181 }
182
183
184 // +z
185 center[0] = 0 * R_E;
186 center[1] = 0 * R_E;
187 center[2] = 2 * R_E;
188
189 if (verbose) {
190 cout << "Cube center at "
191 << center[0] / R_E << ", " << center[1] / R_E << ", " << center[2] / R_E << " (R_E):"
192 << endl;
193 }
194
195 previous_difference = std::numeric_limits<double>::max(),
196 difference = std::numeric_limits<double>::max();
197
198 for (double size = 1 * R_E; size >= 0.01 * R_E; size /= 2) {
199 start[0] = center[0] - size / 2;
200 start[1] = center[1] - size / 2;
201 start[2] = center[2] - size / 2;
202 end[0] = center[0] + size / 2;
203 end[1] = center[1] + size / 2;
204 end[2] = center[2] + size / 2;
205
206 double
207 result_x = 0, result_y = 0, result_z = 0,
208 reference_x = 0, reference_y = 0, reference_z = 0;
209
210 get_dipole_B(center[0], center[1], center[2], reference_x, reference_y, reference_z);
211
212 background_B.BackgroundVolumeAverageFast(start, end, result_x, result_y, result_z);
213
214 const double
215 dx = result_x - reference_x,
216 dy = result_y - reference_y,
217 dz = result_z - reference_z,
218
219 difference = sqrt(dx*dx + dy*dy + dz*dz);
220
221 if (verbose) {
222 cout << "Difference for size " << size / R_E << " (R_E): " << difference << endl;
223 }
224
225 if (difference >= previous_difference) {
226 cerr << __FILE__ << ":" << __LINE__
227 << " New difference (" << difference
228 << ") is larger than previous (" << previous_difference
229 << ")"
230 << endl;
231 abort();
232 }
233
234 previous_difference = difference;
235 }
236
237
238 // +x, -y, -z
239 center[0] = 1 * R_E;
240 center[1] = -2 * R_E;
241 center[2] = 3 * R_E;
242
243 if (verbose) {
244 cout << "Cube center at "
245 << center[0] / R_E << ", " << center[1] / R_E << ", " << center[2] / R_E << " (R_E):"
246 << endl;
247 }
248
249 previous_difference = std::numeric_limits<double>::max(),
250 difference = std::numeric_limits<double>::max();
251
252 for (double size = 1 * R_E; size >= 0.01 * R_E; size /= 2) {
253 start[0] = center[0] - size / 2;
254 start[1] = center[1] - size / 2;
255 start[2] = center[2] - size / 2;
256 end[0] = center[0] + size / 2;
257 end[1] = center[1] + size / 2;
258 end[2] = center[2] + size / 2;
259
260 double
261 result_x = 0, result_y = 0, result_z = 0,
262 reference_x = 0, reference_y = 0, reference_z = 0;
263
264 get_dipole_B(center[0], center[1], center[2], reference_x, reference_y, reference_z);
265
266 background_B.BackgroundVolumeAverageFast(start, end, result_x, result_y, result_z);
267
268 const double
269 dx = result_x - reference_x,
270 dy = result_y - reference_y,
271 dz = result_z - reference_z,
272
273 difference = sqrt(dx*dx + dy*dy + dz*dz);
274
275 if (verbose) {
276 cout << "Difference for size " << size / R_E << " (R_E): " << difference << endl;
277 }
278
279 if (difference >= previous_difference) {
280 cerr << __FILE__ << ":" << __LINE__
281 << " New difference (" << difference
282 << ") is larger than previous (" << previous_difference
283 << ")"
284 << endl;
285 abort();
286 }
287
288 previous_difference = difference;
289 }
290
291
292 // cube overlaps the origin
293 center[0] = 0.26 * R_E;
294 center[1] = 0.26 * R_E;
295 center[2] = 0.26 * R_E;
296
297 if (verbose) {
298 cout << "Cube center at "
299 << center[0] / R_E << ", " << center[1] / R_E << ", " << center[2] / R_E << " (R_E):"
300 << endl;
301 }
302
303 previous_difference = std::numeric_limits<double>::max(),
304 difference = std::numeric_limits<double>::max();
305
306 for (double size = 1 * R_E; size >= 0.01 * R_E; size /= 2) {
307 start[0] = center[0] - size / 2;
308 start[1] = center[1] - size / 2;
309 start[2] = center[2] - size / 2;
310 end[0] = center[0] + size / 2;
311 end[1] = center[1] + size / 2;
312 end[2] = center[2] + size / 2;
313
314 double
315 result_x = 0, result_y = 0, result_z = 0,
316 reference_x = 0, reference_y = 0, reference_z = 0;
317
318 get_dipole_B(center[0], center[1], center[2], reference_x, reference_y, reference_z);
319
320 background_B.BackgroundVolumeAverageFast(start, end, result_x, result_y, result_z);
321
322 const double
323 dx = result_x - reference_x,
324 dy = result_y - reference_y,
325 dz = result_z - reference_z,
326
327 difference = sqrt(dx*dx + dy*dy + dz*dz);
328
329 if (verbose) {
330 cout << "Difference for size " << size / R_E << " (R_E): " << difference << endl;
331 }
332
333 if (difference >= previous_difference) {
334 cerr << __FILE__ << ":" << __LINE__
335 << " New difference (" << difference
336 << ") is larger than previous (" << previous_difference
337 << ")"
338 << endl;
339 abort();
340 }
341
342 previous_difference = difference;
343 }
344
345 cout << "PASSED" << endl;
346
347 return EXIT_SUCCESS;
348}
349
dx
dx
Definition Dispersion.m:38
sqrt
sqrt(1.0+vA *vA/(c *c))) % Ion-acoustic wave cS
R_E
#define R_E
Definition test1.cpp:15
get_dipole_B
void get_dipole_B(const double x, const double y, const double z, double &Bx, double &By, double &Bz)
Definition test1.cpp:19
main
int main()
Definition unit_testing.cpp:368