ASCOT5
Loading...
Searching...
No Matches
step_fo_vpa.c
Go to the documentation of this file.
1
5#include <math.h>
6#include <stdio.h>
7#include "../../offload.h"
8#include "../../ascot5.h"
9#include "../../math.h"
10#include "../../consts.h"
11#include "../../physlib.h"
12#include "../../error.h"
13#include "../../B_field.h"
14#include "../../E_field.h"
15#include "../../boozer.h"
16#include "../../mhd.h"
17#include "../../particle.h"
18#include "step_fo_vpa.h"
19
35void step_fo_vpa(particle_simd_fo* p, real* h, B_field_data* Bdata,
36 E_field_data* Edata) {
37 GPU_DATA_IS_MAPPED(h[0:p->n_mrk])
38 GPU_PARALLEL_LOOP_ALL_LEVELS
39 for(int i = 0; i < p->n_mrk; i++) {
40 if(p->running[i]) {
41 a5err errflag = 0;
42
43 real R0 = p->r[i];
44 real z0 = p->z[i];
45 real t0 = p->time[i];
46 real mass = p->mass[i];
47
48 /* Convert velocity to cartesian coordinates */
49 real prpz[3] = {p->p_r[i], p->p_phi[i], p->p_z[i]};
50 real pxyz[3];
51 math_vec_rpz2xyz(prpz, pxyz, p->phi[i]);
52
53 real posrpz[3] = {p->r[i], p->phi[i], p->z[i]};
54 real posxyz0[3],posxyz[3];
55 math_rpz2xyz(posrpz, posxyz0);
56
57 /* Take a half step and evaluate fields at that position */
58 real gamma = physlib_gamma_pnorm(mass, math_norm(pxyz));
59 posxyz[0] = posxyz0[0] + pxyz[0] * h[i] / (2.0 * gamma * mass);
60 posxyz[1] = posxyz0[1] + pxyz[1] * h[i] / (2.0 * gamma * mass);
61 posxyz[2] = posxyz0[2] + pxyz[2] * h[i] / (2.0 * gamma * mass);
62
63 math_xyz2rpz(posxyz, posrpz);
64
65 real Brpz[3];
66 real Erpz[3];
67 if(!errflag) {
68 errflag = B_field_eval_B(Brpz, posrpz[0], posrpz[1], posrpz[2],
69 t0 + h[i]/2.0, Bdata);
70 }
71 if(!errflag) {
72 errflag = E_field_eval_E(Erpz, posrpz[0], posrpz[1], posrpz[2],
73 t0 + h[i]/2.0, Edata, Bdata);
74 }
75
76 real fposxyz[3]; // final position in cartesian coordinates
77
78 if(!errflag) {
79 /* Electromagnetic fields to cartesian coordinates */
80 real Bxyz[3];
81 real Exyz[3];
82
83 math_vec_rpz2xyz(Brpz, Bxyz, posrpz[1]);
84 math_vec_rpz2xyz(Erpz, Exyz, posrpz[1]);
85
86 /* Evaluate helper variable pminus */
87 real pminus[3];
88 real sigma = p->charge[i]*h[i]/(2*p->mass[i]*CONST_C);
89 pminus[0] = pxyz[0] / (mass * CONST_C) + sigma * Exyz[0];
90 pminus[1] = pxyz[1] / (mass * CONST_C) + sigma * Exyz[1];
91 pminus[2] = pxyz[2] / (mass * CONST_C) + sigma * Exyz[2];
92
93 /* Second helper variable pplus*/
94 real d = (p->charge[i]*h[i]/(2*p->mass[i])) /
95 sqrt( 1 + math_dot(pminus,pminus) );
96 real d2 = d*d;
97
98 real Bhat[9] = { 0, Bxyz[2], -Bxyz[1],
99 -Bxyz[2], 0, Bxyz[0],
100 Bxyz[1], -Bxyz[0], 0};
101 real Bhat2[9];
102 math_matmul(Bhat, Bhat, 3, 3, 3, Bhat2);
103
104 real B2 = Bxyz[0]*Bxyz[0] + Bxyz[1]*Bxyz[1] + Bxyz[2]*Bxyz[2];
105
106 real A[9];
107 for(int j=0; j<9; j++) {
108 A[j] = (Bhat[j] + d*Bhat2[j]) * (2.0*d/(1.0+d2*B2));
109 }
110
111 real pplus[3];
112 math_matmul(pminus, A, 1, 3, 3, pplus);
113
114 /* Take the step */
115 real pfinal[3];
116 pfinal[0] = pminus[0] + pplus[0] + sigma*Exyz[0];
117 pfinal[1] = pminus[1] + pplus[1] + sigma*Exyz[1];
118 pfinal[2] = pminus[2] + pplus[2] + sigma*Exyz[2];
119
120 pxyz[0] = pfinal[0] * mass * CONST_C;
121 pxyz[1] = pfinal[1] * mass * CONST_C;
122 pxyz[2] = pfinal[2] * mass * CONST_C;
123 }
124
125 gamma = physlib_gamma_pnorm(mass, math_norm(pxyz));
126 fposxyz[0] = posxyz[0] + h[i] * pxyz[0] / (2.0 * gamma * mass);
127 fposxyz[1] = posxyz[1] + h[i] * pxyz[1] / (2.0 * gamma * mass);
128 fposxyz[2] = posxyz[2] + h[i] * pxyz[2] / (2.0 * gamma * mass);
129
130 if(!errflag) {
131 /* Back to cylindrical coordinates */
132 p->r[i] = sqrt(fposxyz[0]*fposxyz[0]+fposxyz[1]*fposxyz[1]);
133
134 /* phi is evaluated like this to make sure it is cumulative */
135 p->phi[i] += atan2(
136 posxyz0[0] * fposxyz[1] - posxyz0[1] * fposxyz[0],
137 posxyz0[0] * fposxyz[0] + posxyz0[1] * fposxyz[1] );
138 p->z[i] = fposxyz[2];
139
140 real cosp = cos(p->phi[i]);
141 real sinp = sin(p->phi[i]);
142 p->p_r[i] = pxyz[0] * cosp + pxyz[1] * sinp;
143 p->p_phi[i] = -pxyz[0] * sinp + pxyz[1] * cosp;
144 p->p_z[i] = pxyz[2];
145 }
146
147 /* Evaluate magnetic field (and gradient) and rho at new position */
148 real BdBrpz[15];
149 real psi[1];
150 real rho[2];
151 if(!errflag) {
152 errflag = B_field_eval_B_dB(BdBrpz, p->r[i], p->phi[i], p->z[i],
153 t0 + h[i], Bdata);
154 }
155 if(!errflag) {
156 errflag = B_field_eval_psi(psi, p->r[i], p->phi[i], p->z[i],
157 t0 + h[i], Bdata);
158 }
159 if(!errflag) {
160 errflag = B_field_eval_rho(rho, psi[0], Bdata);
161 }
162
163 if(!errflag) {
164 p->B_r[i] = BdBrpz[0];
165 p->B_r_dr[i] = BdBrpz[1];
166 p->B_r_dphi[i] = BdBrpz[2];
167 p->B_r_dz[i] = BdBrpz[3];
168
169 p->B_phi[i] = BdBrpz[4];
170 p->B_phi_dr[i] = BdBrpz[5];
171 p->B_phi_dphi[i] = BdBrpz[6];
172 p->B_phi_dz[i] = BdBrpz[7];
173
174 p->B_z[i] = BdBrpz[8];
175 p->B_z_dr[i] = BdBrpz[9];
176 p->B_z_dphi[i] = BdBrpz[10];
177 p->B_z_dz[i] = BdBrpz[11];
178
179 p->rho[i] = rho[0];
180
181 /* Evaluate phi and theta angles so that they are cumulative */
182 real axisrz[2];
183 errflag = B_field_get_axis_rz(axisrz, Bdata, p->phi[i]);
184 p->theta[i] += atan2( (R0-axisrz[0]) * (p->z[i]-axisrz[1])
185 - (z0-axisrz[1]) * (p->r[i]-axisrz[0]),
186 (R0-axisrz[0]) * (p->r[i]-axisrz[0])
187 + (z0-axisrz[1]) * (p->z[i]-axisrz[1]) );
188 }
189
190 /* Error handling */
191 if(errflag) {
192 p->err[i] = errflag;
193 p->running[i] = 0;
194 }
195 }
196 }
197}
198
199
212void step_fo_vpa_mhd(particle_simd_fo* p, real* h, B_field_data* Bdata,
213 E_field_data* Edata, boozer_data* boozer, mhd_data* mhd) {
214
215 int i;
216 /* Following loop will be executed simultaneously for all i */
217 #pragma omp simd aligned(h : 64)
218 for(i = 0; i < NSIMD; i++) {
219 if(p->running[i]) {
220 a5err errflag = 0;
221
222 real R0 = p->r[i];
223 real z0 = p->z[i];
224 real t0 = p->time[i];
225 real mass = p->mass[i];
226
227 /* Convert velocity to cartesian coordinates */
228 real prpz[3] = {p->p_r[i], p->p_phi[i], p->p_z[i]};
229 real pxyz[3];
230 math_vec_rpz2xyz(prpz, pxyz, p->phi[i]);
231
232 real posrpz[3] = {p->r[i], p->phi[i], p->z[i]};
233 real posxyz0[3],posxyz[3];
234 math_rpz2xyz(posrpz,posxyz0);
235
236 /* Take a half step and evaluate fields at that position */
237 real gamma = physlib_gamma_pnorm(mass, math_norm(pxyz));
238 posxyz[0] = posxyz0[0] + pxyz[0] * h[i] / (2 * gamma * mass);
239 posxyz[1] = posxyz0[1] + pxyz[1] * h[i] / (2 * gamma * mass);
240 posxyz[2] = posxyz0[2] + pxyz[2] * h[i] / (2 * gamma * mass);
241
242 math_xyz2rpz(posxyz,posrpz);
243
244 real Brpz[3];
245 real Erpz[3];
246 if(!errflag) {
247 errflag = B_field_eval_B(Brpz, posrpz[0], posrpz[1], posrpz[2],
248 t0 + h[i]/2, Bdata);
249 }
250 if(!errflag) {
251 errflag = E_field_eval_E(Erpz, posrpz[0], posrpz[1], posrpz[2],
252 t0 + h[i]/2, Edata, Bdata);
253 }
254
255 real pert[7];
256 int pertonly = 0;
257 if(!errflag) {
258 errflag = mhd_perturbations(
259 pert, posrpz[0], posrpz[1], posrpz[2], t0 + h[i]/2,
260 pertonly, MHD_INCLUDE_ALL, boozer, mhd, Bdata);
261 }
262 Brpz[0] = pert[0];
263 Brpz[1] = pert[1];
264 Brpz[2] = pert[2];
265 Erpz[0] += pert[3];
266 Erpz[1] += pert[4];
267 Erpz[2] += pert[5];
268
269 real fposxyz[3]; // final position in cartesian coordinates
270
271 if(!errflag) {
272 /* Electromagnetic fields to cartesian coordinates */
273 real Bxyz[3];
274 real Exyz[3];
275
276 math_vec_rpz2xyz(Brpz, Bxyz, posrpz[1]);
277 math_vec_rpz2xyz(Erpz, Exyz, posrpz[1]);
278
279 /* Evaluate helper variable pminus */
280 real pminus[3];
281 real sigma = p->charge[i]*h[i]/(2*p->mass[i]*CONST_C);
282 pminus[0] = pxyz[0] / (mass * CONST_C) + sigma * Exyz[0];
283 pminus[1] = pxyz[1] / (mass * CONST_C) + sigma * Exyz[1];
284 pminus[2] = pxyz[2] / (mass * CONST_C) + sigma * Exyz[2];
285
286 /* Second helper variable pplus*/
287 real d = (p->charge[i]*h[i]/(2*p->mass[i])) /
288 sqrt( 1 + math_dot(pminus,pminus) );
289 real d2 = d*d;
290
291 real Bhat[9] = { 0, Bxyz[2], -Bxyz[1],
292 -Bxyz[2], 0, Bxyz[0],
293 Bxyz[1], -Bxyz[0], 0};
294 real Bhat2[9];
295 math_matmul(Bhat, Bhat, 3, 3, 3, Bhat2);
296
297 real B2 = Bxyz[0]*Bxyz[0] + Bxyz[1]*Bxyz[1] + Bxyz[2]*Bxyz[2];
298
299 real A[9];
300 for(int j=0; j<9; j++) {
301 A[j] = (Bhat[j] + d*Bhat2[j]) * (2.0*d/(1+d2*B2));
302 }
303
304 real pplus[3];
305 math_matmul(pminus, A, 1, 3, 3, pplus);
306
307 /* Take the step */
308 real pfinal[3];
309 pfinal[0] = pminus[0] + pplus[0] + sigma*Exyz[0];
310 pfinal[1] = pminus[1] + pplus[1] + sigma*Exyz[1];
311 pfinal[2] = pminus[2] + pplus[2] + sigma*Exyz[2];
312
313 pxyz[0] = pfinal[0] * mass * CONST_C;
314 pxyz[1] = pfinal[1] * mass * CONST_C;
315 pxyz[2] = pfinal[2] * mass * CONST_C;
316 }
317
318 gamma = physlib_gamma_pnorm(mass, math_norm(pxyz));
319 fposxyz[0] = posxyz[0] + h[i] * pxyz[0] / (2 * gamma * mass);
320 fposxyz[1] = posxyz[1] + h[i] * pxyz[1] / (2 * gamma * mass);
321 fposxyz[2] = posxyz[2] + h[i] * pxyz[2] / (2 * gamma * mass);
322
323 if(!errflag) {
324 /* Back to cylindrical coordinates */
325 p->r[i] = sqrt(fposxyz[0]*fposxyz[0]+fposxyz[1]*fposxyz[1]);
326
327 /* phi is evaluated like this to make sure it is cumulative */
328 p->phi[i] += atan2(
329 posxyz0[0] * fposxyz[1] - posxyz0[1] * fposxyz[0],
330 posxyz0[0] * fposxyz[0] + posxyz0[1] * fposxyz[1] );
331 p->z[i] = fposxyz[2];
332
333 real cosp = cos(p->phi[i]);
334 real sinp = sin(p->phi[i]);
335 p->p_r[i] = pxyz[0] * cosp + pxyz[1] * sinp;
336 p->p_phi[i] = -pxyz[0] * sinp + pxyz[1] * cosp;
337 p->p_z[i] = pxyz[2];
338 }
339
340 /* Evaluate magnetic field (and gradient) and rho at new position */
341 real BdBrpz[15];
342 real psi[1];
343 real rho[2];
344 if(!errflag) {
345 errflag = B_field_eval_B_dB(BdBrpz, p->r[i], p->phi[i], p->z[i],
346 t0 + h[i], Bdata);
347 }
348 if(!errflag) {
349 errflag = B_field_eval_psi(psi, p->r[i], p->phi[i], p->z[i],
350 t0 + h[i], Bdata);
351 }
352 if(!errflag) {
353 errflag = B_field_eval_rho(rho, psi[0], Bdata);
354 }
355
356 if(!errflag) {
357 p->B_r[i] = BdBrpz[0];
358 p->B_r_dr[i] = BdBrpz[1];
359 p->B_r_dphi[i] = BdBrpz[2];
360 p->B_r_dz[i] = BdBrpz[3];
361
362 p->B_phi[i] = BdBrpz[4];
363 p->B_phi_dr[i] = BdBrpz[5];
364 p->B_phi_dphi[i] = BdBrpz[6];
365 p->B_phi_dz[i] = BdBrpz[7];
366
367 p->B_z[i] = BdBrpz[8];
368 p->B_z_dr[i] = BdBrpz[9];
369 p->B_z_dphi[i] = BdBrpz[10];
370 p->B_z_dz[i] = BdBrpz[11];
371
372 p->rho[i] = rho[0];
373
374 /* Evaluate phi and theta angles so that they are cumulative */
375 real axisrz[2];
376 errflag = B_field_get_axis_rz(axisrz, Bdata, p->phi[i]);
377 p->theta[i] += atan2( (R0-axisrz[0]) * (p->z[i]-axisrz[1])
378 - (z0-axisrz[1]) * (p->r[i]-axisrz[0]),
379 (R0-axisrz[0]) * (p->r[i]-axisrz[0])
380 + (z0-axisrz[1]) * (p->z[i]-axisrz[1]) );
381 }
382
383 /* Error handling */
384 if(errflag) {
385 p->err[i] = errflag;
386 p->running[i] = 0;
387 }
388 }
389 }
390}
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_eval_B
a5err B_field_eval_B(real B[3], real r, real phi, real z, real t, B_field_data *Bdata)
Evaluate magnetic field.
Definition B_field.c:374
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
error.h
Error module for ASCOT5.
a5err
unsigned long int a5err
Simulation error flag.
Definition error.h:17
math_matmul
void math_matmul(real *matA, real *matB, int d1, int d2, int d3, real *matC)
Matrix multiplication.
Definition math.c:159
math.h
Header file for math.c.
math_dot
#define math_dot(a, b)
Calculate dot product a[3] dot b[3].
Definition math.h:28
math_xyz2rpz
#define math_xyz2rpz(xyz, rpz)
Convert cartesian coordinates xyz to cylindrical coordinates rpz.
Definition math.h:74
math_vec_rpz2xyz
#define math_vec_rpz2xyz(vrpz, vxyz, phi)
Transform vector from cylindrical to cartesian basis: vrpz -> vxyz, phi is the toroidal angle in radi...
Definition math.h:83
math_rpz2xyz
#define math_rpz2xyz(rpz, xyz)
Convert cylindrical coordinates rpz to cartesian coordinates xyz.
Definition math.h:78
math_norm
#define math_norm(a)
Calculate norm of 3D vector a.
Definition math.h:64
mhd_perturbations
a5err mhd_perturbations(real pert_field[7], real r, real phi, real z, real t, int pertonly, int includemode, boozer_data *boozerdata, mhd_data *mhddata, B_field_data *Bdata)
Evaluate perturbed fields Btilde, Etilde and potential Phi explicitly.
Definition mhd.c:147
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.
physlib.h
Methods to evaluate elementary physical quantities.
physlib_gamma_pnorm
#define physlib_gamma_pnorm(m, p)
Evaluate Lorentz factor from momentum norm.
Definition physlib.h:46
step_fo_vpa_mhd
void step_fo_vpa_mhd(particle_simd_fo *p, real *h, B_field_data *Bdata, E_field_data *Edata, boozer_data *boozer, mhd_data *mhd)
Integrate a full orbit step with VPA and MHd modes present.
Definition step_fo_vpa.c:212
step_fo_vpa
void step_fo_vpa(particle_simd_fo *p, real *h, B_field_data *Bdata, E_field_data *Edata)
Integrate a full orbit step for a struct of particles with VPA.
Definition step_fo_vpa.c:35
step_fo_vpa.h
Header file for step_fo_vpa.c.
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_fo
Struct representing NSIMD particle markers.
Definition particle.h:210
particle_simd_fo::p_phi
real * p_phi
Definition particle.h:216
particle_simd_fo::time
real * time
Definition particle.h:220
particle_simd_fo::charge
real * charge
Definition particle.h:219
particle_simd_fo::r
real * r
Definition particle.h:212
particle_simd_fo::B_phi_dphi
real * B_phi_dphi
Definition particle.h:233
particle_simd_fo::B_r_dz
real * B_r_dz
Definition particle.h:235
particle_simd_fo::B_z_dz
real * B_z_dz
Definition particle.h:237
particle_simd_fo::n_mrk
size_t n_mrk
Definition particle.h:256
particle_simd_fo::B_r_dphi
real * B_r_dphi
Definition particle.h:232
particle_simd_fo::rho
real * rho
Definition particle.h:243
particle_simd_fo::B_phi_dr
real * B_phi_dr
Definition particle.h:230
particle_simd_fo::B_phi
real * B_phi
Definition particle.h:226
particle_simd_fo::B_z_dr
real * B_z_dr
Definition particle.h:231
particle_simd_fo::p_z
real * p_z
Definition particle.h:217
particle_simd_fo::B_z
real * B_z
Definition particle.h:227
particle_simd_fo::err
a5err * err
Definition particle.h:254
particle_simd_fo::B_z_dphi
real * B_z_dphi
Definition particle.h:234
particle_simd_fo::B_r
real * B_r
Definition particle.h:225
particle_simd_fo::theta
real * theta
Definition particle.h:244
particle_simd_fo::p_r
real * p_r
Definition particle.h:215
particle_simd_fo::running
integer * running
Definition particle.h:252
particle_simd_fo::mass
real * mass
Definition particle.h:218
particle_simd_fo::B_r_dr
real * B_r_dr
Definition particle.h:229
particle_simd_fo::B_phi_dz
real * B_phi_dz
Definition particle.h:236
particle_simd_fo::phi
real * phi
Definition particle.h:213
particle_simd_fo::z
real * z
Definition particle.h:214
Generated on Wed Mar 19 2025 11:10:23 for ASCOT5 by doxygen 1.13.2