ASCOT5
Loading...
Searching...
No Matches
simulate.c
Go to the documentation of this file.
1
14#include <string.h>
15#include <unistd.h>
16#include "endcond.h"
17#include "particle.h"
18#include "plasma.h"
19#include "wall.h"
20#include "boozer.h"
21#include "mhd.h"
22#include "neutral.h"
23#include "B_field.h"
24#include "E_field.h"
25#include "random.h"
26#include "simulate.h"
27#include "print.h"
32#include "simulate/mccc/mccc.h"
33#include "gctransform.h"
34#include "asigma.h"
35#include "rfof.h"
36
37void sim_monitor(char* filename, volatile int* n, volatile int* finished);
38
84void simulate(int n_particles, particle_state* p, sim_data* sim) {
85
86 // Size = NSIMD on CPU and Size = Total number of particles on GPU
87 int n_queue_size;
88#ifdef GPU
89 n_queue_size = n_particles;
90#else
91 n_queue_size = NSIMD;
92#endif
93 /**************************************************************************/
94 /* 1. Input offload data is unpacked and initialized by calling */
95 /* respective init functions. */
96 /* */
97 /**************************************************************************/
98
99 simulate_init(sim);
100
101#ifdef GPU
102 if(sim->sim_mode != 1) {
103 print_err("Only GO mode ported to GPU. Please set SIM_MODE=1.");
104 exit(1);
105 }
106 if(sim->record_mode) {
107 print_err("RECORD_MODE=1 not ported to GPU. Please disable it.");
108 exit(1);
109 }
110 if(sim->enable_atomic) {
111 print_err("Atomic not yet ported to GPU. Please set ENABLE_ATOMIC=0.");
112 exit(1);
113 }
114 if(sim->enable_mhd) {
115 print_err("MHD not yet ported to GPU. Please set ENABLE_MHD=0.");
116 exit(1);
117 }
118 if(sim->diag_data.diagorb_collect) {
119 print_err(
120 "ENABLE_ORBITWRITE=1 not ported to GPU. Please disable it.");
121 exit(1);
122 }
124 print_err(
125 "ENABLE_TRANSCOEF=1 not ported to GPU. Please disable it.");
126 exit(1);
127 }
128#endif
129
130 diag_init(&sim->diag_data, n_particles);
131 GPU_MAP_TO_DEVICE(sim[0:1])
132 B_field_offload(&sim->B_data);
133 E_field_offload(&sim->E_data);
136 wall_offload(&sim->wall_data);
138 mhd_offload(&sim->mhd_data);
140 diag_offload(&sim->diag_data);
141
142 /**************************************************************************/
143 /* 2. Meta data (e.g. random number generator) is initialized. */
144 /* */
145 /**************************************************************************/
146 random_init(&sim->random_data, 0);
147
148 /**************************************************************************/
149 /* 3. Markers are put into simulation queue. */
150 /* */
151 /**************************************************************************/
153
154 pq.n = 0;
155 for(int i = 0; i < n_particles; i++) {
156 pq.n++;
157 }
158
159 pq.p = (particle_state**) malloc(pq.n * sizeof(particle_state*));
160 pq.finished = 0;
161
162 pq.next = 0;
163 for(int i = 0; i < n_particles; i++) {
164 pq.p[pq.next++] = &p[i];
165
166 }
167 pq.next = 0;
168
169 print_out(VERBOSE_NORMAL, "Simulation begins; %d threads.\n",
170 omp_get_max_threads());
171 fflush(stdout);
172
173 /**************************************************************************/
174 /* 4. Threads are spawned. One thread is dedicated for monitoring */
175 /* progress, if monitoring is active. */
176 /* */
177 /**************************************************************************/
178#ifndef GPU
179 omp_set_max_active_levels(2);
180#endif
181#if !defined(GPU) && VERBOSE > 1
182 #pragma omp parallel sections num_threads(2)
183 {
184 #pragma omp section
185#endif
186 {
187 /******************************************************************/
188 /* 5. Other threads execute marker simulation using the mode the */
189 /* user has chosen. */
190 /* */
191 /******************************************************************/
192 if(pq.n > 0 && (sim->sim_mode == simulate_mode_gc
193 || sim->sim_mode == simulate_mode_hybrid)) {
194 if(sim->enable_ada) {
195 OMP_PARALLEL_CPU_ONLY
196 simulate_gc_adaptive(&pq, sim);
197 }
198 else {
199 OMP_PARALLEL_CPU_ONLY
200 simulate_gc_fixed(&pq, sim);
201 }
202 }
203 else if(pq.n > 0 && sim->sim_mode == simulate_mode_fo) {
204 OMP_PARALLEL_CPU_ONLY
205 simulate_fo_fixed(&pq, sim, n_queue_size);
206 }
207 else if(pq.n > 0 && sim->sim_mode == simulate_mode_ml) {
208 OMP_PARALLEL_CPU_ONLY
209 simulate_ml_adaptive(&pq, sim);
210 }
211 }
212#if !defined(GPU) && VERBOSE > 1
213 #pragma omp section
214 {
215 /* Update progress until simulation is complete. */
216 /* Trim .h5 from filename and replace it with _<QID>.stdout */
217 if(sim->mpi_rank == sim->mpi_root) {
218 char filename[519], outfn[256];
219 strcpy(outfn, sim->hdf5_out);
220 outfn[strlen(outfn)-3] = '\0';
221 sprintf(filename, "%s_%s.stdout", outfn, sim->qid);
222 sim_monitor(filename, &pq.n, &pq.finished);
223 }
224 }
225 }
226#endif
227
228 /**************************************************************************/
229 /* 6. (If hybrid mode is active) Markers with hybrid end condition active */
230 /* are placed on a new queue, and they have their end condition */
231 /* deactivated and they are simulated with simulate_fo_fixed.c until */
232 /* they have met some other end condition. Threads are spawned and */
233 /* progress is monitored as previously. */
234 /* */
235 /**************************************************************************/
236 int n_new = 0;
237 if(sim->sim_mode == simulate_mode_hybrid) {
238
239 /* Determine the number markers that should be run
240 * in fo after previous gc simulation */
241 for(int i = 0; i < pq.n; i++) {
242 if(pq.p[i]->endcond == endcond_hybrid) {
243 /* Check that there was no wall between when moving from
244 gc to fo */
245 real w_coll;
246 int tile = wall_hit_wall(pq.p[i]->r, pq.p[i]->phi, pq.p[i]->z,
247 pq.p[i]->rprt, pq.p[i]->phiprt, pq.p[i]->zprt,
248 &sim->wall_data, &w_coll);
249 if(tile > 0) {
250 pq.p[i]->walltile = tile;
251 pq.p[i]->endcond |= endcond_wall;
252 }
253 else {
254 n_new++;
255 }
256 }
257 }
258 }
259 if(n_new > 0) {
260
261 /* Reset hybrid marker end condition */
262 for(int i = 0; i < pq.n; i++) {
263 if(pq.p[i]->endcond & endcond_hybrid) {
264 pq.p[i]->endcond ^= endcond_hybrid;
265 }
266 }
267 pq.next = 0;
268 pq.finished = 0;
269
270#if !defined(GPU) && VERBOSE > 1
271 #pragma omp parallel sections num_threads(2)
272 {
273 #pragma omp section
274#endif
275 {
276 OMP_PARALLEL_CPU_ONLY
277 simulate_fo_fixed(&pq, sim, n_queue_size);
278 }
279#if !defined(GPU) && VERBOSE > 1
280 #pragma omp section
281 {
282 /* Trim .h5 from filename and replace it with _<qid>.stdout */
283 if(sim->mpi_rank == sim->mpi_root) {
284 char filename[519], outfn[256];
285 strcpy(outfn, sim->hdf5_out);
286 outfn[strlen(outfn)-3] = '\0';
287 sprintf(filename, "%s_%s.stdout", outfn, sim->qid);
288 sim_monitor(filename, &pq.n, &pq.finished);
289 }
290 }
291 }
292#endif
293 }
294
295 /**************************************************************************/
296 /* 7. Simulation data is deallocated. */
297 /**************************************************************************/
298 free(pq.p);
299
300 /**************************************************************************/
301 /* 8. Execution returns to host where this function was called. */
302 /* */
303 /**************************************************************************/
304 print_out(VERBOSE_NORMAL, "Simulation complete.\n");
305}
306
313
316
317 if(sim->disable_gctransform) {
319 }
321}
322
339void sim_monitor(char* filename, volatile int* n, volatile int* finished) {
340 /* Open a file for writing simulation progress */
341 FILE *f = fopen(filename, "w");
342 if (f == NULL) {
344 "Warning. %s could not be opened for progress updates.\n",
345 filename);
346 return;
347 }
348
349 real time_sim_started = A5_WTIME;
350 int stopflag = 1; /* Ensures progress is written one last time at 100% */
351 int n_temp, finished_temp; /* Use these to store volatile variables so that
352 their value does not change during one loop */
353 while(stopflag) {
354 n_temp = *n;
355 finished_temp = *finished;
356 real fracprog = ((real) finished_temp)/n_temp;
357 real timespent = (A5_WTIME)-time_sim_started;
358
359 if(n_temp == finished_temp) {
360 stopflag = 0;
361 }
362
363 if(fracprog == 0) {
364 fprintf(f, "No marker has finished simulation yet. "
365 "Time spent: %.2f h\n", timespent/3600);
366 }
367 else {
368 fprintf(f, "Progress: %d/%d, %.2f %%. Time spent: %.2f h, "
369 "estimated time to finish: %.2f h\n", finished_temp, n_temp,
370 100*fracprog, timespent/3600,
371 (1/fracprog-1)*timespent/3600);
372 }
373 fflush(f);
374 //sleep(A5_PRINTPROGRESSINTERVAL);
375 }
376
377 fprintf(f, "Simulation finished.\n");
378 fclose(f);
379}
void B_field_offload(B_field_data *data)
Offload data to the accelerator.
Definition B_field.c:61
Header file for B_field.c.
void E_field_offload(E_field_data *data)
Offload data to the accelerator.
Definition E_field.c:46
Header file for E_field.c.
double real
Definition ascot5.h:85
#define NSIMD
Number of particles simulated simultaneously in a particle group operations.
Definition ascot5.h:91
#define A5_WTIME
Wall time.
Definition ascot5.h:124
void asigma_extrapolate(int extrapolate)
Toggle extrapolation when evaluating cross sections.
Definition asigma.c:45
void asigma_offload(asigma_data *data)
Offload data to the accelerator.
Definition asigma.c:67
Header file for asigma.c.
void boozer_offload(boozer_data *data)
Offload data to the accelerator.
Definition boozer.c:91
Header file for boozer.c.
void diag_offload(diag_data *data)
Offload data to the accelerator.
Definition diag.c:131
int diag_init(diag_data *data, int Nmrk)
Initializes diagnostics data.
Definition diag.c:37
Header file for endcond.c.
@ endcond_wall
Definition endcond.h:22
@ endcond_hybrid
Definition endcond.h:28
void gctransform_setorder(int order)
Set the order of the transformation.
Definition gctransform.c:62
Header file for gctransform.c.
void mccc_init(mccc_data *mdata, int include_energy, int include_pitch, int include_gcdiff)
Set collision operator data.
Definition mccc.c:17
Header file for mccc package.
void mhd_offload(mhd_data *data)
Offload data to the accelerator.
Definition mhd.c:50
Header file for mhd.c.
void neutral_offload(neutral_data *data)
Offload data to the accelerator.
Definition neutral.c:45
Header file for neutral.c.
Header file for particle.c.
void plasma_offload(plasma_data *data)
Offload data to the accelerator.
Definition plasma.c:57
Header file for plasma.c.
Macros for printing console output.
#define print_out(v,...)
Print to standard output.
Definition print.h:31
@ VERBOSE_NORMAL
Definition print.h:18
@ VERBOSE_DEBUG
Definition print.h:17
#define print_err(...)
Print to standard error.
Definition print.h:42
Header file for random.c.
#define random_init(data, seed)
Definition random.h:94
Contains the functions to be called from the simulation loop when using ICRH.
void simulate(int n_particles, particle_state *p, sim_data *sim)
Execute marker simulation.
Definition simulate.c:84
void simulate_init(sim_data *sim)
Initialize simulation data struct.
Definition simulate.c:312
void sim_monitor(char *filename, volatile int *n, volatile int *finished)
Monitor simulation progress.
Definition simulate.c:339
Header file for simulate.c.
@ simulate_mode_fo
Definition simulate.h:34
@ simulate_mode_ml
Definition simulate.h:45
@ simulate_mode_gc
Definition simulate.h:37
@ simulate_mode_hybrid
Definition simulate.h:42
void simulate_fo_fixed(particle_queue *pq, sim_data *sim, int mrk_array_size)
Simulates particles using fixed time-step.
Header file for simulate_fo_fixed.c.
void simulate_gc_adaptive(particle_queue *pq, sim_data *sim)
Simulates guiding centers using adaptive time-step.
Header file for simulate_gc_adaptive.c.
void simulate_gc_fixed(particle_queue *pq, sim_data *sim)
Simulates guiding centers using fixed time-step.
Header file for simulate_gc_fixed.c.
void simulate_ml_adaptive(particle_queue *pq, sim_data *sim)
Simulates magnetic field-lines using adaptive time-step.
Header file for simulate_ml_adaptive.c.
int diagtrcof_collect
Definition diag.h:28
int diagorb_collect
Definition diag.h:22
Marker queue.
Definition particle.h:154
particle_state ** p
Definition particle.h:156
volatile int next
Definition particle.h:158
volatile int finished
Definition particle.h:159
General representation of a marker.
Definition particle.h:40
integer walltile
Definition particle.h:65
integer endcond
Definition particle.h:64
Simulation data struct.
Definition simulate.h:58
char qid[256]
Definition simulate.h:132
int record_mode
Definition simulate.h:80
int sim_mode
Definition simulate.h:78
int disable_energyccoll
Definition simulate.h:107
plasma_data plasma_data
Definition simulate.h:62
mhd_data mhd_data
Definition simulate.h:66
int enable_atomic
Definition simulate.h:102
int disable_gcdiffccoll
Definition simulate.h:111
E_field_data E_data
Definition simulate.h:61
int mpi_rank
Definition simulate.h:136
int disable_pitchccoll
Definition simulate.h:109
int enable_mhd
Definition simulate.h:101
mccc_data mccc_data
Definition simulate.h:75
random_data random_data
Definition simulate.h:74
neutral_data neutral_data
Definition simulate.h:63
boozer_data boozer_data
Definition simulate.h:65
int enable_ada
Definition simulate.h:79
B_field_data B_data
Definition simulate.h:60
char hdf5_out[256]
Definition simulate.h:131
wall_data wall_data
Definition simulate.h:64
int mpi_root
Definition simulate.h:135
asigma_data asigma_data
Definition simulate.h:67
int disable_gctransform
Definition simulate.h:105
diag_data diag_data
Definition simulate.h:69
void wall_offload(wall_data *data)
Offload data to the accelerator.
Definition wall.c:46
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:80
Header file for wall.c.