diff --git a/src/cdr.c b/src/cdr.c new file mode 100644 index 0000000000000000000000000000000000000000..2a666b15b6cbb44eb73738f40f619f7a64e7ed25 --- /dev/null +++ b/src/cdr.c @@ -0,0 +1,2546 @@ +/** @file cdr.c +* @brief Routines for handling the CDR (convection-diffusion-reaction) equations +* +* Default interpolation methods.\n +* May be changed elsewhere.\n +* \n +* interpol_method_t *cdr_interpol_inside = \&interpol_quadratic;\n +* interpol_method_t *cdr_interpol_bnd = \&interpol_bilin;\n +* \n +*/ + +#include +#include +#include +#include +#include + +#include "cdr.h" +#include "cstream.h" +#include "grid.h" +#include "interpol2.h" +#include "mapper.h" +#include "parameters.h" +#include "poisson.h" +#include "proto.h" +#include "react_table.h" +#include "rz_array.h" +#include "species.h" +#include "tree.h" + +/** Declaration of the functions to interpolate the densities. */ +decl_mapper_funcs(dens); +mapper_t **dens_mappers, **dens_bnd_mappers; + +int cdr_ext_bnd_cond[4]; +double z_cutoff = 10000000.0; +double max_err = 0.0; +react_table mu; +react_table diff; + +extern pois_problem_t *pois_electrostatic; + +static void set_all_bnd (cdr_grid_t *grid, + int sign, int ir, int iz, int itheta, + int dim0, int inout, + int dim1, int dim1_from, int dim1_to, + int dim2, int dim2_from, int dim2_to); +static void set_axis_bnd (cdr_grid_t *grid); +static void max_update (double *x, double cmp); +static double f_ad (rz_array_t *dens, double efield, double dc, double mc, + int ir, int iz, int itheta, int dim); + +static double psi (double theta); +static void prepare_grid (cdr_grid_t *grid); +static double curv_at (cdr_grid_t *grid, rz_array_t *ar, + int ir, int iz, int itheta, double (*f) (double)); +static int needs_refinement (cdr_grid_t *grid, int ir, int iz, int itheta, + int *in_edge); +static int any_needs_refinement (cdr_grid_t *grid, int ir, int iz, + int *in_edge); +static int brick_needs_refinement (cdr_grid_t *grid, int r0, int z0, + int r1, int z1, int *in_edge); +static void refine_in (cdr_grid_t *grid, int cr0, int cz0, int cr1, int cz1, + int contains_edge); + +static void restrict_from (cdr_grid_t *parent, cdr_grid_t *child); +static int match_grids (cdr_grid_t *fro, cdr_grid_t *to); + +double gauss2_xyz (double x, double y, double z); +int curr_seed = 0; +react_table mu; +react_table diff; + +void aux_dump_frames_r (cdr_grid_t *grid, FILE *fp); + +/** @brief Some initialization is needed for this module which has to be + * inside this function. + */ +void +cdr_init (void) +{ + cdr_ext_bnd_cond[BND_BOTTOM] = cdr_bnd_bottom; + cdr_ext_bnd_cond[BND_TOP] = cdr_bnd_top; + cdr_ext_bnd_cond[BND_RIGHT] = cdr_bnd_right; + + dens_mappers = cdr_mappers_a (interpol_methods_index[cdr_interp_in]); + dens_bnd_mappers = cdr_mappers_a (interpol_methods_index[cdr_interp_bnd]); + + react_table_read("input/default.mobility",&mu); + react_table_read("input/default.diffusion",&diff); +} + +/** @brief And the cleaning up has to be here. */ +void +cdr_end (void) +{ + cdr_free_mappers (dens_mappers); + cdr_free_mappers (dens_bnd_mappers); +} + +/** @brief Creates a new 3D cdr grid */ +cdr_grid_t* +cdr_new_3d_a (int r0, int z0, int r1, int z1, int ntheta) +{ + rz_array_t *ez, *er, *etheta, *eabs, *charge, **dens, **d_dens, *photo; + cdr_grid_t *grid; + int i; + REAL *max_dens; + + debug (2, "cdr_new_3d_a (%d, %d, %d, %d)\n", r0, z0, r1, z1); + + ez = rz_new_3d_a (r0, z0, r1, z1, ntheta); + er = rz_new_3d_a (r0, z0, r1, z1, ntheta); + etheta = rz_new_3d_a (r0, z0, r1, z1, ntheta); + eabs = rz_new_3d_a (r0, z0, r1, z1, ntheta); + charge = rz_new_3d_a (r0, z0, r1, z1, ntheta); + photo = has_photoionization? rz_new_3d_a (r0, z0, r1, z1, ntheta): NULL; + + /* Apart from the real species, (electron, ions...) we can create + * a virtual species that enter into the possible reactions. + * For example, the absolute value of the electric field, eabs. + */ + dens = (rz_array_t **) xmalloc (sizeof(rz_array_t*) + * (no_species + N_VIRTUAL_SPECIES)); + d_dens = (rz_array_t **) xmalloc (sizeof(rz_array_t*) * no_species); + + max_dens = (REAL *) xmalloc (sizeof(REAL) * no_species); + + for (i = 0; i < no_species; i++) { + dens[i] = rz_new_3d_a (r0, z0, r1, z1, ntheta); + d_dens[i] = rz_new_3d_a (r0, z0, r1, z1, ntheta); + } + + dens[no_species] = eabs; + + grid = (cdr_grid_t *) xmalloc (sizeof(cdr_grid_t)); + + grid->r0 = r0; + grid->r1 = r1; + grid->z0 = z0; + grid->z1 = z1; + grid->charge = charge; + grid->er = er; + grid->ez = ez; + grid->etheta = etheta; + grid->eabs = eabs; + grid->photo = photo; + + grid->dens = dens; + grid->d_dens = d_dens; + grid->max_dens = max_dens; + + grid->ntheta = ntheta; + + grid->contains_edge = FALSE; + + init_leaf (grid); + + return grid; +} + +/** @brief Creates a guest cdr grid: + * + * one that does not have memory of its own but + * points to the data of another grid (the @a host grid). This is useful + * to easily set boundary conditions. Probably, this function and + * cdr_new_3d_a should be somehow merged. + */ +cdr_grid_t* +cdr_guest (cdr_grid_t *host, int r0, int z0, int r1, int z1) +{ + rz_array_t *ez, *er, *etheta, *eabs, *charge, **dens, **d_dens, *photo; + cdr_grid_t *grid; + int i; + REAL *max_dens; + + debug (2, "cdr_guest (" grid_printf_str ",%d, %d, %d, %d)\n", + grid_printf_args(host), r0, z0, r1, z1); + + ez = rz_guest (host->ez, r0, z0, r1, z1); + er = rz_guest (host->er, r0, z0, r1, z1); + etheta = rz_guest (host->etheta, r0, z0, r1, z1); + eabs = rz_guest (host->eabs, r0, z0, r1, z1); + charge = rz_guest (host->charge, r0, z0, r1, z1); + photo = (has_photoionization? + rz_guest (host->photo, r0, z0, r1, z1): NULL); + + dens = (rz_array_t **) xmalloc (sizeof(rz_array_t*) + * (no_species + N_VIRTUAL_SPECIES)); + d_dens = (rz_array_t **) xmalloc (sizeof(rz_array_t*) * no_species); + + max_dens = (REAL *) xmalloc (sizeof(REAL) * no_species); + + for (i = 0; i < no_species; i++) { + dens[i] = rz_guest (host->dens[i], r0, z0, r1, z1); + d_dens[i] = rz_guest (host->d_dens[i], r0, z0, r1, z1); + } + + dens[no_species] = eabs; + + grid = (cdr_grid_t *) xmalloc (sizeof(cdr_grid_t)); + + grid->r0 = r0; + grid->r1 = r1; + grid->z0 = z0; + grid->z1 = z1; + + grid->charge = charge; + grid->er = er; + grid->ez = ez; + grid->etheta = etheta; + grid->eabs = eabs; + grid->photo = photo; + + grid->dens = dens; + grid->d_dens = d_dens; + grid->max_dens = max_dens; + + grid->ntheta = host->ntheta; + + grid->contains_edge = FALSE; + + set_leaf (grid, host->parent, NULL, NULL, host->level); + + return grid; +} + +/** @brief Returns a grid with the same dimensions as grid and at the same + * level. + * + * Realize that the data is _not_ copied into the new grid. + */ +cdr_grid_t* +cdr_like_a (cdr_grid_t *grid) +{ + cdr_grid_t *n; + + debug (2, "cdr_like_a(" grid_printf_str ")\n", + grid_printf_args(grid)); + + n = cdr_new_3d_a (grid->r0, grid->z0, grid->r1, grid->z1, grid->ntheta); + n->level = grid->level; + n->ext_bound = grid->ext_bound; + n->contains_edge = grid->contains_edge; + + return n; +} + +/** @brief Returns a grid exactly like te one received, with the densities + * copied. + * + * But the descendants are ignored. + */ +cdr_grid_t* +cdr_clone_a (cdr_grid_t *grid) +{ + cdr_grid_t *n; + int itheta, s, nr, nz, nspecs; + + n = cdr_like_a (grid); + + nr = grid->r1 - grid->r0 + 4; + nz = grid->z1 - grid->z0 + 4; + + nspecs = no_species; + + /* When we use the efield as a refinement criterium, we also have to copy + * eabs (which is a @i virtual density). + */ + if (ref_level_eabs >= 0 && ref_threshold_eabs >= 0.0) nspecs++; + +#pragma omp parallel private(s) + { +#pragma omp for + iter_grid_theta_n (grid, itheta, 2) { + for (s = 0; s < nspecs; s++) { + rz_copy_modes (grid->dens[s], grid->r0 - 2, grid->z0 - 2, + n->dens[s], n->r0 - 2, n->z0 - 2, + nr, nz, itheta, itheta); + } + } + } + return n; +} + +/** @brief Impose periodic boundary conditions in \f$\theta\f$ on a grid. */ +void +cdr_set_periodic (cdr_grid_t *grid) +{ + int i; + + rz_set_periodic (grid->etheta); + + for (i = 0; i < no_species; i++) { + /* For negative r, we have the values + * of the densities on an opposite itheta. + */ + + rz_set_periodic (grid->dens[i]); + rz_set_periodic (grid->d_dens[i]); + + } + + /* I am not completely sure that this is actually neccesary */ + rz_set_periodic (grid->er); + rz_set_periodic (grid->ez); +} + +/** @brief Recursive version of cdr_set_periodic. */ +mk_recursive (cdr_set_periodic, cdr_grid_t) + +/** @brief Frees a cdr grid. */ +void +cdr_free (cdr_grid_t *grid) +{ + int i; + + debug (2, "cdr_free\n"); + + rz_free (grid->charge); + rz_free (grid->er); + rz_free (grid->ez); + rz_free (grid->etheta); + rz_free (grid->eabs); + if (has_photoionization) rz_free (grid->photo); + + for (i = 0; i < no_species; i++) { + rz_free (grid->dens[i]); + rz_free (grid->d_dens[i]); + } + + free (grid->dens); + free (grid->d_dens); + free (grid->max_dens); + + free (grid); +} + +/** @brief Recursive version of cdr_free. */ +void +cdr_free_r (cdr_grid_t *grid) +{ + cdr_grid_t *leaf; + + debug (2, "cdr_free_r (" grid_printf_str ")\n", grid_printf_args (grid)); + + free_childs (grid, leaf, cdr_free_r); + + cdr_free (grid); +} + +/** @brief Calculates the charge in a CDR grid. + * + * Important Note: The charge is here divided by ntheta because + * our FFT functions are unnormalized and hence when they are applied + * forth and back they return the initial value multiplied + * by ntheta. I believe that this is the best point to normalize the + * data. + */ +void +cdr_calc_charge (cdr_grid_t *grid) +{ + int s, r, z, t; + double s_charge; + + debug (3, "cdr_calc_charge\n"); + + rz_set_zero (grid->charge); + + for (s = 0; s < no_species; s++) { + s_charge = spec_index[s]->charge; + +#pragma omp parallel private(r, z) + { +#pragma omp for + iter_grid_theta (grid, t) { + for (r = grid->r0 - 2; r < grid->r1 + 2; r++) { + for (z = grid->z0 - 2; z < grid->z1 + 2; z++) { + RZT (grid->charge, r, z, t) += (RZT (grid->dens[s], r, z, t) + * s_charge) / grid->ntheta; + } + } + } + } + } +} + +/** @brief Recursive version of cdr_calc_charge. */ +mk_recursive (cdr_calc_charge, cdr_grid_t) + +/** @brief Fourier-transforms the charge in grid and all its descendants. */ +void +cdr_dft_charge_r (cdr_grid_t *grid, int sign) +{ + cdr_grid_t *leaf; + + debug (2, "cdr_dft_charge_r(" grid_printf_str ")\n", + grid_printf_args(grid)); + + iter_childs (grid, leaf) { + cdr_dft_charge_r (leaf, sign); + } + + dft_transform (grid->charge, grid->charge, sign); +} + +/** @brief Fourier-transforms the electric field in grid and all its descendants. */ +void +cdr_dft_field_r (cdr_grid_t *grid, int sign) +{ + cdr_grid_t *leaf; + + debug (2, "cdr_dft_field_r(" grid_printf_str ")\n", + grid_printf_args(grid)); + + iter_childs (grid, leaf) { + cdr_dft_field_r (leaf, sign); + } + + dft_transform (grid->er, grid->er, sign); + dft_transform (grid->ez, grid->ez, sign); + + dft_diff ((grid_t *) grid, grid->etheta); + dft_transform (grid->etheta, grid->etheta, sign); +} + +/** @brief Creates a grid one level coarser than GRID that covers the same area + * or somewhat more. + * + * Note: The new grid contains information only about the charge, + * not the other functions. + * That's why we don't use cdr_restrict. + */ +cdr_grid_t* +cdr_create_coarser_a (cdr_grid_t *grid) +{ + int r0_half, r1_half, z0_half, z1_half, itheta; + int r, z; + + cdr_grid_t *new_grid; + + debug (2, "cdr_create_coarser_a\n"); + + r0_half = grid->r0 >> 1; + r1_half = (grid->r1 + 1) >> 1; + + z0_half = grid->z0 >> 1; + z1_half = (grid->z1 + 1) >> 1; + + new_grid = cdr_new_3d_a (r0_half, z0_half, r1_half, z1_half, grid->ntheta); + +#pragma omp parallel private (r, z) + { + /** Note that here we may be going one cell out of the aparent boundaries + of the input grids but we allocated enough memory for that. */ +#pragma omp for + iter_grid_3d_n (new_grid, r, z, itheta, 1) { + /* Start cylindrical */ + /* Here we are interpolating the masses, but in many other places + * we interpolate the densities. Why? + */ + RZT (new_grid->charge, r, z, itheta) = + 0.25 + * (cyl_q (r + 0.25) * (RZT (grid->charge, 2 * r, 2 * z, itheta) + + RZT (grid->charge, 2 * r, 2 * z + 1, itheta)) + + cyl_q (r + 0.75) * (RZT (grid->charge, 2 * r + 1, 2 * z, itheta) + + RZT (grid->charge, 2 * r + 1, 2 * z + 1, + itheta))) + / cyl_q (r + 0.5); + /* End cylindrical */ + } + } + return new_grid; +} + +/** @brief Takes a tree of input nodes and adds n coarser grids, returning a + * pointer to the root leaf. + */ +cdr_grid_t* +cdr_add_coarser_grids_a (cdr_grid_t *prev_root, int n) +{ + int i; + cdr_grid_t *grid = NULL; + + debug (2, "cdr_add_coarser_grids_a\n"); + + assert (n > 0); + + for (i = 0; i < n; i++) { + grid = cdr_create_coarser_a (prev_root); + set_leaf (grid, NULL, NULL, prev_root, prev_root->level - 1); + + prev_root->parent = grid; + + prev_root = grid; + } + return grid; +} + +/** @brief Frees the memory allocated by cdr_add_coarser_grids_a. + * + * Note: be careful to pass here only a tree returned by + * cdr_add_coarser_grids_a, since some grids may lose their references and + * stay allocated forever otherwise. + * (see assert below). + */ +void +cdr_free_coarser_grids (cdr_grid_t *prev_root, int n) +{ + int i; + cdr_grid_t *grid; + + debug (2, "cdr_free_coarser_grids_a\n"); + + for (i = 0; i < n; i++) { + assert (NULL == prev_root->next); + grid = prev_root->first_child; + cdr_free (prev_root); + prev_root = grid; + } + + prev_root->parent = NULL; +} + +/** @brief The interface between the cdr and the poisson parts + * + * If return_pois is true, returns the Poisson trees and they have + * to be de-allocated by the calling program. If not, takes care himself + * of deallocation and returns NULL. + */ +pois_grid_t** +cdr_calc_field_r (cdr_grid_t *cdr, int return_pois) +{ + pois_grid_t **pois_modes; + + debug (2, "cdr_calc_field_r (" grid_printf_str ", %d)\n", + grid_printf_args(cdr), return_pois); + + if (cdr->ntheta != 1) + cdr_dft_charge_r (cdr, 1); + + pois_modes = pois_solve_a (cdr, pois_electrostatic); + + if (cdr->ntheta != 1) + cdr_dft_field_r (cdr, -1); + + cdr_add_ext_field_r (cdr); + + /* Formerly, I added the inhomogeneous field here. But this is problematic: + * see pois_add_inhom_phi_r in poisson.c for details. + * + * To return back to adding electric fields and not potentials, uncomment + * these lines: + + * if (pois_inhom) { + * q_factor = pois_inhom_q_factor (pois_modes[0]); + * debug (1, "q_factor = %g\n", q_factor); + * + * cdr_add_inhom_field_r (cdr, q_factor); + * } + */ + + /* Calculates the absolute value of the field. */ + cdr_calc_eabs_r (cdr); + + /* Impose periodic boundary conditions in theta. */ + if (cdr->ntheta != 1) + cdr_set_periodic_r (cdr); + + if (!return_pois) { + int i; + + for (i = 0; i < max_ntheta; i++) { + pois_free_r (pois_modes[i]); + } + free (pois_modes); + pois_modes = NULL; + } + return pois_modes; +} + +/** @brief Adds an external electric field to a cdr @a grid. + * + * Really, it would be more efficient to include that field in some other + * routine and we will avoid one more grid sweep. However, for the sake + * of code clarity, extendability and maintenability I prefer to add it here. + * + * ext_e_r, ext_e_z and ext_e_theta are functions of three real numbers + * (r, z, \f$\theta\f$) that return each of the components of the electric field. + * + * Note that it is the responsibility of e_theta to behave correctly + * (i.e. periodically) when it receives thetas outside \f$[0, 2 \pi]\f$. + */ +void +cdr_add_ext_field (cdr_grid_t *grid) +{ + int ir, iz, itheta; + + debug (2, "cdr_add_ext_field (" grid_printf_str ")\n", + grid_printf_args(grid)); + +#pragma omp parallel private(ir, iz) + { + #pragma omp for + iter_grid_3d_n (grid, ir, iz, itheta, 2) { + RZT(grid->er, ir, iz, itheta) += + ext_e_r (er_r_at (ir, grid->level), + er_z_at (iz, grid->level), + theta_at (itheta)); + + RZT(grid->ez, ir, iz, itheta) += + ext_e_z (ez_r_at (ir, grid->level), + ez_z_at (iz, grid->level), + theta_at (itheta)); + + RZT(grid->etheta, ir, iz, itheta) += + ext_e_theta (r_at (ir, grid->level), + z_at (iz, grid->level), + etheta_theta_at (itheta)); + } + } +} + +/** @brief Recursive version of cdr_add_ext_field. */ +mk_recursive (cdr_add_ext_field, cdr_grid_t) + + +/** @brief Adds to a cdr grid the inhomogeneous field created by a charge + * q located at (r = 0, z = pois_inhom_z). + * + * THIS FUNCTION IS KEPT FOR REFERENCE ONLY. It was superseded by + * pois_add_inhom_phi_r (see poisson.c) + */ +void +cdr_add_inhom_field_r (cdr_grid_t *grid, double q) +{ + int ir, iz, itheta; + cdr_grid_t *child; + + debug (2, "cdr_add_inhom_field_r (" grid_printf_str ", q = %f)\n", + grid_printf_args(grid), q); + + assert (0 != pois_inhom); + +#pragma omp parallel private(ir, iz) + { + #pragma omp for + iter_grid_3d_n (grid, ir, iz, itheta, 2) { + RZT(grid->er, ir, iz, itheta) += + q * pois_inhom_er (er_r_at (ir, grid->level), + er_z_at (iz, grid->level)); + + RZT(grid->ez, ir, iz, itheta) += + q * pois_inhom_ez (ez_r_at (ir, grid->level), + ez_z_at (iz, grid->level)); + } + } + iter_childs (grid, child) { + cdr_add_inhom_field_r (child, q); + } +} + +/** @brief Calculates the absolute value of the electric field by linear + * interpolation of its components. + * + * Of course, it requires that all the components are already present. */ +void +cdr_calc_eabs (cdr_grid_t *grid) +{ + int ir, iz, itheta; + double er, ez, etheta; + +#pragma omp parallel + { +#pragma omp for private(ir, iz, er, ez, etheta) + iter_grid_3d (grid, ir, iz, itheta) { + er = 0.5 * (RZT (grid->er, ir, iz, itheta) + + RZT (grid->er, ir - 1, iz, itheta)); + ez = 0.5 * (RZT (grid->ez, ir, iz, itheta) + + RZT (grid->ez, ir, iz - 1, itheta)); + if (grid->ntheta != 1) { + etheta = 0.5 * (RZT (grid->etheta, ir, iz, itheta) + + RZT (grid->etheta, ir, iz, itheta - 1)); + } else { + etheta = 0; + } + RZT (grid->eabs, ir, iz, itheta) = sqrt (er * er + ez * ez + + etheta * etheta); + + } + } +} + +/** @brief Recursive version of cdr_calc_eabs. */ +mk_recursive (cdr_calc_eabs, cdr_grid_t) + +/******************** + * Time integration * +*********************/ + +/** @brief Makes sure that all species of one grid are non-negative. */ +void +cdr_nonegative (cdr_grid_t *grid) +{ + int s, ir, iz, itheta; + +#pragma omp parallel + { +#pragma omp for private(ir, iz, s) + iter_grid_theta_n (grid, itheta, 2) { + iter_grid_n(grid, ir, iz, 2) { + for (s = 0; s < no_species; s++) + *RZTP (grid->dens[s], ir, iz, itheta) = + MYMAX(0, *RZTP (grid->dens[s], ir, iz, itheta)); + } + } + } +} + +/** @brief Recursive version of cdr_nonegative. */ +mk_recursive (cdr_nonegative, cdr_grid_t) + +/** @brief Sets the boundary values for the root grid and those grids that + have a boundary coinciding with any of the boundaries of the root grid. */ +void +cdr_set_ext_bnd (cdr_grid_t *grid) +{ + int r0, z0, r1, z1, ntheta; + + debug (2, "cdr_set_ext_bnd (" grid_printf_str " [grid->ext_bound = 0x%x])\n", + grid_printf_args(grid), grid->ext_bound); + + r0 = grid->r0; + z0 = grid->z0; + r1 = grid->r1; + z1 = grid->z1; + ntheta = grid->ntheta; + + /* Matching conditions (reduced to Hom. Neumann for ntheta = 1 at + r = r0 */ + if (grid->ext_bound & BND_MASK (BND_LEFT)) { + set_axis_bnd (grid); + } + + /* At r = r1. */ + if (grid->ext_bound & BND_MASK (BND_RIGHT)) + set_all_bnd (grid, cdr_ext_bnd_cond[BND_RIGHT], r1 - 1, z0 - 2, 0, + R_INDX, 1, + Z_INDX, 0, z1 - z0 + 4, + THETA_INDX, 0, ntheta); + + /* At z = z0. */ + if (grid->ext_bound & BND_MASK (BND_BOTTOM)) + set_all_bnd (grid, cdr_ext_bnd_cond[BND_BOTTOM], r0 - 2, z0, 0, + Z_INDX, -1, + R_INDX, 0, r1 - r0 + 4, + THETA_INDX, 0, ntheta); + + /* At z = z1. */ + if (grid->ext_bound & BND_MASK (BND_TOP)) + set_all_bnd (grid, cdr_ext_bnd_cond[BND_TOP], r0 - 2, z1 - 1, 0, + Z_INDX, 1, + R_INDX, 0, r1 - r0 + 4, + THETA_INDX, 0, ntheta); +} + +/** @brief Recursive version of cdr_set_ext_bnd. */ +mk_recursive (cdr_set_ext_bnd, cdr_grid_t) + +/** @brief Sets some boundary conditions for all the relevant variables in a + * grid, to understand the parameters, see rz_array.c:rz_set_bnd(..) + */ +static void +set_all_bnd (cdr_grid_t *grid, + int sign, int ir, int iz, int itheta, int dim0, int inout, + int dim1, int dim1_from, int dim1_to, + int dim2, int dim2_from, int dim2_to) +{ + int s; + + for (s = 0; s < no_species; s++) { + rz_set_bnd (grid->dens[s], sign, RZTP (grid->dens[s], ir, iz, itheta), + dim0, inout, + dim1, dim1_from, dim1_to, + dim2, dim2_from, dim2_to); + } +} + +/** @brief Sets the proper boundary conditions for a grid that has its + * left boundary at r = 0 (e.g. the root grid). + * + * We don't set Hom. Neumann at r = r0 since r = r0 is + * not a real boundary: the proper condition here is matching the + * opposite grid. But note that this condition reduces to Hom. + * Neumann in the case ntheta = 1 (because itheta2 = itheta = 0). + */ +static void +set_axis_bnd (cdr_grid_t *grid) +{ + int s, itheta, itheta2, ntheta; + + debug (2, "set_axis_bnd (" grid_printf_str ")\n", + grid_printf_args(grid)); + + assert (0 == grid->r0); + + ntheta = grid->ntheta; + + for (s = 0; s < no_species; s++) { + iter_grid_theta (grid, itheta) { + itheta2 = (itheta + ntheta / 2) % ntheta; + + rz_copy_bnd (grid->dens[s], grid->dens[s], + /* sign = */ 1, + /* start_from = */ RZTP (grid->dens[s], + grid->r0, grid->z0 - 2, itheta2), + /* start_to = */ RZTP (grid->dens[s], + grid->r0, grid->z0 - 2, itheta), + R_INDX, BND_OUTWARD, BND_INWARD, + Z_INDX, 0, grid->z1 - grid->z0 + 4, + THETA_INDX, 0, 1 /* <- relative to start_xxx */); + } + } +} + +/** @brief Calculates the derivatives of the densities for each point of the + * grid. + */ +void +cdr_calc_d_dens (cdr_grid_t *grid) +{ + react_apply_all (grid); + + cdr_advect_diffu_r (grid); +} + +/* See "An Adaptive Grid Refinement Strategy for the Simulation + * of Negative Streamers", C. Montijn, W. Hundsdorfer and U. Ebert + * for the computations and notation used in this function. + */ + +/** @brief This thing looks really ugly, but what does it mean? + * For example, for stride = dens[electrons]->strides[R_INDX], + * it is the \f$p_{ij}\f$ of the aforementioned paper: + * \f$p_{i,j} = (\sigma_{i,j} - \sigma_{i-1,j}) / (\sigma_{i+1,j} - \sigma_{i,j})\f$ + */ +#define pij(p_, stride_) \ + (((*(p_)) - *((p_) - (stride_))) / (*((p_) + (stride_)) - *(p_))) + +/** @brief This is the F^a + F^d in the paper (times \f$\delta\f${r, z, \f$\theta\f$}). */ +static double +f_ad (rz_array_t *dens, double efield, double dc, double mc, + int ir, int iz, int itheta, int dim) +{ + double sigma, psi_p, sigmasigma; + int ishift; + double *p; + + /* a p is the unmodified density-value in the cell (ir,iz,itheta). + * Normally this should be n_e, as only electrons are mobile. + */ + p = RZTP (dens, ir, iz, itheta); + + ishift = (efield > 0? 0: 1); + + /* sigma is just + * sigma_{i , j} for efield < 0, + * sigma_{i + 1, j} for efield > 0 + */ + sigma = *(p + ishift * dens->strides[dim]); + + /* psi_p has to be + * psi( p_{i , j}) for efield < 0 + * -psi(1/p_{i+1, j}) for efield > 0 + */ + psi_p = pij (p + ishift * dens->strides[dim], dens->strides[dim]); + + /* E > 0 ---> psi_p = -psi( 1 / p_ij ) + * E < 0, p_ij == 0 ---> psi_p = 1 + * E < 0, p_ij != 0 ---> psi_p = psi( 1 / p_ij ) + */ + psi_p = (efield > 0)? -psi (psi_p) : psi_p == 0? 1.0: psi (1.0 / psi_p); + + /* sigmasigma is (sigma_{i+1, j} - sigma_{i, j}) */ + sigmasigma = (*p) - *(p + dens->strides[dim]); + + /* The first term corresponds to the advection term, the second to the + * diffusion term. + * + * See eq. 3.1.5~3.1.8 in the thesis of C. Montijn. Division by cell-size + * is done in the function calling this one. + */ + return efield * mc * (sigma + psi_p * sigmasigma) + dc * sigmasigma; +} + +/** @brief Adds the advection and diffusion terms to the derivatives of the + * densities. + */ +void +cdr_advect_diffu (cdr_grid_t *grid) +{ + int ir, iz, itheta, s; + int ii; + double gdr_inv, gdz_inv, d_iso0, d_iso; + + double field; + int dimd, t, sg; + REAL *d_dens_p; + double r_inv, er_r; + rz_array_t *er,*ez,*etheta; + + debug (2, "cdr_advect_diffu (" grid_printf_str ")\n", + grid_printf_args(grid)); + + gdr_inv = 1.0 / dr[grid->level]; + gdz_inv = 1.0 / dz[grid->level]; + er = grid->er; + ez = grid->ez; + etheta = grid->etheta; + + for (s = 0; s < no_species; s++) { + /* The charge / mass ratio of a species. */ + double cm0, cm; + int *strides; + int ir, rmax, dim2; + rz_array_t *dens, *d_dens, **d_dens_copy; + + if (spec_index[s]->mass <= 0) continue; + + dens = grid->dens[s]; + d_dens = grid->d_dens[s]; + strides = dens->strides; + rmax = strides[1]; dim2 = strides[2]/strides[1]; + + /* Make a copy of d_dens to be able to compare both versions + * of arcos_cdr_advect + int nr = grid->r1 - grid->r0 + 4; + int nz = grid->z1 - grid->z0 + 4; + d_dens_copy = (rz_array_t **) xmalloc (sizeof(rz_array_t*) * no_species); + rz_copy(d_dens,grid->r0,grid->z0,*d_dens_copy,grid->r0,grid->z0,nr,nz); + + strides = dens->strides; + */ + + arcos_cdr_advect_vec(spec_index[s]->mass,spec_index[s]->charge, + dens->data,d_dens->data,er->data,er->len, + ez->data,ez->len,diffusion_coeff, + dr[grid->level],dz[grid->level], + sprite_module,grid->r0,grid->r1,grid->z0,grid->z1); + + /* rz_free(*d_dens_copy); */ + } +} + +mk_recursive (cdr_advect_diffu, cdr_grid_t) + +/** @brief Limiter function, psi(\f$\theta\f$) = max (0, min(1, 1/3 + \f$\theta\f$/6, \f$\theta\f$)) + */ +static double +psi (double theta) +{ + if (theta < 0) return 0; + if (theta < 0.4) return theta; + if (theta < 4) return (1 + 0.5 * theta) / 3.0; + else return 1.0; +} + +/** @brief Returns the minimum tau that satisfies the Courant-Friedrichs-Lewy + restriction (CFL). */ +double +cdr_courant (cdr_grid_t *grid) +{ + int s, ir, iz, itheta, depth, is_3d; + + double dr_min, dz_min, rdtheta_min, emax_r, emax_z, emax_theta, dmax, ddmax; + double tau_a, tau_d, tau_rt, tau_f, max_diff; + + debug (2, "cdr_courant (" grid_printf_str ")\n", grid_printf_args(grid)); + + is_3d = (grid->ntheta == 1? 0: 1); + + /* Should we better keep track of the level of the finest grid? */ + depth = grid_max_depth_r ((grid_t *) grid); + dr_min = dr[depth]; + dz_min = dz[depth]; + rdtheta_min = dtheta * grid_rmin_r ((grid_t *) grid); + + /* This part is not parallel, but if needed it can be parallelized + * easily by creating vectors xmax[0...ntheta-1], where we would store + * the maximum for each ntheta. + */ + emax_r = -1; + emax_z = -1; + emax_theta = -1; + dmax = -1; + ddmax = 0; + + iter_grid_3d (grid, ir, iz, itheta) { + double mu; + if (sprite_module) { + mu = 1.0 / spr_density_at (z_at (iz, grid->level)); + } else { + mu = 1.0; + } + //mu = sprite_module? 1.0 / spr_density_at (z_at (iz, grid->level)): 1.0; + + max_update (&emax_r, mu * fabs (RZT (grid->er, ir, iz, itheta))); + max_update (&emax_z, mu * fabs (RZT (grid->ez, ir, iz, itheta))); + max_update (&emax_theta, mu * fabs (RZT(grid->etheta, ir, iz, itheta))); + + for (s = 0; s < no_species; s++) { + if (spec_index[s]->mass > 0) { + max_update (&dmax, mu * RZT(grid->dens[s], ir, iz, itheta)); + max_update (&ddmax, -(RZT(grid->d_dens[s], ir, iz, itheta) / + RZT(grid->dens[s], ir, iz, itheta))); + } + } + } + + /* This assert fails when the grid has become completely nan. In that + * case, let us not waste cpu cycles any longer. + */ + assert (emax_r >= 0 && emax_z >= 0 && emax_theta >= 0 && dmax >= 0); + + /* Advection. This one usually dominates. */ + tau_a = nu_a / (emax_r / dr_min + emax_z / dz_min + + is_3d * emax_theta / rdtheta_min); + + /* Diffusion. */ + /* */ + if (sprite_module) { + double dens_z1, dens_z0; + /* Since we do not know where the density is higher, we check at both + * extremes of the domain. We are assuming that the density is monotonic + * but nothing else (if it is not an exp. it works still). + */ + dens_z0 = spr_density_at (z_at (grid->z0, grid->level)); + dens_z1 = spr_density_at (z_at (grid->z1, grid->level)); + max_diff = diffusion_coeff / MYMIN(dens_z0, dens_z1); + } else { + max_diff = 0.1; /* diffusion_coeff */ + } + + tau_d = nu_d / (max_diff / dr_min / dr_min + + max_diff / dz_min / dz_min + + is_3d * max_diff / rdtheta_min / rdtheta_min); + + /* */ + + /* Mainly reaction growth / decay time */ + /* tau_f = nu_f / ddmax; */ + + /* If sigma is zero somewhere, tau_f turns inf even if nu_f is very large. + * We correct this here. */ + /*if (tau_f == 0)*/ tau_f = tau_d; /* Just ignore tau_f please. */ + + /* Dielectric relaxation time */ + tau_rt = nu_rt / dmax; + + return MYMIN (tau_a, MYMIN(tau_d, MYMIN(tau_f, tau_rt))); +} + +static +void max_update (double *x, double cmp) +{ + if (*x < cmp) *x = cmp; +} + +/** @brief Writes into the densities of the grid dest the densities of grid orig + * plus its derivatives times h. + * + * If orig = dest, this implements an Euler step. + * But it can also be used as a tool to implement a RK */ +void +cdr_update (cdr_grid_t *orig, cdr_grid_t *dest, double h) +{ + int ir, iz, itheta, s; + + debug (2, "cdr_update (orig = " grid_printf_str ", dest = " + grid_printf_str ", %f)\n", + grid_printf_args(orig), grid_printf_args(dest), h); + + /* We assume that orig and dest have the same shape. */ + assert (orig->r0 == dest->r0 && orig->z0 == dest->z0 && + orig->r1 == dest->r1 && orig->z1 == dest->z1 && + orig->ntheta == dest->ntheta); + +#pragma omp parallel + { +#pragma omp for private(ir, iz, s) + iter_grid_3d_n (orig, ir, iz, itheta, 2) { + for (s = 0; s < no_species; s++) { + RZT (dest->dens[s], ir, iz, itheta) = + RZT (orig->dens[s], ir, iz, itheta) + + h * RZT (orig->d_dens[s], ir, iz, itheta); + } + } + } + cdr_calc_charge (dest); +} + +/** @brief Sets dest->dens = dens_0->dens + h/2 * d_dens_1->d_dens +\n + * h/2 * d_dens_2->d_dens + * + * Useful to implement 2nd order RK time-stepping. */ +void +cdr_rk2_update (cdr_grid_t *dens_0, cdr_grid_t *d_dens_1, + cdr_grid_t *d_dens_2, cdr_grid_t *dest, + double h) +{ + int ir, iz, itheta, s; + + debug (2, "cdr_rk2_update (dens_0 = " grid_printf_str + ",\n\t\td_dens_1 = " grid_printf_str + ",\n\t\td_dens_2 = " grid_printf_str + ",\n\t\tdest = " grid_printf_str ")\n", + grid_printf_args(dens_0), + grid_printf_args(d_dens_1), + grid_printf_args(d_dens_2), + grid_printf_args(dest)); + +#pragma omp parallel + { +#pragma omp for private(ir, iz, s) + iter_grid_3d_n (dest, ir, iz, itheta, 2) { + for (s = 0; s < no_species; s++) { + RZT (dest->dens[s], ir, iz, itheta) = + RZT (dens_0->dens[s], ir, iz, itheta) + + 0.5 * h * RZT (d_dens_1->d_dens[s], ir, iz, itheta) + + 0.5 * h * RZT (d_dens_2->d_dens[s], ir, iz, itheta); + } + } + } + cdr_calc_charge (dest); +} + +/** @brief Recursive version of cdr_rk2_update + * + * Note: we assume here that all the trees are congruent: this is, + * they have the same structure with nodes of the same size. + */ +void +cdr_rk2_update_r (cdr_grid_t *dens_0, cdr_grid_t *d_dens_1, + cdr_grid_t *d_dens_2, cdr_grid_t *dest, + double h) +{ + cdr_grid_t *c_dens_0, *c_d_dens_1, *c_d_dens_2, *c_dest; + + cdr_rk2_update (dens_0, d_dens_1, d_dens_2, dest, h); + + c_dens_0 = dens_0->first_child; + c_d_dens_1 = d_dens_1->first_child; + c_d_dens_2 = d_dens_2->first_child; + + iter_childs (dest, c_dest) { + cdr_rk2_update_r (c_dens_0, c_d_dens_1, c_d_dens_2, c_dest, h); + + c_dens_0 = c_dens_0->next; + c_d_dens_1 = c_d_dens_1->next; + c_d_dens_2 = c_d_dens_2->next; + } +} + +/** @brief Recursive version of the former QQQQ, used to implement Euler + * integration + */ +void +cdr_self_update_r (cdr_grid_t *grid, double h) +{ + cdr_grid_t *child; + cdr_update (grid, grid, h); + + iter_childs (grid, child) { + cdr_self_update_r (child, h); + } +} + +/** @brief Yet another recursive version of cdr_update, used to implement + * higher order Runge-Kutta integrations. + * + * It creates a new tree of cdr grids. + */ +void +cdr_like_update_ar (cdr_grid_t *grid, cdr_grid_t *new_grid, double h) +{ + cdr_grid_t *child, *new_child, *prev = NULL; + + cdr_update (grid, new_grid, h); + + iter_childs (grid, child) { + new_child = cdr_like_a (child); + + /* We can't use add_child here because we read the children of the grid + * in inverse order as we inserted them. + */ + if (NULL == prev) { + new_grid->first_child = new_child; + } else { + prev->next = new_child; + } + set_leaf (new_child, new_grid, NULL, NULL, new_grid->level + 1); + grid_inherit_ext_bound ((grid_t *) new_child); + + prev = new_child; + + cdr_like_update_ar (child, new_child, h); + } +} + +///** @brief Computes the grid error. QQQQ */ +//double +//cdr_grid_error (cdr_grid_t *grid, cdr_grid_t *inter) +//{ +// int ir, iz, itheta; +// int spec_nr; +// double vol, err, val; +// +// spec_nr = find_species_by_name("electrons"); +// +// assert(spec_nr > -1); +// +// vol = dz[grid->level] * dr[grid->level]; +// +// iter_grid_3d(grid, ir, iz, itheta) +// { +// +// val = RZT(grid->dens[spec_nr], ir, iz, itheta); +// err = val > 1e-12 ? (val - RZT(inter->dens[spec_nr], ir, iz, itheta)) / val : 0; +// err = err > 0 ? err : -err; +// if (max_err < err) +// { +// max_err = err; +// } +// } +// return err; +//} +// +///** @brief Recursive version of cdr_grid_error. QQQQ */ +//void +//cdr_grid_error_r (cdr_grid_t *grid, cdr_grid_t *inter) +//{ +// cdr_grid_t *g_child, *i_child; +// +// cdr_grid_error(grid, inter); +// +// i_child = inter->first_child; +// +// for (g_child = grid->first_child; g_child; g_child = g_child->next) +// { +// cdr_grid_error_r(g_child, i_child); +// +// i_child = i_child->next; +// } +//} + +/** @brief Makes a full second order Runge-Kutta step. + * + * Returns the timestep actually performed, which can be smaller that h + * if it does not satisfy the Courant criterium. + */ +double +cdr_rk2 (cdr_grid_t *grid, double h, double t) +{ + /* intermediate step. */ + cdr_grid_t *inter, *child; + double courant_h; + + /* After this number of warnings about small timesteps, the code stops. */ + static int count_min_timestep = 5; + + debug (2, "cdr_rk2 (" grid_printf_str ", %f)\n", + grid_printf_args(grid), h); + + cdr_nonegative_r (grid); + + /* We want to be sure that the boundaries are properly set. */ + prepare_grid (grid); + + /* Experimental: Moving z_max boundary -- DISABLED */ + + cdr_calc_field_r (grid, /*return_pois = */ FALSE); + + /* If we are using the sprite module, we set here the proper magnitudes + corresponding to the head altitude here. */ + if (sprite_module) spr_hook (grid); + + cdr_calc_d_dens (grid); + + /* The actual timestep is limited by the Courant stability criterium. + * + * Note that cdr_courant has to be called AFTER the computation of the fields. + */ + courant_h = cdr_courant (grid); + h = MYMIN (h, courant_h); + + if (h < warn_min_timestep) { + warning ("Time step [h = %g] is rather small. If this is ok, reduce " + "warn_time_step to remove this warning\n", h); + + if (0 == count_min_timestep--) { + fatal ("Too many warnings about small timesteps.\n"); + } + } + + inter = cdr_like_a (grid); + + /* We create a new tree and put there y + dy/dt * h at each level. */ + cdr_like_update_ar (grid, inter, h); + + /* Set the proper boundary conditions also in the intermediate step. */ + prepare_grid (inter); + + cdr_calc_field_r (inter, /* return_pois = */ FALSE); + cdr_calc_d_dens (inter); + + /* This is something like grid = grid + h / 2 * F(grid) + h / 2 F(inter) */ + cdr_rk2_update_r (grid, grid, inter, grid, h); + + /* Also, for the refinement functions, the boundary conditions have to be + * appropriately set. + */ + prepare_grid (grid); + + cdr_free_r (inter); + + return h; +} + +/** @brief Prepares a grid family : + * a) Set the boundaries for all subgrids,\n + * b) Restrict the coarser grids from the fine ones.\n + */ +static void +prepare_grid (cdr_grid_t *grid) +{ + /* Restrict the densities from finer to coarser grids. */ + cdr_restrict_r (grid); + + /* Sets the external boundaries. */ + cdr_set_ext_bnd_r (grid); + + /* Sets the boundaries of children grids interpolating from their parents. */ + cdr_set_bnd_r (grid); + + /* If two grids share a boundary, use the values of one to set the + boundaries of the other. */ + cdr_match_r (grid, NULL); + + /* Sets periodic boundary conditions in theta. */ + cdr_set_periodic_r (grid); +} + +/***************************************************************** + * Refinement and mapping. + *****************************************************************/ + +/** @brief Takes a pointer to a family of cdr grids, creates a new family + * which is stored in that pointer and takes rid of the initial family. + * + * This function is the entry point of the refinement code. +*/ +void +cdr_update_refined (cdr_grid_t **ptree) +{ + cdr_grid_t *old, *new; + + old = *ptree; + + new = cdr_clone_a (old); + cdr_calc_charge (new); + cdr_refine_r (new, old); + + *ptree = new; + + cdr_free_r (old); +} + +/** @brief Updates the values of max_dens and max_charge. + * + * This values are needed for the refinement criterium. + */ +void +cdr_calc_maxs (cdr_grid_t *grid) +{ + int start = TRUE; + int s, ir, iz, itheta; + + debug (2, "cdr_calc_maxs (" grid_printf_str ")\n", + grid_printf_args(grid)); + + /* If we have a parent, we inherit the maxs from him. */ + if (grid->parent) { + grid->max_charge = grid->parent->max_charge; + for (s = 0; s < no_species; s++) { + grid->max_dens[s] = grid->parent->max_dens[s]; + } + return; + } + + iter_grid_3d (grid, ir, iz, itheta) { + if (start) { + start = FALSE; + grid->max_charge = fabs(RZT (grid->charge, ir, iz, itheta)); + for (s = 0; s < no_species; s++) { + grid->max_dens[s] = RZT (grid->dens[s], ir, iz, itheta); + } + } else { + max_update (&grid->max_charge, fabs(RZT (grid->charge, ir, iz, itheta))); + for (s = 0; s < no_species; s++) { + max_update (&grid->max_dens[s], RZT (grid->dens[s], ir, iz, itheta)); + } + } + } + + debug (2, "->max_charge = %g, ->max_dens[0] = %g, ->max_dens[1] = %g\n", + grid->max_charge, grid->max_dens[0], grid->max_dens[1]); +} + +/** @brief Determines the curvature of ar, which has to be a component of + * grid at @a ir, @a iz. + * + * If @a f is nonnull, determines the curvature of \f$f(ar)\f$. \n + * Note that we only consider the curvatures in @a r and @a z (not in + * \f$\theta\f$). + * This is because \f$\theta\f$ is never refined. (can that be changed in some + * future?) + */ +static double +curv_at (cdr_grid_t *grid, rz_array_t *ar, int ir, int iz, int itheta, + double (*f) (double)) +{ + int i, level; + REAL x[3]; + double dur, duz; + + for (i = 0; i < 3; i++) { + x[i] = RZT (ar, ir - 1 + i, iz, itheta); + x[i] = f? f(x[i]): x[i]; + } + + level = grid->level; + /* */ + dur = (cyl_er_r_at (ir, level) * (x[2] - x[1]) + - cyl_er_r_at (ir - 1, level) * (x[1] - x[0])) / cyl_r_at (ir, level); + /* */ + + for (i = 0; i < 3; i++) { + x[i] = RZT (ar, ir, iz - 1 + i, itheta); + x[i] = f? f(x[i]): x[i]; + } + + duz = (x[2] - 2 * x[1] + x[0]); + + return fabs(dur) + fabs(duz); +} + +/** @brief Does the point in grid at @a ir, @a iz, @a itheta require to be + * further refined? + */ +static int +needs_refinement (cdr_grid_t *grid, int ir, int iz, int itheta, int *in_edge) +{ + double curv; + int s; + /* For sprites, the eabs refinement criterium may depend on the location. */ + int s_ref_level_eabs; + double s_ref_threshold_eabs; + + /* If the maximum charge is below this number, we ignore the refinement + * criterium. + */ + const double epsilon_charge = 1e-8; + + if (!sprite_module || ref_level_eabs < 0) { + s_ref_level_eabs = ref_level_eabs; + s_ref_threshold_eabs = ref_threshold_eabs; + } else { + double z, back_dens; + z = z_at (iz, grid->level); + back_dens = spr_density_at (z); + + s_ref_threshold_eabs = ref_threshold_eabs * back_dens; + + /* log2 appears because whenever we increase the background density by a + * factor 2, we should refine up to one more level. + */ + s_ref_level_eabs = ref_level_eabs + (int) floor (log2 (back_dens)); + } + + /* Electric field refinement. If eabs > ref_threshold_eabs, refine up to + * some predefined level. + */ + + if (grid->level < s_ref_level_eabs && + RZT (grid->eabs, ir, iz, itheta) > s_ref_threshold_eabs) { + debug (4, "Refine grid " grid_printf_str + " at ir = %d, iz = %d, itheta = %d\n" + "\t because too high EABS [eabs = %g]\n", + grid_printf_args(grid), ir, iz, itheta, + RZT (grid->eabs, ir, iz, itheta)); + + *in_edge = FALSE; + return TRUE; + } + + curv = curv_at (grid, grid->charge, ir, iz, itheta, NULL); + + if (grid->max_charge > epsilon_charge + && curv / grid->max_charge > ref_threshold_charge) { + debug (4, "Refine grid " grid_printf_str + " at ir = %d, iz = %d, itheta = %d\n" + "\t because too high CHARGE curvature [curv / max_charge = %g]\n", + grid_printf_args(grid), ir, iz, itheta, curv / grid->max_charge); + + *in_edge = FALSE; + return TRUE; + } + + for (s = 0; s < no_species; s++) { + /* We ignore immobile species: */ + if (spec_index[s]->mass <= 0) continue; + + curv = curv_at (grid, grid->dens[s], ir, iz, itheta, NULL); + if (curv / grid->max_dens[s] > ref_threshold_dens) { + debug (4, "Refine grid " grid_printf_str + " at ir = %d, iz = %d, itheta = %d\n" + "\t because too high DENS[%s] curvature [curv / max_dens = %g]\n", + grid_printf_args(grid), ir, iz, itheta, spec_index[s]->name, + curv / grid->max_dens[s]); + + *in_edge = FALSE; + return TRUE; + } + + /* If the grid contains the leading edge, the criterium of + * refinement also takes into account the density. + */ + if (*in_edge && RZT (grid->dens[s], ir, iz, itheta) > ref_threshold_edge) { + debug (4, "Refine grid " grid_printf_str + " at ir = %d, iz = %d, itheta = %d\n" + "\t because too high DENS[%s] in leading edge [dens = %g]\n", + grid_printf_args(grid), ir, iz, itheta, spec_index[s]->name, + RZT (grid->dens[s], ir, iz, itheta)); + return TRUE; + } + } + + /* If we are close to the border, we also check the border itself. + * + * NOTE: Does this have any sense? We are anyway checking for + * the boundaries, aren't we? + */ + if (ir == grid->r0 + 1 && needs_refinement (grid, ir - 1, iz, itheta, + in_edge)) + return TRUE; + + if (iz == grid->z0 + 1 && needs_refinement (grid, ir, iz - 1, itheta, + in_edge)) + return TRUE; + + if (ir == grid->r1 - 2 && needs_refinement (grid, ir + 1, iz, itheta, + in_edge)) + return TRUE; + + if (iz == grid->z1 - 2 && needs_refinement (grid, ir, iz + 1, itheta, + in_edge)) + return TRUE; + + return FALSE; +} + +/** @brief Calls the needs_refinement at all \f$\theta\f$ and returns true if + * the grids needs to be refined in _any_ of the angles. + * + * This is the safest thing to do, but some other criterium could be imagined + * as well. The search is stopped as soon as a cell that requires refinement + * is found, but not if it requires refinement due to the leading edge density + * criterium. This is because we need to find where we have to stop applying + * that criterium. + */ +static int +any_needs_refinement (cdr_grid_t *grid, int ir, int iz, int *in_edge) +{ + int itheta; + int ret = FALSE; + iter_grid_theta (grid, itheta) { + if (needs_refinement (grid, ir, iz, itheta, in_edge)) { + debug (6, grid_printf_str + " needs refinement at ir = %d, iz = %d itheta = %d\n", + grid_printf_args(grid), ir, iz, itheta); + ret = TRUE; + } + if (ret && ! *in_edge) return TRUE; + } + return ret; +} + +/** @brief Determines whether a @a brick needs refinement. + * + * Now we look if we need refinement in a @a brick \n + * i.e. a rectangle in @a r, @a z space + * (usually, @a height = @a cdr_brick_dz, + * @a width = @a cdr_brick_dr), + * but can be smaller if we reach a boundary. + * + * See above, any_needs_refinement to see when do we stop the search. + */ +static int +brick_needs_refinement (cdr_grid_t *grid, int r0, int z0, int r1, int z1, + int *in_edge) +{ + int ir, iz; + int ret = FALSE; + + for (ir = r0; ir < r1; ir ++) { + for (iz = z1 - 1; iz >= z0; iz--) { + ret = (any_needs_refinement (grid, ir, iz, in_edge) || ret); + if (ret && ! *in_edge) return TRUE; + } + } + return ret; +} + +/** @brief Uses the previous routines to refine (if needed) a cdr grid */ +void +cdr_refine (cdr_grid_t *grid) +{ + int ir; + + debug (2, "cdr_refine (" grid_printf_str ")\n", + grid_printf_args(grid)); + + assert (grid->level <= cdr_max_level); + + /* Checks if the maximum refinement level has been reached. */ + if (grid->level == cdr_max_level) return; + + cdr_calc_maxs (grid); + + for (ir = grid->r0; ir < grid->r1; ir += cdr_brick_dr) { + int irmax = MYMIN (grid->r1, ir + cdr_brick_dr); + int tainted, building, z0 = -1, z1 = -1; + int iz; + int in_edge; + int first_edge; + + first_edge = in_edge = grid->contains_edge; + + /* The leading edge extends downwards. */ + for (iz = grid->z0, building = FALSE; + iz < grid->z1; iz += cdr_brick_dz) { + int izmax = MYMIN (grid->z1, iz + cdr_brick_dz); + + tainted = brick_needs_refinement (grid, ir, iz, irmax, izmax, + &first_edge); + if (tainted) { + z1 = izmax; + if (!building) { + z0 = iz; + building = TRUE; + } + } else /* ! tainted */ { + if (building) { + assert (z0 >= 0 && z1 > 0); + /* If a grid does not satisfy the curvature criterium anywhere + * and is here only because of the density threshold (edge) + * criterium, he does not deserve to be refined. + */ + building = FALSE; + if (!first_edge) { + refine_in (grid, ir, z0, irmax, z1, in_edge); + in_edge = FALSE; + } + } + } + } + if (building) { + assert (z0 >= 0 && z1 > 0); + if (!first_edge) refine_in (grid, ir, z0, irmax, z1, in_edge); + } + } +} + +/** @brief Adds a child to a cdr grid. + * + * The coordinates are in the parent's units + */ +static void +refine_in (cdr_grid_t *grid, int cr0, int cz0, int cr1, int cz1, + int contains_edge) +{ + cdr_grid_t *child; + int nr0, nr1, nz0, nz1; + + nr0 = cr0 << 1; + nz0 = cz0 << 1; + nr1 = cr1 << 1; + nz1 = cz1 << 1; + + child = cdr_new_3d_a (nr0, nz0, nr1, nz1, grid->ntheta); + + add_child (grid, child); + grid_inherit_ext_bound ((grid_t*) child); + + child->contains_edge = contains_edge; + + debug (3, "new child created {r0 = %d, z0 = %d, r1 = %d, z1 = %d, " + "level = %d}\n", child->r0, child->z0, child->r1, child->z1, + child->level); +} + +/** @brief Recursively refines a cdr grid and constructs a tree based on it. + * + * if source == NULL\n + * then uses the initial densities\n + * else sets the densities of the newly created grids by interpolating and + * copying from the tree source. + */ +void +cdr_refine_r (cdr_grid_t *grid, cdr_grid_t *source) +{ + cdr_grid_t *child; + int itheta; + + cdr_calc_charge (grid); + + cdr_refine (grid); + + iter_childs (grid, child) { + if (NULL != source) { +#pragma omp parallel + { +#pragma omp for + iter_grid_theta (child, itheta) { + map_grid_r (dens_mappers, (grid_t *) source, (grid_t *) child, + itheta, + /* copy = */ TRUE, + /* interpol = */ TRUE, + /* coarsen = */ FALSE, + /* s_buff = */ 0, /* t_buff = */ 2); + } + } + } else { + cdr_init_dens (child); + } + cdr_refine_r (child, source); + } +} + +/* Matching of the boundaries: + * + * The boundary conditions of a grid that has a common border with another + * one at the same level have to be matched. This is the purpose of this set + * of functions. + */ + +/** @brief Examines pairs of grids and matches them and all their descendants. + * + * If grid2 == NULL, matches only the descendats of grid1. Hence, this + * routine is usually called as cdr_match_r (root, NULL). + */ +void +cdr_match_r (cdr_grid_t *grid1, cdr_grid_t *grid2) +{ + cdr_grid_t *child1, *child2; + int match = TRUE; + + if (grid2) { + debug (2, "cdr_match_r (grid1 = " grid_printf_str + ", grid2 = " grid_printf_str ")\n", + grid_printf_args (grid1), grid_printf_args (grid2)); + } else { + debug (2, "cdr_match_r (grid1 = " grid_printf_str + ", grid2 = NULL)\n", grid_printf_args (grid1)); + } + + if (NULL != grid2) { + match = match_grids (grid1, grid2); + } + + if (!match) return; + + iter_childs (grid1, child1) { + debug (4, "\tchild1 = " grid_printf_str "\n", + grid_printf_args (child1)); + iter_childs (grid1, child2) { + debug (4, "\t\tchild2 = " grid_printf_str "\n", + grid_printf_args (child2)); + + cdr_match_r (child1, (child1 != child2)? child2: NULL); + } + + if (NULL == grid2) + continue; + + iter_childs (grid2, child2) { + cdr_match_r (child1, child2); + } + } +} + +/** @brief Matches the boundaries of grid to reading from grid @a fro. + * + * Returns TRUE if there was a matching, FALSE otherwise. + * Maybe this function should be rewritten to make it more compact and clear. + */ +static int +match_grids (cdr_grid_t *fro, cdr_grid_t *to) +{ + int s, x0, x1; + + debug (2, "match_grids (to = " grid_printf_str + ", fro = " grid_printf_str ")\n", + grid_printf_args(to), grid_printf_args(fro)); + + assert (to->level == fro->level); + + if (to->r1 == fro->r0) { + x0 = MYMAX (to->z0 - 2, fro->z0 - 2); + x1 = MYMIN (to->z1 + 2, fro->z1 + 2); + if (x1 <= x0) return FALSE; + + for (s = 0; s < no_species; s++) { + rz_copy_bnd (fro->dens[s], to->dens[s], 1, + RZTP (fro->dens[s], fro->r0, x0, 0), + RZTP (to->dens[s], to->r1 - 1, x0, 0), + R_INDX, BND_OUTWARD, BND_OUTWARD, + Z_INDX, 0, x1 - x0, + THETA_INDX, 0, fro->ntheta); + } + return TRUE; + } + + if (to->r0 == fro->r1) { + x0 = MYMAX (to->z0 - 2, fro->z0 - 2); + x1 = MYMIN (to->z1 + 2, fro->z1 + 2); + if (x1 <= x0) return FALSE; + + for (s = 0; s < no_species; s++) { + rz_copy_bnd (fro->dens[s], to->dens[s], 1, + RZTP (fro->dens[s], fro->r1 - 1, x0, 0), + RZTP (to->dens[s], to->r0, x0, 0), + R_INDX, BND_INWARD, BND_INWARD, + Z_INDX, 0, x1 - x0, + THETA_INDX, 0, fro->ntheta); + } + return TRUE; + } + + if (to->z1 == fro->z0) { + x0 = MYMAX (to->r0 - 2, fro->r0 - 2); + x1 = MYMIN (to->r1 + 2, fro->r1 + 2); + if (x1 <= x0) return FALSE; + + for (s = 0; s < no_species; s++) { + rz_copy_bnd (fro->dens[s], to->dens[s], 1, + RZTP (fro->dens[s], x0, fro->z0, 0), + RZTP (to->dens[s], x0, to->z1 - 1, 0), + Z_INDX, BND_OUTWARD, BND_OUTWARD, + R_INDX, 0, x1 - x0, + THETA_INDX, 0, fro->ntheta); + } + return TRUE; + } + + if (to->z0 == fro->z1) { + x0 = MYMAX (to->r0 - 2, fro->r0 - 2); + x1 = MYMIN (to->r1 + 2, fro->r1 + 2); + if (x1 <= x0) return FALSE; + + for (s = 0; s < no_species; s++) { + rz_copy_bnd (fro->dens[s], to->dens[s], 1, + RZTP (fro->dens[s], x0, fro->z1 - 1, 0), + RZTP (to->dens[s], x0, to->z0, 0), + Z_INDX, BND_INWARD, BND_INWARD, + R_INDX, 0, x1 - x0, + THETA_INDX, 0, fro->ntheta); + } + return TRUE; + } + return FALSE; +} + +/** @brief Sets the boundary conditions of all children of grid by + * interpolating from grid itself. + */ +void +cdr_set_bnd (cdr_grid_t *grid) +{ + cdr_grid_t *top, *bottom, *left, *right, *parent; + int itheta; + + debug (2, "cdr_set_bnd (" grid_printf_str ")\n", + grid_printf_args(grid)); + + parent = grid->parent; + if (NULL == parent) { + assert (0 == grid->level); + return; + } + + /* Since our mapping functions are easier to implement when we map + * between rectangular regions (grids) we map the boundary conditions + * by creating from grid four "guest" grids that share memory with him + * and that contains each of the four boundaries. + */ + top = cdr_guest (grid, grid->r0 - 2, grid->z1, + grid->r1 + 2, grid->z1 + 2); + bottom = cdr_guest (grid, grid->r0 - 2, grid->z0 - 2, + grid->r1 + 2, grid->z0); + left = cdr_guest (grid, grid->r0 - 2, grid->z0, + grid->r0, grid->z1); + right = cdr_guest (grid, grid->r1, grid->z0, + grid->r1 + 2, grid->z1); + +#pragma omp parallel + { +#pragma omp for + iter_grid_theta (grid, itheta) { + if (0 == (grid->ext_bound & BND_MASK (BND_TOP))) + map_grid (dens_bnd_mappers, (grid_t *) parent, (grid_t *) top, itheta, + FALSE, TRUE, FALSE, 1, 0); + + if (0 == (grid->ext_bound & BND_MASK (BND_BOTTOM))) + map_grid (dens_bnd_mappers, (grid_t *) parent, (grid_t *) bottom, itheta, + FALSE, TRUE, FALSE, 1, 0); + + if (0 == (grid->ext_bound & BND_MASK (BND_LEFT))) { + assert (grid->r0 != 0); + map_grid (dens_bnd_mappers, (grid_t *) parent, (grid_t *) left, itheta, + FALSE, TRUE, FALSE, 1, 0); + + } + + if (0 == (grid->ext_bound & BND_MASK (BND_RIGHT))) + map_grid (dens_bnd_mappers, (grid_t *) parent, (grid_t *) right, itheta, + FALSE, TRUE, FALSE, 1, 0); + } + } + + cdr_free (top); + cdr_free (bottom); + cdr_free (left); + cdr_free (right); +} + +/** @brief Recursive version of cdr_set_bnd. */ +mk_recursive (cdr_set_bnd, cdr_grid_t) + +/** @brief Mapping of the densities. */ +mapper_t** +cdr_mappers_a (interpol_method_t *interp_method) +{ + mapper_t *mapper, **mappers; + int s, nmappers; + + nmappers = no_species; + + /* mapper->extra represent the species whose density we will map. + * There are no_species species, but no_species + 1 represents + * the abs value of the electric field, which is mapped only to use + * as an electric field-based refinement criterium. + */ + if (ref_level_eabs >= 0 && ref_threshold_eabs >= 0.0) nmappers++; + + mappers = (mapper_t **) xmalloc (sizeof(mapper_t*) * (nmappers + 1)); + + for (s = 0; s < nmappers; s++) { + mapper = (mapper_t *) xmalloc (sizeof(mapper_t)); + mapper->extra = s; + mapper->interpol_method = interp_method; + mapper->copy = dens_copy; + mapper->interpol_set = dens_interpol_set; + mapper->interpol = dens_interpol; + mapper->coarsen = NULL; + mapper->shift_r = 0; + mapper->shift_z = 0; + mappers[s] = mapper; + } + + mappers[nmappers] = NULL; + + return mappers; +} + +/** @brief Frees the memory of all mappers QQQQ */ +void +cdr_free_mappers (mapper_t **mappers) +{ + int s; + + for (s = 0; s < no_species; s++) { + free (mappers[s]); + } + free (mappers); +} + +/** @brief dens_copy Copies QQQQ */ +void +dens_copy (mapper_t *mapper, grid_t *source, grid_t *target, + int ir, int iz, int itheta) +{ + cdr_grid_t *cdr_s, *cdr_d; + + cdr_s = (cdr_grid_t*) source; + cdr_d = (cdr_grid_t*) target; + + RZT (cdr_d->dens[mapper->extra], ir, iz, itheta) = + RZT (cdr_s->dens[mapper->extra], ir, iz, itheta); +} + +/** @brief Interpolates a set QQQQ */ +int +dens_interpol_set (mapper_t *mapper, grid_t *source, interpol_t *interpol, + int pr, int pz, int itheta) +{ + cdr_grid_t *cdr; + + cdr = (cdr_grid_t*) source; + + interpol_set_stencil_at (source, interpol, + r_at (pr, cdr->level), + z_at (pz, cdr->level), + cdr->dens[mapper->extra], pr, pz, itheta); + return TRUE; +} + +/** @brief dens_interpol Interpolates QQQQ */ +void +dens_interpol (mapper_t *mapper, grid_t *source, grid_t *target, + interpol_t *interpol, + int ir, int iz, int itheta) +{ + double r, z; + cdr_grid_t *cdr; + + cdr = (cdr_grid_t *) target; + + r = r_at (ir, cdr->level); + z = z_at (iz, cdr->level); + + RZT (cdr->dens[mapper->extra], ir, iz, itheta) = + interpol_apply (interpol, r, z); +} + +/** @brief Restricts the values of the densities from the child grid into the + * parent. + */ +static void +restrict_from (cdr_grid_t *parent, cdr_grid_t *child) +{ + int ir, iz, cr, cz, itheta, s; + REAL sum; + + debug (2, "restrict_from (parent =" grid_printf_str + ", child = " grid_printf_str ")\n", + grid_printf_args(parent), grid_printf_args(child)); + + assert (parent->level == child->level - 1); + +#pragma omp parallel private(sum, ir, iz, cr, cz, s) + { +#pragma omp for + iter_grid_theta (parent, itheta) { + iter_grid_parent (child, ir, iz) { + for (s = 0; s < no_species; s++) { + sum = 0; + for (cr = (ir << 1); cr <= (ir << 1) + 1; cr++) + for (cz = (iz << 1); cz <= (iz << 1) + 1; cz++) + sum += cyl_r_at (cr, child->level) + * RZT (child->dens[s], cr, cz, itheta); + RZT (parent->dens[s], ir, iz, itheta) = + 0.25 * sum / cyl_r_at (ir, parent->level); + } + } + } + } +} + +/** @brief Restricts a grid with all its children. + */ +void +cdr_restrict (cdr_grid_t *grid) +{ + cdr_grid_t *child; + + iter_childs(grid, child) { + restrict_from (grid, child); + } +} + +/** @brief And its recursive version. + * + * Note that this has to be tail recursive, sice we first have to restrict + * the childs. + */ +mk_tail_recursive (cdr_restrict, cdr_grid_t) + + +/** @brief Now the routines to initialize a grid with given densities. + * gauss2 initializes one or two gaussian seeds + */ +double +gauss2_xyz (double x, double y, double z) +{ + double q; + + q = invpi32 * 1 / (seed_index[curr_seed]->sigma_x * seed_index[curr_seed]->sigma_y * + seed_index[curr_seed]->sigma_z) + * exp (- SQ(x - seed_index[curr_seed]->x0) / SQ(seed_index[curr_seed]->sigma_x) + - SQ(y - seed_index[curr_seed]->y0) / SQ(seed_index[curr_seed]->sigma_y) + - SQ(z - seed_index[curr_seed]->z0) / SQ(seed_index[curr_seed]->sigma_z)); + return q; +} + +/** @brief Sets the densities of species species according to the function + * \f$f (x, y, z)\f$. + * + * mode can be one of the following:\n + * SET_DENS_OVERWRITE: Set the dens to f(x, y, x), ignoring what was there + * before.\n + * SET_DENS_ADD: Adds f to the former density.\n + * SET_DENS_SUB: Substracts f from the former density.\n + */ +void +cdr_set_dens (cdr_grid_t *cdr, int species, int mode, double factor, + double (*f) (double, double, double)) +{ + int ir, iz, itheta; + double ctheta, stheta, x, y; + debug (2, "cdr_set_dens (" grid_printf_str ", %s, ...)\n", + grid_printf_args(cdr), spec_index[species]->name); + +#pragma omp parallel private(ir, iz, ctheta, stheta, x, y) + { +#pragma omp for + iter_grid_theta_n (cdr, itheta, 2) { + //sincos (theta_at(itheta), &stheta, &ctheta); + stheta=sin(theta_at(itheta)); + ctheta=cos(theta_at(itheta)); + iter_grid_r_n(cdr, ir, 2) { + x = r_at(ir, cdr->level) * ctheta; + y = r_at(ir, cdr->level) * stheta; + iter_grid_z_n(cdr, iz, 2) { + switch (mode) { + case SET_DENS_OVERWRITE: + *RZTP(cdr->dens[species], ir, iz, itheta) = + factor * f (x, y, z_at(iz, cdr->level)); + break; + case SET_DENS_ADD: + *RZTP(cdr->dens[species], ir, iz, itheta) += + factor * f (x, y, z_at(iz, cdr->level)); + break; + case SET_DENS_SUB: + *RZTP(cdr->dens[species], ir, iz, itheta) -= + factor * f (x, y, z_at(iz, cdr->level)); + break; + } + } + } + } + } +} + +/** @brief Returns 1.0 */ +double f_one(double x, double y, double z) +{ + return 1.0; +} + +/** @brief Inits the densities of all species. */ +void +cdr_init_dens (cdr_grid_t *cdr) +{ + int cnt; + + debug (2, "cdr_init_dens (" grid_printf_str ")\n", + grid_printf_args(cdr)); + + for (cnt = 0; cnt < no_species; cnt++) + { + cdr_set_dens(cdr, cnt, SET_DENS_OVERWRITE, 0.0, f_one); + } + + for (cnt = 0; cnt < no_seed; cnt++) + { + curr_seed = cnt; + if (seed_index[cnt]->type == 0) + cdr_set_dens(cdr, seed_index[cnt]->species, SET_DENS_ADD, seed_index[cnt]->value, gauss2_xyz); + if (seed_index[cnt]->type == 1) + cdr_set_dens(cdr, seed_index[cnt]->species, SET_DENS_ADD, seed_index[cnt]->value, f_one); + } +} + +/** @brief Initializes the densities allocating a new grid family.*/ +cdr_grid_t* +cdr_scratch_init (void) +{ + cdr_grid_t *cdr; + + cdr = cdr_new_3d_a (0, 0, gridpoints_r, gridpoints_z, max_ntheta); + + cdr->level = 0; + cdr->ext_bound = BND_MASK_ALL; + cdr->contains_edge = TRUE; + + cdr_init_dens (cdr); + cdr_refine_r (cdr, NULL); + + if (perturb_epsilon > 0.0 && max_ntheta > 1) { + dft_dens_perturb_r (cdr, electrons, NULL); + + /* Maybe the perturbation has messed up the boundary conditions, so we + * have to repair them. + */ + cdr_set_ext_bnd_r (cdr); + cdr_set_periodic_r (cdr); + } + return cdr; +} + +/** @brief Dumps all the contents of the given grid into filenames given by + * prefix and name + */ +void +cdr_dump (cdr_grid_t *grid, const char *prefix, const char *name) +{ + char *fname; + int s, nt; + + int m = cdr_output_margin; + + /* To make it easier to produce plots in theta, we save one extra + * theta that has to give the same data as the first one. Except if + * we are working in 2D, where anyway we have only one possible theta (0) + */ + nt = grid->ntheta == 1? 1: grid->ntheta + 1; + + +#pragma omp parallel sections private(fname) + { +#pragma omp section + { + asprintf (&fname, "%s/r.%s.tsv", prefix, name); + rz_axis_dump (fname, grid->r0 - m, grid->r1 + m, dr[grid->level]); + free (fname); + } +#pragma omp section + { + asprintf (&fname, "%s/z.%s.tsv", prefix, name); + rz_axis_dump (fname, grid->z0 - m , grid->z1 + m, dz[grid->level]); + free (fname); + } +#pragma omp section + { + /* In all this function we use ntheta + 1 to make matlab's life easier. + * Note that some of these variables have not been made periodic + * and hence the value for ntheta may be undefined. + */ + asprintf (&fname, "%s/theta.%s.tsv", prefix, name); + rz_axis_dump (fname, 0, nt, dtheta); + free (fname); + } +#pragma omp section + { + asprintf (&fname, "%s/charge.%s.tsv", prefix, name); + rz_dump_3d (grid->charge, fname, "w", grid->r0 - m, grid->z0 - m, + grid->r1 + m, grid->z1 + m, nt); + free (fname); + } +#pragma omp section + { + asprintf (&fname, "%s/er.%s.tsv", prefix, name); + rz_dump_3d (grid->er, fname, "w", grid->r0 - m, grid->z0 - m, + grid->r1 + m, grid->z1 + m, nt); + free (fname); + } +#pragma omp section + { + asprintf (&fname, "%s/ez.%s.tsv", prefix, name); + rz_dump_3d (grid->ez, fname, "w", grid->r0 - m, grid->z0 - m, + grid->r1 + m, grid->z1 + m, nt); + free (fname); + } +#pragma omp section + { + asprintf (&fname, "%s/etheta.%s.tsv", prefix, name); + rz_dump_3d (grid->etheta, fname, "w", grid->r0 - m, grid->z0 - m, + grid->r1 + m, grid->z1 + m, nt); + free (fname); + } + +#pragma omp section + { + asprintf (&fname, "%s/eabs.%s.tsv", prefix, name); + rz_dump_3d (grid->eabs, fname, "w", grid->r0 - m, grid->z0 - m, + grid->r1 + m, grid->z1 + m, nt); + free (fname); + } + +#pragma omp section + { + if (has_photoionization) { + asprintf (&fname, "%s/photo.%s.tsv", prefix, name); + rz_dump_3d (grid->photo, fname, "w", grid->r0 - m, grid->z0 - m, + grid->r1 + m, grid->z1 + m, nt); + free (fname); + } + } + +#pragma omp section + { + for (s = 0; s < no_species; s++) { + /* The densities */ + asprintf (&fname, "%s/%s.%s.tsv", prefix, spec_index[s]->name, name); + rz_dump_3d (grid->dens[s], fname, "w", grid->r0 - m, grid->z0 - m, + grid->r1 + m, grid->z1 + m, nt); + free (fname); + + /* ...and their derivatives. */ + asprintf (&fname, "%s/d_%s.%s.tsv", prefix, spec_index[s]->name, name); + rz_dump_3d (grid->d_dens[s], fname, "w", grid->r0 - m, grid->z0 - m, + grid->r1 + m, grid->z1 + m, nt); + free (fname); + } + } + } +} + +/** @brief Dumps all the contents of the given grid into filenames given by + * @a prefix and @a name. + * + * It also writes in tree.NAME.dat the tree structure of grid, in a format + * appropiate for cdr_load_tree to understand. + + @a infp has to be NULL when called from outside. + + The format of a tree.NAME.dat is as follows. Each line is composed + by a "command" and some parameters. the "commands" can be + + * time + parameter: sim_time + Writes the current simulation time, as received in sim_time + + * open + parameters: gridname r0 z0 r1 z1 level ext_bound margin + + Reads the densities from the files given by gridname, with + the corresponding values of the grid parameters. Leaves the grid + open to add childrens to it. Those children will be all the grids + read until the grid is closed. + + * close + parameters: gridname + + Closes the latest opened grid. All subsequent grids will be regarded + as brothers of gridname. + */ +void +cdr_dump_r (cdr_grid_t *grid, const char *prefix, const char *name, + FILE *infp, double sim_time) +{ + cdr_grid_t *child; + char *cname; + int i, nchilds; + char codes[] = "abcdefghijklmnopqrstuvwxyz"; + FILE *fp; + + if (NULL == infp) { + /* Root call */ + asprintf (&cname, "%s/tree.%s.dat", prefix, name); + fp = fopen (cname, "w"); + if (NULL == fp) { + fatal ("Could not open file %s/tree.%s.dat to write\n", + prefix, name); + } + free (cname); + fprintf (fp, "time %g\n", sim_time); + } else { + fp = infp; + } + + /* We want to be sure that the charge that we plot is actually the + * charge and not its Fourier transform. + */ + cdr_calc_charge_r (grid); + + fprintf (fp, "open %s %d %d %d %d %d %d %d\n", name, grid->r0, grid->z0, + grid->r1, grid->z1, grid->level, grid->ext_bound, + cdr_output_margin); + + cdr_dump (grid, prefix, name); + + nchilds = grid_howmany_children ((grid_t *) grid); + + for (i = 0; i < nchilds; i++) { + child = (cdr_grid_t *) grid_get_child ((grid_t*) grid, nchilds - i - 1); + + assert (NULL != child); + + if (i > sizeof(codes) - 2) + asprintf (&cname, "%s{%d}", name, i); + else + asprintf (&cname, "%s%c", name, codes[i]); + + cdr_dump_r (child, prefix, cname, fp, sim_time); + free (cname); + } + + fprintf (fp, "close %s\n", name); + + if (NULL == infp) { + /* Root call */ + fclose (fp); + } +} + +/** @brief Reads the data format produced by cdr_dump_r and creates a complete + * tree from it. + * + * This function is really quick-and-dirty: one should check the fscanfs + * and do everything cleanlier. + */ +cdr_grid_t * +cdr_load_tree_r (const char *prefix, const char *name, FILE *infp) +{ + char command[16], gridname[32], *fname; + int r0, z0, r1, z1, level, ext_bound, margin, nt, s; + int open_close; + cdr_grid_t *grid = NULL, *leaf; + FILE *fp; + + debug (2, "cdr_load_tree_r (prefix = \"%s\", name = \"%s\", ...)\n", + prefix, name); + + nt = max_ntheta == 1? 1: max_ntheta + 1; + + if (NULL == infp) { + asprintf (&fname, "%s/tree.%s.dat", prefix, name); + fp = fopen (fname, "r"); + if (NULL == fp) { + fatal ("Could not open file %s/tree.%s.dat to read\n", + prefix, name); + } + free (fname); + } else { + fp = infp; + } + + do { + fscanf (fp, "%15s %15s", command, gridname); + open_close = TRUE; + + /* We analyze some commands that do not create/close grids, now + * restricted to set time. + */ + if (0 == strcmp (command, "time")) { + double new_time; + sscanf (gridname, "%lf\n", &new_time); + warning ("Starting time [%f] was read from %s/tree.%s.dat\n", + new_time, prefix, name); + + start_t = new_time; + open_close = FALSE; + } + } while (!open_close); + + if (0 == strcmp (command, "open")) { + debug (3, "opening %s\n", gridname); + + fscanf (fp, "%d %d %d %d %d %d %d\n", + &r0, &z0, &r1, &z1, &level, &ext_bound, &margin); + + grid = cdr_new_3d_a (r0, z0, r1, z1, max_ntheta); + grid->ext_bound = ext_bound; + grid->level = level; + + for (s = 0; s < no_species; s++) { + /* The densities */ + asprintf (&fname, "%s/%s.%s.tsv", prefix, spec_index[s]->name, gridname); + debug (3, "Loading %s\n", fname); + rz_dump_3d (grid->dens[s], fname, "r", + grid->r0 - margin, grid->z0 - margin, + grid->r1 + margin, grid->z1 + margin, nt); + free (fname); + } + + do { + leaf = cdr_load_tree_r (prefix, name, fp); + if (NULL != leaf) { + assert (leaf->level == grid->level + 1); + add_child (grid, leaf); + } + + } while (NULL != leaf); + } else if (0 == strcmp (command, "close")) { + debug (3, "closing %s\n", gridname); + grid = NULL; + } + + if (NULL == infp) { + fclose (fp); + } + return grid; +} + +/** @brief Given a grid family, writes in the file named @a fname the + * coordinates of the frames that define the grid and his descendants. + */ +void +cdr_dump_frames (cdr_grid_t *grid, const char *prefix, const char *name) +{ + FILE *fp; + char *fname; + + asprintf (&fname, "%s/frames.%s.tsv", prefix, name); + fp = fopen (fname, "w"); + free (fname); + + if (fp == NULL) { + warning ("Unable to open %s\n", fname); + return; + } + aux_dump_frames_r (grid, fp); + + fclose (fp); +} + +/** @brief aux_dump_frames_r QQQQ */ +void +aux_dump_frames_r (cdr_grid_t *grid, FILE *fp) +{ + cdr_grid_t *child; + int level; + + level = grid->level; + + fprintf (fp, "%g %g %g %g %d\n", + er_r_at (grid->r0 - 1, level), ez_z_at (grid->z0 - 1, level), + er_r_at (grid->r1 - 1, level), ez_z_at (grid->z1 - 1, level), + level); + + iter_childs(grid, child) { + aux_dump_frames_r (child, fp); + } +}