diff --git a/src/reaction.c b/src/reaction.c new file mode 100644 index 0000000000000000000000000000000000000000..ebc92632a5711b1ff424b77f44350333e1f459bb --- /dev/null +++ b/src/reaction.c @@ -0,0 +1,231 @@ +/** @file reaction.c + * @brief Functions to handle the "reaction" part of + * convection-diffusion-reaction equation. + */ +#include +#include +#include + +#include "cdr.h" +#include "grid.h" +#include "parameters.h" +#include "photo.h" +#include "proto.h" +#include "react_table.h" +#include "species.h" +#include "reaction.h" + +static void fill_react_gaps (); + +reaction_t *reactions_list = NULL; +extern double z_cutoff; +double species_written[20]; + +/** @brief Returns the position in the species-array of a given species. + * + * Returns -1 if not found. */ +int +find_species_by_name(const char *spec_name) +{ + int res = -1; + int cnt; + + for (cnt = 0; cnt < no_species; cnt++) + { + if (strcmp(spec_index[cnt]->name, spec_name) == 0) + res = cnt; + } + + if (res == -1) { + printf("Species-lookup-failure: %s!\n",spec_name); + assert(1 == 0); + } + return res; +} + +/** @brief Adds a given reaction to the reaction list */ +void +react_add (reaction_t *react) +{ + react_table *rt; + + debug (3, "react_add (...)\n"); + + printf("Adding reaction with table: %s\n",react->tablefile); + + rt = (react_table *) xmalloc (sizeof(react_table)); + + if (react->tablefile != NULL) { + react_table_read(react->tablefile, rt); + react->rt = rt; + } + + react->next = reactions_list; + reactions_list = react; +} + +/** @brief Applies a reaction to the given @a grid */ +void +react_apply (reaction_t *react, cdr_grid_t *grid, int overwrite) +{ + int i, ir, iz, itheta; + rz_array_t *grid_in[REACTION_MAX_IN], *grid_out[REACTION_MAX_IN + REACTION_MAX_OUT], + *grid_eabs; + double *in = NULL, *out = NULL, eabs; + double test; + double rate; + int pos; + double e1, e2, val1, val2, log_e, res, r_mod; + int cnt, curr_species; + + debug (3, "react_apply (..., " grid_printf_str ")\n", + grid_printf_args (grid)); + + grid_eabs = grid->dens[no_species]; + + for (i = 0; i < react->nin; i++) { + grid_in[i] = grid->dens[react->input[i]]; + grid_out[i] = grid->d_dens[react->input[i]]; + } + + for (i = react->nin; i < react->nin + react->nout; i++) { + grid_out[i] = grid->d_dens[react->output[i - react->nin]]; + } + +#pragma omp parallel private(ir, iz, i, in, out) + { + /* malloc(0) is legal, but I do not want to play with fire. */ + if (react->nin > 0) { + in = (double *) xmalloc (sizeof(double) * react->nin); + } + + /* Do not know what use nout == 0 may have, (perhaps some debugging?) + but we leave it here as theoretically possible. */ + if (react->nout > 0) { + out = (double *) xmalloc (sizeof(double) * (react->nout + react->nin)); + } + + + #pragma omp for + iter_grid_theta(grid, itheta) { //ITER3 + iter_grid_z(grid, iz) { //ITER2 + double back_dens; + if (sprite_module) { + back_dens = spr_density_at (z_at (iz, grid->level)); + } else { + back_dens = 1.0; + } + + + iter_grid_r(grid, ir) { //ITER1 + if ( z_at (iz, grid->level) < z_cutoff ) { //IF1 + eabs = fabs (RZT(grid_eabs, ir, iz, itheta)); + for (i = 0; i < react->nin; i++) { + in[i] = fabs (RZT(grid_in[i], ir, iz, itheta)); + } + + log_e = log10(eabs); + + /* If the supplied fieldstrength falls outside the boundaries of the table, + return predetermined under-/overflow values */ + if (log_e < react->rt->e_min) { + rate = react->rt->underflow; + } else if (log_e > react->rt->e_min + react->rt->e_step * react->rt->steps) { + rate = react->rt->overflow; + } else { + pos = floor((log_e - react->rt->e_min) / react->rt->e_step); + val1 = react->rt->values[pos]; + val2 = react->rt->values[pos+1]; + e1 = pow(10, react->rt->e_min + react->rt->e_step * (double) pos); + e2 = e1 * pow(10, react->rt->e_step); + + rate = val1 + ((val2 - val1) / (e2 - e1)) * (eabs - e1); + } + + for (i = 0; i < react->nin; i++){ rate *= MYMAX(0, in[i]); } + + for (i = 0; i < react->nin + react->nout; i++) { //FOR1 + if (i < react->nin) { + curr_species = react->input[i]; r_mod = -rate; + } else { curr_species = react->output[i - react->nin]; r_mod = rate; } + + if (spec_index[curr_species]->charge != 0.0) { + RZT(grid_out[i], ir, iz, itheta) += r_mod; + } + } //FOR1 + } //IF1 + } //ITER1 + } //ITER2 + } //ITER3 + free (in); + free (out); + } +} + +/** @brief Recursive version of @a react_apply */ +void +react_apply_r (reaction_t *react, cdr_grid_t *grid, int overwrite) +{ + cdr_grid_t *child; + + react_apply (react, grid, overwrite); + + iter_childs (grid, child) { + react_apply_r (react, child, overwrite); + } +} + +/** @brief Sets the d_dens field of a grid to zero */ +void +zero_fill (cdr_grid_t* grid) +{ + int ir, iz, itheta, i; + + iter_grid_3d (grid, ir, iz, itheta) + for (i = 0; i < no_species; i++) { + RZT(grid->d_dens[i], ir, iz, itheta) = 0.0; + } +} + +/** @brief Applies all reactions to the given @a grid and his descendants */ +void +react_apply_all (cdr_grid_t *grid) +{ + reaction_t *react; + int overwrite; + int last = -1; + int cnt; + + zero_fill(grid); + + overwrite = TRUE; + + for (react = reactions_list; react; react = react->next) { + if (react->is_photo) { + photo_calc (photo_terms, grid); + } else { + react_apply_r (react, grid, overwrite); + } + overwrite = FALSE; + } +} + +/** @brief Fill in the gaps that we left in the definitions */ +static void +fill_react_gaps () +{ +} + +/** @brief Initializes the list of reactions. */ +void +react_init () +{ + /* Note that the reactions are applied in inverse order than listed here. */ + /* Rest of the kinetic model: */ + kinetic_init (); + +} + +/** Below this electric field, we do not waste time calculating anything. + Besides, this avoid NaNs for eabs == 0. */ +#define EPS_EABS 1e-6