ASCOT5
Loading...
Searching...
No Matches
endcond.c
Go to the documentation of this file.
1
53#include <math.h>
54#include "endcond.h"
55#include "particle.h"
56#include "simulate.h"
57#include "physlib.h"
58#include "consts.h"
59#include "math.h"
60#include "plasma.h"
61
74 sim_data* sim) {
75
76 /* Note which end conditions are set as active.
77 Only these ones are checked */
78 int active_tlim = sim->endcond_active & endcond_tlim;
79 int active_wall = sim->endcond_active & endcond_wall;
80 int active_emin = sim->endcond_active & endcond_emin;
81 int active_therm = sim->endcond_active & endcond_therm;
82 int active_rhomax = sim->endcond_active & endcond_rhomax;
83 int active_rhomin = sim->endcond_active & endcond_rhomin;
84 int active_polmax = sim->endcond_active & endcond_polmax;
85 int active_tormax = sim->endcond_active & endcond_tormax;
86 int active_cpumax = sim->endcond_active & endcond_cpumax;
87 int active_neutr = sim->endcond_active & endcond_neutr;
88 int active_ioniz = sim->endcond_active & endcond_ioniz;
89
90 GPU_PARALLEL_LOOP_ALL_LEVELS
91 for(int i = 0; i < p_f->n_mrk; i++) {
92 if(p_f->running[i]) {
93
94 /* Update bounces if pitch changed sign */
95 if( (p_i->p_r[i]*p_i->B_r[i] + p_i->p_phi[i]*p_i->B_phi[i]
96 + p_i->p_z[i]*p_i->B_z[i])
97 * (p_f->p_r[i]*p_f->B_r[i] + p_f->p_phi[i]*p_f->B_phi[i]
98 + p_f->p_z[i]*p_f->B_z[i]) < 0 ) {
99 if(p_f->bounces[i] > 0) {
100 /* Half bounce */
101 p_f->bounces[i] *= -1;
102 }
103 else if(p_f->bounces[i] < 0) {
104 /* Bounce complete */
105 p_f->bounces[i] *= -1;
106 p_f->bounces[i] += 1;
107 }
108 else {
109 /* Initial bounce */
110 p_f->bounces[i] += 1;
111 }
112 }
113
114 /* Check if the marker time exceeds simulation time */
115 if(active_tlim) {
116 if(!sim->reverse_time && p_f->time[i] > sim->endcond_lim_simtime) {
117 p_f->endcond[i] |= endcond_tlim;
118 p_f->running[i] = 0;
119 }
120 if(sim->reverse_time && p_f->time[i] < sim->endcond_lim_simtime) {
121 p_f->endcond[i] |= endcond_tlim;
122 p_f->running[i] = 0;
123 }
124 if(p_f->mileage[i] > sim->endcond_max_mileage) {
125 p_f->endcond[i] |= endcond_tlim;
126 p_f->running[i] = 0;
127 }
128 }
129
130 /* Check, using the wall collision module, whether marker hit wall
131 * during this time-step. Store the wall element ID if it did. */
132 if(active_wall) {
133 real w_coll = 0;
134 int tile = wall_hit_wall(
135 p_i->r[i], p_i->phi[i], p_i->z[i],
136 p_f->r[i], p_f->phi[i], p_f->z[i], &sim->wall_data, &w_coll);
137 if(tile > 0) {
138 real w = w_coll;
139 p_f->time[i] = p_i->time[i] + w*(p_f->time[i] - p_i->time[i]);
140 p_f->r[i] = p_i->r[i] + w*(p_f->r[i] - p_i->r[i]);
141 p_f->phi[i] = p_i->phi[i] + w*(p_f->phi[i] - p_i->phi[i]);
142 p_f->z[i] = p_i->z[i] + w*(p_f->z[i] - p_i->z[i]);
143
144 p_f->walltile[i] = tile;
145 p_f->endcond[i] |= endcond_wall;
146 p_f->running[i] = 0;
147 }
148 }
149
150 /* Evaluate marker energy, and check if it is below the minimum
151 * energy limit or local thermal energy limit */
152 if(active_emin || active_therm) {
153 real pnorm = math_normc(
154 p_f->p_r[i], p_f->p_phi[i], p_f->p_z[i]);
155 real ekin = physlib_Ekin_pnorm(p_f->mass[i], pnorm);
156
157 real Ti;
158 a5err errflag =
159 plasma_eval_temp(&Ti, p_f->rho[i], p_f->r[i], p_f->phi[i],
160 p_f->z[i], p_f->time[i], 1,
161 &sim->plasma_data);
162
163 /* Error handling */
164 if(errflag) {
165 p_f->err[i] = errflag;
166 p_f->running[i] = 0;
167 Ti = 0;
168 }
169
170 if( active_emin && (ekin < sim->endcond_min_ekin) ) {
171 p_f->endcond[i] |= endcond_emin;
172 p_f->running[i] = 0;
173 }
174 if( active_therm && (ekin < (sim->endcond_min_thermal * Ti)) ) {
175 p_f->endcond[i] |= endcond_therm;
176 p_f->running[i] = 0;
177 }
178 }
179
180 /* Check if marker is not within the rho limits */
181 if(active_rhomax) {
182 if(p_f->rho[i] > sim->endcond_max_rho) {
183 p_f->endcond[i] |= endcond_rhomax;
184 p_f->running[i] = 0;
185 }
186 }
187 if(active_rhomin) {
188 if(p_f->rho[i] < sim->endcond_min_rho) {
189 p_f->endcond[i] |= endcond_rhomin;
190 p_f->running[i] = 0;
191 }
192 }
193
194 /* Check if marker exceeds toroidal or poloidal limits */
195 int maxorb = 0;
196 if(active_tormax) {
197 if(fabs(p_f->phi[i]) > sim->endcond_max_tororb) {
198 maxorb |= endcond_tormax;
199 }
200 }
201 if(active_polmax) {
202 if(fabs(p_f->theta[i]) > sim->endcond_max_polorb) {
203 maxorb |= endcond_polmax;
204 }
205 else if( p_f->bounces[i] - 1 >=
206 (int)(sim->endcond_max_polorb / CONST_2PI )) {
207 maxorb |= endcond_polmax;
208 }
209 }
210 if( sim->endcond_torandpol &&
211 maxorb & endcond_tormax && maxorb & endcond_polmax ) {
212 p_f->endcond[i] |= maxorb;
213 p_f->running[i] = 0;
214 }
215 else if(!sim->endcond_torandpol && maxorb) {
216 p_f->endcond[i] |= maxorb;
217 p_f->running[i] = 0;
218 }
219 /* Check if the time spent simulating this marker exceeds the
220 * given limit*/
221 if(active_cpumax) {
222 if(p_f->cputime[i] > sim->endcond_max_cputime) {
223 p_f->endcond[i] |= endcond_cpumax;
224 p_f->running[i] = 0;
225 }
226 }
227 /* Check if the particle has been neutralized */
228 if(active_neutr) {
229 if(p_i->charge[i] != 0.0 && p_f->charge[i] == 0.0) {
230 p_f->endcond[i] |= endcond_neutr;
231 p_f->running[i] = 0;
232 }
233 }
234
235 /* Check if the particle has been ionized */
236 if(active_ioniz) {
237 if(p_i->charge[i] == 0.0 && p_f->charge[i] != 0.0) {
238 p_f->endcond[i] |= endcond_ioniz;
239 p_f->running[i] = 0;
240 }
241 }
242
243 /* Zero end condition if error happened in this function */
244 if(p_f->err[i]) {
245 p_f->endcond[i] = 0;
246 }
247 }
248 }
249}
250
263 sim_data* sim) {
264 int i;
265
266 int active_tlim = sim->endcond_active & endcond_tlim;
267 int active_wall = sim->endcond_active & endcond_wall;
268 int active_emin = sim->endcond_active & endcond_emin;
269 int active_therm = sim->endcond_active & endcond_therm;
270 int active_rhomax = sim->endcond_active & endcond_rhomax;
271 int active_rhomin = sim->endcond_active & endcond_rhomin;
272 int active_polmax = sim->endcond_active & endcond_polmax;
273 int active_tormax = sim->endcond_active & endcond_tormax;
274 int active_cpumax = sim->endcond_active & endcond_cpumax;
275
276 #pragma omp simd
277 for(i = 0; i < NSIMD; i++) {
278 if(p_f->running[i]) {
279 /* Update bounces if pitch changed sign */
280 if( p_i->ppar[i] * p_f->ppar[i] < 0 ) {
281 if(p_f->bounces[i] > 0) {
282 /* Half bounce */
283 p_f->bounces[i] *= -1;
284 }
285 else if(p_f->bounces[i] < 0) {
286 /* Bounce complete */
287 p_f->bounces[i] *= -1;
288 p_f->bounces[i] += 1;
289 }
290 else {
291 /* Initial bounce */
292 p_f->bounces[i] += 1;
293 }
294 }
295
296 /* Check if the marker time exceeds simulation time */
297 if(active_tlim) {
298 if(!sim->reverse_time && p_f->time[i] > sim->endcond_lim_simtime) {
299 p_f->endcond[i] |= endcond_tlim;
300 p_f->running[i] = 0;
301 }
302 if(sim->reverse_time && p_f->time[i] < sim->endcond_lim_simtime) {
303 p_f->endcond[i] |= endcond_tlim;
304 p_f->running[i] = 0;
305 }
306 if(p_f->mileage[i] > sim->endcond_max_mileage) {
307 p_f->endcond[i] |= endcond_tlim;
308 p_f->running[i] = 0;
309 }
310 }
311
312 /* Check, using the wall collision module, whether marker hit wall
313 * during this time-step. Store the wall element ID if it did. */
314 if(active_wall) {
315 real w_coll = 0;
316 int tile = wall_hit_wall(p_i->r[i], p_i->phi[i], p_i->z[i],
317 p_f->r[i], p_f->phi[i], p_f->z[i],
318 &sim->wall_data, &w_coll);
319 if(tile > 0) {
320 p_f->walltile[i] = tile;
321 p_f->endcond[i] |= endcond_wall;
322 p_f->running[i] = 0;
323 }
324 }
325
326 /* Evaluate marker energy, and check if it is below the minimum
327 * energy limit or local thermal energy limit */
328 if(active_emin || active_therm) {
329 real Bnorm = math_normc(
330 p_f->B_r[i], p_f->B_phi[i], p_f->B_z[i]);
331 real ekin = physlib_Ekin_ppar(p_f->mass[i], p_f->mu[i],
332 p_f->ppar[i], Bnorm);
333
334
335 real Ti;
336 a5err errflag =
337 plasma_eval_temp(&Ti, p_f->rho[i], p_f->r[i], p_f->phi[i],
338 p_f->z[i], p_f->time[i], 1,
339 &sim->plasma_data);
340
341 /* Error handling */
342 if(errflag) {
343 p_f->err[i] = errflag;
344 p_f->running[i] = 0;
345 Ti = 0;
346 }
347
348 if(active_emin && (ekin < sim->endcond_min_ekin) ) {
349 p_f->endcond[i] |= endcond_emin;
350 p_f->running[i] = 0;
351 }
352 if( active_therm && (ekin < (sim->endcond_min_thermal * Ti)) ) {
353 p_f->endcond[i] |= endcond_therm;
354 p_f->running[i] = 0;
355 }
356 }
357
358 /* Check if marker is not within the rho limits */
359 if(active_rhomax) {
360 if(p_f->rho[i] > sim->endcond_max_rho) {
361 p_f->endcond[i] |= endcond_rhomax;
362 p_f->running[i] = 0;
363 }
364 }
365 if(active_rhomin) {
366 if(p_f->rho[i] < sim->endcond_min_rho) {
367 p_f->endcond[i] |= endcond_rhomin;
368 p_f->running[i] = 0;
369 }
370 }
371
372 /* Check if marker exceeds toroidal or poloidal limits */
373 int maxorb = 0;
374 if(active_tormax) {
375 if(fabs(p_f->phi[i]) > sim->endcond_max_tororb) {
376 maxorb |= endcond_tormax;
377 }
378 }
379 if(active_polmax) {
380 if(fabs(p_f->theta[i]) > sim->endcond_max_polorb) {
381 maxorb |= endcond_polmax;
382 }
383 else if(p_f->bounces[i] - 1 >=
384 (int)(sim->endcond_max_polorb / CONST_2PI)) {
385 maxorb |= endcond_polmax;
386 }
387 }
388 if( sim->endcond_torandpol &&
389 maxorb & endcond_tormax && maxorb & endcond_polmax ) {
390 p_f->endcond[i] |= maxorb;
391 p_f->running[i] = 0;
392 }
393 else if(!sim->endcond_torandpol && maxorb) {
394 p_f->endcond[i] |= maxorb;
395 p_f->running[i] = 0;
396 }
397
398 /* Check if the time spent simulating this marker exceeds the
399 * given limit*/
400 if(active_cpumax) {
401 if(p_f->cputime[i] > sim->endcond_max_cputime) {
402 p_f->endcond[i] |= endcond_cpumax;
403 p_f->running[i] = 0;
404 }
405 }
406
407 /* If hybrid mode is used, check whether this marker meets the hybrid
408 * condition. */
409 if(sim->sim_mode == 3) {
410 if(p_f->rho[i] > sim->endcond_max_rho) {
411 p_f->endcond[i] |= endcond_hybrid;
412 p_f->running[i] = 0;
413 }
414 }
415
416 /* Zero end condition if error happened in this function */
417 if(p_f->err[i]) {
418 p_f->endcond[i] = 0;
419 }
420 }
421 }
422}
423
436 sim_data* sim) {
437 int i;
438
439 int active_tlim = sim->endcond_active & endcond_tlim;
440 int active_wall = sim->endcond_active & endcond_wall;
441 int active_rhomax = sim->endcond_active & endcond_rhomax;
442 int active_rhomin = sim->endcond_active & endcond_rhomin;
443 int active_polmax = sim->endcond_active & endcond_polmax;
444 int active_tormax = sim->endcond_active & endcond_tormax;
445 int active_cpumax = sim->endcond_active & endcond_cpumax;
446
447 #pragma omp simd
448 for(i = 0; i < NSIMD; i++) {
449 if(p_f->running[i]) {
450 /* Check if the marker time exceeds simulation time */
451 if(active_tlim) {
452 if(!sim->reverse_time && p_f->time[i] > sim->endcond_lim_simtime) {
453 p_f->endcond[i] |= endcond_tlim;
454 p_f->running[i] = 0;
455 }
456 if(sim->reverse_time && p_f->time[i] < sim->endcond_lim_simtime) {
457 p_f->endcond[i] |= endcond_tlim;
458 p_f->running[i] = 0;
459 }
460 if(p_f->mileage[i] > sim->endcond_max_mileage) {
461 p_f->endcond[i] |= endcond_tlim;
462 p_f->running[i] = 0;
463 }
464 }
465
466 /* Check, using the wall collision module, whether marker hit wall
467 * during this time-step. Store the wall element ID if it did. */
468 if(active_wall) {
469 real w_coll = 0;
470 int tile = wall_hit_wall(p_i->r[i], p_i->phi[i], p_i->z[i],
471 p_f->r[i], p_f->phi[i], p_f->z[i],
472 &sim->wall_data, &w_coll);
473 if(tile > 0) {
474 p_f->walltile[i] = tile;
475 p_f->endcond[i] |= endcond_wall;
476 p_f->running[i] = 0;
477 }
478 }
479
480 /* Check if marker is not within the rho limits */
481 if(active_rhomax) {
482 if(p_f->rho[i] > sim->endcond_max_rho) {
483 p_f->endcond[i] |= endcond_rhomax;
484 p_f->running[i] = 0;
485 }
486 }
487 if(active_rhomin) {
488 if(p_f->rho[i] < sim->endcond_min_rho) {
489 p_f->endcond[i] |= endcond_rhomin;
490 p_f->running[i] = 0;
491 }
492 }
493
494 /* Check if marker exceeds toroidal or poloidal limits */
495 int maxorb = 0;
496 if(active_tormax) {
497 if(fabs(p_f->phi[i]) > sim->endcond_max_tororb) {
498 maxorb |= endcond_tormax;
499 }
500 }
501 if(active_polmax) {
502 if(fabs(p_f->theta[i]) > sim->endcond_max_polorb) {
503 maxorb |= endcond_polmax;
504 }
505 }
506 if( sim->endcond_torandpol &&
507 maxorb & endcond_tormax && maxorb & endcond_polmax ) {
508 p_f->endcond[i] |= maxorb;
509 p_f->running[i] = 0;
510 }
511 else if(!sim->endcond_torandpol && maxorb) {
512 p_f->endcond[i] |= maxorb;
513 p_f->running[i] = 0;
514 }
515
516 /* Check if the time spent simulating this marker exceeds the
517 * given limit*/
518 if(active_cpumax) {
519 if(p_f->cputime[i] > sim->endcond_max_cputime) {
520 p_f->endcond[i] |= endcond_cpumax;
521 p_f->running[i] = 0;
522 }
523 }
524 }
525 }
526}
527
538void endcond_parse(int endcond, int* endconds) {
539 int i = 0;
540
541 if(endcond & endcond_tlim) {endconds[i++] = 1;};
542 if(endcond & endcond_emin) {endconds[i++] = 2;};
543 if(endcond & endcond_therm) {endconds[i++] = 3;};
544 if(endcond & endcond_wall) {endconds[i++] = 4;};
545 if(endcond & endcond_rhomin) {endconds[i++] = 5;};
546 if(endcond & endcond_rhomax) {endconds[i++] = 6;};
547 if(endcond & endcond_polmax) {endconds[i++] = 7;};
548 if(endcond & endcond_tormax) {endconds[i++] = 8;};
549 if(endcond & endcond_cpumax) {endconds[i++] = 9;};
550 if(endcond & endcond_hybrid) {endconds[i++] = 10;};
551 if(endcond & endcond_neutr) {endconds[i++] = 11;};
552 if(endcond & endcond_ioniz) {endconds[i++] = 12;};
553}
554
564void endcond_parse2str(int endcond, char* str) {
565 int endconds[32];
566 endcond_parse(endcond, endconds);
567
568 switch(endcond) {
569 case 1:
570 sprintf(str, "Sim time limit");
571 break;
572 case 2:
573 sprintf(str, "Min energy");
574 break;
575 case 3:
576 sprintf(str, "Thermalization");
577 break;
578 case 4:
579 sprintf(str, "Wall collision");
580 break;
581 case 5:
582 sprintf(str, "Min rho");
583 break;
584 case 6:
585 sprintf(str, "Max rho");
586 break;
587 case 7:
588 sprintf(str, "Max poloidal orbits");
589 break;
590 case 8:
591 sprintf(str, "Max toroidal orbits");
592 break;
593 case 9:
594 sprintf(str, "CPU time exceeded");
595 break;
596 case 10:
597 sprintf(str, "Hybrid condition");
598 break;
599 case 11:
600 sprintf(str, "Neutralization");
601 break;
602 case 12:
603 sprintf(str, "Ionization");
604 break;
605 }
606}
double real
Definition ascot5.h:85
#define NSIMD
Number of particles simulated simultaneously in a particle group operations.
Definition ascot5.h:91
Header file containing physical and mathematical constants.
#define CONST_2PI
2*pi
Definition consts.h:14
void endcond_parse(int endcond, int *endconds)
Split endcond to an array of end conditions.
Definition endcond.c:538
void endcond_check_gc(particle_simd_gc *p_f, particle_simd_gc *p_i, sim_data *sim)
Check end conditions for GC markers.
Definition endcond.c:262
void endcond_check_fo(particle_simd_fo *p_f, particle_simd_fo *p_i, sim_data *sim)
Check end conditions for FO markers.
Definition endcond.c:73
void endcond_check_ml(particle_simd_ml *p_f, particle_simd_ml *p_i, sim_data *sim)
Check end conditions for ML markers.
Definition endcond.c:435
void endcond_parse2str(int endcond, char *str)
Represent end condition in human-readable format.
Definition endcond.c:564
Header file for endcond.c.
@ endcond_therm
Definition endcond.h:21
@ endcond_tormax
Definition endcond.h:26
@ endcond_emin
Definition endcond.h:20
@ endcond_rhomin
Definition endcond.h:23
@ endcond_tlim
Definition endcond.h:19
@ endcond_neutr
Definition endcond.h:29
@ endcond_cpumax
Definition endcond.h:27
@ endcond_polmax
Definition endcond.h:25
@ endcond_rhomax
Definition endcond.h:24
@ endcond_wall
Definition endcond.h:22
@ endcond_hybrid
Definition endcond.h:28
@ endcond_ioniz
Definition endcond.h:30
unsigned long int a5err
Simulation error flag.
Definition error.h:17
Header file for math.c.
#define math_normc(a1, a2, a3)
Calculate norm of 3D vector from its components a1, a2, a3.
Definition math.h:67
Header file for particle.c.
Methods to evaluate elementary physical quantities.
#define physlib_Ekin_ppar(m, mu, ppar, B)
Evaluate kinetic energy [J] from parallel momentum.
Definition physlib.h:128
#define physlib_Ekin_pnorm(m, p)
Evaluate kinetic energy [J] from momentum norm.
Definition physlib.h:115
a5err plasma_eval_temp(real *temp, real rho, real r, real phi, real z, real t, int species, plasma_data *pls_data)
Evaluate plasma temperature.
Definition plasma.c:184
Header file for plasma.c.
Header file for simulate.c.
Struct representing NSIMD particle markers.
Definition particle.h:210
integer * endcond
Definition particle.h:247
integer * walltile
Definition particle.h:248
integer * running
Definition particle.h:252
Struct representing NSIMD guiding center markers.
Definition particle.h:275
Struct representing NSIMD field line markers.
Definition particle.h:342
Simulation data struct.
Definition simulate.h:154
int endcond_torandpol
Definition simulate.h:220
real endcond_max_rho
Definition simulate.h:214
int sim_mode
Definition simulate.h:173
plasma_data plasma_data
Definition simulate.h:158
real endcond_max_tororb
Definition simulate.h:218
real endcond_min_thermal
Definition simulate.h:216
real endcond_max_mileage
Definition simulate.h:211
int endcond_active
Definition simulate.h:209
real endcond_max_polorb
Definition simulate.h:219
real endcond_max_cputime
Definition simulate.h:212
real endcond_min_rho
Definition simulate.h:213
real endcond_lim_simtime
Definition simulate.h:210
wall_data wall_data
Definition simulate.h:160
int reverse_time
Definition simulate.h:206
int wall_hit_wall(real r1, real phi1, real z1, real r2, real phi2, real z2, wall_data *w, real *w_coll)
Check if a given directed line segment intersects the wall.
Definition wall.c:163