/** @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); } }