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