diff --git a/src/poisson.c b/src/poisson.c new file mode 100644 index 0000000000000000000000000000000000000000..ddde15b19c0ec17a6c1e81150c846b72c11f5710 --- /dev/null +++ b/src/poisson.c @@ -0,0 +1,1564 @@ +/** @file poisson.c + * @brief Poisson/Helmholtz solver, including the routines for + * inhomogeneous (point-plane) configurations. + */ +#include +#include +#include + +#include "cdr.h" +#include "cstream.h" +#include "fishpack.h" +#include "grid.h" +#include "interpol2.h" +#include "mapper.h" +#include "parameters.h" +#include "poisson.h" +#include "proto.h" +#include "rz_array.h" +#include "species.h" + +static void get_bound_cond (pois_grid_t *grid, pois_problem_t *prob, + int *r_bound_cond, int *z_bound_cond, + double lambda); +static int needs_refinement (pois_grid_t *grid, int ir, int iz, + double threshold); +static int refine_in (pois_grid_t *grid, int cr0, int cz0, int cr1, int cz1); + +decl_mapper_funcs(er); +decl_mapper_funcs(ez); +decl_mapper_funcs(etheta); +decl_mapper_funcs(charge); +/*!< Here we set three mapper_t objects that specify how we should do + * the mapping between the potential (once it is calculated) and the three + * components of the electric field. + */ + +/* Global variables for the computations of inhomogeneous fields: + * initialized in pois_inhom_init. */ + +int pois_gridpoints_z, pois_gridpoints_r; +/*!< Number of gridpoints at level 0 of the poisson solver + * + * if !pois_inhom, this has to be equal to the parameter gridpoins_z. + * We also allow a different extension in r, albeit this is not yet used. + */ + +double pois_inhom_L; /*!< Total length of the poisson domain. */ +double pois_inhom_z; /*!< Location of the floating charge. */ +double pois_inhom_dphi; +/*!< Potential difference between the electrodes created by a unit charge + * located at (0, pois_inhom_z) + */ + +double pois_inhom_fixed_q_t; +/*!< Time-dependent fixed charge for the inhomogeneous poisson code + * It is updated in cstream_set_field_at_time acording to the parameter + * rise_time. + */ + +extern double E_x, E_y, E_z; + +static double inhom_phi_term (double r, double z, double b); +static double inhom_e_term (double r, double z, double b); +static double inhom_er_term (double r, double z, double b); +static double inhom_ez_term (double r, double z, double b); +/*!< To calculate the potential and electric fields created by a unit charge + * between two planar electrodes, we use the "mirror charges method" + * the contribution of each mirror charge is given by these functions. + */ + +mapper_t er_mapper = mk_mapper_staggered (er, &interpol_bilin, 0, -1); +mapper_t ez_mapper = mk_mapper_staggered (ez, &interpol_bilin, -1, 0); +mapper_t etheta_mapper = mk_mapper (etheta, &interpol_quadratic); +mapper_t *e_mappers[] = {&er_mapper, &ez_mapper, ðeta_mapper, NULL}; +/*!< These are the mapper structures to translate between cdr and poisson + * grids. See mapper.c and mapper.h for more info. + */ + +mapper_t charge_mapper = mk_mapper_down(charge, &interpol_wackers); +mapper_t *charge_mappers[] = {&charge_mapper, NULL}; + +pois_problem_t *pois_electrostatic; + +/** @brief All initialization of the Poisson solver has to be here. */ +void +pois_init (void) +{ + if (pois_inhom) { + pois_inhom_init (); + } + else { + pois_gridpoints_z = gridpoints_z; + } + pois_gridpoints_r = gridpoints_r; + + pois_electrostatic = xmalloc (sizeof (pois_problem_t)); + + pois_electrostatic->max_level = pois_max_level; + pois_electrostatic->extra_levels = extra_pois_levels; + pois_electrostatic->max_error = pois_max_error; + + pois_electrostatic->bnd_right = pois_bnd_right; + pois_electrostatic->bnd_top = pois_bnd_top; + pois_electrostatic->bnd_bottom = pois_bnd_bottom; +} + +/** @brief Creates a 2D Poisson grid. */ +pois_grid_t* +pois_new_a (int r0, int z0, int r1, int z1) +{ + pois_grid_t *grid; + + grid = pois_new_3d_a (r0, z0, r1, z1, 1); + + return grid; +} + +/** @brief Creates a new 3D Poisson grid. + * + * @a r0, @a z0, @a r1, @a z1 are in GRID units + */ +pois_grid_t* +pois_new_3d_a (int r0, int z0, int r1, int z1, int ntheta) +{ + rz_array_t *phi, *charge, *error; + pois_grid_t *grid; + + debug (2, "pois_new_a (r0 = %d, z0 = %d, r1 = %d, z1 = %d)\n", + r0, z0, r1, z1); + + phi = rz_new_3d_a (r0, z0, r1, z1, ntheta); + error = rz_new_3d_a (r0, z0, r1, z1, ntheta); + +#ifdef F_PHI_SAME_AS_CHARGE + charge = phi; +#else + charge = rz_new_3d_a (r0, z0, r1, z1, ntheta); +#endif /* F_PHI_SAME_AS_CHARGE */ + + grid = (pois_grid_t *) xmalloc (sizeof (pois_grid_t)); + + grid->r0 = r0; + grid->r1 = r1; + grid->z0 = z0; + grid->z1 = z1; + grid->ntheta = ntheta; + + grid->phi = phi; + grid->charge = charge; + grid->error = error; + + /* By default, the array does not have any external boundary. */ + grid->ext_bound = BND_NONE; + + /* Initially, a grid does not have any relatives. */ + init_leaf (grid); + + return grid; +} + +/** @brief Frees the memory allocated by pois_new_a */ +void +pois_free (pois_grid_t *grid) +{ + debug (3, "pois_free (...)\n"); + + /* The decission of whether to use the same memory for the charge and + * the potential can be made elsewhere. + */ + if (grid->phi != grid->charge) + rz_free (grid->charge); + rz_free (grid->phi); + rz_free (grid->error); + + free (grid); +} + +/** @brief Recursively frees a tree of Poisson grids. */ +void +pois_free_r (pois_grid_t *grid) +{ + pois_grid_t *leaf; + + debug (3, "pois_free_r (" grid_printf_str ")\n", grid_printf_args (grid)); + + free_childs (grid, leaf, pois_free_r); + + pois_free (grid); +} + +/** @brief Creates a new Poisson grid. + * + * @a r0, @a r1... are in global units (i.e. independent of the grid level). + */ +pois_grid_t* +pois_new_glob_a (int r0, int z0, int r1, int z1, int level) +{ + debug (3, "pois_new_glob_a (%d, %d, %d, %d)\n", r0, r1, z0, z1); + pois_grid_t *grid; + + if (level >= 0) { + grid = pois_new_a (r0 << level, z0 << level, + r1 << level, z1 << level); + } else { + level = -level; + grid = pois_new_a (r0 >> level, z0 >> level, + r1 >> level, z1 >> level); + } + return grid; +} + +/** @brief Starts the tree of Poisson grids with the two coarsest ones. + * + * Receives the grid dimensions at level 0. + */ +pois_grid_t* +pois_init_tree_a (int r0, int z0, int r1, int z1) +{ + pois_grid_t *coarsest, *sub_coarsest; + int level; + + debug (3, "pois_init_tree_a (%d, %d, %d, %d)\n", r0, z0, r1, z1); + + level = -extra_pois_levels; + + coarsest = pois_new_glob_a (r0, z0, r1, z1, level); + coarsest->ext_bound = BND_MASK_ALL; + coarsest->level = level; + + sub_coarsest = pois_new_glob_a (r0, z0, r1, z1, level + 1); + + add_child (coarsest, sub_coarsest); + grid_inherit_ext_bound ((grid_t *) sub_coarsest); + + return coarsest; +} + +/** @brief Is it possible to use a single function to handle all four + * boundaries? + * + * Maybe this should be rewritten to use the interpol2 module? + */ +REAL* +pois_boundary_a (pois_grid_t *grid, int boundary) +{ + pois_grid_t *parent; + + /* Forget the 0s, they are there to shut up the compiler. */ + int stride = 0, ortho_stride = 0, icells, invert, n; + int rb, zb; + double interp_array[3][3] = {{-4.0, 1.0, 18.0}, + {-2.0, 78.0, 14.0}, + {-6.0, -7.0, 4.0}}; + + const double prefactor = 1 / 96.0; + REAL *bnd, *p, *lstart; + + debug (3, "pois_boundary_a (..., %d)\n", boundary); + + parent = grid->parent; + + /* First we check if the boundary is in r or in z. */ + if (boundary & BND_AT_Z) { + icells = grid->r1 - grid->r0; + zb = 1; rb = 0; + if (NULL != parent) { + stride = parent->phi->strides[R_INDX]; + ortho_stride = parent->phi->strides[Z_INDX]; + } + } else { + zb = 0; rb = 1; + icells = grid->z1 - grid->z0; + if (NULL != parent) { + stride = parent->phi->strides[Z_INDX]; + ortho_stride = parent->phi->strides[R_INDX]; + } + } + + debug (3, " grid->ext_bound = %d\n", grid->ext_bound); + + if (grid->ext_bound & BND_MASK (boundary)) { + /* If the boundary is external, we apply homogeneous Dirichlet/Neumann. */ + bnd = (REAL *) xcalloc (sizeof (REAL), icells); + debug (3, "pois_boundary_a (...) -> [0.0] * %d\n", icells); + + assert (grid->r1 != 0); + return bnd; + } else { + bnd = (REAL *) xmalloc (sizeof (REAL) * icells); + } + /* If the grid is root it should have returned in the previous if */ + assert (NULL != parent); + + /* Now we check if the boundary is at {r,z} minimum or {r,z} maximum */ + if (boundary & BND_MAX) { + lstart = RZP (parent->phi, + (grid->r1 >> 1) - 1 + rb, + (grid->z1 >> 1) - 1 + zb); + stride = -stride; + ortho_stride = -ortho_stride; + invert = TRUE; + } else { + lstart = RZP (parent->phi, (grid->r0 >> 1) - rb, (grid->z0 >> 1) - zb); + invert = FALSE; + } + + for (p = lstart, n = 0; n < icells; p += stride, n += 2) { + /* p points to the cell in the parent grid and p + stride points to + an adjacent cell. */ + double s = 0.0, s2 = 0.0; + int i, j; + for (i = -1; i <= 1; i++) { + for (j = -1; j <= 1; j++) { + s += *(p + i * stride + j * ortho_stride) * interp_array[i + 1][j + 1]; + s2 += *(p - i * stride + j * ortho_stride) * interp_array[i + 1][j + 1]; + } + } + if (!invert) { + bnd[n] = prefactor * s; + if (n < icells - 1) bnd[n + 1] = prefactor * s2; + } else { + bnd[icells - n - 1] = prefactor * s; + if (n < icells - 1) bnd[icells - n - 2] = prefactor * s2; + } + } + return bnd; +} + +/** @brief Solves the Poisson equation, calling the FISHPACK routine. + * + * Assumes that grid->charge is already set. + */ +void +pois_solve_grid (pois_grid_t *grid, pois_problem_t *prob, + double lambda, double s) +{ + REAL *boundaries[4]; + int i, z_bound_cond, r_bound_cond; + double rmin, rmax, zmin, zmax; + + debug (3, "pois_solve_grid (..., lambda = %f, s = %f)\n", lambda, s); + + for (i = 0; i < 4; i++) { + boundaries[i] = pois_boundary_a (grid, i); + debug (3, "boundaries[%d] = %p\n", i, boundaries[i]); + } + +#ifndef F_PHI_SAME_AS_CHARGE + rz_copy (grid->charge, grid->r0, grid->z0, + grid->phi, grid->r0, grid->z0, + grid->r1 - grid->r0, grid->z1 - grid->z0); +#endif + + rmin = grid->r0 * dr[grid->level]; + rmax = grid->r1 * dr[grid->level]; + zmin = grid->z0 * dz[grid->level]; + zmax = grid->z1 * dz[grid->level]; + + debug (3, "HSTCYL/HSTCRT in a grid %d x %d\n", grid->r1 - grid->r0, + grid->z1 - grid->z0); + + get_bound_cond (grid, prob, &r_bound_cond, &z_bound_cond, lambda); + +#ifndef TRUE2D + fish_hstcyl (rmin, rmax, grid->r1 - grid->r0, + r_bound_cond, boundaries[BND_LEFT], boundaries[BND_RIGHT], + zmin, zmax, grid->z1 - grid->z0, + z_bound_cond, boundaries[BND_BOTTOM], boundaries[BND_TOP], + lambda, s, + RZP (grid->phi, grid->r0, grid->z0), + grid->phi->strides[Z_INDX]); +#else + fish_hstcrt (rmin, rmax, grid->r1 - grid->r0, + r_bound_cond, boundaries[BND_LEFT], boundaries[BND_RIGHT], + zmin, zmax, grid->z1 - grid->z0, + z_bound_cond, boundaries[BND_BOTTOM], boundaries[BND_TOP], + lambda, + RZP (grid->phi, grid->r0, grid->z0), + grid->phi->strides[Z_INDX]); +#endif + debug (3, "Finished HSTCYL/HSTCRT in a grid %d x %d\n", grid->r1 - grid->r0, + grid->z1 - grid->z0); + + pois_set_phi_boundaries (grid, boundaries, + (grid->ext_bound & BND_MASK (BND_LEFT)) + && 0. == lambda, + + prob->bnd_right == BND_CND_HNEUMANN + && (grid->ext_bound & BND_MASK (BND_RIGHT)), + + prob->bnd_bottom == BND_CND_HNEUMANN + && (grid->ext_bound & BND_MASK (BND_BOTTOM)), + + prob->bnd_top == BND_CND_HNEUMANN + && (grid->ext_bound & BND_MASK (BND_TOP))); + + for (i = 0; i < 4; i++) { + debug (3, "boundaries[%d] = %p\n", i, boundaries[i]); + free (boundaries[i]); + } + debug (3, " <- pos_solve_grid (...)\n"); +} + +/** @brief Gets the boundary conditions that we pass to FISHPACK. + * + * Note that we use the same conditions for the Poisson equation + * and for the Helmholtz equation used to find the photoionization + * source. From my point of view, it doesn't make sense to use different + * conditions, so that will simply add an unneccesary complication. + */ +static void +get_bound_cond (pois_grid_t *grid, pois_problem_t *prob, + int *r_bound_cond, int *z_bound_cond, double lambda) +{ + static int right_inside[] = {FISH_UNS_DIR, FISH_NEU_DIR, FISH_DIR_DIR}; + static int right_ext_neu[] = {FISH_UNS_NEU, FISH_NEU_NEU, FISH_DIR_NEU}; + static int *right_ext_dir = right_inside; + static int fish_order[4] = { FISH_DIR_DIR, FISH_DIR_NEU, + FISH_NEU_DIR, FISH_NEU_NEU}; + + int bottom, top; + int *right_sel; + + if ((grid->ext_bound & BND_MASK (BND_TOP)) + && prob->bnd_top == BND_CND_HNEUMANN) + top = 1; + else top = 0; + + if ((grid->ext_bound & BND_MASK (BND_BOTTOM)) + && prob->bnd_bottom == BND_CND_HNEUMANN) + bottom = 1; + else bottom = 0; + + *z_bound_cond = fish_order[(bottom << 1) + top]; + + if ((grid->ext_bound & BND_MASK (BND_RIGHT))) { + if (prob->bnd_right == BND_CND_HNEUMANN) + right_sel = right_ext_neu; + else + right_sel = right_ext_dir; + } else { + right_sel = right_inside; + } + +#ifndef TRUE2D + *r_bound_cond = ((grid->ext_bound & BND_MASK (BND_LEFT)) && 0. == lambda)? + right_sel[0]: right_sel[2]; +#else + *r_bound_cond = ((grid->ext_bound & BND_MASK (BND_LEFT))? + right_sel[1]: right_sel[2]); +#endif + +} + +/** @brief Uses the @a boundaries[] vectors to set the boundaries by + * interpolation for the phi array of a grid. + * + * (in the extra space we allocated when we created the array). + * xxx_neu tells us if the xxx boundary has Neumann b.c. + * (or unspecified, as FISHPACK calls them if they are applied in the axis) + * + * TODO: Please rewrite this mess. + */ +void +pois_set_phi_boundaries (pois_grid_t *grid, REAL *boundaries[], + int left_neu, int right_neu, int bottom_neu, + int top_neu) +{ + int i; + + debug (3, "pois_set_phi_boundaries (...)\n"); + + if ((grid->ext_bound & BND_MASK(BND_BOTTOM)) && bottom_neu) { + for (i = grid->r0; i < grid->r1; i++) { + RZ (grid->phi, i, grid->z0 - 1) = RZ (grid->phi, i, grid->z0); + } + } else { + for (i = grid->r0; i < grid->r1; i++) { + RZ (grid->phi, i, grid->z0 - 1) = + 2 * boundaries[BND_BOTTOM][i - grid->r0] - RZ (grid->phi, i, grid->z0); + } + } + + if ((grid->ext_bound & BND_MASK(BND_TOP)) && top_neu) { + for (i = grid->r0; i < grid->r1; i++) { + RZ (grid->phi, i, grid->z1) = RZ (grid->phi, i, grid->z1 - 1); + } + } else { + for (i = grid->r0; i < grid->r1; i++) { + RZ (grid->phi, i, grid->z1) = 2 * boundaries[BND_TOP][i - grid->r0] + - RZ(grid->phi, i, grid->z1 - 1); + } + } + + if (left_neu) { + for (i = grid->z0; i < grid->z1; i++) { + RZ (grid->phi, grid->r0 - 1, i) = RZ(grid->phi, grid->r0, i); + } + } else { + for (i = grid->z0; i < grid->z1; i++) { + RZ(grid->phi, grid->r0 - 1, i) = 2 * boundaries[BND_LEFT][i - grid->z0] + - RZ(grid->phi, grid->r0, i); + } + } + + if (right_neu) { + for (i = grid->z0; i < grid->z1; i++) { + RZ (grid->phi, grid->r1, i) = RZ(grid->phi, grid->r1 - 1, i); + } + } else { + for (i = grid->z0; i < grid->z1; i++) { + RZ(grid->phi, grid->r1, i) = 2 * boundaries[BND_RIGHT][i - grid->z0] + - RZ(grid->phi, grid->r1 - 1, i); + } + } + + /* This is quite nightmarish. Think about rewritting. */ + if ((grid->ext_bound & BND_MASK(BND_RIGHT)) && right_neu) { + RZ(grid->phi, grid->r1, grid->z0 - 1) = RZ(grid->phi, grid->r1 - 1, + grid->z0 - 1); + RZ(grid->phi, grid->r1, grid->z1) = RZ(grid->phi, grid->r1 - 1, + grid->z1); + } else { + RZ(grid->phi, grid->r1, grid->z0 - 1) = + RZ(grid->phi, grid->r1 - 1, grid->z0 - 1) + + RZ(grid->phi, grid->r1, grid->z0) + - 0.5 * (RZ(grid->phi, grid->r1 - 2, grid->z0 - 1) + + RZ(grid->phi, grid->r1, grid->z0 + 1)); + + RZ(grid->phi, grid->r1, grid->z1) = + RZ(grid->phi, grid->r1 - 1, grid->z1) + + RZ(grid->phi, grid->r1, grid->z1 - 1) + - 0.5 * (RZ(grid->phi, grid->r1 - 2, grid->z1) + + RZ(grid->phi, grid->r1, grid->z1 - 2)); + } + + if ((grid->ext_bound & BND_MASK(BND_LEFT)) && left_neu) { + RZ(grid->phi, grid->r0 - 1, grid->z0 - 1) = RZ(grid->phi, grid->r0, + grid->z0 - 1); + RZ(grid->phi, grid->r0 - 1, grid->z1) = RZ(grid->phi, grid->r0, + grid->z1); + } else { + RZ(grid->phi, grid->r0 - 1, grid->z0 - 1) = + RZ(grid->phi, grid->r0, grid->z0 - 1) + + RZ(grid->phi, grid->r0 - 1, grid->z0) + - 0.5 * (RZ(grid->phi, grid->r0 + 1, grid->z0 - 1) + + RZ(grid->phi, grid->r0 - 1, grid->z0 + 1)); + + RZ(grid->phi, grid->r0 - 1, grid->z1) = + RZ(grid->phi, grid->r0, grid->z1) + + RZ(grid->phi, grid->r0 - 1, grid->z1 - 1) + - 0.5 * (RZ(grid->phi, grid->r0 + 1, grid->z1) + + RZ(grid->phi, grid->r0 - 1, grid->z1 - 2)); + } +} + +/** @brief Estimates the error on this grid by comparing with the calculations + * for his parent and performing a Richardson extrapolation. + */ +void +pois_set_error (pois_grid_t *grid) +{ + pois_grid_t *parent; + interpol_t *interpol; + int pr, pz; + + debug (3, "pois_set_error(...)\n"); + + parent = grid->parent; + assert (NULL != parent); + + interpol = interpol_new_a (2.0, 2.0, &interpol_wackers); + + iter_grid_parent (grid, pr, pz) { + int i, j; + + interpol_set_stencil_at ((grid_t *) grid, interpol, + pr * 2.0 + 1, pz * 2.0 + 1, + parent->phi, pr, pz, 0); + for (i = 0; i < 2; i++) + for (j = 0; j < 2; j++) { + int cr, cz; + double v; + + cr = (pr << 1) + i; + cz = (pz << 1) + j; + + v = interpol_apply (interpol, (double) cr + 0.5, (double) cz + 0.5); + + /* The error will only be used as its abs value. Therefore + we calculate it here and forever. */ + RZ (grid->error, cr, cz) = fabs (RZ (grid->phi, cr, cz) - v); + + } + } + interpol_free (interpol); +} + +/** @brief The refinement routine. + * + * Threshold is initially pois_max_error, but if we get too large refinement + * it can be increased (after warning the user that we are not reaching + * the accuracy that he asked for). + */ +int +pois_refine (pois_grid_t *grid, double threshold) +{ + int ir, iz; + int zmin, zmax, rmin, rmax; + /* Be careful: cz1 and cr1 are inclusive: when the child is actually + created, one has to increse them by 1. */ + int cr0, cr1, cz0, cz1; + int all_ok = TRUE; + + int tainted_line = FALSE; + + debug (3, "pois_refine(...)\n"); + + /* Initially we set values out of the grid, so the first comparisons + * always change the values. + */ + cr0 = grid->r1; cr1 = grid->r0 - 1; + cz0 = grid->z1; cz1 = grid->z0 - 1; + + /* The boundary can only coincide with the parent's boundary if this + * is an external boundary. + */ + rmin = (grid->ext_bound & BND_MASK(BND_LEFT))? grid->r0: (grid->r0 + 1); + rmax = (grid->ext_bound & BND_MASK(BND_RIGHT))? grid->r1: (grid->r1 - 1); + zmin = (grid->ext_bound & BND_MASK(BND_BOTTOM))? grid->z0: (grid->z0 + 1); + zmax = (grid->ext_bound & BND_MASK(BND_TOP))? grid->z1: (grid->z1 - 1); + + for (iz = zmin; iz < zmax && all_ok; iz++) { + tainted_line = FALSE; + for (ir = rmin; ir < rmax; ir++) { + if (needs_refinement (grid, ir, iz, threshold)) { + tainted_line = TRUE; + if (ir < cr0) cr0 = ir; + if (ir > cr1) cr1 = ir; + } + } + if (tainted_line) { + if (iz < cz0) cz0 = iz; + if (iz > cz1) cz1 = iz; + } else if (cr1 >= cr0 && cz1 >= cz0) { + all_ok = all_ok && refine_in (grid, cr0, cz0, cr1, cz1); + cr0 = grid->r1; cr1 = grid->r0 - 1; + cz0 = grid->z1; cz1 = grid->z0 - 1; + } + } + + if (cr1 >= cr0 && cz1 >= cz0) { + all_ok = all_ok && refine_in (grid, cr0, cz0, cr1, cz1); + } + + /* If there was some error with some refinement, we delete the created + * sub-grids and go back to the caller, which has to increase the refinement + * threshold. + */ + if (!all_ok) { + pois_grid_t *leaf; + free_childs (grid, leaf, pois_free_r); + set_childless (grid); + } + return all_ok; +} + +#define REFINE_IN_MAX_WARN_CNT 20; + +/** @brief Refines and creates a new grid according to the coordinates given. + * + * (if the FISHPACK limit allows it). Returns TRUE if the refinement + * succeeded, FALSE otherwise (if the FISHPACK limit was exceeded). +*/ +static int +refine_in (pois_grid_t *grid, int cr0, int cz0, int cr1, int cz1) +{ + pois_grid_t *child; + int nr0, nr1, nz0, nz1; + static int warn_cnt = REFINE_IN_MAX_WARN_CNT; + + nr0 = cr0 << 1; + nz0 = cz0 << 1; + nr1 = (++cr1) << 1; + nz1 = (++cz1) << 1; + + if ((nr1 - nr0) > FISH_MAX_GRIDPOINTS || + (nz1 - nz0) > FISH_MAX_GRIDPOINTS) { + if (warn_cnt > 0) { + warning ("FISHPACK limit exceeded. Not refining grid " + grid_printf_str ".\n", nr0, nz0, nr1, nz1, grid->level + 1); + warn_cnt --; + } + if (warn_cnt == 0) { + warning ("No more warnings about FISHPACK limit.\n"); + warn_cnt --; + } + return FALSE; + } + + child = pois_new_a (nr0, nz0, nr1, nz1); + + add_child (grid, child); + grid_inherit_ext_bound ((grid_t *) child); + 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); + + return TRUE; +} + +/** @brief Calculates the error measure of this grid. + * + * For debug/testing purposes + */ +void +pois_error_measures (pois_grid_t *grid, double *L1, double *L2, double *Lmax) +{ + int ir, iz, i; + double err; + *L2 = 0.0; + *L1 = 0.0; + *Lmax = 0.0; + i = 0; + + iter_grid (grid, ir, iz) { + i++; + err = RZ (grid->error, ir, iz); + if (err > *Lmax) *Lmax = err; + *L1 += err; + *L2 += err * err; + } + + *L1 /= i; + *L2 = sqrt (*L2 / i); + +} + +/** @brief Writes some error measures of this grid and its descendants + * into @a fp + */ +void +pois_write_error_r (pois_grid_t *grid, FILE *fp) +{ + double L1, L2, Lmax; + pois_grid_t *child; + + pois_error_measures (grid, &L1, &L2, &Lmax); + + fprintf (fp, "%d %g %g %g %g %g\n", grid->level, dr[grid->level], + dz[grid->level], L1, L2, Lmax); + + iter_childs (grid, child) { + pois_write_error_r (child, fp); + } +} + +/** @brief Takes a @a cdr tree and solves the Poisson equation for it. + * + * Returns a poisson tree. + */ +pois_grid_t ** +pois_solve_a (cdr_grid_t *cdr, pois_problem_t *prob) +{ + /* Call to the generic Poisson/Helmholtz solver */ + return pois_gen_solve_a (cdr, prob, e_mappers, 0.0); +} + +/** @brief Solves a Poisson/Helmholtz equation with FISHPACK and maps + * the result using a given set of mapper. + * + * From this method downwards (in the sense of calling order, not file + * organization), all the code is shared by the Poisson solver and the + * Photoionization (Helmholtz) solver. + */ +pois_grid_t ** +pois_gen_solve_a (cdr_grid_t *cdr, pois_problem_t *prob, + mapper_t **mappers, double es) +{ + cdr_grid_t *tree; + pois_grid_t **pois_modes; + int mode; + + pois_modes = (pois_grid_t**) xmalloc (sizeof(pois_grid_t*) * max_ntheta); + assert (0 == cdr->level); + tree = cdr_add_coarser_grids_a (cdr, prob->extra_levels); + + /* This is the high-level parallelizable loop! */ +#pragma omp parallel + { +#pragma omp for schedule(dynamic) + for (mode = 0; mode < max_ntheta; mode++) { + pois_modes[mode] = pois_solve_mode (tree, cdr, prob, mode, es); + + if (pois_inhom && es == 0.0 && mode == 0) { + /* Add the inhomogeneous laplacian potential to phi. + * It has only to be added to the zero mode. + */ + double q_factor; + q_factor = pois_inhom_q_factor (pois_modes[0]); + + /* Was: + pois_add_inhom_phi_r (pois_modes[0], q_factor); + */ + pois_add_inhom_phi_r (pois_modes[0], -q_factor); + } + + /* For the interpolation, we require an extra buffer cell + * in some boundaries and we check for overflows inside the + * mapper methods. + */ + map_trees_r (mappers, (grid_t*) pois_modes[mode], (grid_t*) cdr, + mode, FALSE, TRUE, FALSE, + /*s_buf = */ 2, + /*t_buf = */ 4); + map_trees_r (mappers, (grid_t*) pois_modes[mode], (grid_t*) cdr, + mode, TRUE, FALSE, TRUE, + /*s_buf = */ 1, + /*t_buf = */ 2); + } + } + cdr_free_coarser_grids (tree, extra_pois_levels); + + return pois_modes; +} + +/** @brief Solves a single Fourier mode. + * + * @a tree is the root of the @a cdr tree with extra_pois_levels added. + * @a cdr is the cdr grid at level 0, contained in tree: + * + * cdr = tree->first_child->firts_child...->first_child + * \_________ extra_pois_levels __________/ + */ +pois_grid_t * +pois_solve_mode (cdr_grid_t *tree, cdr_grid_t *cdr, pois_problem_t *prob, + int mode, double es) +{ + pois_grid_t *pois; + + pois = pois_init_tree_a (0, 0, pois_gridpoints_r, pois_gridpoints_z); + + /* Solves the coarsest Poisson grid. */ + /* Changed Wed Mar 19 13:50:23 2008: was s_buf = 1. */ + map_trees_r (charge_mappers, (grid_t *) tree, (grid_t *) pois, + mode, TRUE, TRUE, FALSE, 0, 2); + pois_solve_grid (pois, prob, -w2k[mode], es); + + /* Wed May 14 11:00:30 2008: + * I divide pois_max_error by abs(wk[mode]) because in the calculation + * of e_theta, the errors get multiplied by wk[mode]. Hence if we + * want to get similar errors in all modes, we have to remove it from the + * threshold. The 1 is because wk[0] = 0, but we do not want to have + * a \inf threshold. + */ + pois_solve_r (pois->first_child, tree, prob, mode, es, + prob->max_error / (1 + abs(wk[mode]))); + + return pois; +} + +#define POIS_SOLVE_R_MAX_WARN 50 + +/** @brief Recursively solves one Poisson grid. */ +void +pois_solve_r (pois_grid_t *pois, cdr_grid_t *cdr, pois_problem_t *prob, + int mode, double es, double threshold) +{ + pois_grid_t *child; + int succeeded = FALSE; + static int warn_cnt = POIS_SOLVE_R_MAX_WARN; + + /* Changed Wed Mar 19 13:51:00 2008: was s_buf = 1 */ + map_grid_r (charge_mappers, (grid_t *) cdr, (grid_t *) pois, + mode, TRUE, TRUE, FALSE, 0, 2); + debug(3, "w2k[mode = %d] = %f\n", mode, w2k[mode]); + + pois_solve_grid (pois, prob, -w2k[mode], es); + + if (pois->level >= prob->max_level) return; + + pois_set_error (pois); + + while (!succeeded) { + succeeded = pois_refine (pois, threshold); + + if (!succeeded) { + threshold *= 2; + if (warn_cnt > 0) { + warning ("Poisson error threshold was too small: trying %g.\ + Consider increasing pois_max_error\n", + threshold); + if (warn_cnt >= 0) warn_cnt --; + } + } + } + iter_childs (pois, child) { + pois_solve_r (child, cdr, prob, mode, es, threshold); + } +} + +/** @brief Returns 1 if the given cell has to be refined. + * + * which means that it satisfies at least one of the following: + * a. The error in the cell is larger than the threshold. + * b. The cell has a neightbour with an error larger than the threshold + */ +static int +needs_refinement (pois_grid_t *grid, int ir, int iz, double threshold) +{ + int i, j; + int r = FALSE; + int irmin, irmax, izmin, izmax; + + irmin = ir > grid->r0 + 3? ir - 1: grid->r0; + irmax = ir < grid->r1 - 4? ir + 1: grid->r1 - 1; + izmin = iz > grid->z0 + 3? iz - 1: grid->z0; + izmax = iz < grid->z1 - 4? iz + 1: grid->z1 - 1; + + for (i = irmin; i <= irmax && !r; i++) + for (j = izmin; j <= izmax; j++) + if (RZ(grid->error, i, j) > threshold) { + r = TRUE; + break; + } + + return r; +} + +/***************************************************************************** + * Functions for the implemetation of a needle-plate discharge. * + * The method used here is that of the "floating charge": * + * A charge is located at pois_inhom_z, which is used to keep the potential * + * at L_z constant. * + ****************************************************************************/ + +/** @brief pois_inhom_init QQQQ */ +void +pois_inhom_init (void) +{ + /* The domain of solution of the Poisson equation will be different + * from that of the CDR equation. The difference between them is + * given by the parameter needle_length, but has to be rounded + * to a multiple of the grid size at level 0. + */ + pois_gridpoints_z = (gridpoints_z + + (int) nearbyint (needle_length / dz[0])); + pois_inhom_L = pois_gridpoints_z * dz[0]; + + pois_inhom_z = L_z + needle_radius; + pois_inhom_dphi = pois_inhom_phi (0, pois_inhom_L) - pois_inhom_phi (0, L_z); + +} + +/* We use the following functions to calculate the fields and potential + * created by a unit charge located at (r = 0, z = b). + * These are created by the method of mirror charges, with a total of + * pois_inhom_reflections charges. + */ + +/* Single term of phi. */ +static double +inhom_phi_term (double r, double z, double b) +{ + double insqrt; + + insqrt = (SQ(z - b) + SQ(r)); + + return invfourpi / sqrt (insqrt); +} + +/** @brief Common decay factor for the electric fields. */ +static double +inhom_e_term (double r, double z, double b) +{ + double insqrt; + + insqrt = (SQ(z - b) + SQ(r)); + insqrt = insqrt * insqrt * insqrt; + + return invfourpi / sqrt (insqrt); +} + +/** @brief Single term of E_r. */ +static double +inhom_er_term (double r, double z, double b) +{ + return r * inhom_e_term (r, z, b); +} + +/** @brief Single term of E_z */ +static double +inhom_ez_term (double r, double z, double b) +{ + return (z - b) * inhom_e_term (r, z, b); +} + +/* Since we have to sum over the same locations of the mirror charges + * but with different "field functions", we use the following macro. + * + * Note that using a single function for inhom_phi, inhom_er, and inhom_ez + * would not be a good idea since they are evaluated at different points. + * Passing that function a pointer to another function would have been somewhat + * slower and, I think, not so readable. + */ +#define sum_reflections(total_, field_, r_, z_) \ + do { \ + int i__; \ + total_ = (field_ (r_, z_, pois_inhom_z) - \ + field_ (r_, z_, -pois_inhom_z)); \ + \ + for (i__ = 2; i__ <= pois_inhom_reflections; i__ += 2) { \ + total_ += (field_ (r_, z_, pois_inhom_z + i__ * pois_inhom_L) + \ + field_ (r_, z_, pois_inhom_z - i__ * pois_inhom_L)); \ + \ + total_ -= (field_ (r_, z_, -pois_inhom_z + i__ * pois_inhom_L) + \ + field_ (r_, z_, -pois_inhom_z - i__ * pois_inhom_L)); \ + } \ + } while (0) + +/** @brief Returns the potential created by a unit charge located at pois_inhom_z + * between an electrode at 0 and a second one at pois_inhom_L. + */ +double +pois_inhom_phi (double r, double z) +{ + double res; + + sum_reflections (res, inhom_phi_term, r, z); + + return res; +} + +/** @brief Radial component of the field. */ +double +pois_inhom_er (double r, double z) +{ + double res; + + sum_reflections (res, inhom_er_term, r, z); + + return res; +} + +/** @brief Axial component of the field. */ +double +pois_inhom_ez (double r, double z) +{ + double res; + + sum_reflections (res, inhom_ez_term, r, z); + + return res; +} + +/** @brief Once we have the electrostatic potential created by the space charges, + * we use this function to compute the multiplying factor of the floating + * charge. + * + * One can also forget about keeping the potential fixed somewhere and impose + * a constant (non-floating, but may depend on time) charge pois_inhom_fixed_q, + * which is used if nonzero. The use of that is to simulate a charged cloud + * close to the earth, which acts as an electrode: usually you would also combine + * that with the sprite module, that allows varing neutral densities. + */ +double +pois_inhom_q_factor (pois_grid_t *pois) +{ + double u; + /* Note that the signs preceding the numerical potentials are inverted since + * in our notation, we use the opposite of the physical one, in order + * to use the charge as the source of the Poisson equation. + */ + + if (pois_inhom_fixed_q != 0.0) { + return pois_inhom_fixed_q_t; + } + + /* Was: + u = - E_z * (pois_inhom_L - L_z) + pois_phi_at (pois, 0.0, L_z, 0.0); + */ + + u = E_z * (L_z - pois_inhom_L) + pois_phi_at (pois, 0.0, L_z, 0.0); + + return -u / pois_inhom_dphi; +} + +/** @brief Adds the potential created by the needle to a Poisson grid and their + * descendants. + * + * Note that \f$\phi\f$ here is in Fourier space and the Laplacian + * potential is always axi-symmetric, so you only have to call this function + * with the zero-mode of the Poisson grids. + * + * The reason that it is better to add the inhomogeneous field to the potential + * and not later to the electric fields is that if we do that, we will + * have to substract two large and very curved functions that give something + * close to zero: the error we make then will be comparable to the result. + */ +void +pois_add_inhom_phi_r (pois_grid_t *grid, double q) +{ + int ir, iz; + pois_grid_t *child; + double res; + + debug (2, "pois_add_inhom_phi_r (" grid_printf_str ", q = %f)\n", + grid_printf_args(grid), q); + + /* Note that there is no factor grid->ntheta coming from thr FFT + un-normalization. The reason is that actually q is sqrt(N) smaller + than it has to be and later, in the inverse FFT we will multiply by + sqrt(N). So everything is correct, I believe. */ + iter_grid_n (grid, ir, iz, 2) { + + sum_reflections (res, inhom_phi_term, r_at(ir, grid->level), z_at(iz, grid->level)); + + + RZ(grid->phi, ir, iz) += q * res; + } + + iter_childs (grid, child) { + pois_add_inhom_phi_r (child, q); + } +} + +/*************************************************************** + * Mapper functions for the components of the electric field. * + * See mapper.h for an explanation of each of these methods. * + ***************************************************************/ + +/** @brief er_copy QQQQ */ +void +er_copy (mapper_t *mapper, grid_t *source, grid_t *target, + int ir, int iz, int itheta) +{ + cdr_grid_t *cdr; + pois_grid_t *pois; + + cdr = (cdr_grid_t*) target; + pois = (pois_grid_t*) source; + + if (ir < pois->r1) + RZT (cdr->er, ir, iz, itheta) = ER_RZ (pois, ir, iz); +} + +/** @brief ez_copy QQQQ */ +void +ez_copy (mapper_t *mapper, grid_t *source, grid_t *target, + int ir, int iz, int itheta) +{ + cdr_grid_t *cdr; + pois_grid_t *pois; + + cdr = (cdr_grid_t*) target; + pois = (pois_grid_t*) source; + + if (iz < pois->z1) + RZT (cdr->ez, ir, iz, itheta) = EZ_RZ (pois, ir, iz); +} + +/** @brief etheta_copy QQQQ */ +void +etheta_copy (mapper_t *mapper, grid_t *source, grid_t *target, + int ir, int iz, int itheta) +{ + cdr_grid_t *cdr; + pois_grid_t *pois; + + cdr = (cdr_grid_t*) target; + pois = (pois_grid_t*) source; + + RZT (cdr->etheta, ir, iz, itheta) = RZ (pois->phi, ir, iz); +} + +/** @brief er_coarsen QQQQ */ +void +er_coarsen (mapper_t *mapper, grid_t *source, grid_t *target, + int ir, int iz, int itheta) +{ + cdr_grid_t *cdr; + pois_grid_t *pois; + int level_diff, er_ir, er_iz; + + cdr = (cdr_grid_t*) target; + pois = (pois_grid_t*) source; + + level_diff = pois->level - cdr->level; + + er_iz = (iz << level_diff) + (1 << (level_diff - 1)); + er_ir = ((ir + 1) << level_diff) - 1; + + if (grid_contains (source, er_ir, er_iz - 1, GRID_INSIDE) && + grid_contains (source, er_ir, er_iz, GRID_INSIDE)) { + + RZT (cdr->er, ir, iz, itheta) = 0.5 * (ER_RZ(pois, er_ir, er_iz) + + ER_RZ(pois, er_ir, er_iz - 1)); + } +} + +/** @brief ez_coarsen QQQQ */ +void +ez_coarsen (mapper_t *mapper, grid_t *source, grid_t *target, + int ir, int iz, int itheta) +{ + cdr_grid_t *cdr; + pois_grid_t *pois; + int level_diff, ez_ir, ez_iz; + + cdr = (cdr_grid_t*) target; + pois = (pois_grid_t*) source; + + level_diff = pois->level - cdr->level; + + ez_iz = ((iz + 1) << level_diff) - 1; + ez_ir = (ir << level_diff) + (1 << (level_diff - 1)); + + if (grid_contains (source, ez_ir - 1, ez_iz, GRID_INSIDE) && + grid_contains (source, ez_ir, ez_iz, GRID_INSIDE)) { + + RZT(cdr->ez, ir, iz, itheta) = 0.5 * (EZ_RZ(pois, ez_ir, ez_iz) + + EZ_RZ(pois, ez_ir - 1, ez_iz)); + } +} + +/** @brief etheta_coarsen QQQQ */ +void +etheta_coarsen (mapper_t *mapper, grid_t *source, grid_t *target, + int ir, int iz, int itheta) +{ + cdr_grid_t *cdr; + pois_grid_t *pois; + int level_diff, ez_ir, er_iz; + + cdr = (cdr_grid_t*) target; + pois = (pois_grid_t*) source; + + level_diff = pois->level - cdr->level; + + er_iz = (iz << level_diff) + (1 << (level_diff - 1)); + ez_ir = (ir << level_diff) + (1 << (level_diff - 1)); + + if (grid_contains (source, ez_ir, er_iz, GRID_INSIDE) && + grid_contains (source, ez_ir - 1, er_iz, GRID_INSIDE) && + grid_contains (source, ez_ir, er_iz - 1, GRID_INSIDE) && + grid_contains (source, ez_ir - 1, er_iz - 1, GRID_INSIDE)) { + + RZT (cdr->etheta, ir, iz, itheta) = 0.25 * + (RZ (pois->phi, ez_ir, er_iz) + + RZ (pois->phi, ez_ir - 1, er_iz) + + RZ (pois->phi, ez_ir, er_iz - 1) + + RZ (pois->phi, ez_ir - 1, er_iz - 1)); + } +} + +/** @brief er_interpol QQQQ */ +int +er_interpol_set (mapper_t *mapper, grid_t *source, interpol_t *interpol, + int pr, int pz, int itheta) +{ + pois_grid_t *pois; + + pois = (pois_grid_t*) source; + + /* When we set Neumann b.c. in one boundary and Dirichlet in the opposite + * a larger-than-usual error appears in the boundaries because we are + * matching potetntials in one side and fields in the other. The only way + * I see of avoiding this is to forget about the extrapolated values of + * the grids. But we still need those values on the external boundaries. + */ + if (pr < pois->r0 || pz < pois->z0 + || pr > pois->r1 - 1 || pz > pois->z1) + return FALSE; + + interpol_set_stencil (interpol, + er_r_at (pr - 1, pois->level), + er_z_at (pz - 1, pois->level), + ER_RZ (pois, pr - 1, pz - 1), + ER_RZ (pois, pr - 1, pz), + ER_RZ (pois, pr, pz - 1), + ER_RZ (pois, pr, pz)); + + return TRUE; +} + +/** @brief ez_interpol_set QQQQ */ +int +ez_interpol_set (mapper_t *mapper, grid_t *source, interpol_t *interpol, + int pr, int pz, int itheta) +{ + pois_grid_t *pois; + + pois = (pois_grid_t*) source; + + if (pr < pois->r0 || pz < pois->z0 + || pr > pois->r1 || pz > pois->z1 - 1) + return FALSE; + + interpol_set_stencil (interpol, + ez_r_at (pr - 1, pois->level), + ez_z_at (pz - 1, pois->level), + EZ_RZ (pois, pr - 1, pz - 1), + EZ_RZ (pois, pr - 1, pz), + EZ_RZ (pois, pr, pz - 1), + EZ_RZ (pois, pr, pz)); + return TRUE; +} + +/** @brief etheta_interpol_set QQQQ */ +int +etheta_interpol_set (mapper_t *mapper, grid_t *source, interpol_t *interpol, + int pr, int pz, int itheta) +{ + pois_grid_t *pois; + + pois = (pois_grid_t*) source; + + if (pr < pois->r0 || pz < pois->z0 + || pr > pois->r1 - 1|| pz > pois->z1 - 1) + return FALSE; + + interpol_set_stencil_at (source, interpol, + r_at (pr, pois->level), + z_at (pz, pois->level), + pois->phi, pr, pz, 0); + return TRUE; +} + +/** @brief er_interpol QQQQ */ +void +er_interpol (mapper_t *mapper, grid_t *source, grid_t *target, + interpol_t *interpol, int ir, int iz, int itheta) +{ + double er_r, er_z; + cdr_grid_t *cdr; + pois_grid_t *pois; + int level_diff; + + cdr = (cdr_grid_t *) target; + pois = (pois_grid_t *) source; + + level_diff = cdr->level - pois->level; + + /* Note that when pr = pois->r1, the electric field is calculated + * from the extrapolation of phi and it is hence less accurate. + * The same applies below to pz and pois->z1. + */ + if ((ir >> level_diff) >= pois->r1) return; + + er_r = er_r_at (ir, cdr->level); + er_z = er_z_at (iz, cdr->level); + + RZT (cdr->er, ir, iz, itheta) = interpol_apply (interpol, er_r, er_z); +} + +/** @brief ez_interpol QQQQ */ +void +ez_interpol (mapper_t *mapper, grid_t *source, grid_t *target, + interpol_t *interpol, int ir, int iz, int itheta) +{ + double ez_r, ez_z; + cdr_grid_t *cdr; + pois_grid_t *pois; + int level_diff; + + cdr = (cdr_grid_t *) target; + pois = (pois_grid_t *) source; + + level_diff = cdr->level - pois->level; + + if ((iz >> level_diff) >= pois->z1) return; + + ez_r = ez_r_at (ir, cdr->level); + ez_z = ez_z_at (iz, cdr->level); + + RZT (cdr->ez, ir, iz, itheta) = interpol_apply (interpol, ez_r, ez_z); +} + +/** @brief etheta_interpol QQQQ */ +void +etheta_interpol (mapper_t *mapper, grid_t *source, grid_t *target, + interpol_t *interpol, int ir, int iz, int itheta) +{ + double eth_r, eth_z; + cdr_grid_t *cdr; + + cdr = (cdr_grid_t *) target; + + eth_r = r_at (ir, cdr->level); + eth_z = z_at (iz, cdr->level); + + RZT (cdr->etheta, ir, iz, itheta) = + interpol_apply (interpol, eth_r, eth_z); +} + +/** @brief Mapping of the charge. */ +void +charge_copy (mapper_t *mapper, grid_t *source, grid_t *target, + int ir, int iz, int itheta) +{ + cdr_grid_t *cdr; + pois_grid_t *pois; + + pois = (pois_grid_t*) target; + cdr = (cdr_grid_t*) source; + + RZ (pois->charge, ir, iz) = RZT (cdr->charge, ir, iz, itheta); +} + +/** @brief charge_interpol_set QQQQ. */ +int +charge_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; + + if (pr < cdr->r0 || pz < cdr->z0 + || pr > cdr->r1 - 1|| pz > cdr->z1 - 1) + return FALSE; + + interpol_set_stencil_at (source, interpol, + r_at (pr, cdr->level), + z_at (pz, cdr->level), + cdr->charge, pr, pz, itheta); + return TRUE; +} + +/** @brief charge_interpol QQQQ. */ +void +charge_interpol (mapper_t *mapper, grid_t *source, grid_t *target, + interpol_t *interpol, + int ir, int iz, int itheta) +{ + double r, z; + pois_grid_t *pois; + + pois = (pois_grid_t *) target; + + r = r_at (ir, pois->level); + z = z_at (iz, pois->level); + + RZT (pois->charge, ir, iz, 0) = interpol_apply (interpol, r, z); +} + +/** @brief Finds the best-possible approximation for the potential at a given + * point (@a r, @a z, @a theta) by interpolating from the finest grid that covers + * that point. + */ +double +pois_phi_at (pois_grid_t *grid, double r, double z, double theta) +{ + int ir, iz, itheta, it, itn; + REAL phi_it[2]; + interpol_t *phi_at_interpol = NULL; + + debug (3, "pois_phi_at (" grid_printf_str ", r = %g, z = %g, theta = %g)\n", + grid_printf_args (grid), r, z, theta); + + grid = (pois_grid_t *) grid_finest_containing_r ((grid_t*) grid, r, z); + + if (NULL == grid) { + fatal ("In pois_phi_at: potential outside the grid"); + } + + phi_at_interpol = interpol_new_a (dr[grid->level], dz[grid->level], + &interpol_quadratic); + + ir = (int) (r / dr[grid->level]); + iz = (int) (z / dz[grid->level]); + + /* If r == 0, it doesn't matter which theta we use, so we use 0. + * also if we are in a 2D simulation. + */ + if (grid->ntheta > 1 && r > 0.0) { + itheta = (int) (theta / dtheta); + itn = 2; + } else { + itheta = 0; + itn = 1; + } + + for (it = 0; it < itn; it++) { + interpol_set_stencil_at ((grid_t *) grid, phi_at_interpol, + r_at (ir, grid->level), + z_at (iz, grid->level), + grid->phi, ir, iz, itheta + it); + + phi_it[it] = interpol_apply (phi_at_interpol, r, z); + } + + interpol_free (phi_at_interpol); + + if (grid->ntheta > 1 && r > 0.0) { + return (phi_it[0] * (theta - theta_at (itheta)) + + phi_it[1] * (theta_at (itheta + 1) - theta)) / dtheta; + } else { + return phi_it[0]; + } +} + +/** @brief Dumps all the contents of the given grid into filenames given by + * prefix and name + */ +void +pois_dump (pois_grid_t *grid, const char *prefix, const char *name) +{ + char *fname; + int m = pois_output_margin; + + asprintf (&fname, "%s/r.%s.tsv", prefix, name); + rz_axis_dump (fname, grid->r0 - m, grid->r1 + m, dr[grid->level]); + free (fname); + + asprintf (&fname, "%s/z.%s.tsv", prefix, name); + rz_axis_dump (fname, grid->z0 - m, grid->z1 + m, dz[grid->level]); + free (fname); + + asprintf (&fname, "%s/phi.%s.tsv", prefix, name); + rz_dump (grid->phi, fname, "w", grid->r0 - m, grid->z0 - m, + grid->r1 + m, grid->z1 + m); + free (fname); + + asprintf (&fname, "%s/charge.%s.tsv", prefix, name); + rz_dump (grid->charge, fname, "w", grid->r0 - m, grid->z0 - m, + grid->r1 + m, grid->z1 + m); + free (fname); + + asprintf (&fname, "%s/error.%s.tsv", prefix, name); + rz_dump (grid->error, fname, "w", grid->r0 - m, grid->z0 - m, + grid->r1 + m, grid->z1 + m); + free (fname); +} + +/** @brief Dumps all the contents of the given grid into filenames given by + * prefix and name + */ +void +pois_dump_r (pois_grid_t *grid, const char *prefix, const char *name) +{ + pois_grid_t *child; + char *cname; + int i = 0, nchilds; + char codes[] = "abcdefghijklmnopqrstuvwxyz"; + + pois_dump (grid, prefix, name); + + nchilds = grid_howmany_children ((grid_t *) grid); + + for (i = 0; i < nchilds; i++) { + child = (pois_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]); + + pois_dump_r (child, prefix, cname); + free (cname); + } +}