ASCOT5
Loading...
Searching...
No Matches
hist.c
Go to the documentation of this file.
1
9#include <stdlib.h>
10#include "../ascot5.h"
11#include "../math.h"
12#include "../consts.h"
13#include "../offload.h"
14#include "../particle.h"
15#include "hist.h"
16
20int hist_init(histogram* data, int dimensions, hist_coordinate* coordinates,
21 real* binmin, real* binmax, size_t* nbin) {
22 data->axes[0].name = R;
23 data->axes[1].name = PHI;
24 data->axes[2].name = Z;
25 data->axes[3].name = RHO;
26 data->axes[4].name = THETA;
27 data->axes[5].name = PPAR;
28 data->axes[6].name = PPERP;
29 data->axes[7].name = PR;
30 data->axes[8].name = PPHI;
31 data->axes[9].name = PZ;
32 data->axes[10].name = EKIN;
33 data->axes[11].name = XI;
34 data->axes[12].name = MU;
35 data->axes[13].name = PTOR;
36 data->axes[14].name = TIME;
37 data->axes[15].name = CHARGE;
38
39 data->nbin = 1;
40 for(int i = HIST_ALLDIM-1; i >= 0; i--) {
41 data->axes[i].n = 0;
42 data->axes[i].min = 0;
43 data->axes[i].max = 1;
44 for(int k = 0; k < dimensions; k++) {
45 if(coordinates[k] == data->axes[i].name) {
46 data->axes[i].min = binmin[k];
47 data->axes[i].max = binmax[k];
48 data->axes[i].n = nbin[k];
49 data->nbin *= nbin[k];
50 }
51 }
52 if(i < HIST_ALLDIM-1) {
53 data->strides[i] = data->axes[i+1].n ? data->axes[i+1].n : 1;
54 if( i < HIST_ALLDIM-2) {
55 data->strides[i] *= data->strides[i+1];
56 }
57 }
58 }
59 data->bins = (real*) calloc(data->nbin, sizeof(real));
60 return 0;
61}
62
66void hist_free(histogram* data) {
67 free(data->bins);
68 data->nbin = 0;
69}
70
75 GPU_MAP_TO_DEVICE(
76 data->axes[0:data->dimensions], data->bins[0:data->ns],
77 )
78}
79
84 particle_simd_fo* p_i) {
85 GPU_PARALLEL_LOOP_ALL_LEVELS
86 for(int i = 0; i < p_f->n_mrk; i++) {
87 hist_axis* axis;
88
89 real phi = fmod(p_f->phi[i], 2*CONST_PI);
90 phi += (phi<0) * 2*CONST_PI;
91
92 real theta = fmod(p_f->theta[i], 2*CONST_PI);
93 theta += (theta<0) * 2*CONST_PI;
94
95 real ppar = ( p_f->p_r[i] * p_f->B_r[i]
96 + p_f->p_phi[i] * p_f->B_phi[i]
97 + p_f->p_z[i] * p_f->B_z[i])
98 / sqrt( p_f->B_r[i] * p_f->B_r[i]
99 + p_f->B_phi[i]* p_f->B_phi[i]
100 + p_f->B_z[i] * p_f->B_z[i]);
101
102 real pperp = sqrt(
103 p_f->p_r[i] * p_f->p_r[i]
104 + p_f->p_phi[i] * p_f->p_phi[i]
105 + p_f->p_z[i] * p_f->p_z[i]
106 - ppar * ppar);
107
108 axis = &hist->axes[15];
109 size_t nq = axis->n;
110 size_t i15 = 0;
111
112 axis = &hist->axes[14];
113 size_t i14 = math_bin_index(
114 p_f->time[i], axis->n, axis->min, axis->max);
115
116 axis = &hist->axes[13];
117 size_t nptor = axis->n;
118 size_t i13 = 0;
119
120 axis = &hist->axes[12];
121 size_t nmu = axis->n;
122 size_t i12 = 0;
123
124 axis = &hist->axes[11];
125 size_t nxi = axis->n;
126 size_t i11 = 0;
127
128 axis = &hist->axes[10];
129 size_t nekin = axis->n;
130 size_t i10 = 0;
131
132 axis = &hist->axes[9];
133 size_t i9 = math_bin_index(p_f->p_z[i], axis->n, axis->min, axis->max);
134
135 axis = &hist->axes[8];
136 size_t i8 = math_bin_index(
137 p_f->p_phi[i], axis->n, axis->min, axis->max);
138
139 axis = &hist->axes[7];
140 size_t i7 = math_bin_index(p_f->p_r[i], axis->n, axis->min, axis->max);
141
142 axis = &hist->axes[6];
143 size_t i6 = math_bin_index(pperp, axis->n, axis->min, axis->max);
144
145 axis = &hist->axes[5];
146 size_t i5 = math_bin_index(ppar, axis->n, axis->min, axis->max);
147
148 axis = &hist->axes[4];
149 size_t i4 = math_bin_index(theta, axis->n, axis->min, axis->max);
150
151 axis = &hist->axes[3];
152 size_t i3 = math_bin_index(p_f->rho[i], axis->n, axis->min, axis->max);
153
154 axis = &hist->axes[2];
155 size_t i2 = math_bin_index(p_f->z[i], axis->n, axis->min, axis->max);
156
157 axis = &hist->axes[1];
158 size_t i1 = math_bin_index(phi, axis->n, axis->min, axis->max);
159
160 axis = &hist->axes[0];
161 size_t i0 = math_bin_index(p_f->r[i], axis->n, axis->min, axis->max);
162
163 size_t* n = hist->strides;
164 size_t index = i0*n[0] + i1*n[1] + i2*n[2] + i3*n[3] + i4*n[4]
165 + i5*n[5] + i6*n[6] + i7*n[7] + i8*n[8] + i9*n[9]
166 + i10*n[10] + i11*n[11] + i12*n[12] + i13*n[13] + i14*n[14]
167 + i15;
168 GPU_ATOMIC
169 hist->bins[index] += p_f->weight[i];
170 }
171}
172
177 particle_simd_gc* p_i) {
178 //TODO
179}
Main header file for ASCOT5.
double real
Definition ascot5.h:85
Header file containing physical and mathematical constants.
#define CONST_PI
pi
Definition consts.h:11
void hist_update_fo(histogram *hist, particle_simd_fo *p_f, particle_simd_fo *p_i)
Update the histogram in the particle picture.
Definition hist.c:83
void hist_offload(histogram *data)
Offload the data to the accelerator.
Definition hist.c:74
int hist_init(histogram *data, int dimensions, hist_coordinate *coordinates, real *binmin, real *binmax, size_t *nbin)
Initialize the histogram object.
Definition hist.c:20
void hist_free(histogram *data)
Free allocated resources.
Definition hist.c:66
void hist_update_gc(histogram *hist, particle_simd_gc *p_f, particle_simd_gc *p_i)
Update the histogram in the GC picture.
Definition hist.c:176
Header file for hist.c.
hist_coordinate
Quantities that can be used as histogram axis coordinates.
Definition hist.h:17
@ R
Definition hist.h:18
@ XI
Definition hist.h:29
@ THETA
Definition hist.h:22
@ PZ
Definition hist.h:27
@ PPHI
Definition hist.h:26
@ RHO
Definition hist.h:21
@ PHI
Definition hist.h:19
@ PPERP
Definition hist.h:24
@ PR
Definition hist.h:25
@ MU
Definition hist.h:30
@ PPAR
Definition hist.h:23
@ Z
Definition hist.h:20
@ EKIN
Definition hist.h:28
@ PTOR
Definition hist.h:31
@ TIME
Definition hist.h:32
@ CHARGE
Definition hist.h:33
real fmod(real x, real y)
Compute the modulus of two real numbers.
Definition math.c:22
Header file for math.c.
#define math_bin_index(x, nx, xmin, xmax)
Find the bin index on a uniform grid.
Definition math.h:21
Header file for particle.c.
Coordinate axis for the histogram.
Definition hist.h:39
real min
Definition hist.h:41
real max
Definition hist.h:42
hist_coordinate name
Definition hist.h:40
size_t n
Definition hist.h:43
Histogram parameters.
Definition hist.h:52
hist_axis axes[HIST_ALLDIM]
Definition hist.h:53
size_t nbin
Definition hist.h:55
size_t strides[HIST_ALLDIM-1]
Definition hist.h:54
real * bins
Definition hist.h:56
Struct representing NSIMD particle markers.
Definition particle.h:210
Struct representing NSIMD guiding center markers.
Definition particle.h:275