ASCOT5
Loading...
Searching...
No Matches
dist_5D.c
Go to the documentation of this file.
1
5#include <stdio.h>
6#include <stdlib.h>
7#include <math.h>
8#include "../ascot5.h"
9#include "../consts.h"
10#include "../physlib.h"
11#include "dist_5D.h"
12#include "../particle.h"
13
17size_t dist_5D_index(int i_r, int i_phi, int i_z, int i_ppara, int i_pperp,
18 int i_time, int i_q, size_t step_6, size_t step_5,
19 size_t step_4, size_t step_3, size_t step_2,
20 size_t step_1) {
21 return (size_t)(i_r) * step_6
22 + (size_t)(i_phi) * step_5
23 + (size_t)(i_z) * step_4
24 + (size_t)(i_ppara) * step_3
25 + (size_t)(i_pperp) * step_2
26 + (size_t)(i_time) * step_1
27 + (size_t)(i_q);
28}
29
36
37 size_t n_q = (size_t)(data->n_q);
38 size_t n_time = (size_t)(data->n_time);
39 size_t n_pperp = (size_t)(data->n_pperp);
40 size_t n_ppara = (size_t)(data->n_ppara);
41 size_t n_z = (size_t)(data->n_z);
42 size_t n_phi = (size_t)(data->n_phi);
43 data->step_6 = n_q * n_time * n_pperp * n_ppara * n_z * n_phi;
44 data->step_5 = n_q * n_time * n_pperp * n_ppara * n_z;
45 data->step_4 = n_q * n_time * n_pperp * n_ppara;
46 data->step_3 = n_q * n_time * n_pperp;
47 data->step_2 = n_q * n_time;
48 data->step_1 = n_q;
49
50 data->histogram = calloc(data->step_6 * (size_t)data->n_r, sizeof(real));
51 return data->histogram == NULL;
52}
53
60 free(data->histogram);
61}
62
69 GPU_MAP_TO_DEVICE(
70 data->histogram[0:data->n_r*data->n_phi*data->n_z*data->n_ppara*data->n_pperp*data->n_time*data->n_q]
71 )
72}
73
86 particle_simd_fo* p_i) {
87
88#ifdef GPU
89 size_t index;
90 real weight;
91#else
92 size_t index[NSIMD];
93 real weight[NSIMD];
94#endif
95
96 GPU_PARALLEL_LOOP_ALL_LEVELS
97 for(int i = 0; i < p_f->n_mrk; i++) {
98 if(p_f->running[i]) {
99 real i_r = floor((p_f->r[i] - dist->min_r)
100 / ((dist->max_r - dist->min_r)/dist->n_r));
101
102 real phi = fmod(p_f->phi[i], 2*CONST_PI);
103 if(phi < 0) {
104 phi += 2*CONST_PI;
105 }
106 int i_phi = floor((phi - dist->min_phi)
107 / ((dist->max_phi - dist->min_phi)/dist->n_phi));
108
109 int i_z = floor((p_f->z[i] - dist->min_z)
110 / ((dist->max_z - dist->min_z) / dist->n_z));
111
112 real ppara = ( p_f->p_r[i] * p_f->B_r[i]
113 + p_f->p_phi[i] * p_f->B_phi[i]
114 + p_f->p_z[i] * p_f->B_z[i])
115 / sqrt( p_f->B_r[i] * p_f->B_r[i]
116 + p_f->B_phi[i]* p_f->B_phi[i]
117 + p_f->B_z[i] * p_f->B_z[i]);
118 int i_ppara = floor((ppara - dist->min_ppara)
119 / ((dist->max_ppara - dist->min_ppara) / dist->n_ppara));
120
121 real pperp = sqrt(
122 p_f->p_r[i] * p_f->p_r[i]
123 + p_f->p_phi[i] * p_f->p_phi[i]
124 + p_f->p_z[i] * p_f->p_z[i]
125 - ppara * ppara);
126 int i_pperp = floor((pperp - dist->min_pperp)
127 / ((dist->max_pperp - dist->min_pperp) / dist->n_pperp));
128
129 int i_time = floor((p_f->time[i] - dist->min_time)
130 / ((dist->max_time - dist->min_time) / dist->n_time));
131
132 int i_q = floor((p_f->charge[i]/CONST_E - dist->min_q)
133 / ((dist->max_q - dist->min_q) / dist->n_q));
134
135 if(i_r >= 0 && i_r <= dist->n_r - 1 &&
136 i_phi >= 0 && i_phi <= dist->n_phi - 1 &&
137 i_z >= 0 && i_z <= dist->n_z - 1 &&
138 i_ppara >= 0 && i_ppara <= dist->n_ppara - 1 &&
139 i_pperp >= 0 && i_pperp <= dist->n_pperp - 1 &&
140 i_time >= 0 && i_time <= dist->n_time - 1 &&
141 i_q >= 0 && i_q <= dist->n_q - 1 ) {
142#ifdef GPU
143 index = dist_5D_index(
144 i_r, i_phi, i_z, i_ppara, i_pperp, i_time,
145 i_q, dist->step_6, dist->step_5, dist->step_4,
146 dist->step_3, dist->step_2, dist->step_1);
147 weight = p_f->weight[i] * (p_f->time[i] - p_i->time[i]);
148 GPU_ATOMIC
149 dist->histogram[index] += weight;
150#else
151 index[i] = dist_5D_index(
152 i_r, i_phi, i_z, i_ppara, i_pperp, i_time,
153 i_q, dist->step_6, dist->step_5, dist->step_4,
154 dist->step_3, dist->step_2, dist->step_1);
155 weight[i] = p_f->weight[i] * (p_f->time[i] - p_i->time[i]);
156#endif
157 }
158 }
159 }
160
161#ifndef GPU
162 for(int i = 0; i < p_f->n_mrk; i++) {
163 if(p_f->running[i] && index[i] >= 0 &&
164 index[i] < dist->step_6 * dist->n_r) {
165 GPU_ATOMIC
166 dist->histogram[index[i]] += weight[i];
167 }
168 }
169#endif
170}
171
184 particle_simd_gc* p_i) {
185 real phi[NSIMD];
186 real pperp[NSIMD];
187
188 int i_r[NSIMD];
189 int i_phi[NSIMD];
190 int i_z[NSIMD];
191 int i_ppara[NSIMD];
192 int i_pperp[NSIMD];
193 int i_time[NSIMD];
194 int i_q[NSIMD];
195
196 int ok[NSIMD];
197 real weight[NSIMD];
198
199 #pragma omp simd
200 for(int i = 0; i < NSIMD; i++) {
201 if(p_f->running[i]) {
202 i_r[i] = floor((p_f->r[i] - dist->min_r)
203 / ((dist->max_r - dist->min_r)/dist->n_r));
204
205 phi[i] = fmod(p_f->phi[i], 2*CONST_PI);
206 if(phi[i] < 0) {
207 phi[i] = phi[i] + 2*CONST_PI;
208 }
209 i_phi[i] = floor((phi[i] - dist->min_phi)
210 / ((dist->max_phi - dist->min_phi)/dist->n_phi));
211
212 i_z[i] = floor((p_f->z[i] - dist->min_z)
213 / ((dist->max_z - dist->min_z) / dist->n_z));
214
215 i_ppara[i] = floor((p_f->ppar[i] - dist->min_ppara)
216 / ((dist->max_ppara - dist->min_ppara) / dist->n_ppara));
217
218 pperp[i] = sqrt(2 * sqrt( p_f->B_r[i] * p_f->B_r[i]
219 + p_f->B_phi[i] * p_f->B_phi[i]
220 + p_f->B_z[i] * p_f->B_z[i] )
221 * p_f->mu[i] * p_f->mass[i]);
222 i_pperp[i] = floor((pperp[i] - dist->min_pperp)
223 / ((dist->max_pperp - dist->min_pperp) / dist->n_pperp));
224
225 i_time[i] = floor((p_f->time[i] - dist->min_time)
226 / ((dist->max_time - dist->min_time) / dist->n_time));
227
228 i_q[i] = floor((p_f->charge[i]/CONST_E - dist->min_q)
229 / ((dist->max_q - dist->min_q) / dist->n_q));
230
231 if(i_r[i] >= 0 && i_r[i] <= dist->n_r - 1 &&
232 i_phi[i] >= 0 && i_phi[i] <= dist->n_phi - 1 &&
233 i_z[i] >= 0 && i_z[i] <= dist->n_z - 1 &&
234 i_ppara[i] >= 0 && i_ppara[i] <= dist->n_ppara - 1 &&
235 i_pperp[i] >= 0 && i_pperp[i] <= dist->n_pperp - 1 &&
236 i_time[i] >= 0 && i_time[i] <= dist->n_time - 1 &&
237 i_q[i] >= 0 && i_q[i] <= dist->n_q - 1 ) {
238 ok[i] = 1;
239 weight[i] = p_f->weight[i] * (p_f->time[i] - p_i->time[i]);
240 }
241 else {
242 ok[i] = 0;
243 }
244 }
245 }
246
247 for(int i = 0; i < NSIMD; i++) {
248 if(p_f->running[i] && ok[i]) {
249 size_t index = dist_5D_index(
250 i_r[i], i_phi[i], i_z[i], i_ppara[i], i_pperp[i], i_time[i],
251 i_q[i], dist->step_6, dist->step_5, dist->step_4,
252 dist->step_3, dist->step_2, dist->step_1);
253 #pragma omp atomic
254 dist->histogram[index] += weight[i];
255 }
256 }
257}
Main header file for ASCOT5.
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_PI
pi
Definition consts.h:11
#define CONST_E
Elementary charge [C].
Definition consts.h:35
void dist_5D_update_gc(dist_5D_data *dist, particle_simd_gc *p_f, particle_simd_gc *p_i)
Update the histogram from guiding center markers.
Definition dist_5D.c:183
size_t dist_5D_index(int i_r, int i_phi, int i_z, int i_ppara, int i_pperp, int i_time, int i_q, size_t step_6, size_t step_5, size_t step_4, size_t step_3, size_t step_2, size_t step_1)
Function for calculating the index in the histogram array.
Definition dist_5D.c:17
void dist_5D_update_fo(dist_5D_data *dist, particle_simd_fo *p_f, particle_simd_fo *p_i)
Update the histogram from full-orbit particles.
Definition dist_5D.c:85
void dist_5D_free(dist_5D_data *data)
Free allocated resources.
Definition dist_5D.c:59
int dist_5D_init(dist_5D_data *data)
Initializes distribution from offload data.
Definition dist_5D.c:35
void dist_5D_offload(dist_5D_data *data)
Offload data to the accelerator.
Definition dist_5D.c:68
Header file for dist_5D.c.
real fmod(real x, real y)
Compute the modulus of two real numbers.
Definition math.c:22
Header file for math.c.
Header file for particle.c.
Methods to evaluate elementary physical quantities.
Histogram parameters.
Definition dist_5D.h:15
real min_time
Definition dist_5D.h:37
real max_r
Definition dist_5D.h:18
real max_pperp
Definition dist_5D.h:34
size_t step_5
Definition dist_5D.h:48
real min_phi
Definition dist_5D.h:21
size_t step_3
Definition dist_5D.h:46
real max_phi
Definition dist_5D.h:22
real min_ppara
Definition dist_5D.h:29
size_t step_1
Definition dist_5D.h:44
real max_z
Definition dist_5D.h:26
real max_time
Definition dist_5D.h:38
size_t step_2
Definition dist_5D.h:45
int n_ppara
Definition dist_5D.h:28
size_t step_6
Definition dist_5D.h:49
int n_pperp
Definition dist_5D.h:32
real min_z
Definition dist_5D.h:25
real max_ppara
Definition dist_5D.h:30
real min_r
Definition dist_5D.h:17
real max_q
Definition dist_5D.h:42
real * histogram
Definition dist_5D.h:51
real min_pperp
Definition dist_5D.h:33
size_t step_4
Definition dist_5D.h:47
real min_q
Definition dist_5D.h:41
Struct representing NSIMD particle markers.
Definition particle.h:210
integer * running
Definition particle.h:252
Struct representing NSIMD guiding center markers.
Definition particle.h:275