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 data->simmode[idx]= DIAG_ORB_FO;
402 }
403 }
404 }
405 }
406 }
407}
408
418void diag_orb_update_gc(diag_orb_data* data, particle_simd_gc* p_f,
419 particle_simd_gc* p_i) {
420
421 if(data->mode == DIAG_ORB_INTERVAL) {
422 #pragma omp simd
423 for(int i= 0; i < NSIMD; i++) {
424
425 /* Mask dummy markers */
426 if(p_f->id[i] > 0) {
427 integer imrk = p_f->index[i];
428 integer ipoint = data->mrk_pnt[imrk];
429 integer idx = imrk * data->Npnt + ipoint;
430
431 /* If this is the first time-step, record marker position. */
432 if( data->id[imrk * data->Npnt] == 0 ) {
433 data->id[idx] = (real)(p_i->id[i]);
434 data->mileage[idx]= p_i->mileage[i];
435 data->r[idx] = p_i->r[i];
436 data->phi[idx] = p_i->phi[i];
437 data->z[idx] = p_i->z[i];
438 data->ppar[idx] = p_i->ppar[i];
439 data->mu[idx] = p_i->mu[i];
440 data->zeta[idx] = p_i->zeta[i];
441 data->weight[idx] = p_i->weight[i];
442 data->charge[idx] = p_i->charge[i];
443 data->rho[idx] = p_i->rho[i];
444 data->theta[idx] = p_i->theta[i];
445 data->B_r[idx] = p_i->B_r[i];
446 data->B_phi[idx] = p_i->B_phi[i];
447 data->B_z[idx] = p_i->B_z[i];
448 data->simmode[idx]= DIAG_ORB_GC;
449
450 ipoint++;
451 if(ipoint == data->Npnt) {
452 ipoint = 0;
453 }
454 data->mrk_pnt[imrk] = ipoint;
455 data->mrk_recorded[imrk] = p_i->mileage[i];
456 }
457
458 /* Record marker if enough time has passed from last record, or
459 if marker has met some end condition. */
460 real dt = data->mrk_recorded[imrk] + data->writeInterval
461 - p_f->mileage[i];
462
463 if( dt <= 0 || p_f->endcond[i] > 0 ) {
464 idx = imrk * data->Npnt + ipoint;
465
466 data->id[idx] = (real)p_f->id[i];
467 data->mileage[idx]= p_f->mileage[i];
468 data->r[idx] = p_f->r[i];
469 data->phi[idx] = p_f->phi[i];
470 data->z[idx] = p_f->z[i];
471 data->ppar[idx] = p_f->ppar[i];
472 data->mu[idx] = p_f->mu[i];
473 data->zeta[idx] = p_f->zeta[i];
474 data->weight[idx] = p_f->weight[i];
475 data->charge[idx] = p_f->charge[i];
476 data->rho[idx] = p_f->rho[i];
477 data->theta[idx] = p_f->theta[i];
478 data->B_r[idx] = p_f->B_r[i];
479 data->B_phi[idx] = p_f->B_phi[i];
480 data->B_z[idx] = p_f->B_z[i];
481 data->simmode[idx]= DIAG_ORB_GC;
482
483 ipoint++;
484 if(ipoint == data->Npnt) {
485 ipoint = 0;
486 }
487 data->mrk_pnt[imrk] = ipoint;
488 data->mrk_recorded[imrk] = p_f->mileage[i];
489 }
490 }
491 }
492 }
493 else if(data->mode == DIAG_ORB_POINCARE) {
494 #pragma omp simd
495 for(int i= 0; i < NSIMD; i++) {
496 /* Mask dummy markers and those whose time-step was rejected. */
497 if( p_f->id[i] > 0 && (p_f->mileage[i] != p_i->mileage[i]) ) {
498
499 real k;
500 integer imrk = p_f->index[i];
501 integer ipoint = data->mrk_pnt[imrk];
502 integer idx = imrk * data->Npnt + ipoint;
503
504 /* Check and store toroidal crossings. */
505 for(int j=0; j < data->ntoroidalplots; j++) {
506 k = diag_orb_check_plane_crossing(p_f->phi[i], p_i->phi[i],
507 data->toroidalangles[j]);
508 if(k) {
509 real d = 1-k;
510 idx = imrk * data->Npnt + ipoint;
511 data->id[idx] = (real)p_f->id[i];
512 data->mileage[idx]= k*p_f->mileage[i]+ d*p_i->mileage[i];
513 data->r[idx] = k*p_f->r[i] + d*p_i->r[i];
514 data->phi[idx] = k*p_f->phi[i] + d*p_i->phi[i];
515 data->z[idx] = k*p_f->z[i] + d*p_i->z[i];
516 data->ppar[idx] = k*p_f->ppar[i] + d*p_i->ppar[i];
517 data->mu[idx] = k*p_f->mu[i] + d*p_i->mu[i];
518 data->zeta[idx] = k*p_f->zeta[i] + d*p_i->zeta[i];
519 data->weight[idx] = k*p_f->weight[i] + d*p_i->weight[i];
520 data->charge[idx] = p_i->charge[i];
521 data->rho[idx] = k*p_f->rho[i] + d*p_i->rho[i];
522 data->theta[idx] = k*p_f->theta[i] + d*p_i->theta[i];
523 data->B_r[idx] = k*p_f->B_r[i] + d*p_i->B_r[i];
524 data->B_phi[idx] = k*p_f->B_phi[i] + d*p_i->B_phi[i];
525 data->B_z[idx] = k*p_f->B_z[i] + d*p_i->B_z[i];
526 data->pncrid[idx] = j;
527 data->pncrdi[idx] = 1 - 2 * (p_f->phi[i] < p_i->phi[i]);
528 data->simmode[idx]= DIAG_ORB_GC;
529
530 ipoint++;
531 if(ipoint == data->Npnt) {
532 ipoint = 0;
533 }
534 data->mrk_pnt[imrk] = ipoint;
535 data->mrk_recorded[imrk] = p_f->mileage[i];
536 }
537 }
538
539 /* Check and store poloidal crossings. */
540 for(int j=0; j < data->npoloidalplots; j++) {
541 k = diag_orb_check_plane_crossing(p_f->theta[i],
542 p_i->theta[i],
543 data->poloidalangles[j]);
544 if(k) {
545 real d = 1-k;
546 idx = imrk * data->Npnt + ipoint;
547 data->id[idx] = (real)p_f->id[i];
548 data->mileage[idx]= k*p_f->mileage[i]+ d*p_i->mileage[i];
549 data->r[idx] = k*p_f->r[i] + d*p_i->r[i];
550 data->phi[idx] = k*p_f->phi[i] + d*p_i->phi[i];
551 data->z[idx] = k*p_f->z[i] + d*p_i->z[i];
552 data->ppar[idx] = k*p_f->ppar[i] + d*p_i->ppar[i];
553 data->mu[idx] = k*p_f->mu[i] + d*p_i->mu[i];
554 data->zeta[idx] = k*p_f->zeta[i] + d*p_i->zeta[i];
555 data->weight[idx] = k*p_f->weight[i] + d*p_i->weight[i];
556 data->charge[idx] = p_i->charge[i];
557 data->rho[idx] = k*p_f->rho[i] + d*p_i->rho[i];
558 data->theta[idx] = k*p_f->theta[i] + d*p_i->theta[i];
559 data->B_r[idx] = k*p_f->B_r[i] + d*p_i->B_r[i];
560 data->B_phi[idx] = k*p_f->B_phi[i] + d*p_i->B_phi[i];
561 data->B_z[idx] = k*p_f->B_z[i] + d*p_i->B_z[i];
562 data->pncrid[idx] = j + data->ntoroidalplots;
563 data->pncrdi[idx] = 1 - 2 * (p_f->theta[i] < p_i->theta[i]);
564 data->simmode[idx]= DIAG_ORB_GC;
565
566 ipoint++;
567 if(ipoint == data->Npnt) {
568 ipoint = 0;
569 }
570 data->mrk_pnt[imrk] = ipoint;
571 data->mrk_recorded[imrk] = p_f->mileage[i];
572 }
573 }
574
575
576 /* Check and store radial crossings. */
577 for(int j=0; j < data->nradialplots; j++) {
578 k = diag_orb_check_radial_crossing(p_f->rho[i],
579 p_i->rho[i],
580 data->radialdistances[j]);
581 if(k) {
582 real d = 1-k;
583 idx = imrk * data->Npnt + ipoint;
584 data->id[idx] = (real)p_f->id[i];
585 data->mileage[idx]= k*p_f->mileage[i]+ d*p_i->mileage[i];
586 data->r[idx] = k*p_f->r[i] + d*p_i->r[i];
587 data->phi[idx] = k*p_f->phi[i] + d*p_i->phi[i];
588 data->z[idx] = k*p_f->z[i] + d*p_i->z[i];
589 data->ppar[idx] = k*p_f->ppar[i] + d*p_i->ppar[i];
590 data->mu[idx] = k*p_f->mu[i] + d*p_i->mu[i];
591 data->zeta[idx] = k*p_f->zeta[i] + d*p_i->zeta[i];
592 data->weight[idx] = k*p_f->weight[i] + d*p_i->weight[i];
593 data->charge[idx] = p_i->charge[i];
594 data->rho[idx] = k*p_f->rho[i] + d*p_i->rho[i];
595 data->theta[idx] = k*p_f->theta[i] + d*p_i->theta[i];
596 data->B_r[idx] = k*p_f->B_r[i] + d*p_i->B_r[i];
597 data->B_phi[idx] = k*p_f->B_phi[i] + d*p_i->B_phi[i];
598 data->B_z[idx] = k*p_f->B_z[i] + d*p_i->B_z[i];
599 data->pncrid[idx] =
600 j + data->ntoroidalplots + data->npoloidalplots;
601 data->pncrdi[idx] = 1 - 2 * (p_f->rho[i] < p_i->rho[i]);
602
603 ipoint++;
604 if(ipoint == data->Npnt) {
605 ipoint = 0;
606 }
607 data->mrk_pnt[imrk] = ipoint;
608 data->mrk_recorded[imrk] = p_f->mileage[i];
609 data->simmode[idx]= DIAG_ORB_GC;
610 }
611 }
612 }
613 }
614 }
615}
616
626void diag_orb_update_ml(diag_orb_data* data, particle_simd_ml* p_f,
627 particle_simd_ml* p_i) {
628
629 if(data->mode == DIAG_ORB_INTERVAL) {
630
631 #pragma omp simd
632 for(int i= 0; i < NSIMD; i++) {
633
634 /* Mask dummy markers */
635 if(p_f->id[i] > 0) {
636 integer imrk = p_f->index[i];
637 integer ipoint = data->mrk_pnt[imrk];
638 integer idx = imrk * data->Npnt + ipoint;
639
640 /* If this is the first time-step, record marker position. */
641 if( data->id[imrk * data->Npnt] == 0 ) {
642 data->id[idx] = (real)p_i->id[i];
643 data->mileage[idx] = p_i->mileage[i];
644 data->r[idx] = p_i->r[i];
645 data->phi[idx] = p_i->phi[i];
646 data->z[idx] = p_i->z[i];
647 data->rho[idx] = p_i->rho[i];
648 data->theta[idx] = p_i->theta[i];
649 data->B_r[idx] = p_i->B_r[i];
650 data->B_phi[idx] = p_i->B_phi[i];
651 data->B_z[idx] = p_i->B_z[i];
652 data->simmode[idx] = DIAG_ORB_ML;
653
654 ipoint++;
655 if(ipoint == data->Npnt) {
656 ipoint = 0;
657 }
658 data->mrk_pnt[imrk] = ipoint;
659 data->mrk_recorded[imrk] = p_i->mileage[i];
660 }
661
662 /* Record marker if enough time has passed from last record, or
663 if marker has met some end condition. */
664 real dt = data->mrk_recorded[imrk] + data->writeInterval
665 - p_f->mileage[i];
666 if( dt <= 0 || p_f->endcond[i] > 0 ) {
667 idx = imrk * data->Npnt + ipoint;
668 data->id[idx] = (real)p_f->id[i];
669 data->mileage[idx] = p_f->mileage[i];
670 data->r[idx] = p_f->r[i];
671 data->phi[idx] = p_f->phi[i];
672 data->z[idx] = p_f->z[i];
673 data->rho[idx] = p_f->rho[i];
674 data->theta[idx] = p_f->theta[i];
675 data->B_r[idx] = p_f->B_r[i];
676 data->B_phi[idx] = p_f->B_phi[i];
677 data->B_z[idx] = p_f->B_z[i];
678 data->simmode[idx] = DIAG_ORB_ML;
679
680 ipoint++;
681 if(ipoint == data->Npnt) {
682 ipoint = 0;
683 }
684 data->mrk_pnt[imrk] = ipoint;
685 data->mrk_recorded[imrk] = p_f->mileage[i];
686 }
687 }
688 }
689 }
690 else if(data->mode == DIAG_ORB_POINCARE) {
691 #pragma omp simd
692 for(int i= 0; i < NSIMD; i++) {
693 /* Mask dummy markers and thosw whose time-step was rejected. */
694 if( p_f->id[i] > 0 && (p_f->mileage[i] != p_i->mileage[i]) ) {
695
696 real k;
697 integer imrk = p_f->index[i];
698 integer ipoint = data->mrk_pnt[imrk];
699 integer idx = imrk * data->Npnt + ipoint;
700
701 /* Check and store toroidal crossings. */
702 for(int j=0; j < data->ntoroidalplots; j++) {
703 k = diag_orb_check_plane_crossing(p_f->phi[i], p_i->phi[i],
704 data->toroidalangles[j]);
705 if(k) {
706 real d = 1-k;
707 idx = imrk * data->Npnt + ipoint;
708 data->id[idx] = (real)p_f->id[i];
709 data->mileage[idx]= k*p_f->mileage[i]+ d*p_i->mileage[i];
710 data->r[idx] = k*p_f->r[i] + d*p_i->r[i];
711 data->phi[idx] = k*p_f->phi[i] + d*p_i->phi[i];
712 data->z[idx] = k*p_f->z[i] + d*p_i->z[i];
713 data->rho[idx] = k*p_f->rho[i] + d*p_i->rho[i];
714 data->theta[idx] = k*p_f->theta[i] + d*p_i->theta[i];
715 data->B_r[idx] = k*p_f->B_r[i] + d*p_i->B_r[i];
716 data->B_phi[idx] = k*p_f->B_phi[i] + d*p_i->B_phi[i];
717 data->B_z[idx] = k*p_f->B_z[i] + d*p_i->B_z[i];
718 data->pncrid[idx] = j;
719 data->pncrdi[idx] = 1 - 2 * (p_f->phi[i] < p_i->phi[i]);
720 data->simmode[idx]= DIAG_ORB_ML;
721
722 ipoint++;
723 if(ipoint == data->Npnt) {
724 ipoint = 0;
725 }
726 data->mrk_pnt[imrk] = ipoint;
727 data->mrk_recorded[imrk] = p_f->mileage[i];
728 }
729 }
730
731 /* Check and store poloidal crossings. */
732 for(int j=0; j < data->npoloidalplots; j++) {
733 k = diag_orb_check_plane_crossing(p_f->theta[i],
734 p_i->theta[i],
735 data->poloidalangles[j]);
736 if(k) {
737 real d = 1-k;
738 idx = imrk * data->Npnt + ipoint;
739 data->id[idx] = (real)p_f->id[i];
740 data->mileage[idx]= k*p_f->mileage[i] + d*p_i->mileage[i];
741 data->r[idx] = k*p_f->r[i] + d*p_i->r[i];
742 data->phi[idx] = k*p_f->phi[i] + d*p_i->phi[i];
743 data->z[idx] = k*p_f->z[i] + d*p_i->z[i];
744 data->rho[idx] = k*p_f->rho[i] + d*p_i->rho[i];
745 data->theta[idx] = k*p_f->theta[i] + d*p_i->theta[i];
746 data->B_r[idx] = k*p_f->B_r[i] + d*p_i->B_r[i];
747 data->B_phi[idx] = k*p_f->B_phi[i] + d*p_i->B_phi[i];
748 data->B_z[idx] = k*p_f->B_z[i] + d*p_i->B_z[i];
749 data->pncrid[idx] = j + data->ntoroidalplots;
750 data->pncrdi[idx] = 1 - 2 * (p_f->theta[i] < p_i->theta[i]);
751 data->simmode[idx]= DIAG_ORB_ML;
752
753 ipoint++;
754 if(ipoint == data->Npnt) {
755 ipoint = 0;
756 }
757 data->mrk_pnt[imrk] = ipoint;
758 data->mrk_recorded[imrk] = p_f->mileage[i];
759 }
760 }
761
762 /* Check and store radial crossings. */
763 for(int j=0; j < data->nradialplots; j++) {
764 k = diag_orb_check_radial_crossing(p_f->rho[i],
765 p_i->rho[i],
766 data->radialdistances[j]);
767 if(k) {
768 real d = 1-k;
769 idx = imrk * data->Npnt + ipoint;
770 data->id[idx] = (real)p_f->id[i];
771 data->mileage[idx]= k*p_f->mileage[i] + d*p_i->mileage[i];
772 data->r[idx] = k*p_f->r[i] + d*p_i->r[i];
773 data->phi[idx] = k*p_f->phi[i] + d*p_i->phi[i];
774 data->z[idx] = k*p_f->z[i] + d*p_i->z[i];
775 data->rho[idx] = k*p_f->rho[i] + d*p_i->rho[i];
776 data->theta[idx] = k*p_f->theta[i] + d*p_i->theta[i];
777 data->B_r[idx] = k*p_f->B_r[i] + d*p_i->B_r[i];
778 data->B_phi[idx] = k*p_f->B_phi[i] + d*p_i->B_phi[i];
779 data->B_z[idx] = k*p_f->B_z[i] + d*p_i->B_z[i];
780 data->pncrid[idx] =
781 j + data->ntoroidalplots + data->npoloidalplots;
782
783 ipoint++;
784 if(ipoint == data->Npnt) {
785 ipoint = 0;
786 }
787 data->simmode[idx]= DIAG_ORB_ML;
788 data->mrk_pnt[imrk] = ipoint;
789 data->mrk_recorded[imrk] = p_f->mileage[i];
790 }
791 }
792 }
793 }
794 }
795}
796
810real diag_orb_check_plane_crossing(real fang, real iang, real ang0) {
811
812 real k = 0;
813 /* Check whether nag0 is between iang and fang. Note that this *
814 * implementation works only because iang and fang are cumulative. */
815 if( floor( (fang - ang0)/CONST_2PI ) != floor( (iang - ang0)/CONST_2PI ) ) {
816
817 /* Move iang to interval [0, 2pi] */
818 real a = fmod(iang, CONST_2PI);
819 if(a < 0){
820 a = CONST_2PI + a;
821 }
822
823 a = fabs(ang0 - a);
824 if(a > CONST_PI) {
825 a = CONST_2PI - a;
826 }
827 k = fabs(a / (fang - iang));
828 }
829 return k;
830}
831
844real diag_orb_check_radial_crossing(real frho, real irho, real rho0){
845
846 real k = 0;
847 if((frho <= rho0 && irho > rho0) || (irho <= rho0 && frho > rho0)){
848 k = fabs((irho - rho0) / (frho - irho));
849 }
850 return k;
851}
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:810
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:626
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:418
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:844
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:34
simulate_mode_ml
@ simulate_mode_ml
Definition simulate.h:45
simulate_mode_gc
@ simulate_mode_gc
Definition simulate.h:37
simulate_mode_hybrid
@ simulate_mode_hybrid
Definition simulate.h:42
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_gc::weight
real * weight
Definition particle.h:307
particle_simd_gc::z
real * z
Definition particle.h:279
particle_simd_gc::ppar
real * ppar
Definition particle.h:280
particle_simd_gc::mu
real * mu
Definition particle.h:281
particle_simd_gc::zeta
real * zeta
Definition particle.h:282
particle_simd_gc::rho
real * rho
Definition particle.h:309
particle_simd_gc::id
integer * id
Definition particle.h:312
particle_simd_gc::phi
real * phi
Definition particle.h:278
particle_simd_gc::index
integer * index
Definition particle.h:323
particle_simd_gc::charge
real * charge
Definition particle.h:284
particle_simd_gc::r
real * r
Definition particle.h:277
particle_simd_gc::B_phi
real * B_phi
Definition particle.h:290
particle_simd_gc::theta
real * theta
Definition particle.h:310
particle_simd_gc::mileage
real * mileage
Definition particle.h:318
particle_simd_gc::B_z
real * B_z
Definition particle.h:292
particle_simd_gc::B_r
real * B_r
Definition particle.h:288
particle_simd_ml
Struct representing NSIMD field line markers.
Definition particle.h:343
Generated on for ASCOT5 by doxygen 1.18.0