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