ASCOT5
Loading...
Searching...
No Matches
step_gc_cashkarp.c
Go to the documentation of this file.
1
5#include <stdlib.h>
6#include <stdio.h>
7#include <math.h>
8#include <float.h>
9#include "../../ascot5.h"
10#include "../../B_field.h"
11#include "../../math.h"
12#include "../../consts.h"
13#include "../../particle.h"
14#include "../../error.h"
15#include "step_gc_cashkarp.h"
16#include "step_gceom.h"
17#include "step_gceom_mhd.h"
18
35void step_gc_cashkarp(particle_simd_gc* p, real* h, real* hnext, real tol,
36 B_field_data* Bdata, E_field_data* Edata) {
37
38 int i;
39 /* Following loop will be executed simultaneously for all i */
40 #pragma omp simd aligned(h, hnext : 64)
41 for(i = 0; i < NSIMD; i++) {
42 if(p->running[i]) {
43 a5err errflag = 0;
44
45 real k1[6], k2[6], k3[6], k4[6], k5[6], k6[6];
46 real tempy[6];
47 real yprev[6];
48
49 real mass = p->mass[i];;
50 real charge = p->charge[i];
51
52 real B_dB[15];
53 real E[3];
54
55 real R0 = p->r[i];
56 real z0 = p->z[i];
57 real t0 = p->time[i];
58
59 /* Coordinates are copied from the struct into an array to make
60 * passing parameters easier */
61 yprev[0] = p->r[i];
62 yprev[1] = p->phi[i];
63 yprev[2] = p->z[i];
64 yprev[3] = p->ppar[i];
65 yprev[4] = p->mu[i];
66 yprev[5] = p->zeta[i];
67
68 /* Magnetic field at initial position already known */
69 B_dB[0] = p->B_r[i];
70 B_dB[1] = p->B_r_dr[i];
71 B_dB[2] = p->B_r_dphi[i];
72 B_dB[3] = p->B_r_dz[i];
73
74 B_dB[4] = p->B_phi[i];
75 B_dB[5] = p->B_phi_dr[i];
76 B_dB[6] = p->B_phi_dphi[i];
77 B_dB[7] = p->B_phi_dz[i];
78
79 B_dB[8] = p->B_z[i];
80 B_dB[9] = p->B_z_dr[i];
81 B_dB[10] = p->B_z_dphi[i];
82 B_dB[11] = p->B_z_dz[i];
83
84 if(!errflag) {
85 errflag = E_field_eval_E(E, yprev[0], yprev[1], yprev[2],
86 t0, Edata, Bdata);
87 }
88 if(!errflag) {
89 step_gceom(k1, yprev, mass, charge, B_dB, E);
90 }
91 for(int j = 0; j < 6; j++) {
92 tempy[j] = yprev[j]
93 + h[i]*(
94 (1.0/5) * k1[j] );
95 }
96
97
98 if(!errflag) {
99 errflag = B_field_eval_B_dB(B_dB, tempy[0], tempy[1], tempy[2],
100 t0 + (1.0/5)*h[i], Bdata);
101 }
102 if(!errflag) {
103 errflag = E_field_eval_E(E, tempy[0], tempy[1], tempy[2],
104 t0 + (1.0/5)*h[i], Edata, Bdata);
105 }
106 if(!errflag) {
107 step_gceom(k2, tempy, mass, charge, B_dB, E);
108 }
109 for(int j = 0; j < 6; j++) {
110 tempy[j] = yprev[j]
111 + h[i]*(
112 (3.0/40) * k1[j]
113 + (9.0/40) * k2[j] );
114 }
115
116
117 if(!errflag) {
118 errflag = B_field_eval_B_dB(B_dB, tempy[0], tempy[1], tempy[2],
119 t0 + (3.0/10)*h[i], Bdata);
120 }
121 if(!errflag) {
122 errflag = E_field_eval_E(E, tempy[0], tempy[1], tempy[2],
123 t0 + (3.0/10)*h[i], Edata, Bdata);
124 }
125 if(!errflag) {
126 step_gceom(k3, tempy, mass, charge, B_dB, E);
127 }
128 for(int j = 0; j < 6; j++) {
129 tempy[j] = yprev[j]
130 + h[i]*(
131 ( 3.0/10) * k1[j]
132 + (-9.0/10) * k2[j]
133 + ( 6.0/5 ) * k3[j] );
134 }
135
136
137 if(!errflag) {
138 errflag = B_field_eval_B_dB(B_dB, tempy[0], tempy[1], tempy[2],
139 t0 + (3.0/5)*h[i], Bdata);
140 }
141 if(!errflag) {
142 errflag = E_field_eval_E(E, tempy[0], tempy[1], tempy[2],
143 t0 + (3.0/5)*h[i], Edata, Bdata);
144 }
145 if(!errflag) {
146 step_gceom(k4, tempy, mass, charge, B_dB, E);
147 }
148 for(int j = 0; j < 6; j++) {
149 tempy[j] = yprev[j]
150 + h[i]*(
151 (-11.0/54) * k1[j]
152 + ( 5.0/2 ) * k2[j]
153 + (-70.0/27) * k3[j]
154 + ( 35.0/27) * k4[j] );
155 }
156
157
158 if(!errflag) {
159 errflag = B_field_eval_B_dB(B_dB, tempy[0], tempy[1], tempy[2],
160 t0 + h[i], Bdata);
161 }
162 if(!errflag) {
163 errflag = E_field_eval_E(E, tempy[0], tempy[1], tempy[2],
164 t0 + h[i], Edata, Bdata);
165 }
166 if(!errflag) {
167 step_gceom(k5, tempy, mass, charge, B_dB, E);
168 }
169 for(int j = 0; j < 6; j++) {
170 tempy[j] = yprev[j]
171 + h[i]*(
172 ( 1631.0/55296 ) * k1[j]
173 + ( 175.0/512 ) * k2[j]
174 + ( 575.0/13824 ) * k3[j]
175 + (44275.0/110592) * k4[j]
176 + ( 253.0/4096 ) * k5[j] );
177 }
178
179
180 if(!errflag) {
181 errflag = B_field_eval_B_dB(B_dB, tempy[0], tempy[1], tempy[2],
182 t0 + (7.0/8)*h[i], Bdata);
183 }
184 if(!errflag) {
185 errflag = E_field_eval_E(E, tempy[0], tempy[1], tempy[2],
186 t0 + (7.0/8)*h[i], Edata, Bdata);
187 }
188 if(!errflag) {
189 step_gceom(k6, tempy, mass, charge, B_dB, E);
190 }
191
192 /* Error estimate is a difference between RK4 and RK5 solutions. If
193 * time-step is accepted, the RK5 solution will be used to advance
194 * marker. */
195 real rk5[6], rk4[6];
196 if(!errflag) {
197 real err = 0.0;
198 for(int j = 0; j < 6; j++) {
199 rk5[j] = yprev[j]
200 + h[i]*(
201 ( 37.0/378 ) * k1[j]
202 + (250.0/621 ) * k3[j]
203 + (125.0/594 ) * k4[j]
204 + (512.0/1771) * k6[j] );
205
206 rk4[j] = yprev[j] +
207 h[i]*(
208 ( 2825.0/27648) * k1[j]
209 + (18575.0/48384) * k3[j]
210 + (13525.0/55296) * k4[j]
211 + ( 277.0/14336) * k5[j]
212 + ( 1.0/4 ) * k6[j] );
213 if(j==3) {
214 real yerr = fabs(rk5[j] - rk4[j]);
215 real ytol = fabs(yprev[j]) + fabs(k1[j]*h[i])
216 + DBL_EPSILON;
217 err = fmax( err, yerr/ytol );
218 }
219 else if(j==2) {
220 real rk1[3] = {k1[0]*h[i], k1[1]*h[i], k1[2]*h[i]};
221 real yerr =
222 rk5[0] * rk5[0] + rk4[0] * rk4[0]
223 - 2 * rk5[0] * rk4[0] * cos(rk5[1] - rk4[1])
224 + ( rk5[2] - rk4[2] ) * ( rk5[2] - rk4[2] );
225 real ytol =
226 yprev[0] * yprev[0] + rk1[0] * rk1[0]
227 - 2 * yprev[0] * rk1[0] * cos(yprev[1] - rk1[1])
228 + ( yprev[2] - rk1[2] ) * ( yprev[2] - rk1[2] )
229 + DBL_EPSILON;
230 err = fmax( err, sqrt(yerr/ytol) );
231 }
232 }
233
234 err = err/tol;
235 if(err <= 1){
236 /* Time step accepted */
237 hnext[i] = 0.85*h[i]*pow(err,-0.2);
238
239 /* Make sure we don't make a huge jump */
240 if(hnext[i] > 1.5*h[i]) {
241 hnext[i] = 1.5*h[i];
242 }
243 }
244 else{
245 /* Time step rejected */
246 hnext[i] = -0.85*h[i]*pow(err,-0.25);
247 }
248 }
249
250 /* Test that results are physical */
251 if(!errflag && fabs(hnext[i]) < A5_EXTREMELY_SMALL_TIMESTEP) {
252 errflag = error_raise(
253 ERR_INVALID_TIMESTEP, __LINE__, EF_STEP_GC_CASHKARP);
254 }
255 else if(!errflag && rk5[0] <= 0) {
256 errflag = error_raise(
257 ERR_INTEGRATION, __LINE__, EF_STEP_GC_CASHKARP);
258 }
259 else if(!errflag && rk5[4] < 0) {
260 errflag = error_raise(
261 ERR_INTEGRATION, __LINE__, EF_STEP_GC_CASHKARP);
262 }
263
264 /* Update gc phase space position */
265 if(!errflag) {
266 p->r[i] = rk5[0];
267 p->phi[i] = rk5[1];
268 p->z[i] = rk5[2];
269 p->ppar[i] = rk5[3];
270 p->mu[i] = rk5[4];
271 p->zeta[i] = fmod( rk5[5], CONST_2PI );
272 if(p->zeta[i]<0) {
273 p->zeta[i] = CONST_2PI + p->zeta[i];
274 }
275 }
276
277 /* Evaluate magnetic field (and gradient) and rho at new position */
278 real psi[1];
279 real rho[2];
280 if(!errflag) {
281 errflag = B_field_eval_B_dB(B_dB, p->r[i], p->phi[i], p->z[i],
282 p->time[i] + h[i], Bdata);
283 }
284 if(!errflag) {
285 errflag = B_field_eval_psi(psi, p->r[i], p->phi[i], p->z[i],
286 p->time[i] + h[i], Bdata);
287 }
288 if(!errflag) {
289 errflag = B_field_eval_rho(rho, psi[0], Bdata);
290 }
291
292 if(!errflag) {
293 p->B_r[i] = B_dB[0];
294 p->B_r_dr[i] = B_dB[1];
295 p->B_r_dphi[i] = B_dB[2];
296 p->B_r_dz[i] = B_dB[3];
297
298 p->B_phi[i] = B_dB[4];
299 p->B_phi_dr[i] = B_dB[5];
300 p->B_phi_dphi[i] = B_dB[6];
301 p->B_phi_dz[i] = B_dB[7];
302
303 p->B_z[i] = B_dB[8];
304 p->B_z_dr[i] = B_dB[9];
305 p->B_z_dphi[i] = B_dB[10];
306 p->B_z_dz[i] = B_dB[11];
307 p->rho[i] = rho[0];
308
309 /* Evaluate theta angle so that it is cumulative */
310 real axisrz[2];
311 errflag = B_field_get_axis_rz(axisrz, Bdata, p->phi[i]);
312 p->theta[i] += atan2( (R0-axisrz[0]) * (p->z[i]-axisrz[1])
313 - (z0-axisrz[1]) * (p->r[i]-axisrz[0]),
314 (R0-axisrz[0]) * (p->r[i]-axisrz[0])
315 + (z0-axisrz[1]) * (p->z[i]-axisrz[1]) );
316 }
317
318 /* Error handling */
319 if(errflag) {
320 p->err[i] = errflag;
321 p->running[i] = 0;
322 hnext[i] = h[i];
323 }
324 }
325 }
326}
327
344void step_gc_cashkarp_mhd(particle_simd_gc* p, real* h, real* hnext, real tol,
345 B_field_data* Bdata, E_field_data* Edata,
346 boozer_data* boozer, mhd_data* mhd) {
347
348 int i;
349 /* Following loop will be executed simultaneously for all i */
350#pragma omp simd aligned(h, hnext : 64)
351 for(i = 0; i < NSIMD; i++) {
352 if(p->running[i]) {
353 a5err errflag = 0;
354
355 real k1[6], k2[6], k3[6], k4[6], k5[6], k6[6];
356 real tempy[6];
357 real yprev[6];
358
359 real mass = p->mass[i];;
360 real charge = p->charge[i];
361
362 real B_dB[15];
363 real E[3];
364 real mhd_dmhd[10];
365
366 real R0 = p->r[i];
367 real z0 = p->z[i];
368 real t0 = p->time[i];
369
370 /* Coordinates are copied from the struct into an array to make
371 * passing parameters easier */
372 yprev[0] = p->r[i];
373 yprev[1] = p->phi[i];
374 yprev[2] = p->z[i];
375 yprev[3] = p->ppar[i];
376 yprev[4] = p->mu[i];
377 yprev[5] = p->zeta[i];
378
379 /* Magnetic field at initial position already known */
380 B_dB[0] = p->B_r[i];
381 B_dB[1] = p->B_r_dr[i];
382 B_dB[2] = p->B_r_dphi[i];
383 B_dB[3] = p->B_r_dz[i];
384
385 B_dB[4] = p->B_phi[i];
386 B_dB[5] = p->B_phi_dr[i];
387 B_dB[6] = p->B_phi_dphi[i];
388 B_dB[7] = p->B_phi_dz[i];
389
390 B_dB[8] = p->B_z[i];
391 B_dB[9] = p->B_z_dr[i];
392 B_dB[10] = p->B_z_dphi[i];
393 B_dB[11] = p->B_z_dz[i];
394
395 if(!errflag) {
396 errflag = E_field_eval_E(E, yprev[0], yprev[1], yprev[2],
397 t0, Edata, Bdata);
398 }
399 if(!errflag) {
400 errflag = mhd_eval(mhd_dmhd, yprev[0], yprev[1], yprev[2],
401 t0, MHD_INCLUDE_ALL, boozer, mhd, Bdata);
402 }
403 if(!errflag) {
404 step_gceom_mhd(k1, yprev, mass, charge, B_dB, E, mhd_dmhd);
405 }
406 for(int j = 0; j < 6; j++) {
407 tempy[j] = yprev[j]
408 + h[i]*(
409 (1.0/5) * k1[j] );
410 }
411
412
413 if(!errflag) {
414 errflag = B_field_eval_B_dB(B_dB, tempy[0], tempy[1], tempy[2],
415 t0 + (1.0/5)*h[i], Bdata);
416 }
417 if(!errflag) {
418 errflag = E_field_eval_E(E, tempy[0], tempy[1], tempy[2],
419 t0 + (1.0/5)*h[i], Edata, Bdata);
420 }
421 if(!errflag) {
422 errflag = mhd_eval(mhd_dmhd, tempy[0], tempy[1], tempy[2],
423 t0 + (1.0/5)*h[i], MHD_INCLUDE_ALL, boozer,
424 mhd, Bdata);
425 }
426 if(!errflag) {
427 step_gceom_mhd(k2, tempy, mass, charge, B_dB, E, mhd_dmhd);
428 }
429 for(int j = 0; j < 6; j++) {
430 tempy[j] = yprev[j]
431 + h[i]*(
432 (3.0/40) * k1[j]
433 + (9.0/40) * k2[j] );
434 }
435
436
437 if(!errflag) {
438 errflag = B_field_eval_B_dB(B_dB, tempy[0], tempy[1], tempy[2],
439 t0 + (3.0/10)*h[i], Bdata);
440 }
441 if(!errflag) {
442 errflag = E_field_eval_E(E, tempy[0], tempy[1], tempy[2],
443 t0 + (3.0/10)*h[i], Edata, Bdata);
444 }
445 if(!errflag) {
446 errflag = mhd_eval(mhd_dmhd, tempy[0], tempy[1], tempy[2],
447 t0 + (3.0/10)*h[i], MHD_INCLUDE_ALL, boozer,
448 mhd, Bdata);
449 }
450 if(!errflag) {
451 step_gceom_mhd(k3, tempy, mass, charge, B_dB, E, mhd_dmhd);
452 }
453 for(int j = 0; j < 6; j++) {
454 tempy[j] = yprev[j]
455 + h[i]*(
456 ( 3.0/10) * k1[j]
457 + (-9.0/10) * k2[j]
458 + ( 6.0/5 ) * k3[j] );
459 }
460
461
462 if(!errflag) {
463 errflag = B_field_eval_B_dB(B_dB, tempy[0], tempy[1], tempy[2],
464 t0 + (3.0/5)*h[i], Bdata);
465 }
466 if(!errflag) {
467 errflag = E_field_eval_E(E, tempy[0], tempy[1], tempy[2],
468 t0 + (3.0/5)*h[i], Edata, Bdata);
469 }
470 if(!errflag) {
471 errflag = mhd_eval(mhd_dmhd, tempy[0], tempy[1], tempy[2],
472 t0 + (3.0/5)*h[i], MHD_INCLUDE_ALL, boozer,
473 mhd, Bdata);
474 }
475 if(!errflag) {
476 step_gceom_mhd(k4, tempy, mass, charge, B_dB, E, mhd_dmhd);
477 }
478 for(int j = 0; j < 6; j++) {
479 tempy[j] = yprev[j]
480 + h[i]*(
481 (-11.0/54) * k1[j]
482 + ( 5.0/2 ) * k2[j]
483 + (-70.0/27) * k3[j]
484 + ( 35.0/27) * k4[j] );
485 }
486
487
488 if(!errflag) {
489 errflag = B_field_eval_B_dB(B_dB, tempy[0], tempy[1], tempy[2],
490 t0 + h[i], Bdata);
491 }
492 if(!errflag) {
493 errflag = E_field_eval_E(E, tempy[0], tempy[1], tempy[2],
494 t0 + h[i], Edata, Bdata);
495 }
496 if(!errflag) {
497 errflag = mhd_eval(mhd_dmhd, tempy[0], tempy[1], tempy[2],
498 t0 + h[i], MHD_INCLUDE_ALL, boozer, mhd,
499 Bdata);
500 }
501 if(!errflag) {
502 step_gceom_mhd(k5, tempy, mass, charge, B_dB, E, mhd_dmhd);
503 }
504 for(int j = 0; j < 6; j++) {
505 tempy[j] = yprev[j]
506 + h[i]*(
507 ( 1631.0/55296 ) * k1[j]
508 + ( 175.0/512 ) * k2[j]
509 + ( 575.0/13824 ) * k3[j]
510 + (44275.0/110592) * k4[j]
511 + ( 253.0/4096 ) * k5[j] );
512 }
513
514
515 if(!errflag) {
516 errflag = B_field_eval_B_dB(B_dB, tempy[0], tempy[1], tempy[2],
517 t0 + (7.0/8)*h[i], Bdata);
518 }
519 if(!errflag) {
520 errflag = E_field_eval_E(E, tempy[0], tempy[1], tempy[2],
521 t0 + (7.0/8)*h[i], Edata, Bdata);
522 }
523 if(!errflag) {
524 errflag = mhd_eval(mhd_dmhd, tempy[0], tempy[1], tempy[2],
525 t0 + (7.0/8)*h[i], MHD_INCLUDE_ALL, boozer,
526 mhd, Bdata);
527 }
528 if(!errflag) {
529 step_gceom_mhd(k6, tempy, mass, charge, B_dB, E, mhd_dmhd);
530 }
531
532 /* Error estimate is a difference between RK4 and RK5 solutions. If
533 * time-step is accepted, the RK5 solution will be used to advance
534 * marker. */
535 real rk5[6], rk4[6];
536 if(!errflag) {
537 real err = 0.0;
538 for(int j = 0; j < 5; j++) {
539 rk5[j] = yprev[j]
540 + h[i]*(
541 ( 37.0/378 ) * k1[j]
542 + (250.0/621 ) * k3[j]
543 + (125.0/594 ) * k4[j]
544 + (512.0/1771) * k6[j] );
545
546 rk4[j] = yprev[j] +
547 h[i]*(
548 ( 2825.0/27648) * k1[j]
549 + (18575.0/48384) * k3[j]
550 + (13525.0/55296) * k4[j]
551 + ( 277.0/14336) * k5[j]
552 + ( 1.0/4 ) * k6[j] );
553 if(j==3) {
554 real yerr = fabs(rk5[j] - rk4[j]);
555 real ytol = fabs(yprev[j]) + fabs(k1[j]*h[i])
556 + DBL_EPSILON;
557 err = fmax( err, yerr/ytol );
558 }
559 else if(j==2) {
560 real rk1[3] = {k1[0]*h[i], k1[1]*h[i], k1[2]*h[i]};
561 real yerr =
562 rk5[0] * rk5[0] + rk4[0] * rk4[0]
563 - 2 * rk5[0] * rk4[0] * cos(rk5[1] - rk4[1])
564 + ( rk5[2] - rk4[2] ) * ( rk5[2] - rk4[2] );
565 real ytol =
566 yprev[0] * yprev[0] + rk1[0] * rk1[0]
567 - 2 * yprev[0] * rk1[0] * cos(yprev[1] - rk1[1])
568 + ( yprev[2] - rk1[2] ) * ( yprev[2] - rk1[2] )
569 + DBL_EPSILON;
570 err = fmax( err, sqrt(yerr/ytol) );
571 }
572 }
573
574 err = err/tol;
575 if(err <= 1){
576 /* Time step accepted */
577 hnext[i] = 0.85*h[i]*pow(err,-0.2);
578
579 /* Make sure we don't make a huge jump */
580 if(hnext[i] > 1.5*h[i]) {
581 hnext[i] = 1.5*h[i];
582 }
583 }
584 else{
585 /* Time step rejected */
586 hnext[i] = -0.85*h[i]*pow(err,-0.25);
587 }
588 }
589
590 /* Update gc phase space position */
591 if(!errflag) {
592 p->r[i] = rk5[0];
593 p->phi[i] = rk5[1];
594 p->z[i] = rk5[2];
595 p->ppar[i] = rk5[3];
596 p->mu[i] = rk5[4];
597 p->zeta[i] = fmod( rk5[5], CONST_2PI );
598 if(p->zeta[i]<0) {
599 p->zeta[i] = CONST_2PI + p->zeta[i];
600 }
601 }
602
603 /* Evaluate magnetic field (and gradient) and rho at new position */
604 real psi[1];
605 real rho[2];
606 if(!errflag) {
607 errflag = B_field_eval_B_dB(B_dB, p->r[i], p->phi[i], p->z[i],
608 p->time[i] + h[i], Bdata);
609 }
610 if(!errflag) {
611 errflag = B_field_eval_psi(psi, p->r[i], p->phi[i], p->z[i],
612 p->time[i] + h[i], Bdata);
613 }
614 if(!errflag) {
615 errflag = B_field_eval_rho(rho, psi[0], Bdata);
616 }
617
618 if(!errflag) {
619 p->B_r[i] = B_dB[0];
620 p->B_r_dr[i] = B_dB[1];
621 p->B_r_dphi[i] = B_dB[2];
622 p->B_r_dz[i] = B_dB[3];
623
624 p->B_phi[i] = B_dB[4];
625 p->B_phi_dr[i] = B_dB[5];
626 p->B_phi_dphi[i] = B_dB[6];
627 p->B_phi_dz[i] = B_dB[7];
628
629 p->B_z[i] = B_dB[8];
630 p->B_z_dr[i] = B_dB[9];
631 p->B_z_dphi[i] = B_dB[10];
632 p->B_z_dz[i] = B_dB[11];
633 p->rho[i] = rho[0];
634
635 /* Evaluate theta angle so that it is cumulative */
636 real axisrz[2];
637 errflag = B_field_get_axis_rz(axisrz, Bdata, p->phi[i]);
638 p->theta[i] += atan2( (R0-axisrz[0]) * (p->z[i]-axisrz[1])
639 - (z0-axisrz[1]) * (p->r[i]-axisrz[0]),
640 (R0-axisrz[0]) * (p->r[i]-axisrz[0])
641 + (z0-axisrz[1]) * (p->z[i]-axisrz[1]) );
642 }
643
644 /* Error handling */
645 if(errflag) {
646 p->err[i] = errflag;
647 p->running[i] = 0;
648 hnext[i] = h[i];
649 }
650 }
651 }
652}
B_field_eval_rho
a5err B_field_eval_rho(real rho[2], real psi, B_field_data *Bdata)
Evaluate normalized poloidal flux rho and its psi derivative.
Definition B_field.c:228
B_field_eval_psi
a5err B_field_eval_psi(real *psi, real r, real phi, real z, real t, B_field_data *Bdata)
Evaluate poloidal flux psi.
Definition B_field.c:102
B_field_eval_B_dB
a5err B_field_eval_B_dB(real B_dB[15], real r, real phi, real z, real t, B_field_data *Bdata)
Evaluate magnetic field and its derivatives.
Definition B_field.c:449
B_field_get_axis_rz
a5err B_field_get_axis_rz(real rz[2], B_field_data *Bdata, real phi)
Return magnetic axis Rz-coordinates.
Definition B_field.c:501
B_field.h
Header file for B_field.c.
E_field_eval_E
a5err E_field_eval_E(real E[3], real r, real phi, real z, real t, E_field_data *Edata, B_field_data *Bdata)
Evaluate electric field.
Definition E_field.c:82
ascot5.h
Main header file for ASCOT5.
real
double real
Definition ascot5.h:85
NSIMD
#define NSIMD
Number of particles simulated simultaneously in a particle group operations.
Definition ascot5.h:91
A5_EXTREMELY_SMALL_TIMESTEP
#define A5_EXTREMELY_SMALL_TIMESTEP
If adaptive time step falls below this value, produce an error.
Definition ascot5.h:118
consts.h
Header file containing physical and mathematical constants.
CONST_2PI
#define CONST_2PI
2*pi
Definition consts.h:14
error.h
Error module for ASCOT5.
a5err
unsigned long int a5err
Simulation error flag.
Definition error.h:17
EF_STEP_GC_CASHKARP
@ EF_STEP_GC_CASHKARP
Definition error.h:31
ERR_INVALID_TIMESTEP
@ ERR_INVALID_TIMESTEP
Definition error.h:67
ERR_INTEGRATION
@ ERR_INTEGRATION
Definition error.h:69
error_raise
static DECLARE_TARGET_SIMD a5err error_raise(error_type type, int line, error_file file)
Raise a new error.
Definition error.h:86
fmod
real fmod(real x, real y)
Compute the modulus of two real numbers.
Definition math.c:22
math.h
Header file for math.c.
mhd_eval
a5err mhd_eval(real mhd_dmhd[10], real r, real phi, real z, real t, int includemode, boozer_data *boozerdata, mhd_data *mhddata, B_field_data *Bdata)
Evaluate the needed quantities from MHD mode for orbit following.
Definition mhd.c:91
MHD_INCLUDE_ALL
#define MHD_INCLUDE_ALL
includemode parameter to include all modes (default)
Definition mhd.h:19
particle.h
Header file for particle.c.
step_gc_cashkarp_mhd
void step_gc_cashkarp_mhd(particle_simd_gc *p, real *h, real *hnext, real tol, B_field_data *Bdata, E_field_data *Edata, boozer_data *boozer, mhd_data *mhd)
Integrate a guiding center step for a struct of markers with MHD.
Definition step_gc_cashkarp.c:344
step_gc_cashkarp
void step_gc_cashkarp(particle_simd_gc *p, real *h, real *hnext, real tol, B_field_data *Bdata, E_field_data *Edata)
Integrate a guiding center step for a struct of markers.
Definition step_gc_cashkarp.c:35
step_gceom.h
Guiding center equations of motion.
step_gceom
static DECLARE_TARGET_SIMD void step_gceom(real *ydot, real *y, real mass, real charge, real *B_dB, real *E)
Calculate guiding center equations of motion for a single particle.
Definition step_gceom.h:26
step_gceom_mhd.h
Guiding center equations of motion with MHD activity.
step_gceom_mhd
static DECLARE_TARGET_SIMD void step_gceom_mhd(real *ydot, real *y, real mass, real charge, real *B_dB, real *E, real *mhd_dmhd)
Calculate guiding center equations of motion for a single particle.
Definition step_gceom_mhd.h:28
B_field_data
Magnetic field simulation data.
Definition B_field.h:41
E_field_data
Electric field simulation data.
Definition E_field.h:36
boozer_data
Data for mapping between the cylindrical and Boozer coordinates.
Definition boozer.h:16
mhd_data
MHD simulation data.
Definition mhd.h:35
particle_simd_gc
Struct representing NSIMD guiding center markers.
Definition particle.h:275
Generated on Wed Mar 19 2025 11:10:23 for ASCOT5 by doxygen 1.13.2