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
36void step_fo_vpa(particle_simd_fo* p, real* h, B_field_data* Bdata,
37 E_field_data* Edata, int aldforce) {
38 GPU_DATA_IS_MAPPED(h[0:p->n_mrk])
39 GPU_PARALLEL_LOOP_ALL_LEVELS
40 for(int i = 0; i < p->n_mrk; i++) {
41 if(p->running[i]) {
42 a5err errflag = 0;
43
44 real R0 = p->r[i];
45 real z0 = p->z[i];
46 real t0 = p->time[i];
47 real mass = p->mass[i];
48
49 /* Convert velocity to cartesian coordinates */
50 real prpz[3] = {p->p_r[i], p->p_phi[i], p->p_z[i]};
51 real pxyz[3];
52 math_vec_rpz2xyz(prpz, pxyz, p->phi[i]);
53
54 real posrpz[3] = {p->r[i], p->phi[i], p->z[i]};
55 real posxyz0[3],posxyz[3];
56 math_rpz2xyz(posrpz, posxyz0);
57
58 /* Take a half step and evaluate fields at that position */
59 real gamma = physlib_gamma_pnorm(mass, math_norm(pxyz));
60 posxyz[0] = posxyz0[0] + pxyz[0] * h[i] / (2.0 * gamma * mass);
61 posxyz[1] = posxyz0[1] + pxyz[1] * h[i] / (2.0 * gamma * mass);
62 posxyz[2] = posxyz0[2] + pxyz[2] * h[i] / (2.0 * gamma * mass);
63
64 math_xyz2rpz(posxyz, posrpz);
65
66 real Brpz[3];
67 real Erpz[3];
68 if(!errflag) {
69 errflag = B_field_eval_B(Brpz, posrpz[0], posrpz[1], posrpz[2],
70 t0 + h[i]/2.0, Bdata);
71 }
72 if(!errflag) {
73 errflag = E_field_eval_E(Erpz, posrpz[0], posrpz[1], posrpz[2],
74 t0 + h[i]/2.0, Edata, Bdata);
75 }
76
77 real fposxyz[3]; // final position in cartesian coordinates
78
79 if(!errflag) {
80 /* Electromagnetic fields to cartesian coordinates */
81 real Bxyz[3];
82 real Exyz[3];
83
84 math_vec_rpz2xyz(Brpz, Bxyz, posrpz[1]);
85 math_vec_rpz2xyz(Erpz, Exyz, posrpz[1]);
86
87 /* Evaluate helper variable pminus */
88 real pminus[3];
89 real sigma = p->charge[i]*h[i]/(2*p->mass[i]*CONST_C);
90 pminus[0] = pxyz[0] / (mass * CONST_C) + sigma * Exyz[0];
91 pminus[1] = pxyz[1] / (mass * CONST_C) + sigma * Exyz[1];
92 pminus[2] = pxyz[2] / (mass * CONST_C) + sigma * Exyz[2];
93
94 /* Second helper variable pplus*/
95 real d = (p->charge[i]*h[i]/(2*p->mass[i])) /
96 sqrt( 1 + math_dot(pminus,pminus) );
97 real d2 = d*d;
98
99 real Bhat[9] = { 0, Bxyz[2], -Bxyz[1],
100 -Bxyz[2], 0, Bxyz[0],
101 Bxyz[1], -Bxyz[0], 0};
102 real Bhat2[9];
103 math_matmul(Bhat, Bhat, 3, 3, 3, Bhat2);
104
105 real B2 = Bxyz[0]*Bxyz[0] + Bxyz[1]*Bxyz[1] + Bxyz[2]*Bxyz[2];
106
107 real A[9];
108 for(int j=0; j<9; j++) {
109 A[j] = (Bhat[j] + d*Bhat2[j]) * (2.0*d/(1.0+d2*B2));
110 }
111
112 real pplus[3];
113 math_matmul(pminus, A, 1, 3, 3, pplus);
114
115 /* Take the step */
116 real pfinal[3];
117 pfinal[0] = pminus[0] + pplus[0] + sigma*Exyz[0];
118 pfinal[1] = pminus[1] + pplus[1] + sigma*Exyz[1];
119 pfinal[2] = pminus[2] + pplus[2] + sigma*Exyz[2];
120
121 pxyz[0] = pfinal[0] * mass * CONST_C;
122 pxyz[1] = pfinal[1] * mass * CONST_C;
123 pxyz[2] = pfinal[2] * mass * CONST_C;
124 }
125
126 gamma = physlib_gamma_pnorm(mass, math_norm(pxyz));
127 fposxyz[0] = posxyz[0] + h[i] * pxyz[0] / (2.0 * gamma * mass);
128 fposxyz[1] = posxyz[1] + h[i] * pxyz[1] / (2.0 * gamma * mass);
129 fposxyz[2] = posxyz[2] + h[i] * pxyz[2] / (2.0 * gamma * mass);
130
131 if(!errflag) {
132 /* Back to cylindrical coordinates */
133 p->r[i] = sqrt(fposxyz[0]*fposxyz[0]+fposxyz[1]*fposxyz[1]);
134
135 /* phi is evaluated like this to make sure it is cumulative */
136 p->phi[i] += atan2(
137 posxyz0[0] * fposxyz[1] - posxyz0[1] * fposxyz[0],
138 posxyz0[0] * fposxyz[0] + posxyz0[1] * fposxyz[1] );
139 p->z[i] = fposxyz[2];
140
141 real cosp = cos(p->phi[i]);
142 real sinp = sin(p->phi[i]);
143 p->p_r[i] = pxyz[0] * cosp + pxyz[1] * sinp;
144 p->p_phi[i] = -pxyz[0] * sinp + pxyz[1] * cosp;
145 p->p_z[i] = pxyz[2];
146 }
147
148 /* Evaluate magnetic field (and gradient) and rho at new position */
149 real BdBrpz[15];
150 real psi[1];
151 real rho[2];
152 if(!errflag) {
153 errflag = B_field_eval_B_dB(BdBrpz, p->r[i], p->phi[i], p->z[i],
154 t0 + h[i], Bdata);
155 }
156 if(!errflag) {
157 errflag = B_field_eval_psi(psi, p->r[i], p->phi[i], p->z[i],
158 t0 + h[i], Bdata);
159 }
160 if(!errflag) {
161 errflag = B_field_eval_rho(rho, psi[0], Bdata);
162 }
163
164 if(!errflag) {
165 p->B_r[i] = BdBrpz[0];
166 p->B_r_dr[i] = BdBrpz[1];
167 p->B_r_dphi[i] = BdBrpz[2];
168 p->B_r_dz[i] = BdBrpz[3];
169
170 p->B_phi[i] = BdBrpz[4];
171 p->B_phi_dr[i] = BdBrpz[5];
172 p->B_phi_dphi[i] = BdBrpz[6];
173 p->B_phi_dz[i] = BdBrpz[7];
174
175 p->B_z[i] = BdBrpz[8];
176 p->B_z_dr[i] = BdBrpz[9];
177 p->B_z_dphi[i] = BdBrpz[10];
178 p->B_z_dz[i] = BdBrpz[11];
179
180 p->rho[i] = rho[0];
181
182 /* Evaluate phi and theta angles so that they are cumulative */
183 real axisrz[2];
184 errflag = B_field_get_axis_rz(axisrz, Bdata, p->phi[i]);
185 p->theta[i] += atan2( (R0-axisrz[0]) * (p->z[i]-axisrz[1])
186 - (z0-axisrz[1]) * (p->r[i]-axisrz[0]),
187 (R0-axisrz[0]) * (p->r[i]-axisrz[0])
188 + (z0-axisrz[1]) * (p->z[i]-axisrz[1]) );
189 }
190
191 /* Evaluate Abraham-Lorentz-Dirac force (if enabled) is evaluated
192 * separately using the Euler method */
193 real Bnorm = math_normc(p->B_r[i], p->B_phi[i], p->B_z[i]);
194 real pnorm = math_normc(p->p_r[i], p->p_phi[i], p->p_z[i]);
195 real t_ald = phys_ald_force_chartime(
196 p->charge[i], p->mass[i], Bnorm, gamma) * aldforce;
197 real pparbhatperB = (
198 p->p_r[i]*p->B_r[i] + p->p_phi[i]*p->B_phi[i]
199 + p->p_z[i]*p->B_z[i] ) / ( Bnorm * Bnorm * pnorm );
200 real pperpvec[3] = {
201 p->p_r[i] - pparbhatperB * p->B_r[i],
202 p->p_phi[i] - pparbhatperB * p->B_phi[i],
203 p->p_z[i] - pparbhatperB * p->B_z[i] };
204 real C = ( pperpvec[0]*pperpvec[0] + pperpvec[1]*pperpvec[1]
205 + pperpvec[2]*pperpvec[2] )
206 / ( p->mass[i]*p->mass[i] * CONST_C2 );
207 p->p_r[i] -= t_ald * ( pperpvec[0] + C * p->p_r[i] );
208 p->p_phi[i] -= t_ald * ( pperpvec[1] + C * p->p_phi[i] );
209 p->p_z[i] -= t_ald * ( pperpvec[2] + C * p->p_z[i] );
210
211 /* Error handling */
212 if(errflag) {
213 p->err[i] = errflag;
214 p->running[i] = 0;
215 }
216 }
217 }
218}
219
220
234void step_fo_vpa_mhd(
235 particle_simd_fo* 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 R0 = p->r[i];
246 real z0 = p->z[i];
247 real t0 = p->time[i];
248 real mass = p->mass[i];
249
250 /* Convert velocity to cartesian coordinates */
251 real prpz[3] = {p->p_r[i], p->p_phi[i], p->p_z[i]};
252 real pxyz[3];
253 math_vec_rpz2xyz(prpz, pxyz, p->phi[i]);
254
255 real posrpz[3] = {p->r[i], p->phi[i], p->z[i]};
256 real posxyz0[3],posxyz[3];
257 math_rpz2xyz(posrpz,posxyz0);
258
259 /* Take a half step and evaluate fields at that position */
260 real gamma = physlib_gamma_pnorm(mass, math_norm(pxyz));
261 posxyz[0] = posxyz0[0] + pxyz[0] * h[i] / (2 * gamma * mass);
262 posxyz[1] = posxyz0[1] + pxyz[1] * h[i] / (2 * gamma * mass);
263 posxyz[2] = posxyz0[2] + pxyz[2] * h[i] / (2 * gamma * mass);
264
265 math_xyz2rpz(posxyz,posrpz);
266
267 real Brpz[3];
268 real Erpz[3];
269 if(!errflag) {
270 errflag = B_field_eval_B(Brpz, posrpz[0], posrpz[1], posrpz[2],
271 t0 + h[i]/2, Bdata);
272 }
273 if(!errflag) {
274 errflag = E_field_eval_E(Erpz, posrpz[0], posrpz[1], posrpz[2],
275 t0 + h[i]/2, Edata, Bdata);
276 }
277
278 real pert[7];
279 int pertonly = 0;
280 if(!errflag) {
281 errflag = mhd_perturbations(
282 pert, posrpz[0], posrpz[1], posrpz[2], t0 + h[i]/2,
283 pertonly, MHD_INCLUDE_ALL, boozer, mhd, Bdata);
284 }
285 Brpz[0] = pert[0];
286 Brpz[1] = pert[1];
287 Brpz[2] = pert[2];
288 Erpz[0] += pert[3];
289 Erpz[1] += pert[4];
290 Erpz[2] += pert[5];
291
292 real fposxyz[3]; // final position in cartesian coordinates
293
294 if(!errflag) {
295 /* Electromagnetic fields to cartesian coordinates */
296 real Bxyz[3];
297 real Exyz[3];
298
299 math_vec_rpz2xyz(Brpz, Bxyz, posrpz[1]);
300 math_vec_rpz2xyz(Erpz, Exyz, posrpz[1]);
301
302 /* Evaluate helper variable pminus */
303 real pminus[3];
304 real sigma = p->charge[i]*h[i]/(2*p->mass[i]*CONST_C);
305 pminus[0] = pxyz[0] / (mass * CONST_C) + sigma * Exyz[0];
306 pminus[1] = pxyz[1] / (mass * CONST_C) + sigma * Exyz[1];
307 pminus[2] = pxyz[2] / (mass * CONST_C) + sigma * Exyz[2];
308
309 /* Second helper variable pplus*/
310 real d = (p->charge[i]*h[i]/(2*p->mass[i])) /
311 sqrt( 1 + math_dot(pminus,pminus) );
312 real d2 = d*d;
313
314 real Bhat[9] = { 0, Bxyz[2], -Bxyz[1],
315 -Bxyz[2], 0, Bxyz[0],
316 Bxyz[1], -Bxyz[0], 0};
317 real Bhat2[9];
318 math_matmul(Bhat, Bhat, 3, 3, 3, Bhat2);
319
320 real B2 = Bxyz[0]*Bxyz[0] + Bxyz[1]*Bxyz[1] + Bxyz[2]*Bxyz[2];
321
322 real A[9];
323 for(int j=0; j<9; j++) {
324 A[j] = (Bhat[j] + d*Bhat2[j]) * (2.0*d/(1+d2*B2));
325 }
326
327 real pplus[3];
328 math_matmul(pminus, A, 1, 3, 3, pplus);
329
330 /* Take the step */
331 real pfinal[3];
332 pfinal[0] = pminus[0] + pplus[0] + sigma*Exyz[0];
333 pfinal[1] = pminus[1] + pplus[1] + sigma*Exyz[1];
334 pfinal[2] = pminus[2] + pplus[2] + sigma*Exyz[2];
335
336 pxyz[0] = pfinal[0] * mass * CONST_C;
337 pxyz[1] = pfinal[1] * mass * CONST_C;
338 pxyz[2] = pfinal[2] * mass * CONST_C;
339 }
340
341 gamma = physlib_gamma_pnorm(mass, math_norm(pxyz));
342 fposxyz[0] = posxyz[0] + h[i] * pxyz[0] / (2 * gamma * mass);
343 fposxyz[1] = posxyz[1] + h[i] * pxyz[1] / (2 * gamma * mass);
344 fposxyz[2] = posxyz[2] + h[i] * pxyz[2] / (2 * gamma * mass);
345
346 if(!errflag) {
347 /* Back to cylindrical coordinates */
348 p->r[i] = sqrt(fposxyz[0]*fposxyz[0]+fposxyz[1]*fposxyz[1]);
349
350 /* phi is evaluated like this to make sure it is cumulative */
351 p->phi[i] += atan2(
352 posxyz0[0] * fposxyz[1] - posxyz0[1] * fposxyz[0],
353 posxyz0[0] * fposxyz[0] + posxyz0[1] * fposxyz[1] );
354 p->z[i] = fposxyz[2];
355
356 real cosp = cos(p->phi[i]);
357 real sinp = sin(p->phi[i]);
358 p->p_r[i] = pxyz[0] * cosp + pxyz[1] * sinp;
359 p->p_phi[i] = -pxyz[0] * sinp + pxyz[1] * cosp;
360 p->p_z[i] = pxyz[2];
361 }
362
363 /* Evaluate magnetic field (and gradient) and rho at new position */
364 real BdBrpz[15];
365 real psi[1];
366 real rho[2];
367 if(!errflag) {
368 errflag = B_field_eval_B_dB(BdBrpz, p->r[i], p->phi[i], p->z[i],
369 t0 + h[i], Bdata);
370 }
371 if(!errflag) {
372 errflag = B_field_eval_psi(psi, p->r[i], p->phi[i], p->z[i],
373 t0 + h[i], Bdata);
374 }
375 if(!errflag) {
376 errflag = B_field_eval_rho(rho, psi[0], Bdata);
377 }
378
379 if(!errflag) {
380 p->B_r[i] = BdBrpz[0];
381 p->B_r_dr[i] = BdBrpz[1];
382 p->B_r_dphi[i] = BdBrpz[2];
383 p->B_r_dz[i] = BdBrpz[3];
384
385 p->B_phi[i] = BdBrpz[4];
386 p->B_phi_dr[i] = BdBrpz[5];
387 p->B_phi_dphi[i] = BdBrpz[6];
388 p->B_phi_dz[i] = BdBrpz[7];
389
390 p->B_z[i] = BdBrpz[8];
391 p->B_z_dr[i] = BdBrpz[9];
392 p->B_z_dphi[i] = BdBrpz[10];
393 p->B_z_dz[i] = BdBrpz[11];
394
395 p->rho[i] = rho[0];
396
397 /* Evaluate phi and theta angles so that they are cumulative */
398 real axisrz[2];
399 errflag = B_field_get_axis_rz(axisrz, Bdata, p->phi[i]);
400 p->theta[i] += atan2( (R0-axisrz[0]) * (p->z[i]-axisrz[1])
401 - (z0-axisrz[1]) * (p->r[i]-axisrz[0]),
402 (R0-axisrz[0]) * (p->r[i]-axisrz[0])
403 + (z0-axisrz[1]) * (p->z[i]-axisrz[1]) );
404 }
405
406 /* Error handling */
407 if(errflag) {
408 p->err[i] = errflag;
409 p->running[i] = 0;
410 }
411 }
412 }
413}
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
CONST_C2
#define CONST_C2
Speed of light squared [m^2/s^2].
Definition consts.h:26
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:31
math_xyz2rpz
#define math_xyz2rpz(xyz, rpz)
Convert cartesian coordinates xyz to cylindrical coordinates rpz.
Definition math.h:77
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:86
math_rpz2xyz
#define math_rpz2xyz(rpz, xyz)
Convert cylindrical coordinates rpz to cartesian coordinates xyz.
Definition math.h:81
math_normc
#define math_normc(a1, a2, a3)
Calculate norm of 3D vector from its components a1, a2, a3.
Definition math.h:70
math_norm
#define math_norm(a)
Calculate norm of 3D vector a.
Definition math.h:67
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
void step_fo_vpa(particle_simd_fo *p, real *h, B_field_data *Bdata, E_field_data *Edata, int aldforce)
Integrate a full orbit step for a struct of particles with VPA.
Definition step_fo_vpa.c:36
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, int aldforce)
Integrate a full orbit step with VPA and MHd modes present.
Definition step_fo_vpa.c:234
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 Mon Apr 14 2025 13:15:19 for ASCOT5 by doxygen 1.13.2