ASCOT5
Loading...
Searching...
No Matches
asigma_loc.c
Go to the documentation of this file.
1
9#include <stdio.h>
10#include <stdlib.h>
11#include "../ascot5.h"
12#include "../print.h"
13#include "../error.h"
14#include "../spline/interp.h"
15#include "../consts.h"
16#include "../math.h"
17#include "../physlib.h"
18#include "../suzuki.h"
19#include "../asigma.h"
20#include "asigma_loc.h"
21
53 real** offload_array) {
54 int err = 0;
55 int N_reac = offload_data->N_reac;
56
57 /* Find how much space is needed for the output array */
58 int temp_arr_length = 6 * N_reac;
59 for(int i_reac = 0; i_reac < N_reac; i_reac++) {
60 int N_E = offload_data->N_E[i_reac];
61 int N_n = offload_data->N_n[i_reac];
62 int N_T = offload_data->N_T[i_reac];
63 int dim = (N_E > 1) + (N_n > 1) + (N_T > 1);
64 switch(dim) {
65 case 1:
66 temp_arr_length += N_E * NSIZE_COMP1D;
67 break;
68
69 case 2:
70 temp_arr_length += N_E * N_T * NSIZE_COMP2D;
71 break;
72
73 case 3:
74 temp_arr_length += N_E * N_n * N_T * NSIZE_COMP3D;
75 break;
76
77 default:
78 /* Unrecognized dimensionality. Produce error. */
79 print_err("Error: Unrecognized abscissa dimensionality\n");
80 err = 1;
81 break;
82 }
83 }
84 real* temp_array = (real*) malloc(temp_arr_length*sizeof(real));
85
86 /* Helper pointers to keep track of the current memory positions in the
87 reaction data parts of the arrays */
88 real * temp_arr_pos = temp_array + 6 * N_reac;
89 real * offload_arr_pos = *offload_array + 6 * N_reac;
90
91 /* Copy over reaction identifiers and abscissae parameters, and evaluate
92 spline coefficients according to dimensionality of reaction data */
93 for(int i_reac = 0; i_reac < N_reac; i_reac++) {
94 int N_E = offload_data->N_E[i_reac];
95 int N_n = offload_data->N_n[i_reac];
96 int N_T = offload_data->N_T[i_reac];
97 int dim = (N_E > 1) + (N_n > 1) + (N_T > 1);
98
99 real E_min = (*offload_array)[0*N_reac+i_reac];
100 real E_max = (*offload_array)[1*N_reac+i_reac];
101 real n_min = (*offload_array)[2*N_reac+i_reac];
102 real n_max = (*offload_array)[3*N_reac+i_reac];
103 real T_min = (*offload_array)[4*N_reac+i_reac];
104 real T_max = (*offload_array)[5*N_reac+i_reac];
105 temp_array[ 0*N_reac+i_reac] = E_min;
106 temp_array[ 1*N_reac+i_reac] = E_max;
107 temp_array[ 2*N_reac+i_reac] = n_min;
108 temp_array[ 3*N_reac+i_reac] = n_max;
109 temp_array[ 4*N_reac+i_reac] = T_min;
110 temp_array[ 5*N_reac+i_reac] = T_max;
111 switch(dim) {
112 case 1:
114 temp_arr_pos, offload_arr_pos,
115 N_E,
116 NATURALBC,
117 E_min, E_max);
118 temp_arr_pos += N_E * NSIZE_COMP1D;
119 offload_arr_pos += N_E;
120 break;
121
122 case 2:
124 temp_arr_pos, offload_arr_pos,
125 N_E, N_T,
127 E_min, E_max,
128 T_min, T_max);
129 temp_arr_pos += N_E * N_T * NSIZE_COMP2D;
130 offload_arr_pos += N_E * N_T;
131 break;
132
133 case 3:
135 temp_arr_pos, offload_arr_pos,
136 N_E, N_n, N_T,
138 E_min, E_max,
139 n_min, n_max,
140 T_min, T_max);
141 temp_arr_pos += N_E * N_n * N_T * NSIZE_COMP3D;
142 offload_arr_pos += N_E * N_n * N_T;
143 break;
144
145 default:
146 /* Unrecognized dimensionality. Produce error. */
147 print_err("Error: Unrecognized abscissa dimensionality\n");
148 err = 1;
149 break;
150 }
151 }
152
153 free(*offload_array);
154 *offload_array = temp_array;
155 offload_data->offload_array_length = temp_arr_length;
156
157 if(err) {
158 return err;
159 }
160
161 print_out(VERBOSE_IO,"\nFound data for %d atomic reactions:\n", N_reac);
162 for(int i_reac = 0; i_reac < N_reac; i_reac++) {
164 "Reaction number / Total number of reactions = %d / %d\n"
165 " Reactant species Z_1 / A_1, Z_2 / A_2 = %d / %d, %d / %d\n"
166 " Min/Max energy = %1.2le / %1.2le\n"
167 " Min/Max density = %1.2le / %1.2le\n"
168 " Min/Max temperature = %1.2le / %1.2le\n"
169 " Number of energy grid points = %d\n"
170 " Number of density grid points = %d\n"
171 " Number of temperature grid points = %d\n",
172 i_reac+1, N_reac,
173 offload_data->z_1[i_reac], offload_data->a_1[i_reac],
174 offload_data->z_2[i_reac], offload_data->a_2[i_reac],
175 (*offload_array)[0 * N_reac + i_reac],
176 (*offload_array)[1 * N_reac + i_reac],
177 (*offload_array)[2 * N_reac + i_reac],
178 (*offload_array)[3 * N_reac + i_reac],
179 (*offload_array)[4 * N_reac + i_reac],
180 (*offload_array)[5 * N_reac + i_reac],
181 offload_data->N_E[i_reac],
182 offload_data->N_n[i_reac],
183 offload_data->N_T[i_reac]);
184 }
185
186 return 0;
187}
188
198 real** offload_array) {
199 free(*offload_array);
200 *offload_array = NULL;
201}
202
218 real* offload_array) {
219 /* Copy over number of reactions and store it in a helper variable */
220 int N_reac = offload_data->N_reac;
221 asigma_data->N_reac = N_reac;
222
223 /* Helper pointer to keep track of position in offload array */
224 real* offload_arr_pos = offload_array + 6 * N_reac;
225
226 /* Copy data from offload array to atomic sigma struct,
227 initialize spline structs and determine reaction availability */
228 for(int i_reac = 0; i_reac < N_reac; i_reac++) {
229
230 /* Reaction identifiers */
231 asigma_data->z_1[i_reac] = offload_data->z_1[i_reac];
232 asigma_data->a_1[i_reac] = offload_data->a_1[i_reac];
233 asigma_data->z_2[i_reac] = offload_data->z_2[i_reac];
234 asigma_data->a_2[i_reac] = offload_data->a_2[i_reac];
235 asigma_data->reac_type[i_reac] = offload_data->reac_type[i_reac];
236
237 /* Initialize spline struct according to dimensionality of
238 reaction data (and mark reaction availability) */
239 int N_E = offload_data->N_E[i_reac];
240 real E_min = offload_array[0*N_reac+i_reac];
241 real E_max = offload_array[1*N_reac+i_reac];
242 int N_n = offload_data->N_n[i_reac];
243 real n_min = offload_array[2*N_reac+i_reac];
244 real n_max = offload_array[3*N_reac+i_reac];
245 int N_T = offload_data->N_T[i_reac];
246 real T_min = offload_array[4*N_reac+i_reac];
247 real T_max = offload_array[5*N_reac+i_reac];
248 int dim = (N_E > 1) + (N_n > 1) + (N_T > 1);
249 switch(dim) {
250 case 1:
252 &(asigma_data->sigma[i_reac]), offload_arr_pos,
253 N_E,
254 NATURALBC,
255 E_min, E_max);
256 offload_arr_pos += N_E *NSIZE_COMP1D;
257 break;
258
259 case 2:
261 &(asigma_data->sigmav[i_reac]), offload_arr_pos,
262 N_E, N_T,
264 E_min, E_max,
265 T_min, T_max);
266 offload_arr_pos += N_E * N_T * NSIZE_COMP2D;
267 break;
268
269 case 3:
271 &(asigma_data->BMSsigmav[i_reac]), offload_arr_pos,
272 N_E, N_n, N_T,
274 E_min, E_max,
275 n_min, n_max,
276 T_min, T_max);
277 offload_arr_pos += N_E * N_n * N_T * NSIZE_COMP3D;
278 break;
279
280 default:
281 /* Unrecognized dimensionality. Produce error. */
282 print_err("Error: Unrecognized abscissa dimensionality\n");
283 break;
284 }
285 }
286}
287
310 real* sigma, int z_1, int a_1, int z_2, int a_2, real E_coll_per_amu,
311 int reac_type, int extrapolate, asigma_loc_data* asigma_data) {
312 a5err err = 0;
313
314 /* We look for a match of the reaction identifiers in asigma_data to
315 determine if the reaction of interest has been initialized */
316 int reac_found = -1, i_reac;
317 for(i_reac = 0; i_reac < asigma_data->N_reac; i_reac++) {
318 if(z_1 == asigma_data->z_1[i_reac] &&
319 a_1 == asigma_data->a_1[i_reac] &&
320 z_2 == asigma_data->z_2[i_reac] &&
321 a_2 == asigma_data->a_2[i_reac] &&
322 reac_type == asigma_data->reac_type[i_reac]) {
323 reac_found = i_reac;
324 }
325 }
326 i_reac = reac_found;
327
328 /* The cross-section is evaluated if reaction data was found,
329 is available, and its interpolation implemented. Otherwise,
330 the cross-section is set to zero to avoid further problems. */
331 if(reac_found < 0) {
332 /* Reaction not found. Raise error. */
334 } else {
335 if(asigma_data->reac_type[i_reac] == sigma_ioniz ||
336 asigma_data->reac_type[i_reac] == sigma_recomb ||
337 asigma_data->reac_type[i_reac] == sigma_CX) {
338 int interperr = 0;
339 interperr += interp1Dcomp_eval_f(sigma,
340 &asigma_data->sigma[i_reac],
341 E_coll_per_amu);
342 if(interperr) {
343 /* Energy is outside spline domain */
344 if(extrapolate) {
345 *sigma = 0.0;
346 } else {
347 err = error_raise( ERR_INPUT_EVALUATION, __LINE__,
349 }
350 }
351 } else {
352 /* Interpolation of cross-section not implemented. Raise error. */
354 }
355 }
356 return err;
357}
358
385 real* sigmav, int z_1, int a_1, real m_1, int z_2, int a_2,
386 real E, real T_e, real T_0, real n_i, int reac_type, int extrapolate,
388 a5err err = 0;
389
390 /* Convert Joule to eV */
391 E /= CONST_E;
392 T_e /= CONST_E;
393 T_0 /= CONST_E;
394
395 /* Find the matching reaction. Note that BMS data is same for all
396 * isotopes, so we don't compare anums */
397 int reac_found = -1, i_reac;
398 for(i_reac = 0; i_reac < asigma_data->N_reac; i_reac++) {
399 if(reac_type == sigmav_BMS &&
400 z_1 == asigma_data->z_1[i_reac] &&
401 z_2 == asigma_data->z_2[i_reac] &&
402 reac_type == asigma_data->reac_type[i_reac]) {
403 reac_found = i_reac;
404 } else if(z_1 == asigma_data->z_1[i_reac] &&
405 a_1 == asigma_data->a_1[i_reac] &&
406 z_2 == asigma_data->z_2[i_reac] &&
407 a_2 == asigma_data->a_2[i_reac] &&
408 reac_type == asigma_data->reac_type[i_reac]) {
409 reac_found = i_reac;
410 }
411 }
412 i_reac = reac_found;
413
414 if(reac_found < 0) {
415 /* Reaction not found. Raise error. */
417 } else {
418 /* Interpolation error means the data has to be extrapolated */
419 if(reac_type == sigmav_ioniz ||
420 reac_type == sigmav_recomb ||
421 reac_type == sigmav_CX) {
422 int interperr = interp2Dcomp_eval_f(
423 sigmav, &asigma_data->sigmav[i_reac], E, T_0);
424 if(interperr) {
425 if(extrapolate) {
426 *sigmav = 0.0;
427 } else {
428 err = error_raise( ERR_INPUT_EVALUATION, __LINE__,
430 }
431 }
432 } else if(reac_type == sigmav_BMS) {
433 int interperr = interp3Dcomp_eval_f(
434 sigmav, &asigma_data->BMSsigmav[i_reac], E/a_2,
435 z_2*n_i, T_e);
436 if(interperr) {
437 if(extrapolate) {
438 *sigmav = 0.0;
439 } else {
440 err = error_raise( ERR_INPUT_EVALUATION, __LINE__,
442 }
443 }
444 } else {
445 /* Interpolation of rate coefficient not implemented.
446 Raise error. */
448 }
449 }
450
451 return err;
452}
453
478 real* ratecoeff, int z_1, int a_1, real E, real mass, int nspec,
479 const int* znum, const int* anum, real T_0, real* n_0, int extrapolate,
481 a5err err = 0;
482
483 /* Convert Joule to eV */
484 E /= CONST_E;
485 T_0 /= CONST_E;
486 *ratecoeff = 0;
487 for(int i_spec = 0; i_spec < nspec; i_spec++) {
488
489 /* Find the matching reaction */
490 int reac_found = -1, i_reac;
491 for(i_reac = 0; i_reac < asigma_data->N_reac; i_reac++) {
492 if(asigma_data->z_1[i_reac] == z_1 &&
493 asigma_data->a_1[i_reac] == a_1 &&
494 asigma_data->z_2[i_reac] == znum[i_spec] &&
495 asigma_data->a_2[i_reac] == anum[i_spec] &&
496 asigma_data->reac_type[i_reac] == sigmav_CX) {
497 reac_found = i_reac;
498 }
499 }
500 i_reac = reac_found;
501
502 if(reac_found < 0) {
503 /* Reaction not found. Raise error. */
505 } else {
506 real sigmav;
507 int interperr = interp2Dcomp_eval_f(
508 &sigmav, &asigma_data->sigmav[i_reac], E, T_0);
509
510 /* Interpolation error means the data has to be extrapolated */
511 if(interperr) {
512 if(extrapolate) {
513 sigmav = 0.0;
514 } else {
515 err = error_raise( ERR_INPUT_EVALUATION, __LINE__,
517 }
518 }
519 *ratecoeff += sigmav*n_0[i_spec];
520 }
521 }
522
523 return err;
524}
525
550 real* ratecoeff, int z_1, int a_1, real E, real mass, int nion,
551 const int* znum, const int* anum, real T_e, real* n_i, int extrapolate,
553 a5err err = 0;
554
555 /* Convert Joule to eV */
556 real E_eV = E / CONST_E;
557 T_e /= CONST_E;
558
559 /* Find the matching reaction. Note that BMS data is same for all
560 * isotopes, so we don't compare anums */
561 int reac_found = -1; real n_e = 0; *ratecoeff = 0;
562 for(int i_spec = 0; i_spec < nion; i_spec++) {
563 n_e += znum[i_spec] * n_i[i_spec];
564 for(int i_reac = 0; i_reac < asigma_data->N_reac; i_reac++) {
565 if(asigma_data->z_1[i_reac] == z_1 &&
566 asigma_data->z_2[i_reac] == znum[i_spec] &&
567 asigma_data->reac_type[i_reac] == sigmav_BMS) {
568 reac_found = i_reac;
569 real sigmav;
570 int interperr = \
571 interp3Dcomp_eval_f(
572 &sigmav, &asigma_data->BMSsigmav[i_reac],
573 E_eV/anum[i_spec], znum[i_spec] * n_i[i_spec], T_e);
574
575 /* Interpolation error means the data has to be extrapolated */
576 if(interperr) {
577 if(extrapolate) {
578 sigmav = 0.0;
579 } else {
580 err = error_raise( ERR_INPUT_EVALUATION, __LINE__,
582 }
583 }
584 *ratecoeff += sigmav * ( znum[i_spec] * n_i[i_spec]);
585 }
586 }
587 }
588 *ratecoeff /= n_e;
589
590 if(reac_found < 0) {
591 /* Reaction not found. Try Suzuki before throwing error. */
592 T_e *= CONST_E;
593 real gamma = physlib_gamma_Ekin(mass, E);
594 real vnorm = physlib_vnorm_gamma(gamma);
595 if(suzuki_sigmav(ratecoeff, E/a_1, vnorm, n_e, T_e, nion, n_i, anum,
596 znum)) {
598 }
599 }
600
601 return err;
602}
Main header file for ASCOT5.
double real
Definition ascot5.h:85
Header file for asigma.c.
void asigma_loc_free_offload(asigma_loc_offload_data *offload_data, real **offload_array)
Free offload array and reset parameters.
Definition asigma_loc.c:197
a5err asigma_loc_eval_cx(real *ratecoeff, int z_1, int a_1, real E, real mass, int nspec, const int *znum, const int *anum, real T_0, real *n_0, int extrapolate, asigma_loc_data *asigma_data)
Evaluate atomic reaction rate coefficient.
Definition asigma_loc.c:477
a5err asigma_loc_eval_sigma(real *sigma, int z_1, int a_1, int z_2, int a_2, real E_coll_per_amu, int reac_type, int extrapolate, asigma_loc_data *asigma_data)
Evaluate atomic reaction cross-section.
Definition asigma_loc.c:309
void asigma_loc_init(asigma_loc_data *asigma_data, asigma_loc_offload_data *offload_data, real *offload_array)
Initialize atomic reaction data struct on target.
Definition asigma_loc.c:216
a5err asigma_loc_eval_sigmav(real *sigmav, int z_1, int a_1, real m_1, int z_2, int a_2, real E, real T_e, real T_0, real n_i, int reac_type, int extrapolate, asigma_loc_data *asigma_data)
Evaluate atomic reaction rate coefficient.
Definition asigma_loc.c:384
int asigma_loc_init_offload(asigma_loc_offload_data *offload_data, real **offload_array)
Initialize local file atomic data and check inputs.
Definition asigma_loc.c:52
a5err asigma_loc_eval_bms(real *ratecoeff, int z_1, int a_1, real E, real mass, int nion, const int *znum, const int *anum, real T_e, real *n_i, int extrapolate, asigma_loc_data *asigma_data)
Evaluate beam stopping rate coefficient.
Definition asigma_loc.c:549
Header file for asigma_loc.c.
Header file containing physical and mathematical constants.
#define CONST_E
Elementary charge [C]
Definition consts.h:32
Error module for ASCOT5.
unsigned long int a5err
Simulation error flag.
Definition error.h:17
@ EF_ASIGMA_LOC
Definition error.h:52
@ ERR_INPUT_EVALUATION
Definition error.h:63
static DECLARE_TARGET_SIMD a5err error_raise(error_type type, int line, error_file file)
Raise a new error.
Definition error.h:86
Spline interpolation library.
DECLARE_TARGET_END a5err interp3Dcomp_eval_f(real *f, interp3D_data *str, real x, real y, real z)
Evaluate interpolated value of 3D scalar field.
int interp3Dcomp_init_coeff(real *c, real *f, int n_x, int n_y, int n_z, int bc_x, int bc_y, int bc_z, real x_min, real x_max, real y_min, real y_max, real z_min, real z_max)
Calculate tricubic spline interpolation coefficients for 3D data.
@ NATURALBC
Definition interp.h:37
int interp1Dcomp_init_coeff(real *c, real *f, int n_x, int bc_x, real x_min, real x_max)
Calculate cubic spline interpolation coefficients for scalar 1D data.
int interp2Dcomp_init_coeff(real *c, real *f, int n_x, int n_y, int bc_x, int bc_y, real x_min, real x_max, real y_min, real y_max)
Calculate bicubic spline interpolation coefficients for scalar 2D data.
a5err interp1Dcomp_eval_f(real *f, interp1D_data *str, real x)
Evaluate interpolated value of 1D scalar field.
void interp3Dcomp_init_spline(interp3D_data *str, real *c, int n_x, int n_y, int n_z, int bc_x, int bc_y, int bc_z, real x_min, real x_max, real y_min, real y_max, real z_min, real z_max)
Initialize a tricubic spline.
void interp2Dcomp_init_spline(interp2D_data *str, real *c, int n_x, int n_y, int bc_x, int bc_y, real x_min, real x_max, real y_min, real y_max)
Initialize a bicubic spline.
void interp1Dcomp_init_spline(interp1D_data *str, real *c, int n_x, int bc_x, real x_min, real x_max)
Initialize a cubic spline.
DECLARE_TARGET_END a5err interp2Dcomp_eval_f(real *f, interp2D_data *str, real x, real y)
Evaluate interpolated value of a 2D field.
Header file for math.c.
Methods to evaluate elementary physical quantities.
#define physlib_vnorm_gamma(gamma)
Evaluate velocity norm from Lorentz factor.
Definition physlib.h:33
#define physlib_gamma_Ekin(m, ekin)
Evaluate Lorentz factor from kinetic energy [J].
Definition physlib.h:103
Macros for printing console output.
#define print_out(v,...)
Print to standard output.
Definition print.h:31
@ VERBOSE_IO
Definition print.h:20
#define print_err(...)
Print to standard error.
Definition print.h:42
Atomic reaction simulation data.
Definition asigma.h:66
Local-files atomic reaction simulation data.
Definition asigma_loc.h:31
Local-files atomic reaction offload data.
Definition asigma_loc.h:15
int a_2[MAX_ATOMIC]
Definition asigma_loc.h:20
int z_1[MAX_ATOMIC]
Definition asigma_loc.h:17
int reac_type[MAX_ATOMIC]
Definition asigma_loc.h:24
int N_n[MAX_ATOMIC]
Definition asigma_loc.h:22
int a_1[MAX_ATOMIC]
Definition asigma_loc.h:18
int N_E[MAX_ATOMIC]
Definition asigma_loc.h:21
int z_2[MAX_ATOMIC]
Definition asigma_loc.h:19
int N_T[MAX_ATOMIC]
Definition asigma_loc.h:23
a5err suzuki_sigmav(real *sigmav, real EperAmu, real vnorm, real ne, real te, integer nion, real *ni, const int *anum, const int *znum)
Calculate beam-stopping cross-section according to Suzuki model.
Definition suzuki.c:139
Header file for suzuki.c.