ASCOT5
Loading...
Searching...
No Matches
step_gc_rk4.c
Go to the documentation of this file.
1
5#include <stdio.h>
6#include <math.h>
7#include "../../ascot5.h"
8#include "../../consts.h"
9#include "../../B_field.h"
10#include "../../E_field.h"
11#include "../../boozer.h"
12#include "../../mhd.h"
13#include "../../math.h"
14#include "../../particle.h"
15#include "../../error.h"
16#include "step_gceom.h"
17#include "step_gceom_mhd.h"
18#include "step_gc_rk4.h"
19
34void step_gc_rk4(particle_simd_gc* p, real* h, B_field_data* Bdata,
35 E_field_data* Edata, int aldforce) {
36
37 int i;
38 /* Following loop will be executed simultaneously for all i */
39 GPU_DATA_IS_MAPPED(h[0:p->n_mrk])
40 GPU_PARALLEL_LOOP_ALL_LEVELS
41 for(i = 0; i < p->n_mrk; i++) {
42 if(p->running[i]) {
43 a5err errflag = 0;
44
45 real k1[6], k2[6], k3[6], k4[6];
46 real tempy[6];
47 real yprev[6];
48 real y[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
93
94 /* particle coordinates for the subsequent ydot evaluations are
95 * stored in tempy */
96 for(int j = 0; j < 6; j++) {
97 tempy[j] = yprev[j] + h[i]*k1[j]/2.0;
98 }
99
100
101 if(!errflag) {
102 errflag = B_field_eval_B_dB(B_dB, tempy[0], tempy[1], tempy[2],
103 t0 + h[i]/2.0, Bdata);
104 }
105 if(!errflag) {
106 errflag = E_field_eval_E(E, tempy[0], tempy[1], tempy[2],
107 t0 + h[i]/2.0, Edata, Bdata);
108 }
109 if(!errflag) {
110 step_gceom(k2, tempy, mass, charge, B_dB, E, aldforce);
111 }
112 for(int j = 0; j < 6; j++) {
113 tempy[j] = yprev[j] + h[i]*k2[j]/2.0;
114 }
115
116
117 if(!errflag) {
118 errflag = B_field_eval_B_dB(B_dB, tempy[0], tempy[1], tempy[2],
119 t0 + h[i]/2.0, Bdata);
120 }
121 if(!errflag) {
122 errflag = E_field_eval_E(E, tempy[0], tempy[1], tempy[2],
123 t0 + h[i]/2.0, Edata, Bdata);
124 }
125 if(!errflag) {
126 step_gceom(k3, tempy, mass, charge, B_dB, E, aldforce);
127 }
128 for(int j = 0; j < 6; j++) {
129 tempy[j] = yprev[j] + h[i]*k3[j];
130 }
131
132
133 if(!errflag) {
134 errflag = B_field_eval_B_dB(B_dB, tempy[0], tempy[1], tempy[2],
135 t0 + h[i], Bdata);
136 }
137 if(!errflag) {
138 errflag = E_field_eval_E(E, tempy[0], tempy[1], tempy[2],
139 t0 + h[i], Edata, Bdata);}
140 if(!errflag) {
141 step_gceom(k4, tempy, mass, charge, B_dB, E, aldforce);
142 }
143 for(int j = 0; j < 6; j++) {
144 y[j] = yprev[j]
145 + h[i]/6.0 * (k1[j] + 2*k2[j] + 2*k3[j] + k4[j]);
146 }
147
148
149 /* Test that results are physical */
150 if(!errflag && y[0] <= 0) {
151 errflag = error_raise(ERR_INTEGRATION, __LINE__, EF_STEP_GC_RK4);
152 }
153 if(!errflag && y[4] < 0) {
154 errflag = error_raise(ERR_INTEGRATION, __LINE__, EF_STEP_GC_RK4);
155 }
156
157 /* Update gc phase space position */
158 if(!errflag) {
159 p->r[i] = y[0];
160 p->phi[i] = y[1];
161 p->z[i] = y[2];
162 p->ppar[i] = y[3];
163 p->mu[i] = y[4];
164 p->zeta[i] = fmod(y[5],CONST_2PI);
165 if(p->zeta[i]<0) {
166 p->zeta[i] = CONST_2PI + p->zeta[i];
167 }
168 }
169
170 /* Evaluate magnetic field (and gradient) and rho at new position */
171 real psi[1];
172 real rho[2];
173 if(!errflag) {
174 errflag = B_field_eval_B_dB(B_dB, p->r[i], p->phi[i], p->z[i],
175 t0 + h[i], Bdata);
176 }
177 if(!errflag) {
178 errflag = B_field_eval_psi(psi, p->r[i], p->phi[i], p->z[i],
179 t0 + h[i], Bdata);
180 }
181 if(!errflag) {
182 errflag = B_field_eval_rho(rho, psi[0], Bdata);
183 }
184
185 if(!errflag) {
186 p->B_r[i] = B_dB[0];
187 p->B_r_dr[i] = B_dB[1];
188 p->B_r_dphi[i] = B_dB[2];
189 p->B_r_dz[i] = B_dB[3];
190
191 p->B_phi[i] = B_dB[4];
192 p->B_phi_dr[i] = B_dB[5];
193 p->B_phi_dphi[i] = B_dB[6];
194 p->B_phi_dz[i] = B_dB[7];
195
196 p->B_z[i] = B_dB[8];
197 p->B_z_dr[i] = B_dB[9];
198 p->B_z_dphi[i] = B_dB[10];
199 p->B_z_dz[i] = B_dB[11];
200
201 p->rho[i] = rho[0];
202
203 /* Evaluate theta angle so that it is cumulative */
204 real axisrz[2];
205 errflag = B_field_get_axis_rz(axisrz, Bdata, p->phi[i]);
206 p->theta[i] += atan2( (R0-axisrz[0]) * (p->z[i]-axisrz[1])
207 - (z0-axisrz[1]) * (p->r[i]-axisrz[0]),
208 (R0-axisrz[0]) * (p->r[i]-axisrz[0])
209 + (z0-axisrz[1]) * (p->z[i]-axisrz[1]) );
210 }
211
212 /* Error handling */
213 if(errflag) {
214 p->err[i] = errflag;
215 p->running[i] = 0;
216 }
217 }
218 }
219}
220
234void step_gc_rk4_mhd(
235 particle_simd_gc* p, real* h, B_field_data* Bdata, E_field_data* Edata,
236 boozer_data* boozer, mhd_data* mhd, int aldforce) {
237
238 int i;
239 /* Following loop will be executed simultaneously for all i */
240 #pragma omp simd aligned(h : 64)
241 for(i = 0; i < NSIMD; i++) {
242 if(p->running[i]) {
243 a5err errflag = 0;
244
245 real k1[6], k2[6], k3[6], k4[6];
246 real tempy[6];
247 real yprev[6];
248 real y[6];
249
250 real mass = p->mass[i];
251 real charge = p->charge[i];
252
253 real B_dB[15];
254 real E[3];
255 real mhd_dmhd[10];
256
257 real R0 = p->r[i];
258 real z0 = p->z[i];
259 real t0 = p->time[i];
260
261 /* Coordinates are copied from the struct into an array to make
262 * passing parameters easier */
263 yprev[0] = p->r[i];
264 yprev[1] = p->phi[i];
265 yprev[2] = p->z[i];
266 yprev[3] = p->ppar[i];
267 yprev[4] = p->mu[i];
268 yprev[5] = p->zeta[i];
269
270 /* Magnetic field at initial position already known */
271 B_dB[0] = p->B_r[i];
272 B_dB[1] = p->B_r_dr[i];
273 B_dB[2] = p->B_r_dphi[i];
274 B_dB[3] = p->B_r_dz[i];
275
276 B_dB[4] = p->B_phi[i];
277 B_dB[5] = p->B_phi_dr[i];
278 B_dB[6] = p->B_phi_dphi[i];
279 B_dB[7] = p->B_phi_dz[i];
280
281 B_dB[8] = p->B_z[i];
282 B_dB[9] = p->B_z_dr[i];
283 B_dB[10] = p->B_z_dphi[i];
284 B_dB[11] = p->B_z_dz[i];
285
286 if(!errflag) {
287 errflag = E_field_eval_E(E, yprev[0], yprev[1], yprev[2], t0,
288 Edata, Bdata);
289 }
290 if(!errflag) {
291 errflag = mhd_eval(mhd_dmhd, yprev[0], yprev[1], yprev[2], t0,
292 MHD_INCLUDE_ALL, boozer, mhd, Bdata);
293 }
294 if(!errflag) {
295 step_gceom_mhd(k1, yprev, mass, charge, B_dB, E, mhd_dmhd,
296 aldforce);
297 }
298
299 /* particle coordinates for the subsequent ydot evaluations are
300 * stored in tempy */
301 for(int j = 0; j < 6; j++) {
302 tempy[j] = yprev[j] + h[i]/2.0*k1[j];
303 }
304
305 if(!errflag) {
306 errflag = B_field_eval_B_dB(B_dB, tempy[0], tempy[1], tempy[2],
307 t0 + h[i]/2.0, Bdata);
308 }
309 if(!errflag) {
310 errflag = E_field_eval_E(E, tempy[0], tempy[1], tempy[2],
311 t0 + h[i]/2.0, Edata, Bdata);
312 }
313 if(!errflag) {
314 errflag = mhd_eval(mhd_dmhd, tempy[0], tempy[1], tempy[2],
315 t0 + h[i]/2.0, MHD_INCLUDE_ALL, boozer,
316 mhd, Bdata);
317 }
318 if(!errflag) {
319 step_gceom_mhd(k2, tempy, mass, charge, B_dB, E, mhd_dmhd,
320 aldforce);
321 }
322 for(int j = 0; j < 6; j++) {
323 tempy[j] = yprev[j] + h[i]/2.0*k2[j];
324 }
325
326 if(!errflag) {
327 errflag = B_field_eval_B_dB(B_dB, tempy[0], tempy[1], tempy[2],
328 t0 + h[i]/2.0, Bdata);
329 }
330 if(!errflag) {
331 errflag = E_field_eval_E(E, tempy[0], tempy[1], tempy[2],
332 t0 + h[i]/2.0, Edata, Bdata);
333 }
334 if(!errflag) {
335 errflag = mhd_eval(mhd_dmhd, tempy[0], tempy[1], tempy[2],
336 t0 + h[i]/2.0, MHD_INCLUDE_ALL, boozer,
337 mhd, Bdata);
338 }
339 if(!errflag) {
340 step_gceom_mhd(k3, tempy, mass, charge, B_dB, E, mhd_dmhd,
341 aldforce);
342 }
343 for(int j = 0; j < 6; j++) {
344 tempy[j] = yprev[j] + h[i]*k3[j];
345 }
346
347 if(!errflag) {
348 errflag = B_field_eval_B_dB(B_dB, tempy[0], tempy[1], tempy[2],
349 t0 + h[i], Bdata);
350 }
351 if(!errflag) {
352 errflag = E_field_eval_E(E, tempy[0], tempy[1], tempy[2],
353 t0 + h[i], Edata, Bdata);
354 }
355 if(!errflag) {
356 errflag = mhd_eval(mhd_dmhd, tempy[0], tempy[1], tempy[2],
357 t0 + h[i], MHD_INCLUDE_ALL, boozer,
358 mhd, Bdata);
359 }
360 if(!errflag) {
361 step_gceom_mhd(k4, tempy, mass, charge, B_dB, E, mhd_dmhd,
362 aldforce);
363 }
364 for(int j = 0; j < 6; j++) {
365 y[j] = yprev[j]
366 + h[i]/6.0 * (k1[j] + 2*k2[j] + 2*k3[j] + k4[j]);
367 }
368
369 /* Test that results are physical */
370 if(!errflag && y[0] <= 0) {
371 errflag = error_raise(ERR_INTEGRATION, __LINE__, EF_STEP_GC_RK4);
372 }
373 else if(!errflag && fabs(y[4]) >= CONST_C) {
374 errflag = error_raise(ERR_INTEGRATION, __LINE__, EF_STEP_GC_RK4);
375 }
376 else if(!errflag && y[4] < 0) {
377 errflag = error_raise(ERR_INTEGRATION, __LINE__, EF_STEP_GC_RK4);
378 }
379
380 /* Update gc phase space position */
381 if(!errflag) {
382 p->r[i] = y[0];
383 p->phi[i] = y[1];
384 p->z[i] = y[2];
385 p->ppar[i] = y[3];
386 p->mu[i] = y[4];
387 p->zeta[i] = fmod(y[5],CONST_2PI);
388 if(p->zeta[i]<0){
389 p->zeta[i] = CONST_2PI + p->zeta[i];
390 }
391 }
392
393 /* Evaluate magnetic field (and gradient) and rho at new position */
394 real psi[1];
395 real rho[2];
396 if(!errflag) {
397 errflag = B_field_eval_B_dB(B_dB, p->r[i], p->phi[i], p->z[i],
398 t0 + h[i], Bdata);
399 }
400 if(!errflag) {
401 errflag = B_field_eval_psi(psi, p->r[i], p->phi[i], p->z[i],
402 t0 + h[i], Bdata);
403 }
404 if(!errflag) {
405 errflag = B_field_eval_rho(rho, psi[0], Bdata);
406 }
407
408 if(!errflag) {
409 p->B_r[i] = B_dB[0];
410 p->B_r_dr[i] = B_dB[1];
411 p->B_r_dphi[i] = B_dB[2];
412 p->B_r_dz[i] = B_dB[3];
413
414 p->B_phi[i] = B_dB[4];
415 p->B_phi_dr[i] = B_dB[5];
416 p->B_phi_dphi[i] = B_dB[6];
417 p->B_phi_dz[i] = B_dB[7];
418
419 p->B_z[i] = B_dB[8];
420 p->B_z_dr[i] = B_dB[9];
421 p->B_z_dphi[i] = B_dB[10];
422 p->B_z_dz[i] = B_dB[11];
423
424 p->rho[i] = rho[0];
425
426 /* Evaluate pol angle so that it is cumulative */
427 real axisrz[2];
428 errflag = B_field_get_axis_rz(axisrz, Bdata, p->phi[i]);
429 p->theta[i] += atan2( (R0-axisrz[0]) * (p->z[i]-axisrz[1])
430 - (z0-axisrz[1]) * (p->r[i]-axisrz[0]),
431 (R0-axisrz[0]) * (p->r[i]-axisrz[0])
432 + (z0-axisrz[1]) * (p->z[i]-axisrz[1]) );
433 }
434
435 /* Error handling */
436 if(errflag) {
437 p->err[i] = errflag;
438 p->running[i] = 0;
439 }
440 }
441 }
442
443}
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
E_field.h
Header file for E_field.c.
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
boozer.h
Header file for boozer.c.
consts.h
Header file containing physical and mathematical constants.
CONST_C
#define CONST_C
Speed of light [m/s].
Definition consts.h:23
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_RK4
@ EF_STEP_GC_RK4
Definition error.h:32
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.h
Header file for mhd.c.
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_rk4
void step_gc_rk4(particle_simd_gc *p, real *h, B_field_data *Bdata, E_field_data *Edata, int aldforce)
Integrate a guiding center step for a struct of markers with RK4.
Definition step_gc_rk4.c:34
step_gc_rk4_mhd
void step_gc_rk4_mhd(particle_simd_gc *p, real *h, B_field_data *Bdata, E_field_data *Edata, boozer_data *boozer, mhd_data *mhd, int aldforce)
Integrate a guiding center step with RK4 with MHD modes present.
Definition step_gc_rk4.c:234
step_gc_rk4.h
Header file for step_gc_rk4.c.
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