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