ASCOT5
Loading...
Searching...
No Matches
diag_orb.c
Go to the documentation of this file.
1
5#include "diag_orb.h"
6
7#include <math.h>
8#include <stdlib.h>
9#include <string.h>
10#include "../ascot5.h"
11#include "../consts.h"
12#include "../simulate.h"
13
26void diag_orb_init(diag_orb_data* data) {
27
28 int step = data->Nmrk*data->Npnt;
29 data->id = (real*) calloc( step, sizeof(real) );
30 if(data->mode == DIAG_ORB_POINCARE) {
31 data->pncrid = (real*) calloc( step, sizeof(real) );
32 data->pncrdi = (real*) calloc( step, sizeof(real) );
33 }
34 switch(data->record_mode) {
35
36 case simulate_mode_fo:
37 data->mileage = (real*) calloc( step, sizeof(real) );
38 data->r = (real*) calloc( step, sizeof(real) );
39 data->phi = (real*) calloc( step, sizeof(real) );
40 data->z = (real*) calloc( step, sizeof(real) );
41 data->p_r = (real*) calloc( step, sizeof(real) );
42 data->p_phi = (real*) calloc( step, sizeof(real) );
43 data->p_z = (real*) calloc( step, sizeof(real) );
44 data->weight = (real*) calloc( step, sizeof(real) );
45 data->charge = (real*) calloc( step, sizeof(real) );
46 data->rho = (real*) calloc( step, sizeof(real) );
47 data->theta = (real*) calloc( step, sizeof(real) );
48 data->B_r = (real*) calloc( step, sizeof(real) );
49 data->B_phi = (real*) calloc( step, sizeof(real) );
50 data->B_z = (real*) calloc( step, sizeof(real) );
51 data->simmode = (real*) calloc( step, sizeof(real) );
52 break;
53
54 case simulate_mode_gc:
55 data->mileage = (real*) calloc( step, sizeof(real) );
56 data->r = (real*) calloc( step, sizeof(real) );
57 data->phi = (real*) calloc( step, sizeof(real) );
58 data->z = (real*) calloc( step, sizeof(real) );
59 data->ppar = (real*) calloc( step, sizeof(real) );
60 data->mu = (real*) calloc( step, sizeof(real) );
61 data->zeta = (real*) calloc( step, sizeof(real) );
62 data->weight = (real*) calloc( step, sizeof(real) );
63 data->charge = (real*) calloc( step, sizeof(real) );
64 data->rho = (real*) calloc( step, sizeof(real) );
65 data->theta = (real*) calloc( step, sizeof(real) );
66 data->B_r = (real*) calloc( step, sizeof(real) );
67 data->B_phi = (real*) calloc( step, sizeof(real) );
68 data->B_z = (real*) calloc( step, sizeof(real) );
69 data->simmode = (real*) calloc( step, sizeof(real) );
70 break;
71
72 case simulate_mode_ml:
73 data->mileage = (real*) calloc( step, sizeof(real) );
74 data->r = (real*) calloc( step, sizeof(real) );
75 data->phi = (real*) calloc( step, sizeof(real) );
76 data->z = (real*) calloc( step, sizeof(real) );
77 data->rho = (real*) calloc( step, sizeof(real) );
78 data->theta = (real*) calloc( step, sizeof(real) );
79 data->B_r = (real*) calloc( step, sizeof(real) );
80 data->B_phi = (real*) calloc( step, sizeof(real) );
81 data->B_z = (real*) calloc( step, sizeof(real) );
82 data->simmode = (real*) calloc( step, sizeof(real) );
83 break;
84
85 case simulate_mode_hybrid:
86 data->mileage = (real*) calloc( step, sizeof(real) );
87 data->r = (real*) calloc( step, sizeof(real) );
88 data->phi = (real*) calloc( step, sizeof(real) );
89 data->z = (real*) calloc( step, sizeof(real) );
90 data->p_r = (real*) calloc( step, sizeof(real) );
91 data->p_phi = (real*) calloc( step, sizeof(real) );
92 data->p_z = (real*) calloc( step, sizeof(real) );
93 data->ppar = (real*) calloc( step, sizeof(real) );
94 data->mu = (real*) calloc( step, sizeof(real) );
95 data->zeta = (real*) calloc( step, sizeof(real) );
96 data->weight = (real*) calloc( step, sizeof(real) );
97 data->charge = (real*) calloc( step, sizeof(real) );
98 data->rho = (real*) calloc( step, sizeof(real) );
99 data->theta = (real*) calloc( step, sizeof(real) );
100 data->B_r = (real*) calloc( step, sizeof(real) );
101 data->B_phi = (real*) calloc( step, sizeof(real) );
102 data->B_z = (real*) calloc( step, sizeof(real) );
103 data->simmode = (real*) calloc( step, sizeof(real) );
104 break;
105 }
106
107 data->mrk_pnt = (integer*) malloc( data->Nmrk*sizeof(integer) );
108 data->mrk_recorded = (real*) malloc( data->Nmrk*sizeof(real) );
109
110 memset(data->mrk_pnt, 0, data->Nmrk*sizeof(integer));
111 memset(data->mrk_recorded, 0, data->Nmrk*sizeof(real));
112}
113
119void diag_orb_free(diag_orb_data* data){
120 free(data->mrk_pnt);
121 free(data->mrk_recorded);
122
123 free(data->id);
124 if(data->mode == DIAG_ORB_POINCARE) {
125 free(data->pncrid);
126 free(data->pncrdi);
127 }
128 switch(data->record_mode) {
129
130 case simulate_mode_fo:
131 free(data->mileage);
132 free(data->r);
133 free(data->phi);
134 free(data->z);
135 free(data->p_r);
136 free(data->p_phi);
137 free(data->p_z);
138 free(data->weight);
139 free(data->charge );
140 free(data->rho);
141 free(data->theta);
142 free(data->B_r);
143 free(data->B_phi);
144 free(data->B_z);
145 free(data->simmode);
146 break;
147
148 case simulate_mode_gc:
149 free(data->mileage);
150 free(data->r);
151 free(data->phi);
152 free(data->z);
153 free(data->ppar);
154 free(data->mu);
155 free(data->zeta);
156 free(data->weight);
157 free(data->charge );
158 free(data->rho);
159 free(data->theta);
160 free(data->B_r);
161 free(data->B_phi);
162 free(data->B_z);
163 free(data->simmode);
164 break;
165
166 case simulate_mode_ml:
167 free(data->mileage);
168 free(data->r);
169 free(data->phi);
170 free(data->z);
171 free(data->rho);
172 free(data->theta);
173 free(data->B_r);
174 free(data->B_phi);
175 free(data->B_z);
176 free(data->simmode);
177 break;
178
179 case simulate_mode_hybrid:
180 free(data->mileage);
181 free(data->r);
182 free(data->phi);
183 free(data->z);
184 free(data->p_r);
185 free(data->p_phi);
186 free(data->p_z);
187 free(data->ppar);
188 free(data->mu);
189 free(data->zeta);
190 free(data->weight);
191 free(data->charge );
192 free(data->rho);
193 free(data->theta);
194 free(data->B_r);
195 free(data->B_phi);
196 free(data->B_z);
197 free(data->simmode);
198 break;
199 }
200}
201
211void diag_orb_update_fo(diag_orb_data* data, particle_simd_fo* p_f,
212 particle_simd_fo* p_i) {
213
214 if(data->mode == DIAG_ORB_INTERVAL) {
215
216 #pragma omp simd
217 for(int i= 0; i < NSIMD; i++) {
218
219 /* Mask dummy markers */
220 if(p_f->id[i] > 0) {
221
222 integer imrk = p_f->index[i];
223 integer ipoint = data->mrk_pnt[imrk];
224 integer idx = imrk * data->Npnt + ipoint;
225
226 /* If this is the first time-step, record marker position. */
227 if( data->id[imrk * data->Npnt] == 0 ) {
228 data->id[idx] = (real)p_i->id[i];
229 data->mileage[idx]= p_i->mileage[i];
230 data->r[idx] = p_i->r[i];
231 data->phi[idx] = p_i->phi[i];
232 data->z[idx] = p_i->z[i];
233 data->p_r[idx] = p_i->p_r[i];
234 data->p_phi[idx] = p_i->p_phi[i];
235 data->p_z[idx] = p_i->p_z[i];
236 data->weight[idx] = p_i->weight[i];
237 data->charge[idx] = p_i->charge[i];
238 data->rho[idx] = p_i->rho[i];
239 data->theta[idx] = p_i->theta[i];
240 data->B_r[idx] = p_i->B_r[i];
241 data->B_phi[idx] = p_i->B_phi[i];
242 data->B_z[idx] = p_i->B_z[i];
243 data->simmode[idx]= DIAG_ORB_FO;
244
245 ipoint++;
246 if(ipoint == data->Npnt) {
247 ipoint = 0;
248 }
249 data->mrk_pnt[imrk] = ipoint;
250 data->mrk_recorded[imrk] = p_i->mileage[i];
251 }
252
253 /* Record marker if enough time has passed from last record, or
254 if marker has met some end condition. */
255 real dt = data->mrk_recorded[imrk] + data->writeInterval
256 - p_f->mileage[i];
257 if( dt <= 0 || p_f->endcond[i] > 0 ) {
258 idx = imrk * data->Npnt + ipoint;
259
260 data->id[idx] = (real)p_f->id[i];
261 data->mileage[idx]= p_f->mileage[i];
262 data->r[idx] = p_f->r[i];
263 data->phi[idx] = p_f->phi[i];
264 data->z[idx] = p_f->z[i];
265 data->p_r[idx] = p_f->p_r[i];
266 data->p_phi[idx] = p_f->p_phi[i];
267 data->p_z[idx] = p_f->p_z[i];
268 data->weight[idx] = p_f->weight[i];
269 data->charge[idx] = p_f->charge[i];
270 data->rho[idx] = p_f->rho[i];
271 data->theta[idx] = p_f->theta[i];
272 data->B_r[idx] = p_f->B_r[i];
273 data->B_phi[idx] = p_f->B_phi[i];
274 data->B_z[idx] = p_f->B_z[i];
275 data->simmode[idx]= DIAG_ORB_FO;
276
277 ipoint++;
278 if(ipoint == data->Npnt) {
279 ipoint = 0;
280 }
281 data->mrk_pnt[imrk] = ipoint;
282 data->mrk_recorded[imrk] = p_f->mileage[i];
283 }
284 }
285 }
286 }
287 else if(data->mode == DIAG_ORB_POINCARE) {
288
289 #pragma omp simd
290 for(int i= 0; i < NSIMD; i++) {
291 /* Mask dummy markers and those whose time-step was rejected. */
292 if( p_f->id[i] > 0 && (p_f->mileage[i] != p_i->mileage[i]) ) {
293
294 real k;
295 integer imrk = p_f->index[i];
296 integer ipoint = data->mrk_pnt[imrk];
297 integer idx = imrk * data->Npnt + ipoint;
298
299 /* Check and store toroidal crossings. */
300 for(int j=0; j < data->ntoroidalplots; j++) {
301 k = diag_orb_check_plane_crossing(p_f->phi[i], p_i->phi[i],
302 data->toroidalangles[j]);
303 if(k) {
304 real d = 1-k;
305 idx = imrk * data->Npnt + ipoint;
306 data->id[idx] = (real)p_f->id[i];
307 data->mileage[idx]= k*p_f->mileage[i]+ d*p_i->mileage[i];
308 data->r[idx] = k*p_f->r[i] + d*p_i->r[i];
309 data->phi[idx] = k*p_f->phi[i] + d*p_i->phi[i];
310 data->z[idx] = k*p_f->z[i] + d*p_i->z[i];
311 data->p_r[idx] = k*p_f->p_r[i] + d*p_i->p_r[i];
312 data->p_phi[idx] = k*p_f->p_phi[i] + d*p_i->p_phi[i];
313 data->p_z[idx] = k*p_f->p_z[i] + d*p_i->p_z[i];
314 data->weight[idx] = k*p_f->weight[i] + d*p_i->weight[i];
315 data->charge[idx] = p_i->charge[i];
316 data->rho[idx] = k*p_f->rho[i] + d*p_i->rho[i];
317 data->theta[idx] = k*p_f->theta[i] + d*p_i->theta[i];
318 data->B_r[idx] = k*p_f->B_r[i] + d*p_i->B_r[i];
319 data->B_phi[idx] = k*p_f->B_phi[i] + d*p_i->B_phi[i];
320 data->B_z[idx] = k*p_f->B_z[i] + d*p_i->B_z[i];
321 data->pncrid[idx] = j;
322 data->pncrdi[idx] = 1 - 2 * (p_f->phi[i] < p_i->phi[i]);
323 data->simmode[idx]= DIAG_ORB_FO;
324
325 ipoint++;
326 if(ipoint == data->Npnt) {
327 ipoint = 0;
328 }
329 data->mrk_pnt[imrk] = ipoint;
330 data->mrk_recorded[imrk] = p_f->mileage[i];
331 }
332 }
333
334 /* Check and store poloidal crossings. */
335 for(int j=0; j < data->npoloidalplots; j++) {
336 k = diag_orb_check_plane_crossing(p_f->theta[i],
337 p_i->theta[i],
338 data->poloidalangles[j]);
339 if(k) {
340 real d = 1-k;
341 idx = imrk * data->Npnt + ipoint;
342 data->id[idx] = (real)p_f->id[i];
343 data->mileage[idx]= k*p_f->mileage[i]+ d*p_i->mileage[i];
344 data->r[idx] = k*p_f->r[i] + d*p_i->r[i];
345 data->phi[idx] = k*p_f->phi[i] + d*p_i->phi[i];
346 data->z[idx] = k*p_f->z[i] + d*p_i->z[i];
347 data->p_r[idx] = k*p_f->p_r[i] + d*p_i->p_r[i];
348 data->p_phi[idx] = k*p_f->p_phi[i] + d*p_i->p_phi[i];
349 data->p_z[idx] = k*p_f->p_z[i] + d*p_i->p_z[i];
350 data->weight[idx] = k*p_f->weight[i] + d*p_i->weight[i];
351 data->charge[idx] = p_i->charge[i];
352 data->rho[idx] = k*p_f->rho[i] + d*p_i->rho[i];
353 data->theta[idx] = k*p_f->theta[i] + d*p_i->theta[i];
354 data->B_r[idx] = k*p_f->B_r[i] + d*p_i->B_r[i];
355 data->B_phi[idx] = k*p_f->B_phi[i] + d*p_i->B_phi[i];
356 data->B_z[idx] = k*p_f->B_z[i] + d*p_i->B_z[i];
357 data->pncrid[idx] = j + data->ntoroidalplots;
358 data->pncrdi[idx] = 1 - 2 * (p_f->theta[i] < p_i->theta[i]);
359 data->simmode[idx]= DIAG_ORB_FO;
360
361 ipoint++;
362 if(ipoint == data->Npnt) {
363 ipoint = 0;
364 }
365 data->mrk_pnt[imrk] = ipoint;
366 data->mrk_recorded[imrk] = p_f->mileage[i];
367 }
368 }
369
370 /* Check and store radial crossings. */
371 for(int j=0; j < data->nradialplots; j++) {
372 k = diag_orb_check_radial_crossing(p_f->rho[i],p_i->rho[i],
373 data->radialdistances[j]);
374 if(k) {
375 real d = k;
376 k = 1-d;
377 idx = imrk * data->Npnt + ipoint;
378 data->id[idx] = (real)p_f->id[i];
379 data->mileage[idx]= k*p_f->mileage[i]+ d*p_i->mileage[i];
380 data->r[idx] = k*p_f->r[i] + d*p_i->r[i];
381 data->phi[idx] = k*p_f->phi[i] + d*p_i->phi[i];
382 data->z[idx] = k*p_f->z[i] + d*p_i->z[i];
383 data->p_r[idx] = k*p_f->p_r[i] + d*p_i->p_r[i];
384 data->p_phi[idx] = k*p_f->p_phi[i] + d*p_i->p_phi[i];
385 data->p_z[idx] = k*p_f->p_z[i] + d*p_i->p_z[i];
386 data->weight[idx] = k*p_f->weight[i] + d*p_i->weight[i];
387 data->charge[idx] = p_i->charge[i];
388 data->rho[idx] = k*p_f->rho[i] + d*p_i->rho[i];
389 data->theta[idx] = k*p_f->theta[i] + d*p_i->theta[i];
390 data->B_r[idx] = k*p_f->B_r[i] + d*p_i->B_r[i];
391 data->B_phi[idx] = k*p_f->B_phi[i] + d*p_i->B_phi[i];
392 data->B_z[idx] = k*p_f->B_z[i] + d*p_i->B_z[i];
393 data->pncrid[idx] = j + data->ntoroidalplots + data->npoloidalplots;
394 data->pncrdi[idx] = 1 - 2 * (p_f->rho[i] < p_i->rho[i]);
395 ipoint++;
396 if(ipoint == data->Npnt) {
397 ipoint = 0;
398 }
399 data->mrk_pnt[imrk] = ipoint;
400 data->mrk_recorded[imrk] = p_f->mileage[i];
401 }
402 }
403 }
404 }
405 }
406}
407
417void diag_orb_update_gc(diag_orb_data* data, particle_simd_gc* p_f,
418 particle_simd_gc* p_i) {
419
420 if(data->mode == DIAG_ORB_INTERVAL) {
421 #pragma omp simd
422 for(int i= 0; i < NSIMD; i++) {
423
424 /* Mask dummy markers */
425 if(p_f->id[i] > 0) {
426 integer imrk = p_f->index[i];
427 integer ipoint = data->mrk_pnt[imrk];
428 integer idx = imrk * data->Npnt + ipoint;
429
430 /* If this is the first time-step, record marker position. */
431 if( data->id[imrk * data->Npnt] == 0 ) {
432 data->id[idx] = (real)(p_i->id[i]);
433 data->mileage[idx]= p_i->mileage[i];
434 data->r[idx] = p_i->r[i];
435 data->phi[idx] = p_i->phi[i];
436 data->z[idx] = p_i->z[i];
437 data->ppar[idx] = p_i->ppar[i];
438 data->mu[idx] = p_i->mu[i];
439 data->zeta[idx] = p_i->zeta[i];
440 data->weight[idx] = p_i->weight[i];
441 data->charge[idx] = p_i->charge[i];
442 data->rho[idx] = p_i->rho[i];
443 data->theta[idx] = p_i->theta[i];
444 data->B_r[idx] = p_i->B_r[i];
445 data->B_phi[idx] = p_i->B_phi[i];
446 data->B_z[idx] = p_i->B_z[i];
447 data->simmode[idx]= DIAG_ORB_GC;
448
449 ipoint++;
450 if(ipoint == data->Npnt) {
451 ipoint = 0;
452 }
453 data->mrk_pnt[imrk] = ipoint;
454 data->mrk_recorded[imrk] = p_i->mileage[i];
455 }
456
457 /* Record marker if enough time has passed from last record, or
458 if marker has met some end condition. */
459 real dt = data->mrk_recorded[imrk] + data->writeInterval
460 - p_f->mileage[i];
461
462 if( dt <= 0 || p_f->endcond[i] > 0 ) {
463 idx = imrk * data->Npnt + ipoint;
464
465 data->id[idx] = (real)p_f->id[i];
466 data->mileage[idx]= p_f->mileage[i];
467 data->r[idx] = p_f->r[i];
468 data->phi[idx] = p_f->phi[i];
469 data->z[idx] = p_f->z[i];
470 data->ppar[idx] = p_f->ppar[i];
471 data->mu[idx] = p_f->mu[i];
472 data->zeta[idx] = p_f->zeta[i];
473 data->weight[idx] = p_f->weight[i];
474 data->charge[idx] = p_f->charge[i];
475 data->rho[idx] = p_f->rho[i];
476 data->theta[idx] = p_f->theta[i];
477 data->B_r[idx] = p_f->B_r[i];
478 data->B_phi[idx] = p_f->B_phi[i];
479 data->B_z[idx] = p_f->B_z[i];
480 data->simmode[idx]= DIAG_ORB_GC;
481
482 ipoint++;
483 if(ipoint == data->Npnt) {
484 ipoint = 0;
485 }
486 data->mrk_pnt[imrk] = ipoint;
487 data->mrk_recorded[imrk] = p_f->mileage[i];
488 }
489 }
490 }
491 }
492 else if(data->mode == DIAG_ORB_POINCARE) {
493 #pragma omp simd
494 for(int i= 0; i < NSIMD; i++) {
495 /* Mask dummy markers and those whose time-step was rejected. */
496 if( p_f->id[i] > 0 && (p_f->mileage[i] != p_i->mileage[i]) ) {
497
498 real k;
499 integer imrk = p_f->index[i];
500 integer ipoint = data->mrk_pnt[imrk];
501 integer idx = imrk * data->Npnt + ipoint;
502
503 /* Check and store toroidal crossings. */
504 for(int j=0; j < data->ntoroidalplots; j++) {
505 k = diag_orb_check_plane_crossing(p_f->phi[i], p_i->phi[i],
506 data->toroidalangles[j]);
507 if(k) {
508 real d = 1-k;
509 idx = imrk * data->Npnt + ipoint;
510 data->id[idx] = (real)p_f->id[i];
511 data->mileage[idx]= k*p_f->mileage[i]+ d*p_i->mileage[i];
512 data->r[idx] = k*p_f->r[i] + d*p_i->r[i];
513 data->phi[idx] = k*p_f->phi[i] + d*p_i->phi[i];
514 data->z[idx] = k*p_f->z[i] + d*p_i->z[i];
515 data->ppar[idx] = k*p_f->ppar[i] + d*p_i->ppar[i];
516 data->mu[idx] = k*p_f->mu[i] + d*p_i->mu[i];
517 data->zeta[idx] = k*p_f->zeta[i] + d*p_i->zeta[i];
518 data->weight[idx] = k*p_f->weight[i] + d*p_i->weight[i];
519 data->charge[idx] = p_i->charge[i];
520 data->rho[idx] = k*p_f->rho[i] + d*p_i->rho[i];
521 data->theta[idx] = k*p_f->theta[i] + d*p_i->theta[i];
522 data->B_r[idx] = k*p_f->B_r[i] + d*p_i->B_r[i];
523 data->B_phi[idx] = k*p_f->B_phi[i] + d*p_i->B_phi[i];
524 data->B_z[idx] = k*p_f->B_z[i] + d*p_i->B_z[i];
525 data->pncrid[idx] = j;
526 data->pncrdi[idx] = 1 - 2 * (p_f->phi[i] < p_i->phi[i]);
527 data->simmode[idx]= DIAG_ORB_GC;
528
529 ipoint++;
530 if(ipoint == data->Npnt) {
531 ipoint = 0;
532 }
533 data->mrk_pnt[imrk] = ipoint;
534 data->mrk_recorded[imrk] = p_f->mileage[i];
535 }
536 }
537
538 /* Check and store poloidal crossings. */
539 for(int j=0; j < data->npoloidalplots; j++) {
540 k = diag_orb_check_plane_crossing(p_f->theta[i],
541 p_i->theta[i],
542 data->poloidalangles[j]);
543 if(k) {
544 real d = 1-k;
545 idx = imrk * data->Npnt + ipoint;
546 data->id[idx] = (real)p_f->id[i];
547 data->mileage[idx]= k*p_f->mileage[i]+ d*p_i->mileage[i];
548 data->r[idx] = k*p_f->r[i] + d*p_i->r[i];
549 data->phi[idx] = k*p_f->phi[i] + d*p_i->phi[i];
550 data->z[idx] = k*p_f->z[i] + d*p_i->z[i];
551 data->ppar[idx] = k*p_f->ppar[i] + d*p_i->ppar[i];
552 data->mu[idx] = k*p_f->mu[i] + d*p_i->mu[i];
553 data->zeta[idx] = k*p_f->zeta[i] + d*p_i->zeta[i];
554 data->weight[idx] = k*p_f->weight[i] + d*p_i->weight[i];
555 data->charge[idx] = p_i->charge[i];
556 data->rho[idx] = k*p_f->rho[i] + d*p_i->rho[i];
557 data->theta[idx] = k*p_f->theta[i] + d*p_i->theta[i];
558 data->B_r[idx] = k*p_f->B_r[i] + d*p_i->B_r[i];
559 data->B_phi[idx] = k*p_f->B_phi[i] + d*p_i->B_phi[i];
560 data->B_z[idx] = k*p_f->B_z[i] + d*p_i->B_z[i];
561 data->pncrid[idx] = j + data->ntoroidalplots;
562 data->pncrdi[idx] = 1 - 2 * (p_f->theta[i] < p_i->theta[i]);
563 data->simmode[idx]= DIAG_ORB_GC;
564
565 ipoint++;
566 if(ipoint == data->Npnt) {
567 ipoint = 0;
568 }
569 data->mrk_pnt[imrk] = ipoint;
570 data->mrk_recorded[imrk] = p_f->mileage[i];
571 }
572 }
573
574
575 /* Check and store radial crossings. */
576 for(int j=0; j < data->nradialplots; j++) {
577 k = diag_orb_check_radial_crossing(p_f->rho[i],
578 p_i->rho[i],
579 data->radialdistances[j]);
580 if(k) {
581 real d = 1-k;
582 idx = imrk * data->Npnt + ipoint;
583 data->id[idx] = (real)p_f->id[i];
584 data->mileage[idx]= k*p_f->mileage[i]+ d*p_i->mileage[i];
585 data->r[idx] = k*p_f->r[i] + d*p_i->r[i];
586 data->phi[idx] = k*p_f->phi[i] + d*p_i->phi[i];
587 data->z[idx] = k*p_f->z[i] + d*p_i->z[i];
588 data->ppar[idx] = k*p_f->ppar[i] + d*p_i->ppar[i];
589 data->mu[idx] = k*p_f->mu[i] + d*p_i->mu[i];
590 data->zeta[idx] = k*p_f->zeta[i] + d*p_i->zeta[i];
591 data->weight[idx] = k*p_f->weight[i] + d*p_i->weight[i];
592 data->charge[idx] = p_i->charge[i];
593 data->rho[idx] = k*p_f->rho[i] + d*p_i->rho[i];
594 data->theta[idx] = k*p_f->theta[i] + d*p_i->theta[i];
595 data->B_r[idx] = k*p_f->B_r[i] + d*p_i->B_r[i];
596 data->B_phi[idx] = k*p_f->B_phi[i] + d*p_i->B_phi[i];
597 data->B_z[idx] = k*p_f->B_z[i] + d*p_i->B_z[i];
598 data->pncrid[idx] =
599 j + data->ntoroidalplots + data->npoloidalplots;
600 data->pncrdi[idx] = 1 - 2 * (p_f->rho[i] < p_i->rho[i]);
601
602 ipoint++;
603 if(ipoint == data->Npnt) {
604 ipoint = 0;
605 }
606 data->mrk_pnt[imrk] = ipoint;
607 data->mrk_recorded[imrk] = p_f->mileage[i];
608 }
609 }
610 }
611 }
612 }
613}
614
624void diag_orb_update_ml(diag_orb_data* data, particle_simd_ml* p_f,
625 particle_simd_ml* p_i) {
626
627 if(data->mode == DIAG_ORB_INTERVAL) {
628
629 #pragma omp simd
630 for(int i= 0; i < NSIMD; i++) {
631
632 /* Mask dummy markers */
633 if(p_f->id[i] > 0) {
634 integer imrk = p_f->index[i];
635 integer ipoint = data->mrk_pnt[imrk];
636 integer idx = imrk * data->Npnt + ipoint;
637
638 /* If this is the first time-step, record marker position. */
639 if( data->id[imrk * data->Npnt] == 0 ) {
640 data->id[idx] = (real)p_i->id[i];
641 data->mileage[idx] = p_i->mileage[i];
642 data->r[idx] = p_i->r[i];
643 data->phi[idx] = p_i->phi[i];
644 data->z[idx] = p_i->z[i];
645 data->rho[idx] = p_i->rho[i];
646 data->theta[idx] = p_i->theta[i];
647 data->B_r[idx] = p_i->B_r[i];
648 data->B_phi[idx] = p_i->B_phi[i];
649 data->B_z[idx] = p_i->B_z[i];
650 data->simmode[idx] = DIAG_ORB_ML;
651
652 ipoint++;
653 if(ipoint == data->Npnt) {
654 ipoint = 0;
655 }
656 data->mrk_pnt[imrk] = ipoint;
657 data->mrk_recorded[imrk] = p_i->mileage[i];
658 }
659
660 /* Record marker if enough time has passed from last record, or
661 if marker has met some end condition. */
662 real dt = data->mrk_recorded[imrk] + data->writeInterval
663 - p_f->mileage[i];
664 if( dt <= 0 || p_f->endcond[i] > 0 ) {
665 idx = imrk * data->Npnt + ipoint;
666 data->id[idx] = (real)p_f->id[i];
667 data->mileage[idx] = p_f->mileage[i];
668 data->r[idx] = p_f->r[i];
669 data->phi[idx] = p_f->phi[i];
670 data->z[idx] = p_f->z[i];
671 data->rho[idx] = p_f->rho[i];
672 data->theta[idx] = p_f->theta[i];
673 data->B_r[idx] = p_f->B_r[i];
674 data->B_phi[idx] = p_f->B_phi[i];
675 data->B_z[idx] = p_f->B_z[i];
676 data->simmode[idx] = DIAG_ORB_ML;
677
678 ipoint++;
679 if(ipoint == data->Npnt) {
680 ipoint = 0;
681 }
682 data->mrk_pnt[imrk] = ipoint;
683 data->mrk_recorded[imrk] = p_f->mileage[i];
684 }
685 }
686 }
687 }
688 else if(data->mode == DIAG_ORB_POINCARE) {
689 #pragma omp simd
690 for(int i= 0; i < NSIMD; i++) {
691 /* Mask dummy markers and thosw whose time-step was rejected. */
692 if( p_f->id[i] > 0 && (p_f->mileage[i] != p_i->mileage[i]) ) {
693
694 real k;
695 integer imrk = p_f->index[i];
696 integer ipoint = data->mrk_pnt[imrk];
697 integer idx = imrk * data->Npnt + ipoint;
698
699 /* Check and store toroidal crossings. */
700 for(int j=0; j < data->ntoroidalplots; j++) {
701 k = diag_orb_check_plane_crossing(p_f->phi[i], p_i->phi[i],
702 data->toroidalangles[j]);
703 if(k) {
704 real d = 1-k;
705 idx = imrk * data->Npnt + ipoint;
706 data->id[idx] = (real)p_f->id[i];
707 data->mileage[idx]= k*p_f->mileage[i]+ d*p_i->mileage[i];
708 data->r[idx] = k*p_f->r[i] + d*p_i->r[i];
709 data->phi[idx] = k*p_f->phi[i] + d*p_i->phi[i];
710 data->z[idx] = k*p_f->z[i] + d*p_i->z[i];
711 data->rho[idx] = k*p_f->rho[i] + d*p_i->rho[i];
712 data->theta[idx] = k*p_f->theta[i] + d*p_i->theta[i];
713 data->B_r[idx] = k*p_f->B_r[i] + d*p_i->B_r[i];
714 data->B_phi[idx] = k*p_f->B_phi[i] + d*p_i->B_phi[i];
715 data->B_z[idx] = k*p_f->B_z[i] + d*p_i->B_z[i];
716 data->pncrid[idx] = j;
717 data->pncrdi[idx] = 1 - 2 * (p_f->phi[i] < p_i->phi[i]);
718 data->simmode[idx]= DIAG_ORB_ML;
719
720 ipoint++;
721 if(ipoint == data->Npnt) {
722 ipoint = 0;
723 }
724 data->mrk_pnt[imrk] = ipoint;
725 data->mrk_recorded[imrk] = p_f->mileage[i];
726 }
727 }
728
729 /* Check and store poloidal crossings. */
730 for(int j=0; j < data->npoloidalplots; j++) {
731 k = diag_orb_check_plane_crossing(p_f->theta[i],
732 p_i->theta[i],
733 data->poloidalangles[j]);
734 if(k) {
735 real d = 1-k;
736 idx = imrk * data->Npnt + ipoint;
737 data->id[idx] = (real)p_f->id[i];
738 data->mileage[idx]= k*p_f->mileage[i] + d*p_i->mileage[i];
739 data->r[idx] = k*p_f->r[i] + d*p_i->r[i];
740 data->phi[idx] = k*p_f->phi[i] + d*p_i->phi[i];
741 data->z[idx] = k*p_f->z[i] + d*p_i->z[i];
742 data->rho[idx] = k*p_f->rho[i] + d*p_i->rho[i];
743 data->theta[idx] = k*p_f->theta[i] + d*p_i->theta[i];
744 data->B_r[idx] = k*p_f->B_r[i] + d*p_i->B_r[i];
745 data->B_phi[idx] = k*p_f->B_phi[i] + d*p_i->B_phi[i];
746 data->B_z[idx] = k*p_f->B_z[i] + d*p_i->B_z[i];
747 data->pncrid[idx] = j + data->ntoroidalplots;
748 data->pncrdi[idx] = 1 - 2 * (p_f->theta[i] < p_i->theta[i]);
749 data->simmode[idx]= DIAG_ORB_ML;
750
751 ipoint++;
752 if(ipoint == data->Npnt) {
753 ipoint = 0;
754 }
755 data->mrk_pnt[imrk] = ipoint;
756 data->mrk_recorded[imrk] = p_f->mileage[i];
757 }
758 }
759
760 /* Check and store radial crossings. */
761 for(int j=0; j < data->nradialplots; j++) {
762 k = diag_orb_check_radial_crossing(p_f->rho[i],
763 p_i->rho[i],
764 data->radialdistances[j]);
765 if(k) {
766 real d = 1-k;
767 idx = imrk * data->Npnt + ipoint;
768 data->id[idx] = (real)p_f->id[i];
769 data->mileage[idx]= k*p_f->mileage[i] + d*p_i->mileage[i];
770 data->r[idx] = k*p_f->r[i] + d*p_i->r[i];
771 data->phi[idx] = k*p_f->phi[i] + d*p_i->phi[i];
772 data->z[idx] = k*p_f->z[i] + d*p_i->z[i];
773 data->rho[idx] = k*p_f->rho[i] + d*p_i->rho[i];
774 data->theta[idx] = k*p_f->theta[i] + d*p_i->theta[i];
775 data->B_r[idx] = k*p_f->B_r[i] + d*p_i->B_r[i];
776 data->B_phi[idx] = k*p_f->B_phi[i] + d*p_i->B_phi[i];
777 data->B_z[idx] = k*p_f->B_z[i] + d*p_i->B_z[i];
778 data->pncrid[idx] =
779 j + data->ntoroidalplots + data->npoloidalplots;
780
781 ipoint++;
782 if(ipoint == data->Npnt) {
783 ipoint = 0;
784 }
785
786 data->mrk_pnt[imrk] = ipoint;
787 data->mrk_recorded[imrk] = p_f->mileage[i];
788 }
789 }
790 }
791 }
792 }
793}
794
808real diag_orb_check_plane_crossing(real fang, real iang, real ang0) {
809
810 real k = 0;
811 /* Check whether nag0 is between iang and fang. Note that this *
812 * implementation works only because iang and fang are cumulative. */
813 if( floor( (fang - ang0)/CONST_2PI ) != floor( (iang - ang0)/CONST_2PI ) ) {
814
815 /* Move iang to interval [0, 2pi] */
816 real a = fmod(iang, CONST_2PI);
817 if(a < 0){
818 a = CONST_2PI + a;
819 }
820
821 a = fabs(ang0 - a);
822 if(a > CONST_PI) {
823 a = CONST_2PI - a;
824 }
825 k = fabs(a / (fang - iang));
826 }
827 return k;
828}
829
842real diag_orb_check_radial_crossing(real frho, real irho, real rho0){
843
844 real k = 0;
845 if((frho <= rho0 && irho > rho0) || (irho <= rho0 && frho > rho0)){
846 k = fabs((irho - rho0) / (frho - irho));
847 }
848 return k;
849}
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
integer
long integer
Definition ascot5.h:84
consts.h
Header file containing physical and mathematical constants.
CONST_PI
#define CONST_PI
pi
Definition consts.h:11
CONST_2PI
#define CONST_2PI
2*pi
Definition consts.h:14
diag_orb_update_fo
void diag_orb_update_fo(diag_orb_data *data, particle_simd_fo *p_f, particle_simd_fo *p_i)
Collects orbit diagnostics when marker represents a particle.
Definition diag_orb.c:211
diag_orb_check_plane_crossing
real diag_orb_check_plane_crossing(real fang, real iang, real ang0)
Check if marker has crossed a plane.
Definition diag_orb.c:808
diag_orb_update_ml
void diag_orb_update_ml(diag_orb_data *data, particle_simd_ml *p_f, particle_simd_ml *p_i)
Collects orbit diagnostics when marker represents a field line.
Definition diag_orb.c:624
diag_orb_update_gc
void diag_orb_update_gc(diag_orb_data *data, particle_simd_gc *p_f, particle_simd_gc *p_i)
Collects orbit diagnostics when marker represents a guiding center.
Definition diag_orb.c:417
diag_orb_free
void diag_orb_free(diag_orb_data *data)
Free allocated resources.
Definition diag_orb.c:119
diag_orb_check_radial_crossing
real diag_orb_check_radial_crossing(real frho, real irho, real rho0)
Check if marker has crossed given rho.
Definition diag_orb.c:842
diag_orb_init
void diag_orb_init(diag_orb_data *data)
Initializes orbit diagnostics data.
Definition diag_orb.c:26
diag_orb.h
Header file for diag_orb.c.
DIAG_ORB_POINCARE
#define DIAG_ORB_POINCARE
Definition diag_orb.h:13
DIAG_ORB_ML
#define DIAG_ORB_ML
Definition diag_orb.h:24
DIAG_ORB_INTERVAL
#define DIAG_ORB_INTERVAL
Definition diag_orb.h:14
DIAG_ORB_GC
#define DIAG_ORB_GC
Definition diag_orb.h:23
DIAG_ORB_FO
#define DIAG_ORB_FO
Definition diag_orb.h:22
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.
simulate.h
Header file for simulate.c.
simulate_mode_fo
@ simulate_mode_fo
Definition simulate.h:33
simulate_mode_ml
@ simulate_mode_ml
Definition simulate.h:44
simulate_mode_gc
@ simulate_mode_gc
Definition simulate.h:36
simulate_mode_hybrid
@ simulate_mode_hybrid
Definition simulate.h:41
diag_orb_data
Orbit diagnostics data struct.
Definition diag_orb.h:42
diag_orb_data::mileage
real * mileage
Definition diag_orb.h:44
diag_orb_data::poloidalangles
real poloidalangles[DIAG_ORB_MAXPOINCARES]
Definition diag_orb.h:78
diag_orb_data::ppar
real * ppar
Definition diag_orb.h:51
diag_orb_data::simmode
real * simmode
Definition diag_orb.h:61
diag_orb_data::B_r
real * B_r
Definition diag_orb.h:58
diag_orb_data::Nmrk
int Nmrk
Definition diag_orb.h:70
diag_orb_data::ntoroidalplots
int ntoroidalplots
Definition diag_orb.h:73
diag_orb_data::Npnt
int Npnt
Definition diag_orb.h:69
diag_orb_data::toroidalangles
real toroidalangles[DIAG_ORB_MAXPOINCARES]
Definition diag_orb.h:77
diag_orb_data::pncrdi
real * pncrdi
Definition diag_orb.h:63
diag_orb_data::writeInterval
real writeInterval
Definition diag_orb.h:72
diag_orb_data::r
real * r
Definition diag_orb.h:45
diag_orb_data::zeta
real * zeta
Definition diag_orb.h:53
diag_orb_data::p_phi
real * p_phi
Definition diag_orb.h:49
diag_orb_data::theta
real * theta
Definition diag_orb.h:57
diag_orb_data::z
real * z
Definition diag_orb.h:47
diag_orb_data::radialdistances
real radialdistances[DIAG_ORB_MAXPOINCARES]
Definition diag_orb.h:79
diag_orb_data::rho
real * rho
Definition diag_orb.h:56
diag_orb_data::pncrid
real * pncrid
Definition diag_orb.h:62
diag_orb_data::weight
real * weight
Definition diag_orb.h:54
diag_orb_data::mrk_recorded
real * mrk_recorded
Definition diag_orb.h:66
diag_orb_data::B_z
real * B_z
Definition diag_orb.h:60
diag_orb_data::phi
real * phi
Definition diag_orb.h:46
diag_orb_data::mrk_pnt
integer * mrk_pnt
Definition diag_orb.h:65
diag_orb_data::nradialplots
int nradialplots
Definition diag_orb.h:75
diag_orb_data::mu
real * mu
Definition diag_orb.h:52
diag_orb_data::p_z
real * p_z
Definition diag_orb.h:50
diag_orb_data::mode
int mode
Definition diag_orb.h:68
diag_orb_data::record_mode
int record_mode
Definition diag_orb.h:76
diag_orb_data::charge
real * charge
Definition diag_orb.h:55
diag_orb_data::id
real * id
Definition diag_orb.h:43
diag_orb_data::p_r
real * p_r
Definition diag_orb.h:48
diag_orb_data::npoloidalplots
int npoloidalplots
Definition diag_orb.h:74
diag_orb_data::B_phi
real * B_phi
Definition diag_orb.h:59
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::charge
real * charge
Definition particle.h:219
particle_simd_fo::r
real * r
Definition particle.h:212
particle_simd_fo::id
integer * id
Definition particle.h:246
particle_simd_fo::weight
real * weight
Definition particle.h:241
particle_simd_fo::rho
real * rho
Definition particle.h:243
particle_simd_fo::mileage
real * mileage
Definition particle.h:251
particle_simd_fo::B_phi
real * B_phi
Definition particle.h:226
particle_simd_fo::index
integer * index
Definition particle.h:255
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::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::phi
real * phi
Definition particle.h:213
particle_simd_fo::z
real * z
Definition particle.h:214
particle_simd_gc
Struct representing NSIMD guiding center markers.
Definition particle.h:275
particle_simd_ml
Struct representing NSIMD field line markers.
Definition particle.h:342
Generated on Wed Mar 19 2025 11:10:23 for ASCOT5 by doxygen 1.13.2