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