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