ASCOT5
Loading...
Searching...
No Matches
step_ml_cashkarp.c
Go to the documentation of this file.
1
5#include <stdlib.h>
6#include <float.h>
7#include <math.h>
8#include "../../ascot5.h"
9#include "../../B_field.h"
10#include "../../math.h"
11#include "../../particle.h"
12#include "../../boozer.h"
13#include "../../mhd.h"
14#include "step_ml_cashkarp.h"
15
32void step_ml_cashkarp(particle_simd_ml* p, real* h, real* hnext, real tol,
33 B_field_data* Bdata) {
34
35 int i;
36 /* Following loop will be executed simultaneously for all i */
37 #pragma omp simd
38 for(i = 0; i < NSIMD; i++) {
39 if(p->running[i]) {
40 a5err errflag = 0;
41
42 real k1[3], k2[3], k3[3], k4[3], k5[3], k6[3];
43 real tempy[3];
44 real yprev[3];
45
46 real normB;
47
48 real R0 = p->r[i];
49 real z0 = p->z[i];
50 real t0 = p->time[i];
51
52 /* Direction */
53 int direction = 1;
54 if(p->pitch[i] < 0) {
55 direction = -1;
56 }
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
64 /* Magnetic field at initial position already known */
65 k1[0] = p->B_r[i];
66 k1[1] = p->B_phi[i];
67 k1[2] = p->B_z[i];
68
69 normB = (math_normc(k1[0], k1[1], k1[2])) * direction;
70 k1[0] /= normB;
71 k1[1] /= normB;
72 k1[2] /= normB;
73 k1[1] /= yprev[0];
74
75 for(int j = 0; j < 3; j++) {
76 tempy[j] = yprev[j]
77 + h[i] * ( (1.0/5) * k1[j] );
78 }
79 if(!errflag) {
80 errflag = B_field_eval_B(k2, tempy[0], tempy[1], tempy[2],
81 t0 + (1.0/5)*h[i], Bdata);
82 }
83 normB = (math_normc(k2[0], k2[1], k2[2])) * direction;
84 k2[0] /= normB;
85 k2[1] /= normB;
86 k2[2] /= normB;
87 k2[1] /= tempy[0];
88
89 for(int j = 0; j < 3; j++) {
90 tempy[j] = yprev[j]
91 + h[i] * (
92 (3.0/40) * k1[j]
93 + (9.0/40) * k2[j] );
94 }
95 if(!errflag) {
96 errflag = B_field_eval_B(k3, tempy[0], tempy[1], tempy[2],
97 t0 + (3.0/10)*h[i], Bdata);
98 }
99 normB = (math_normc(k3[0], k3[1], k3[2])) * direction;
100 k3[0] /= normB;
101 k3[1] /= normB;
102 k3[2] /= normB;
103 k3[1] /= tempy[0];
104
105 for(int j = 0; j < 3; j++) {
106 tempy[j] = yprev[j]
107 + h[i] * (
108 ( 3.0/10) * k1[j]
109 + (-9.0/10) * k2[j]
110 + ( 6.0/5 ) * k3[j] );
111 }
112 if(!errflag) {
113 errflag = B_field_eval_B(k4, tempy[0], tempy[1], tempy[2],
114 t0 + (3.0/5)*h[i], Bdata);
115 }
116 normB = (math_normc(k4[0], k4[1], k4[2])) * direction;
117 k4[0] /= normB;
118 k4[1] /= normB;
119 k4[2] /= normB;
120 k4[1] /= tempy[0];
121
122 for(int j = 0; j < 3; j++) {
123 tempy[j] = yprev[j]
124 + h[i] * (
125 (-11.0/54) * k1[j]
126 + ( 5.0/2 ) * k2[j]
127 + (-70.0/27) * k3[j]
128 + ( 35.0/27) * k4[j] );
129 }
130 if(!errflag) {
131 errflag = B_field_eval_B(k5, tempy[0], tempy[1], tempy[2],
132 t0 + h[i], Bdata);
133 }
134 normB = (math_normc(k5[0], k5[1], k5[2])) * direction;
135 k5[0] /= normB;
136 k5[1] /= normB;
137 k5[2] /= normB;
138 k5[1] /= tempy[0];
139
140 for(int j = 0; j < 3; j++) {
141 tempy[j] = yprev[j]
142 + h[i] * (
143 ( 1631.0/55296 ) * k1[j]
144 + ( 175.0/512 ) * k2[j]
145 + ( 575.0/13824 ) * k3[j]
146 + (44275.0/110592) * k4[j]
147 + ( 253.0/4096 ) * k5[j] );
148 }
149 if(!errflag) {
150 errflag = B_field_eval_B(k6, tempy[0], tempy[1], tempy[2],
151 t0 + (7.0/8)*h[i], Bdata);
152 }
153 normB = (math_normc(k6[0], k6[1], k6[2])) * direction;
154 k6[0] /= normB;
155 k6[1] /= normB;
156 k6[2] /= normB;
157 k6[1] /= tempy[0];
158
159 /* Error estimate is a difference between RK4 and RK5 solutions. If
160 * time-step is accepted, the RK5 solution will be used to advance
161 * marker. */
162 real rk5[3];
163 real err = 0.0;
164 for(int j = 0; j < 3; j++) {
165 rk5[j] = yprev[j]
166 + h[i] * (
167 ( 37.0/378 ) * k1[j]
168 + (250.0/621 ) * k3[j]
169 + (125.0/594 ) * k4[j]
170 + (512.0/1771) * k6[j] );
171
172 real rk4 = yprev[j] +
173 h[i] * (
174 ( 2825.0/27648) * k1[j]
175 + (18575.0/48384) * k3[j]
176 + (13525.0/55296) * k4[j]
177 + ( 277.0/14336) * k5[j]
178 + ( 1.0/4 ) * k6[j] );
179
180 real yerr = fabs(rk5[j] - rk4);
181 real ytol = fabs(yprev[j]) + fabs( k1[j] * h[i] ) + DBL_EPSILON;
182 err = fmax( err, yerr/ytol );
183 }
184
185 err = err/tol;
186 if(err <= 1){
187 /* Time step accepted */
188 hnext[i] = 0.85*h[i]*pow(err,-0.2);
189
190 }
191 else{
192 /* Time step rejected */
193 hnext[i] = -0.85*h[i]*pow(err,-0.25);
194
195 }
196
197 p->r[i] = rk5[0];
198 p->phi[i] = rk5[1];
199 p->z[i] = rk5[2];
200
201 /* Evaluate magnetic field (and gradient) and rho at new position */
202 real B_dB[15];
203 if(!errflag) {
204 errflag = B_field_eval_B_dB(B_dB, p->r[i], p->phi[i], p->z[i],
205 p->time[i] + h[i], Bdata);
206 }
207 p->B_r[i] = B_dB[0];
208 p->B_r_dr[i] = B_dB[1];
209 p->B_r_dphi[i] = B_dB[2];
210 p->B_r_dz[i] = B_dB[3];
211
212 p->B_phi[i] = B_dB[4];
213 p->B_phi_dr[i] = B_dB[5];
214 p->B_phi_dphi[i] = B_dB[6];
215 p->B_phi_dz[i] = B_dB[7];
216
217 p->B_z[i] = B_dB[8];
218 p->B_z_dr[i] = B_dB[9];
219 p->B_z_dphi[i] = B_dB[10];
220 p->B_z_dz[i] = B_dB[11];
221
222
223 real psi[1];
224 real rho[2];
225 if(!errflag) {
226 errflag = B_field_eval_psi(psi, p->r[i], p->phi[i], p->z[i],
227 p->time[i] + h[i], Bdata);
228 }
229 if(!errflag) {
230 errflag = B_field_eval_rho(rho, psi[0], Bdata);
231 }
232 p->rho[i] = rho[0];
233
234
235 /* Evaluate theta angle so that it is cumulative */
236 real axisrz[2];
237 errflag = B_field_get_axis_rz(axisrz, Bdata, p->phi[i]);
238 p->theta[i] += atan2( (R0-axisrz[0]) * (p->z[i]-axisrz[1])
239 - (z0-axisrz[1]) * (p->r[i]-axisrz[0]),
240 (R0-axisrz[0]) * (p->r[i]-axisrz[0])
241 + (z0-axisrz[1]) * (p->z[i]-axisrz[1]) );
242
243 }
244 }
245}
246
265void step_ml_cashkarp_mhd(particle_simd_ml* p, real* h, real* hnext, real tol,
266 B_field_data* Bdata, boozer_data* boozerdata,
267 mhd_data* mhddata) {
268
269 int i;
270 /* Following loop will be executed simultaneously for all i */
271 #pragma omp simd
272 for(i = 0; i < NSIMD; i++) {
273 if(p->running[i]) {
274 a5err errflag = 0;
275
276 real k1[3], k2[3], k3[3], k4[3], k5[3], k6[3];
277 real tempy[3];
278 real yprev[3];
279
280 real normB;
281 real pert_field[7];
282
283 real R0 = p->r[i];
284 real z0 = p->z[i];
285 real t0 = p->time[i];
286
287 /* Direction */
288 int direction = 1;
289 if(p->pitch[i] < 0) {
290 direction = -1;
291 }
292
293 int pertonly = 0;
294
295 /* Coordinates are copied from the struct into an array to make
296 * passing parameters easier */
297 yprev[0] = p->r[i];
298 yprev[1] = p->phi[i];
299 yprev[2] = p->z[i];
300
301 if(!errflag) {
302 errflag = mhd_perturbations(pert_field, p->r[i], p->phi[i],
303 p->z[i], p->time[i], pertonly,
304 MHD_INCLUDE_ALL, boozerdata,
305 mhddata, Bdata);
306 }
307 k1[0] = pert_field[0];
308 k1[1] = pert_field[1];
309 k1[2] = pert_field[2];
310
311 normB = (math_normc(k1[0], k1[1], k1[2])) * direction;
312 k1[0] /= normB;
313 k1[1] /= normB;
314 k1[2] /= normB;
315 k1[1] /= yprev[0];
316
317 for(int j = 0; j < 3; j++) {
318 tempy[j] = yprev[j]
319 + h[i] * ( (1.0/5) * k1[j] );
320 }
321 if(!errflag) {
322 errflag = mhd_perturbations(pert_field, tempy[0], tempy[1],
323 tempy[2], t0, pertonly,
324 MHD_INCLUDE_ALL, boozerdata,
325 mhddata, Bdata);
326 }
327 k2[0] = pert_field[0];
328 k2[1] = pert_field[1];
329 k2[2] = pert_field[2];
330
331 normB = (math_normc(k2[0], k2[1], k2[2])) * direction;
332 k2[0] /= normB;
333 k2[1] /= normB;
334 k2[2] /= normB;
335 k2[1] /= tempy[0];
336
337 for(int j = 0; j < 3; j++) {
338 tempy[j] = yprev[j]
339 + h[i] * (
340 (3.0/40) * k1[j]
341 + (9.0/40) * k2[j] );
342 }
343 if(!errflag) {
344 errflag = mhd_perturbations(pert_field, tempy[0], tempy[1],
345 tempy[2], t0, pertonly,
346 MHD_INCLUDE_ALL, boozerdata,
347 mhddata, Bdata);
348 }
349 k3[0] = pert_field[0];
350 k3[1] = pert_field[1];
351 k3[2] = pert_field[2];
352
353 normB = (math_normc(k3[0], k3[1], k3[2])) * direction;
354 k3[0] /= normB;
355 k3[1] /= normB;
356 k3[2] /= normB;
357 k3[1] /= tempy[0];
358
359 for(int j = 0; j < 3; j++) {
360 tempy[j] = yprev[j]
361 + h[i] * (
362 ( 3.0/10) * k1[j]
363 + (-9.0/10) * k2[j]
364 + ( 6.0/5 ) * k3[j] );
365 }
366 if(!errflag) {
367 errflag = mhd_perturbations(pert_field, tempy[0], tempy[1],
368 tempy[2], t0, pertonly,
369 MHD_INCLUDE_ALL, boozerdata,
370 mhddata, Bdata);
371 }
372 k4[0] = pert_field[0];
373 k4[1] = pert_field[1];
374 k4[2] = pert_field[2];
375
376 normB = (math_normc(k4[0], k4[1], k4[2])) * direction;
377 k4[0] /= normB;
378 k4[1] /= normB;
379 k4[2] /= normB;
380 k4[1] /= tempy[0];
381
382 for(int j = 0; j < 3; j++) {
383 tempy[j] = yprev[j]
384 + h[i] * (
385 (-11.0/54) * k1[j]
386 + ( 5.0/2 ) * k2[j]
387 + (-70.0/27) * k3[j]
388 + ( 35.0/27) * k4[j] );
389 }
390 if(!errflag) {
391 errflag = mhd_perturbations(pert_field, tempy[0], tempy[1],
392 tempy[2], t0, pertonly,
393 MHD_INCLUDE_ALL, boozerdata,
394 mhddata, Bdata);
395 }
396 k5[0] = pert_field[0];
397 k5[1] = pert_field[1];
398 k5[2] = pert_field[2];
399
400 normB = (math_normc(k5[0], k5[1], k5[2])) * direction;
401 k5[0] /= normB;
402 k5[1] /= normB;
403 k5[2] /= normB;
404 k5[1] /= tempy[0];
405
406 for(int j = 0; j < 3; j++) {
407 tempy[j] = yprev[j]
408 + h[i] * (
409 ( 1631.0/55296 ) * k1[j]
410 + ( 175.0/512 ) * k2[j]
411 + ( 575.0/13824 ) * k3[j]
412 + (44275.0/110592) * k4[j]
413 + ( 253.0/4096 ) * k5[j] );
414 }
415 if(!errflag) {
416 errflag = mhd_perturbations(pert_field, tempy[0], tempy[1],
417 tempy[2], t0, pertonly,
418 MHD_INCLUDE_ALL, boozerdata,
419 mhddata, Bdata);
420 }
421 k6[0] = pert_field[0];
422 k6[1] = pert_field[1];
423 k6[2] = pert_field[2];
424
425 normB = (math_normc(k6[0], k6[1], k6[2])) * direction;
426 k6[0] /= normB;
427 k6[1] /= normB;
428 k6[2] /= normB;
429 k6[1] /= tempy[0];
430
431 /* Error estimate is a difference between RK4 and RK5 solutions. If
432 * time-step is accepted, the RK5 solution will be used to advance
433 * marker. */
434 real rk5[3];
435 real err = 0.0;
436 for(int j = 0; j < 3; j++) {
437 rk5[j] = yprev[j]
438 + h[i] * (
439 ( 37.0/378 ) * k1[j]
440 + (250.0/621 ) * k3[j]
441 + (125.0/594 ) * k4[j]
442 + (512.0/1771) * k6[j] );
443
444 real rk4 = yprev[j] +
445 h[i] * (
446 ( 2825.0/27648) * k1[j]
447 + (18575.0/48384) * k3[j]
448 + (13525.0/55296) * k4[j]
449 + ( 277.0/14336) * k5[j]
450 + ( 1.0/4 ) * k6[j] );
451
452 real yerr = fabs(rk5[j] - rk4);
453 real ytol = fabs(yprev[j]) + fabs( k1[j] * h[i] ) + DBL_EPSILON;
454 err = fmax( err, yerr/ytol );
455 }
456
457 err = err/tol;
458 if(err <= 1){
459 /* Time step accepted */
460 hnext[i] = 0.85*h[i]*pow(err,-0.2);
461
462 }
463 else{
464 /* Time step rejected */
465 hnext[i] = -0.85*h[i]*pow(err,-0.25);
466
467 }
468
469 p->r[i] = rk5[0];
470 p->phi[i] = rk5[1];
471 p->z[i] = rk5[2];
472
473 /* Evaluate magnetic field (and gradient) and rho at new position */
474 real B_dB[15];
475 if(!errflag) {
476 errflag = B_field_eval_B_dB(B_dB, p->r[i], p->phi[i], p->z[i],
477 p->time[i] + h[i], Bdata);
478 }
479 p->B_r[i] = B_dB[0];
480 p->B_r_dr[i] = B_dB[1];
481 p->B_r_dphi[i] = B_dB[2];
482 p->B_r_dz[i] = B_dB[3];
483
484 p->B_phi[i] = B_dB[4];
485 p->B_phi_dr[i] = B_dB[5];
486 p->B_phi_dphi[i] = B_dB[6];
487 p->B_phi_dz[i] = B_dB[7];
488
489 p->B_z[i] = B_dB[8];
490 p->B_z_dr[i] = B_dB[9];
491 p->B_z_dphi[i] = B_dB[10];
492 p->B_z_dz[i] = B_dB[11];
493
494
495 real psi[1];
496 real rho[2];
497 if(!errflag) {
498 errflag = B_field_eval_psi(psi, p->r[i], p->phi[i], p->z[i],
499 p->time[i] + h[i], Bdata);
500 }
501 if(!errflag) {
502 errflag = B_field_eval_rho(rho, psi[0], Bdata);
503 }
504 p->rho[i] = rho[0];
505
506
507 /* Evaluate theta angle so that it is cumulative */
508 real axisrz[2];
509 errflag = B_field_get_axis_rz(axisrz, Bdata, p->phi[i]);
510 p->theta[i] += atan2( (R0-axisrz[0]) * (p->z[i]-axisrz[1])
511 - (z0-axisrz[1]) * (p->r[i]-axisrz[0]),
512 (R0-axisrz[0]) * (p->r[i]-axisrz[0])
513 + (z0-axisrz[1]) * (p->z[i]-axisrz[1]) );
514
515 }
516 }
517}
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:327
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:201
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:547
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:599
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:472
B_field.h
Header file for B_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.
a5err
unsigned long int a5err
Simulation error flag.
Definition error.h:17
math.h
Header file for math.c.
math_normc
#define math_normc(a1, a2, a3)
Calculate norm of 3D vector from its components a1, a2, a3.
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:234
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_ml_cashkarp_mhd
void step_ml_cashkarp_mhd(particle_simd_ml *p, real *h, real *hnext, real tol, B_field_data *Bdata, boozer_data *boozerdata, mhd_data *mhddata)
Integrate a magnetic field line step for a struct of markers.
Definition step_ml_cashkarp.c:265
step_ml_cashkarp
void step_ml_cashkarp(particle_simd_ml *p, real *h, real *hnext, real tol, B_field_data *Bdata)
Integrate a magnetic field line step for a struct of markers.
Definition step_ml_cashkarp.c:32
step_ml_cashkarp.h
Header file for step_ml_cashkarp.c.
B_field_data
Magnetic field simulation data.
Definition B_field.h:63
boozer_data
Boozer parameters on the target.
Definition boozer.h:29
mhd_data
MHD simulation data.
Definition mhd.h:57
particle_simd_ml
Struct representing NSIMD field line markers.
Definition particle.h:342
Generated on Wed Sep 25 2024 07:10:18 for ASCOT5 by doxygen 1.10.0