ASCOT5
Loading...
Searching...
No Matches
dist_rho5D.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_rho5D.h"
12#include "../particle.h"
13
17size_t dist_rho5D_index(int i_rho, int i_theta, int i_phi, int i_ppara,
18 int i_pperp, int i_time, int i_q, size_t step_6,
19 size_t step_5, size_t step_4, size_t step_3,
20 size_t step_2, size_t step_1) {
21 return (size_t)(i_rho) * step_6
22 + (size_t)(i_theta) * step_5
23 + (size_t)(i_phi) * 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_phi = (size_t)(data->n_phi);
42 size_t n_theta = (size_t)(data->n_theta);
43 data->step_6 = n_q * n_time * n_pperp * n_ppara * n_phi * n_theta;
44 data->step_5 = n_q * n_time * n_pperp * n_ppara * n_phi;
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_rho, 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_rho*data->n_theta*data->n_phi*data->n_ppara*data->n_pperp*data->n_time*data->n_q]
71 )
72}
73
86 particle_simd_fo* p_i) {
87
88 GPU_PARALLEL_LOOP_ALL_LEVELS
89 for(int i = 0; i < p_f->n_mrk; i++) {
90 if(p_f->running[i]) {
91
92 int i_rho = floor((p_f->rho[i] - dist->min_rho)
93 / ((dist->max_rho - dist->min_rho)/dist->n_rho));
94
95 real phi = fmod(p_f->phi[i], 2*CONST_PI);
96 if(phi < 0) {
97 phi += 2*CONST_PI;
98 }
99 int i_phi = floor((phi - dist->min_phi)
100 / ((dist->max_phi - dist->min_phi)/dist->n_phi));
101
102 real theta = fmod(p_f->theta[i], 2*CONST_PI);
103 if(theta < 0) {
104 theta += 2*CONST_PI;
105 }
106 int i_theta = floor((theta - dist->min_theta)
107 / ( (dist->max_theta - dist->min_theta)
108 / dist->n_theta) );
109
110 real ppara = ( p_f->p_r[i] * p_f->B_r[i]
111 + p_f->p_phi[i] * p_f->B_phi[i]
112 + p_f->p_z[i] * p_f->B_z[i])
113 / sqrt( p_f->B_r[i] * p_f->B_r[i]
114 + p_f->B_phi[i]* p_f->B_phi[i]
115 + p_f->B_z[i] * p_f->B_z[i]);
116 int i_ppara = floor((ppara - dist->min_ppara)
117 / ((dist->max_ppara - dist->min_ppara)
118 / dist->n_ppara));
119
120 real pperp = sqrt(
121 p_f->p_r[i] * p_f->p_r[i]
122 + p_f->p_phi[i] * p_f->p_phi[i]
123 + p_f->p_z[i] * p_f->p_z[i]
124 - ppara * ppara);
125 int i_pperp = floor((pperp - dist->min_pperp)
126 / ((dist->max_pperp - dist->min_pperp)
127 / 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_rho >= 0 && i_rho <= dist->n_rho - 1 &&
136 i_phi >= 0 && i_phi <= dist->n_phi - 1 &&
137 i_theta >= 0 && i_theta <= dist->n_theta - 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 real weight = p_f->weight[i] * (p_f->time[i] - p_i->time[i]);
143 size_t index = dist_rho5D_index(
144 i_rho, i_theta, i_phi, i_ppara, i_pperp,
145 i_time, i_q, dist->step_6, dist->step_5, dist->step_4,
146 dist->step_3, dist->step_2, dist->step_1);
147 GPU_ATOMIC
148 dist->histogram[index] += weight;
149 }
150 }
151 }
152}
153
166 particle_simd_gc* p_i) {
167 real phi[NSIMD];
168 real theta[NSIMD];
169 real pperp[NSIMD];
170
171 int i_rho[NSIMD];
172 int i_phi[NSIMD];
173 int i_theta[NSIMD];
174 int i_ppara[NSIMD];
175 int i_pperp[NSIMD];
176 int i_time[NSIMD];
177 int i_q[NSIMD];
178
179 int ok[NSIMD];
180 real weight[NSIMD];
181
182 #pragma omp simd
183 for(int i = 0; i < NSIMD; i++) {
184 if(p_f->running[i]) {
185
186 i_rho[i] = floor((p_f->rho[i] - dist->min_rho)
187 / ((dist->max_rho - dist->min_rho)/dist->n_rho));
188
189 phi[i] = fmod(p_f->phi[i], 2*CONST_PI);
190 if(phi[i] < 0) {
191 phi[i] = phi[i] + 2*CONST_PI;
192 }
193 i_phi[i] = floor((phi[i] - dist->min_phi)
194 / ((dist->max_phi - dist->min_phi)/dist->n_phi));
195
196 theta[i] = fmod(p_f->theta[i], 2*CONST_PI);
197 if(theta[i] < 0) {
198 theta[i] = theta[i] + 2*CONST_PI;
199 }
200 i_theta[i] = floor((theta[i] - dist->min_theta)
201 / ((dist->max_theta - dist->min_theta)
202 / dist->n_theta));
203
204 i_ppara[i] = floor((p_f->ppar[i] - dist->min_ppara)
205 / ((dist->max_ppara - dist->min_ppara) / dist->n_ppara));
206
207 pperp[i] = sqrt(2 * sqrt( p_f->B_r[i] * p_f->B_r[i]
208 + p_f->B_phi[i] * p_f->B_phi[i]
209 + p_f->B_z[i] * p_f->B_z[i] )
210 * p_f->mu[i] * p_f->mass[i]);
211 i_pperp[i] = floor((pperp[i] - dist->min_pperp)
212 / ((dist->max_pperp - dist->min_pperp)
213 / dist->n_pperp));
214
215 i_time[i] = floor((p_f->time[i] - dist->min_time)
216 / ((dist->max_time - dist->min_time) / dist->n_time));
217
218 i_q[i] = floor((p_f->charge[i]/CONST_E - dist->min_q)
219 / ((dist->max_q - dist->min_q) / dist->n_q));
220
221 if(i_rho[i] >= 0 && i_rho[i] <= dist->n_rho - 1 &&
222 i_phi[i] >= 0 && i_phi[i] <= dist->n_phi - 1 &&
223 i_theta[i] >= 0 && i_theta[i] <= dist->n_theta - 1 &&
224 i_ppara[i] >= 0 && i_ppara[i] <= dist->n_ppara - 1 &&
225 i_pperp[i] >= 0 && i_pperp[i] <= dist->n_pperp - 1 &&
226 i_time[i] >= 0 && i_time[i] <= dist->n_time - 1 &&
227 i_q[i] >= 0 && i_q[i] <= dist->n_q - 1 ) {
228 ok[i] = 1;
229 weight[i] = p_f->weight[i] * (p_f->time[i] - p_i->time[i]);
230 }
231 else {
232 ok[i] = 0;
233 }
234 }
235 }
236
237 for(int i = 0; i < NSIMD; i++) {
238 if(p_f->running[i] && ok[i]) {
239 size_t index = dist_rho5D_index(
240 i_rho[i], i_theta[i], i_phi[i], i_ppara[i], i_pperp[i],
241 i_time[i], i_q[i], dist->step_6, dist->step_5, dist->step_4,
242 dist->step_3, dist->step_2, dist->step_1);
243 #pragma omp atomic
244 dist->histogram[index] += weight[i];
245 }
246 }
247}
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:32
int dist_rho5D_init(dist_rho5D_data *data)
Initializes distribution data.
Definition dist_rho5D.c:35
void dist_rho5D_free(dist_rho5D_data *data)
Free the allocated resources.
Definition dist_rho5D.c:59
void dist_rho5D_offload(dist_rho5D_data *data)
Offload data to the accelerator.
Definition dist_rho5D.c:68
void dist_rho5D_update_gc(dist_rho5D_data *dist, particle_simd_gc *p_f, particle_simd_gc *p_i)
Update the histogram from guiding center markers.
Definition dist_rho5D.c:165
size_t dist_rho5D_index(int i_rho, int i_theta, int i_phi, 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)
Internal function calculating the index in the histogram array.
Definition dist_rho5D.c:17
void dist_rho5D_update_fo(dist_rho5D_data *dist, particle_simd_fo *p_f, particle_simd_fo *p_i)
Update the histogram from full-orbit particles.
Definition dist_rho5D.c:85
Header file for dist_rho5D.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_rho5D.h:15
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