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 = malloc( data->Nmrk*sizeof(integer) );
108 data->mrk_recorded = 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
133void diag_orb_update_fo(diag_orb_data* data, particle_simd_fo* p_f,
134 particle_simd_fo* p_i) {
135
136 if(data->mode == DIAG_ORB_INTERVAL) {
137
138 #pragma omp simd
139 for(int i= 0; i < NSIMD; i++) {
140
141 /* Mask dummy markers */
142 if(p_f->id[i] > 0) {
143
144 integer imrk = p_f->index[i];
145 integer ipoint = data->mrk_pnt[imrk];
146 integer idx = imrk * data->Npnt + ipoint;
147
148 /* If this is the first time-step, record marker position. */
149 if( data->id[imrk * data->Npnt] == 0 ) {
150 data->id[idx] = (real)p_i->id[i];
151 data->mileage[idx]= p_i->mileage[i];
152 data->r[idx] = p_i->r[i];
153 data->phi[idx] = p_i->phi[i];
154 data->z[idx] = p_i->z[i];
155 data->p_r[idx] = p_i->p_r[i];
156 data->p_phi[idx] = p_i->p_phi[i];
157 data->p_z[idx] = p_i->p_z[i];
158 data->weight[idx] = p_i->weight[i];
159 data->charge[idx] = p_i->charge[i];
160 data->rho[idx] = p_i->rho[i];
161 data->theta[idx] = p_i->theta[i];
162 data->B_r[idx] = p_i->B_r[i];
163 data->B_phi[idx] = p_i->B_phi[i];
164 data->B_z[idx] = p_i->B_z[i];
165 data->simmode[idx]= DIAG_ORB_FO;
166
167 ipoint++;
168 if(ipoint == data->Npnt) {
169 ipoint = 0;
170 }
171 data->mrk_pnt[imrk] = ipoint;
172 data->mrk_recorded[imrk] = p_i->mileage[i];
173 }
174
175 /* Record marker if enough time has passed from last record, or
176 if marker has met some end condition. */
177 real dt = data->mrk_recorded[imrk] + data->writeInterval
178 - p_f->mileage[i];
179 if( dt <= 0 || p_f->endcond[i] > 0 ) {
180 idx = imrk * data->Npnt + ipoint;
181
182 data->id[idx] = (real)p_f->id[i];
183 data->mileage[idx]= p_f->mileage[i];
184 data->r[idx] = p_f->r[i];
185 data->phi[idx] = p_f->phi[i];
186 data->z[idx] = p_f->z[i];
187 data->p_r[idx] = p_f->p_r[i];
188 data->p_phi[idx] = p_f->p_phi[i];
189 data->p_z[idx] = p_f->p_z[i];
190 data->weight[idx] = p_f->weight[i];
191 data->charge[idx] = p_f->charge[i];
192 data->rho[idx] = p_f->rho[i];
193 data->theta[idx] = p_f->theta[i];
194 data->B_r[idx] = p_f->B_r[i];
195 data->B_phi[idx] = p_f->B_phi[i];
196 data->B_z[idx] = p_f->B_z[i];
197 data->simmode[idx]= DIAG_ORB_FO;
198
199 ipoint++;
200 if(ipoint == data->Npnt) {
201 ipoint = 0;
202 }
203 data->mrk_pnt[imrk] = ipoint;
204 data->mrk_recorded[imrk] = p_f->mileage[i];
205 }
206 }
207 }
208 }
209 else if(data->mode == DIAG_ORB_POINCARE) {
210
211 #pragma omp simd
212 for(int i= 0; i < NSIMD; i++) {
213 /* Mask dummy markers and those whose time-step was rejected. */
214 if( p_f->id[i] > 0 && (p_f->mileage[i] != p_i->mileage[i]) ) {
215
216 real k;
217 integer imrk = p_f->index[i];
218 integer ipoint = data->mrk_pnt[imrk];
219 integer idx = imrk * data->Npnt + ipoint;
220
221 /* Check and store toroidal crossings. */
222 for(int j=0; j < data->ntoroidalplots; j++) {
223 k = diag_orb_check_plane_crossing(p_f->phi[i], p_i->phi[i],
224 data->toroidalangles[j]);
225 if(k) {
226 real d = 1-k;
227 idx = imrk * data->Npnt + ipoint;
228 data->id[idx] = (real)p_f->id[i];
229 data->mileage[idx]= k*p_f->mileage[i]+ d*p_i->mileage[i];
230 data->r[idx] = k*p_f->r[i] + d*p_i->r[i];
231 data->phi[idx] = k*p_f->phi[i] + d*p_i->phi[i];
232 data->z[idx] = k*p_f->z[i] + d*p_i->z[i];
233 data->p_r[idx] = k*p_f->p_r[i] + d*p_i->p_r[i];
234 data->p_phi[idx] = k*p_f->p_phi[i] + d*p_i->p_phi[i];
235 data->p_z[idx] = k*p_f->p_z[i] + d*p_i->p_z[i];
236 data->weight[idx] = k*p_f->weight[i] + d*p_i->weight[i];
237 data->charge[idx] = p_i->charge[i];
238 data->rho[idx] = k*p_f->rho[i] + d*p_i->rho[i];
239 data->theta[idx] = k*p_f->theta[i] + d*p_i->theta[i];
240 data->B_r[idx] = k*p_f->B_r[i] + d*p_i->B_r[i];
241 data->B_phi[idx] = k*p_f->B_phi[i] + d*p_i->B_phi[i];
242 data->B_z[idx] = k*p_f->B_z[i] + d*p_i->B_z[i];
243 data->pncrid[idx] = j;
244 data->pncrdi[idx] = 1 - 2 * (p_f->phi[i] < p_i->phi[i]);
245 data->simmode[idx]= DIAG_ORB_FO;
246
247 ipoint++;
248 if(ipoint == data->Npnt) {
249 ipoint = 0;
250 }
251 data->mrk_pnt[imrk] = ipoint;
252 data->mrk_recorded[imrk] = p_f->mileage[i];
253 }
254 }
255
256 /* Check and store poloidal crossings. */
257 for(int j=0; j < data->npoloidalplots; j++) {
258 k = diag_orb_check_plane_crossing(p_f->theta[i],
259 p_i->theta[i],
260 data->poloidalangles[j]);
261 if(k) {
262 real d = 1-k;
263 idx = imrk * data->Npnt + ipoint;
264 data->id[idx] = (real)p_f->id[i];
265 data->mileage[idx]= k*p_f->mileage[i]+ d*p_i->mileage[i];
266 data->r[idx] = k*p_f->r[i] + d*p_i->r[i];
267 data->phi[idx] = k*p_f->phi[i] + d*p_i->phi[i];
268 data->z[idx] = k*p_f->z[i] + d*p_i->z[i];
269 data->p_r[idx] = k*p_f->p_r[i] + d*p_i->p_r[i];
270 data->p_phi[idx] = k*p_f->p_phi[i] + d*p_i->p_phi[i];
271 data->p_z[idx] = k*p_f->p_z[i] + d*p_i->p_z[i];
272 data->weight[idx] = k*p_f->weight[i] + d*p_i->weight[i];
273 data->charge[idx] = p_i->charge[i];
274 data->rho[idx] = k*p_f->rho[i] + d*p_i->rho[i];
275 data->theta[idx] = k*p_f->theta[i] + d*p_i->theta[i];
276 data->B_r[idx] = k*p_f->B_r[i] + d*p_i->B_r[i];
277 data->B_phi[idx] = k*p_f->B_phi[i] + d*p_i->B_phi[i];
278 data->B_z[idx] = k*p_f->B_z[i] + d*p_i->B_z[i];
279 data->pncrid[idx] = j + data->ntoroidalplots;
280 data->pncrdi[idx] = 1 - 2 * (p_f->theta[i] < p_i->theta[i]);
281 data->simmode[idx]= DIAG_ORB_FO;
282
283 ipoint++;
284 if(ipoint == data->Npnt) {
285 ipoint = 0;
286 }
287 data->mrk_pnt[imrk] = ipoint;
288 data->mrk_recorded[imrk] = p_f->mileage[i];
289 }
290 }
291
292 /* Check and store radial crossings. */
293 for(int j=0; j < data->nradialplots; j++) {
294 k = diag_orb_check_radial_crossing(p_f->rho[i],p_i->rho[i],
295 data->radialdistances[j]);
296 if(k) {
297 real d = k;
298 k = 1-d;
299 idx = imrk * data->Npnt + ipoint;
300 data->id[idx] = (real)p_f->id[i];
301 data->mileage[idx]= k*p_f->mileage[i]+ d*p_i->mileage[i];
302 data->r[idx] = k*p_f->r[i] + d*p_i->r[i];
303 data->phi[idx] = k*p_f->phi[i] + d*p_i->phi[i];
304 data->z[idx] = k*p_f->z[i] + d*p_i->z[i];
305 data->p_r[idx] = k*p_f->p_r[i] + d*p_i->p_r[i];
306 data->p_phi[idx] = k*p_f->p_phi[i] + d*p_i->p_phi[i];
307 data->p_z[idx] = k*p_f->p_z[i] + d*p_i->p_z[i];
308 data->weight[idx] = k*p_f->weight[i] + d*p_i->weight[i];
309 data->charge[idx] = p_i->charge[i];
310 data->rho[idx] = k*p_f->rho[i] + d*p_i->rho[i];
311 data->theta[idx] = k*p_f->theta[i] + d*p_i->theta[i];
312 data->B_r[idx] = k*p_f->B_r[i] + d*p_i->B_r[i];
313 data->B_phi[idx] = k*p_f->B_phi[i] + d*p_i->B_phi[i];
314 data->B_z[idx] = k*p_f->B_z[i] + d*p_i->B_z[i];
315 data->pncrid[idx] = j + data->ntoroidalplots + data->npoloidalplots;
316 data->pncrdi[idx] = 1 - 2 * (p_f->rho[i] < p_i->rho[i]);
317 ipoint++;
318 if(ipoint == data->Npnt) {
319 ipoint = 0;
320 }
321 data->mrk_pnt[imrk] = ipoint;
322 data->mrk_recorded[imrk] = p_f->mileage[i];
323 }
324 }
325 }
326 }
327 }
328}
329
339void diag_orb_update_gc(diag_orb_data* data, particle_simd_gc* p_f,
340 particle_simd_gc* p_i) {
341
342 if(data->mode == DIAG_ORB_INTERVAL) {
343 #pragma omp simd
344 for(int i= 0; i < NSIMD; i++) {
345
346 /* Mask dummy markers */
347 if(p_f->id[i] > 0) {
348 integer imrk = p_f->index[i];
349 integer ipoint = data->mrk_pnt[imrk];
350 integer idx = imrk * data->Npnt + ipoint;
351
352 /* If this is the first time-step, record marker position. */
353 if( data->id[imrk * data->Npnt] == 0 ) {
354 data->id[idx] = (real)(p_i->id[i]);
355 data->mileage[idx]= p_i->mileage[i];
356 data->r[idx] = p_i->r[i];
357 data->phi[idx] = p_i->phi[i];
358 data->z[idx] = p_i->z[i];
359 data->ppar[idx] = p_i->ppar[i];
360 data->mu[idx] = p_i->mu[i];
361 data->zeta[idx] = p_i->zeta[i];
362 data->weight[idx] = p_i->weight[i];
363 data->charge[idx] = p_i->charge[i];
364 data->rho[idx] = p_i->rho[i];
365 data->theta[idx] = p_i->theta[i];
366 data->B_r[idx] = p_i->B_r[i];
367 data->B_phi[idx] = p_i->B_phi[i];
368 data->B_z[idx] = p_i->B_z[i];
369 data->simmode[idx]= DIAG_ORB_GC;
370
371 ipoint++;
372 if(ipoint == data->Npnt) {
373 ipoint = 0;
374 }
375 data->mrk_pnt[imrk] = ipoint;
376 data->mrk_recorded[imrk] = p_i->mileage[i];
377 }
378
379 /* Record marker if enough time has passed from last record, or
380 if marker has met some end condition. */
381 real dt = data->mrk_recorded[imrk] + data->writeInterval
382 - p_f->mileage[i];
383
384 if( dt <= 0 || p_f->endcond[i] > 0 ) {
385 idx = imrk * data->Npnt + ipoint;
386
387 data->id[idx] = (real)p_f->id[i];
388 data->mileage[idx]= p_f->mileage[i];
389 data->r[idx] = p_f->r[i];
390 data->phi[idx] = p_f->phi[i];
391 data->z[idx] = p_f->z[i];
392 data->ppar[idx] = p_f->ppar[i];
393 data->mu[idx] = p_f->mu[i];
394 data->zeta[idx] = p_f->zeta[i];
395 data->weight[idx] = p_f->weight[i];
396 data->charge[idx] = p_f->charge[i];
397 data->rho[idx] = p_f->rho[i];
398 data->theta[idx] = p_f->theta[i];
399 data->B_r[idx] = p_f->B_r[i];
400 data->B_phi[idx] = p_f->B_phi[i];
401 data->B_z[idx] = p_f->B_z[i];
402 data->simmode[idx]= DIAG_ORB_GC;
403
404 ipoint++;
405 if(ipoint == data->Npnt) {
406 ipoint = 0;
407 }
408 data->mrk_pnt[imrk] = ipoint;
409 data->mrk_recorded[imrk] = p_f->mileage[i];
410 }
411 }
412 }
413 }
414 else if(data->mode == DIAG_ORB_POINCARE) {
415 #pragma omp simd
416 for(int i= 0; i < NSIMD; i++) {
417 /* Mask dummy markers and those whose time-step was rejected. */
418 if( p_f->id[i] > 0 && (p_f->mileage[i] != p_i->mileage[i]) ) {
419
420 real k;
421 integer imrk = p_f->index[i];
422 integer ipoint = data->mrk_pnt[imrk];
423 integer idx = imrk * data->Npnt + ipoint;
424
425 /* Check and store toroidal crossings. */
426 for(int j=0; j < data->ntoroidalplots; j++) {
427 k = diag_orb_check_plane_crossing(p_f->phi[i], p_i->phi[i],
428 data->toroidalangles[j]);
429 if(k) {
430 real d = 1-k;
431 idx = imrk * data->Npnt + ipoint;
432 data->id[idx] = (real)p_f->id[i];
433 data->mileage[idx]= k*p_f->mileage[i]+ d*p_i->mileage[i];
434 data->r[idx] = k*p_f->r[i] + d*p_i->r[i];
435 data->phi[idx] = k*p_f->phi[i] + d*p_i->phi[i];
436 data->z[idx] = k*p_f->z[i] + d*p_i->z[i];
437 data->ppar[idx] = k*p_f->ppar[i] + d*p_i->ppar[i];
438 data->mu[idx] = k*p_f->mu[i] + d*p_i->mu[i];
439 data->zeta[idx] = k*p_f->zeta[i] + d*p_i->zeta[i];
440 data->weight[idx] = k*p_f->weight[i] + d*p_i->weight[i];
441 data->charge[idx] = p_i->charge[i];
442 data->rho[idx] = k*p_f->rho[i] + d*p_i->rho[i];
443 data->theta[idx] = k*p_f->theta[i] + d*p_i->theta[i];
444 data->B_r[idx] = k*p_f->B_r[i] + d*p_i->B_r[i];
445 data->B_phi[idx] = k*p_f->B_phi[i] + d*p_i->B_phi[i];
446 data->B_z[idx] = k*p_f->B_z[i] + d*p_i->B_z[i];
447 data->pncrid[idx] = j;
448 data->pncrdi[idx] = 1 - 2 * (p_f->phi[i] < p_i->phi[i]);
449 data->simmode[idx]= DIAG_ORB_GC;
450
451 ipoint++;
452 if(ipoint == data->Npnt) {
453 ipoint = 0;
454 }
455 data->mrk_pnt[imrk] = ipoint;
456 data->mrk_recorded[imrk] = p_f->mileage[i];
457 }
458 }
459
460 /* Check and store poloidal crossings. */
461 for(int j=0; j < data->npoloidalplots; j++) {
462 k = diag_orb_check_plane_crossing(p_f->theta[i],
463 p_i->theta[i],
464 data->poloidalangles[j]);
465 if(k) {
466 real d = 1-k;
467 idx = imrk * data->Npnt + ipoint;
468 data->id[idx] = (real)p_f->id[i];
469 data->mileage[idx]= k*p_f->mileage[i]+ d*p_i->mileage[i];
470 data->r[idx] = k*p_f->r[i] + d*p_i->r[i];
471 data->phi[idx] = k*p_f->phi[i] + d*p_i->phi[i];
472 data->z[idx] = k*p_f->z[i] + d*p_i->z[i];
473 data->ppar[idx] = k*p_f->ppar[i] + d*p_i->ppar[i];
474 data->mu[idx] = k*p_f->mu[i] + d*p_i->mu[i];
475 data->zeta[idx] = k*p_f->zeta[i] + d*p_i->zeta[i];
476 data->weight[idx] = k*p_f->weight[i] + d*p_i->weight[i];
477 data->charge[idx] = p_i->charge[i];
478 data->rho[idx] = k*p_f->rho[i] + d*p_i->rho[i];
479 data->theta[idx] = k*p_f->theta[i] + d*p_i->theta[i];
480 data->B_r[idx] = k*p_f->B_r[i] + d*p_i->B_r[i];
481 data->B_phi[idx] = k*p_f->B_phi[i] + d*p_i->B_phi[i];
482 data->B_z[idx] = k*p_f->B_z[i] + d*p_i->B_z[i];
483 data->pncrid[idx] = j + data->ntoroidalplots;
484 data->pncrdi[idx] = 1 - 2 * (p_f->theta[i] < p_i->theta[i]);
485 data->simmode[idx]= DIAG_ORB_GC;
486
487 ipoint++;
488 if(ipoint == data->Npnt) {
489 ipoint = 0;
490 }
491 data->mrk_pnt[imrk] = ipoint;
492 data->mrk_recorded[imrk] = p_f->mileage[i];
493 }
494 }
495
496
497 /* Check and store radial crossings. */
498 for(int j=0; j < data->nradialplots; j++) {
499 k = diag_orb_check_radial_crossing(p_f->rho[i],
500 p_i->rho[i],
501 data->radialdistances[j]);
502 if(k) {
503 real d = 1-k;
504 idx = imrk * data->Npnt + ipoint;
505 data->id[idx] = (real)p_f->id[i];
506 data->mileage[idx]= k*p_f->mileage[i]+ d*p_i->mileage[i];
507 data->r[idx] = k*p_f->r[i] + d*p_i->r[i];
508 data->phi[idx] = k*p_f->phi[i] + d*p_i->phi[i];
509 data->z[idx] = k*p_f->z[i] + d*p_i->z[i];
510 data->ppar[idx] = k*p_f->ppar[i] + d*p_i->ppar[i];
511 data->mu[idx] = k*p_f->mu[i] + d*p_i->mu[i];
512 data->zeta[idx] = k*p_f->zeta[i] + d*p_i->zeta[i];
513 data->weight[idx] = k*p_f->weight[i] + d*p_i->weight[i];
514 data->charge[idx] = p_i->charge[i];
515 data->rho[idx] = k*p_f->rho[i] + d*p_i->rho[i];
516 data->theta[idx] = k*p_f->theta[i] + d*p_i->theta[i];
517 data->B_r[idx] = k*p_f->B_r[i] + d*p_i->B_r[i];
518 data->B_phi[idx] = k*p_f->B_phi[i] + d*p_i->B_phi[i];
519 data->B_z[idx] = k*p_f->B_z[i] + d*p_i->B_z[i];
520 data->pncrid[idx] =
521 j + data->ntoroidalplots + data->npoloidalplots;
522 data->pncrdi[idx] = 1 - 2 * (p_f->rho[i] < p_i->rho[i]);
523
524 ipoint++;
525 if(ipoint == data->Npnt) {
526 ipoint = 0;
527 }
528 data->mrk_pnt[imrk] = ipoint;
529 data->mrk_recorded[imrk] = p_f->mileage[i];
530 }
531 }
532 }
533 }
534 }
535}
536
546void diag_orb_update_ml(diag_orb_data* data, particle_simd_ml* p_f,
547 particle_simd_ml* p_i) {
548
549 if(data->mode == DIAG_ORB_INTERVAL) {
550
551 #pragma omp simd
552 for(int i= 0; i < NSIMD; i++) {
553
554 /* Mask dummy markers */
555 if(p_f->id[i] > 0) {
556 integer imrk = p_f->index[i];
557 integer ipoint = data->mrk_pnt[imrk];
558 integer idx = imrk * data->Npnt + ipoint;
559
560 /* If this is the first time-step, record marker position. */
561 if( data->id[imrk * data->Npnt] == 0 ) {
562 data->id[idx] = (real)p_i->id[i];
563 data->mileage[idx] = p_i->mileage[i];
564 data->r[idx] = p_i->r[i];
565 data->phi[idx] = p_i->phi[i];
566 data->z[idx] = p_i->z[i];
567 data->rho[idx] = p_i->rho[i];
568 data->theta[idx] = p_i->theta[i];
569 data->B_r[idx] = p_i->B_r[i];
570 data->B_phi[idx] = p_i->B_phi[i];
571 data->B_z[idx] = p_i->B_z[i];
572 data->simmode[idx] = DIAG_ORB_ML;
573
574 ipoint++;
575 if(ipoint == data->Npnt) {
576 ipoint = 0;
577 }
578 data->mrk_pnt[imrk] = ipoint;
579 data->mrk_recorded[imrk] = p_i->mileage[i];
580 }
581
582 /* Record marker if enough time has passed from last record, or
583 if marker has met some end condition. */
584 real dt = data->mrk_recorded[imrk] + data->writeInterval
585 - p_f->mileage[i];
586 if( dt <= 0 || p_f->endcond[i] > 0 ) {
587 idx = imrk * data->Npnt + ipoint;
588 data->id[idx] = (real)p_f->id[i];
589 data->mileage[idx] = p_f->mileage[i];
590 data->r[idx] = p_f->r[i];
591 data->phi[idx] = p_f->phi[i];
592 data->z[idx] = p_f->z[i];
593 data->rho[idx] = p_f->rho[i];
594 data->theta[idx] = p_f->theta[i];
595 data->B_r[idx] = p_f->B_r[i];
596 data->B_phi[idx] = p_f->B_phi[i];
597 data->B_z[idx] = p_f->B_z[i];
598 data->simmode[idx] = DIAG_ORB_ML;
599
600 ipoint++;
601 if(ipoint == data->Npnt) {
602 ipoint = 0;
603 }
604 data->mrk_pnt[imrk] = ipoint;
605 data->mrk_recorded[imrk] = p_f->mileage[i];
606 }
607 }
608 }
609 }
610 else if(data->mode == DIAG_ORB_POINCARE) {
611 #pragma omp simd
612 for(int i= 0; i < NSIMD; i++) {
613 /* Mask dummy markers and thosw whose time-step was rejected. */
614 if( p_f->id[i] > 0 && (p_f->mileage[i] != p_i->mileage[i]) ) {
615
616 real k;
617 integer imrk = p_f->index[i];
618 integer ipoint = data->mrk_pnt[imrk];
619 integer idx = imrk * data->Npnt + ipoint;
620
621 /* Check and store toroidal crossings. */
622 for(int j=0; j < data->ntoroidalplots; j++) {
623 k = diag_orb_check_plane_crossing(p_f->phi[i], p_i->phi[i],
624 data->toroidalangles[j]);
625 if(k) {
626 real d = 1-k;
627 idx = imrk * data->Npnt + ipoint;
628 data->id[idx] = (real)p_f->id[i];
629 data->mileage[idx]= k*p_f->mileage[i]+ d*p_i->mileage[i];
630 data->r[idx] = k*p_f->r[i] + d*p_i->r[i];
631 data->phi[idx] = k*p_f->phi[i] + d*p_i->phi[i];
632 data->z[idx] = k*p_f->z[i] + d*p_i->z[i];
633 data->rho[idx] = k*p_f->rho[i] + d*p_i->rho[i];
634 data->theta[idx] = k*p_f->theta[i] + d*p_i->theta[i];
635 data->B_r[idx] = k*p_f->B_r[i] + d*p_i->B_r[i];
636 data->B_phi[idx] = k*p_f->B_phi[i] + d*p_i->B_phi[i];
637 data->B_z[idx] = k*p_f->B_z[i] + d*p_i->B_z[i];
638 data->pncrid[idx] = j;
639 data->pncrdi[idx] = 1 - 2 * (p_f->phi[i] < p_i->phi[i]);
640 data->simmode[idx]= DIAG_ORB_ML;
641
642 ipoint++;
643 if(ipoint == data->Npnt) {
644 ipoint = 0;
645 }
646 data->mrk_pnt[imrk] = ipoint;
647 data->mrk_recorded[imrk] = p_f->mileage[i];
648 }
649 }
650
651 /* Check and store poloidal crossings. */
652 for(int j=0; j < data->npoloidalplots; j++) {
653 k = diag_orb_check_plane_crossing(p_f->theta[i],
654 p_i->theta[i],
655 data->poloidalangles[j]);
656 if(k) {
657 real d = 1-k;
658 idx = imrk * data->Npnt + ipoint;
659 data->id[idx] = (real)p_f->id[i];
660 data->mileage[idx]= k*p_f->mileage[i] + d*p_i->mileage[i];
661 data->r[idx] = k*p_f->r[i] + d*p_i->r[i];
662 data->phi[idx] = k*p_f->phi[i] + d*p_i->phi[i];
663 data->z[idx] = k*p_f->z[i] + d*p_i->z[i];
664 data->rho[idx] = k*p_f->rho[i] + d*p_i->rho[i];
665 data->theta[idx] = k*p_f->theta[i] + d*p_i->theta[i];
666 data->B_r[idx] = k*p_f->B_r[i] + d*p_i->B_r[i];
667 data->B_phi[idx] = k*p_f->B_phi[i] + d*p_i->B_phi[i];
668 data->B_z[idx] = k*p_f->B_z[i] + d*p_i->B_z[i];
669 data->pncrid[idx] = j + data->ntoroidalplots;
670 data->pncrdi[idx] = 1 - 2 * (p_f->theta[i] < p_i->theta[i]);
671 data->simmode[idx]= DIAG_ORB_ML;
672
673 ipoint++;
674 if(ipoint == data->Npnt) {
675 ipoint = 0;
676 }
677 data->mrk_pnt[imrk] = ipoint;
678 data->mrk_recorded[imrk] = p_f->mileage[i];
679 }
680 }
681
682 /* Check and store radial crossings. */
683 for(int j=0; j < data->nradialplots; j++) {
684 k = diag_orb_check_radial_crossing(p_f->rho[i],
685 p_i->rho[i],
686 data->radialdistances[j]);
687 if(k) {
688 real d = 1-k;
689 idx = imrk * data->Npnt + ipoint;
690 data->id[idx] = (real)p_f->id[i];
691 data->mileage[idx]= k*p_f->mileage[i] + d*p_i->mileage[i];
692 data->r[idx] = k*p_f->r[i] + d*p_i->r[i];
693 data->phi[idx] = k*p_f->phi[i] + d*p_i->phi[i];
694 data->z[idx] = k*p_f->z[i] + d*p_i->z[i];
695 data->rho[idx] = k*p_f->rho[i] + d*p_i->rho[i];
696 data->theta[idx] = k*p_f->theta[i] + d*p_i->theta[i];
697 data->B_r[idx] = k*p_f->B_r[i] + d*p_i->B_r[i];
698 data->B_phi[idx] = k*p_f->B_phi[i] + d*p_i->B_phi[i];
699 data->B_z[idx] = k*p_f->B_z[i] + d*p_i->B_z[i];
700 data->pncrid[idx] =
701 j + data->ntoroidalplots + data->npoloidalplots;
702
703 ipoint++;
704 if(ipoint == data->Npnt) {
705 ipoint = 0;
706 }
707
708 data->mrk_pnt[imrk] = ipoint;
709 data->mrk_recorded[imrk] = p_f->mileage[i];
710 }
711 }
712 }
713 }
714 }
715}
716
730real diag_orb_check_plane_crossing(real fang, real iang, real ang0) {
731
732 real k = 0;
733 /* Check whether nag0 is between iang and fang. Note that this *
734 * implementation works only because iang and fang are cumulative. */
735 if( floor( (fang - ang0)/CONST_2PI ) != floor( (iang - ang0)/CONST_2PI ) ) {
736
737 /* Move iang to interval [0, 2pi] */
738 real a = fmod(iang, CONST_2PI);
739 if(a < 0){
740 a = CONST_2PI + a;
741 }
742
743 a = fabs(ang0 - a);
744 if(a > CONST_PI) {
745 a = CONST_2PI - a;
746 }
747 k = fabs(a / (fang - iang));
748 }
749 return k;
750}
751
764real diag_orb_check_radial_crossing(real frho, real irho, real rho0){
765
766 real k = 0;
767 if((frho <= rho0 && irho > rho0) || (irho <= rho0 && frho > rho0)){
768 k = fabs((irho - rho0) / (frho - irho));
769 }
770 return k;
771}
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:133
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:730
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:546
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:339
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:764
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 Tue Dec 3 2024 11:57:29 for ASCOT5 by doxygen 1.11.0