Poisson/Helmholtz solver, including the routines for inhomogeneous (point-plane) configurations. More...

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#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"

Include dependency graph for poisson.c:

Go to the source code of this file.

Macros
#define	REFINE_IN_MAX_WARN_CNT 20;

#define	POIS_SOLVE_R_MAX_WARN 50

#define	sum_reflections(total_, field_, r_, z_)

Functions
static void	get_bound_cond (pois_grid_t grid, pois_problem_t prob, int r_bound_cond, int z_bound_cond, double lambda)
	Gets the boundary conditions that we pass to FISHPACK.

static int	needs_refinement (pois_grid_t *grid, int ir, int iz, double threshold)
	Returns 1 if the given cell has to be refined.

static int	refine_in (pois_grid_t *grid, int cr0, int cz0, int cr1, int cz1)
	Refines and creates a new grid according to the coordinates given.

	decl_mapper_funcs (er)

	decl_mapper_funcs (ez)

	decl_mapper_funcs (etheta)

	decl_mapper_funcs (charge)

static double	inhom_phi_term (double r, double z, double b)

static double	inhom_e_term (double r, double z, double b)
	Common decay factor for the electric fields.

static double	inhom_er_term (double r, double z, double b)
	Single term of E_r.

static double	inhom_ez_term (double r, double z, double b)
	Single term of E_z.

void	pois_init (void)
	All initialization of the Poisson solver has to be here.

pois_grid_t *	pois_new_a (int r0, int z0, int r1, int z1)
	Creates a 2D Poisson grid.

pois_grid_t *	pois_new_3d_a (int r0, int z0, int r1, int z1, int ntheta)
	Creates a new 3D Poisson grid.

void	pois_free (pois_grid_t *grid)
	Frees the memory allocated by pois_new_a.

void	pois_free_r (pois_grid_t *grid)
	Recursively frees a tree of Poisson grids.

pois_grid_t *	pois_new_glob_a (int r0, int z0, int r1, int z1, int level)
	Creates a new Poisson grid.

pois_grid_t *	pois_init_tree_a (int r0, int z0, int r1, int z1)
	Starts the tree of Poisson grids with the two coarsest ones.

REAL *	pois_boundary_a (pois_grid_t *grid, int boundary)
	Is it possible to use a single function to handle all four boundaries?

void	pois_solve_grid (pois_grid_t grid, pois_problem_t prob, double lambda, double s)
	Solves the Poisson equation, calling the FISHPACK routine.

void	pois_set_phi_boundaries (pois_grid_t grid, REAL boundaries[], int left_neu, int right_neu, int bottom_neu, int top_neu)
	Uses the boundaries[] vectors to set the boundaries by interpolation for the phi array of a grid.

void	pois_set_error (pois_grid_t *grid)
	Estimates the error on this grid by comparing with the calculations for his parent and performing a Richardson extrapolation.

int	pois_refine (pois_grid_t *grid, double threshold)
	The refinement routine.

void	pois_error_measures (pois_grid_t grid, double L1, double L2, double Lmax)
	Calculates the error measure of this grid.

void	pois_write_error_r (pois_grid_t grid, FILE fp)
	Writes some error measures of this grid and its descendants into fp.

pois_grid_t **	pois_solve_a (cdr_grid_t cdr, pois_problem_t prob)
	Takes a cdr tree and solves the Poisson equation for it.

pois_grid_t **	pois_gen_solve_a (cdr_grid_t cdr, pois_problem_t prob, mapper_t **mappers, double es)
	Solves a Poisson/Helmholtz equation with FISHPACK and maps the result using a given set of mapper.

pois_grid_t *	pois_solve_mode (cdr_grid_t tree, cdr_grid_t cdr, pois_problem_t *prob, int mode, double es)
	Solves a single Fourier mode.

void	pois_solve_r (pois_grid_t pois, cdr_grid_t cdr, pois_problem_t *prob, int mode, double es, double threshold)
	Recursively solves one Poisson grid.

void	pois_inhom_init (void)
	pois_inhom_init QQQQ

double	pois_inhom_phi (double r, double z)
	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_er (double r, double z)
	Radial component of the field.

double	pois_inhom_ez (double r, double z)
	Axial component of the field.

double	pois_inhom_q_factor (pois_grid_t *pois)
	Once we have the electrostatic potential created by the space charges, we use this function to compute the multiplying factor of the floating charge.

void	pois_add_inhom_phi_r (pois_grid_t *grid, double q)
	Adds the potential created by the needle to a Poisson grid and their descendants.

void	er_copy (mapper_t mapper, grid_t source, grid_t *target, int ir, int iz, int itheta)
	er_copy QQQQ

void	ez_copy (mapper_t mapper, grid_t source, grid_t *target, int ir, int iz, int itheta)
	ez_copy QQQQ

void	etheta_copy (mapper_t mapper, grid_t source, grid_t *target, int ir, int iz, int itheta)
	etheta_copy QQQQ

void	er_coarsen (mapper_t mapper, grid_t source, grid_t *target, int ir, int iz, int itheta)
	er_coarsen QQQQ

void	ez_coarsen (mapper_t mapper, grid_t source, grid_t *target, int ir, int iz, int itheta)
	ez_coarsen QQQQ

void	etheta_coarsen (mapper_t mapper, grid_t source, grid_t *target, int ir, int iz, int itheta)
	etheta_coarsen QQQQ

int	er_interpol_set (mapper_t mapper, grid_t source, interpol_t *interpol, int pr, int pz, int itheta)
	er_interpol QQQQ

int	ez_interpol_set (mapper_t mapper, grid_t source, interpol_t *interpol, int pr, int pz, int itheta)
	ez_interpol_set QQQQ

int	etheta_interpol_set (mapper_t mapper, grid_t source, interpol_t *interpol, int pr, int pz, int itheta)
	etheta_interpol_set QQQQ

void	er_interpol (mapper_t mapper, grid_t source, grid_t target, interpol_t interpol, int ir, int iz, int itheta)
	er_interpol QQQQ

void	ez_interpol (mapper_t mapper, grid_t source, grid_t target, interpol_t interpol, int ir, int iz, int itheta)
	ez_interpol QQQQ

void	etheta_interpol (mapper_t mapper, grid_t source, grid_t target, interpol_t interpol, int ir, int iz, int itheta)
	etheta_interpol QQQQ

void	charge_copy (mapper_t mapper, grid_t source, grid_t *target, int ir, int iz, int itheta)
	Mapping of the charge.

int	charge_interpol_set (mapper_t mapper, grid_t source, interpol_t *interpol, int pr, int pz, int itheta)
	charge_interpol_set QQQQ.

void	charge_interpol (mapper_t mapper, grid_t source, grid_t target, interpol_t interpol, int ir, int iz, int itheta)
	charge_interpol QQQQ.

double	pois_phi_at (pois_grid_t *grid, double r, double z, double theta)
	Finds the best-possible approximation for the potential at a given point (r, z, theta) by interpolating from the finest grid that covers that point.

void	pois_dump (pois_grid_t grid, const char prefix, const char *name)
	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)
	Dumps all the contents of the given grid into filenames given by prefix and name.

Variables
int	pois_gridpoints_z

int	pois_gridpoints_r

double	pois_inhom_L

double	pois_inhom_z

double	pois_inhom_dphi

double	pois_inhom_fixed_q_t

double	E_x

double	E_y

double	E_z

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, &etheta_mapper, NULL}

mapper_t	charge_mapper = mk_mapper_down(charge, &interpol_wackers)

mapper_t *	charge_mappers [] = {&charge_mapper, NULL}

pois_problem_t *	pois_electrostatic

Detailed Description

Poisson/Helmholtz solver, including the routines for inhomogeneous (point-plane) configurations.

Definition in file poisson.c.

Macro Definition Documentation

#define REFINE_IN_MAX_WARN_CNT 20;

Definition at line 668 of file poisson.c.

#define POIS_SOLVE_R_MAX_WARN 50

Definition at line 862 of file poisson.c.

#define sum_reflections	(	total_,
		field_,
		r_,
		z_
	)

Value:

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)

Definition at line 1006 of file poisson.c.

Function Documentation

static void get_bound_cond	(	pois_grid_t *	grid,
		pois_problem_t *	prob,
		int *	r_bound_cond,
		int *	z_bound_cond,
		double	lambda
	)

static

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.

Definition at line 412 of file poisson.c.

Here is the caller graph for this function:

static int needs_refinement	(	pois_grid_t *	grid,
		int	ir,
		int	iz,
		double	threshold
	)

static

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

Definition at line 909 of file poisson.c.

Here is the caller graph for this function:

static int refine_in	(	pois_grid_t *	grid,
		int	cr0,
		int	cz0,
		int	cr1,
		int	cz1
	)

static

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).

Definition at line 676 of file poisson.c.

Here is the call graph for this function:

Here is the caller graph for this function:

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.

static double inhom_phi_term	(	double	r,
		double	z,
		double	b
	)

static

Definition at line 963 of file poisson.c.

Here is the caller graph for this function:

static double inhom_e_term	(	double	r,
		double	z,
		double	b
	)

static

Common decay factor for the electric fields.

Definition at line 974 of file poisson.c.

Here is the caller graph for this function:

static double inhom_er_term	(	double	r,
		double	z,
		double	b
	)

static

Single term of E_r.

Definition at line 986 of file poisson.c.

Here is the call graph for this function:

Here is the caller graph for this function:

static double inhom_ez_term	(	double	r,
		double	z,
		double	b
	)

static

Single term of E_z.

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.

Definition at line 993 of file poisson.c.

Here is the call graph for this function:

Here is the caller graph for this function:

void pois_init ( void )

All initialization of the Poisson solver has to be here.

Definition at line 86 of file poisson.c.

Here is the call graph for this function:

Here is the caller graph for this function:

pois_grid_t* pois_new_a	(	int	r0,
		int	z0,
		int	r1,
		int	z1
	)

Creates a 2D Poisson grid.

Definition at line 109 of file poisson.c.

Here is the call graph for this function:

Here is the caller graph for this function:

pois_grid_t* pois_new_3d_a	(	int	r0,
		int	z0,
		int	r1,
		int	z1,
		int	ntheta
	)

Creates a new 3D Poisson grid.

r0, z0, r1, z1 are in GRID units

Definition at line 123 of file poisson.c.

Here is the call graph for this function:

Here is the caller graph for this function:

void pois_free ( pois_grid_t * grid )

Frees the memory allocated by pois_new_a.

Definition at line 163 of file poisson.c.

Here is the call graph for this function:

Here is the caller graph for this function:

void pois_free_r ( pois_grid_t * grid )

Recursively frees a tree of Poisson grids.

Definition at line 180 of file poisson.c.

Here is the call graph for this function:

Here is the caller graph for this function:

pois_grid_t* pois_new_glob_a	(	int	r0,
		int	z0,
		int	r1,
		int	z1,
		int	level
	)

Creates a new Poisson grid.

r0, r1... are in global units (i.e. independent of the grid level).

Definition at line 196 of file poisson.c.

Here is the call graph for this function:

Here is the caller graph for this function:

pois_grid_t* pois_init_tree_a	(	int	r0,
		int	z0,
		int	r1,
		int	z1
	)

Starts the tree of Poisson grids with the two coarsest ones.

Receives the grid dimensions at level 0.

Definition at line 217 of file poisson.c.

Here is the call graph for this function:

Here is the caller graph for this function:

REAL* pois_boundary_a	(	pois_grid_t *	grid,
		int	boundary
	)

Is it possible to use a single function to handle all four boundaries?

Maybe this should be rewritten to use the interpol2 module?

Definition at line 244 of file poisson.c.

Here is the call graph for this function:

Here is the caller graph for this function:

void pois_solve_grid	(	pois_grid_t *	grid,
		pois_problem_t *	prob,
		double	lambda,
		double	s
	)

Solves the Poisson equation, calling the FISHPACK routine.

Assumes that grid->charge is already set.

Definition at line 334 of file poisson.c.

Here is the call graph for this function:

Here is the caller graph for this function:

void pois_set_phi_boundaries	(	pois_grid_t *	grid,
		REAL *	boundaries[],
		int	left_neu,
		int	right_neu,
		int	bottom_neu,
		int	top_neu
	)

Uses the 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.

Definition at line 465 of file poisson.c.

Here is the caller graph for this function:

void pois_set_error ( pois_grid_t * grid )

Estimates the error on this grid by comparing with the calculations for his parent and performing a Richardson extrapolation.

Definition at line 561 of file poisson.c.

Here is the call graph for this function:

Here is the caller graph for this function:

int pois_refine	(	pois_grid_t *	grid,
		double	threshold
	)

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).

Definition at line 606 of file poisson.c.

Here is the call graph for this function:

Here is the caller graph for this function:

void pois_error_measures	(	pois_grid_t *	grid,
		double *	L1,
		double *	L2,
		double *	Lmax
	)

Calculates the error measure of this grid.

For debug/testing purposes

Definition at line 717 of file poisson.c.

Here is the caller graph for this function:

void pois_write_error_r	(	pois_grid_t *	grid,
		FILE *	fp
	)

Writes some error measures of this grid and its descendants into fp.

Definition at line 743 of file poisson.c.

Here is the call graph for this function:

Here is the caller graph for this function:

pois_grid_t** pois_solve_a	(	cdr_grid_t *	cdr,
		pois_problem_t *	prob
	)

Takes a cdr tree and solves the Poisson equation for it.

Returns a poisson tree.

Definition at line 763 of file poisson.c.

Here is the call graph for this function:

Here is the caller graph for this function:

pois_grid_t** pois_gen_solve_a	(	cdr_grid_t *	cdr,
		pois_problem_t *	prob,
		mapper_t **	mappers,
		double	es
	)

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.

Definition at line 777 of file poisson.c.

Here is the call graph for this function:

Here is the caller graph for this function:

pois_grid_t* pois_solve_mode	(	cdr_grid_t *	tree,
		cdr_grid_t *	cdr,
		pois_problem_t *	prob,
		int	mode,
		double	es
	)

Solves a single Fourier mode.

tree is the root of the cdr tree with extra_pois_levels added. cdr is the cdr grid at level 0, contained in tree:

cdr = tree->first_child->firts_child...->first_child _________ extra_pois_levels __________/

Definition at line 836 of file poisson.c.

Here is the call graph for this function:

Here is the caller graph for this function:

void pois_solve_r	(	pois_grid_t *	pois,
		cdr_grid_t *	cdr,
		pois_problem_t *	prob,
		int	mode,
		double	es,
		double	threshold
	)

Recursively solves one Poisson grid.

Definition at line 866 of file poisson.c.

Here is the call graph for this function:

Here is the caller graph for this function:

void pois_inhom_init ( void )

pois_inhom_init QQQQ

Definition at line 939 of file poisson.c.

Here is the call graph for this function:

Here is the caller graph for this function:

double pois_inhom_phi	(	double	r,
		double	z
	)

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.

Definition at line 1025 of file poisson.c.

Here is the call graph for this function:

Here is the caller graph for this function:

double pois_inhom_er	(	double	r,
		double	z
	)

Radial component of the field.

Definition at line 1036 of file poisson.c.

Here is the call graph for this function:

Here is the caller graph for this function:

double pois_inhom_ez	(	double	r,
		double	z
	)

Axial component of the field.

Definition at line 1047 of file poisson.c.

Here is the call graph for this function:

Here is the caller graph for this function:

double pois_inhom_q_factor ( pois_grid_t * pois )

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.

Definition at line 1067 of file poisson.c.

Here is the call graph for this function:

Here is the caller graph for this function:

void pois_add_inhom_phi_r	(	pois_grid_t *	grid,
		double	q
	)

Adds the potential created by the needle to a Poisson grid and their descendants.

Note that $\phi$ 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.

Definition at line 1101 of file poisson.c.

Here is the call graph for this function:

Here is the caller graph for this function:

void er_copy	(	mapper_t *	mapper,
		grid_t *	source,
		grid_t *	target,
		int	ir,
		int	iz,
		int	itheta
	)

er_copy QQQQ

Definition at line 1134 of file poisson.c.

void ez_copy	(	mapper_t *	mapper,
		grid_t *	source,
		grid_t *	target,
		int	ir,
		int	iz,
		int	itheta
	)

ez_copy QQQQ

Definition at line 1149 of file poisson.c.

void etheta_copy	(	mapper_t *	mapper,
		grid_t *	source,
		grid_t *	target,
		int	ir,
		int	iz,
		int	itheta
	)

etheta_copy QQQQ

Definition at line 1164 of file poisson.c.

void er_coarsen	(	mapper_t *	mapper,
		grid_t *	source,
		grid_t *	target,
		int	ir,
		int	iz,
		int	itheta
	)

er_coarsen QQQQ

Definition at line 1178 of file poisson.c.

Here is the call graph for this function:

void ez_coarsen	(	mapper_t *	mapper,
		grid_t *	source,
		grid_t *	target,
		int	ir,
		int	iz,
		int	itheta
	)

ez_coarsen QQQQ

Definition at line 1203 of file poisson.c.

Here is the call graph for this function:

void etheta_coarsen	(	mapper_t *	mapper,
		grid_t *	source,
		grid_t *	target,
		int	ir,
		int	iz,
		int	itheta
	)

etheta_coarsen QQQQ

Definition at line 1228 of file poisson.c.

Here is the call graph for this function:

int er_interpol_set	(	mapper_t *	mapper,
		grid_t *	source,
		interpol_t *	interpol,
		int	pr,
		int	pz,
		int	itheta
	)

er_interpol QQQQ

Definition at line 1258 of file poisson.c.

Here is the call graph for this function:

int ez_interpol_set	(	mapper_t *	mapper,
		grid_t *	source,
		interpol_t *	interpol,
		int	pr,
		int	pz,
		int	itheta
	)

ez_interpol_set QQQQ

Definition at line 1288 of file poisson.c.

Here is the call graph for this function:

int etheta_interpol_set	(	mapper_t *	mapper,
		grid_t *	source,
		interpol_t *	interpol,
		int	pr,
		int	pz,
		int	itheta
	)

etheta_interpol_set QQQQ

Definition at line 1311 of file poisson.c.

Here is the call graph for this function:

void er_interpol	(	mapper_t *	mapper,
		grid_t *	source,
		grid_t *	target,
		interpol_t *	interpol,
		int	ir,
		int	iz,
		int	itheta
	)

er_interpol QQQQ

Definition at line 1331 of file poisson.c.

Here is the call graph for this function:

void ez_interpol	(	mapper_t *	mapper,
		grid_t *	source,
		grid_t *	target,
		interpol_t *	interpol,
		int	ir,
		int	iz,
		int	itheta
	)

ez_interpol QQQQ

Definition at line 1358 of file poisson.c.

Here is the call graph for this function:

void etheta_interpol	(	mapper_t *	mapper,
		grid_t *	source,
		grid_t *	target,
		interpol_t *	interpol,
		int	ir,
		int	iz,
		int	itheta
	)

etheta_interpol QQQQ

Definition at line 1381 of file poisson.c.

Here is the call graph for this function:

void charge_copy	(	mapper_t *	mapper,
		grid_t *	source,
		grid_t *	target,
		int	ir,
		int	iz,
		int	itheta
	)

Mapping of the charge.

Definition at line 1398 of file poisson.c.

int charge_interpol_set	(	mapper_t *	mapper,
		grid_t *	source,
		interpol_t *	interpol,
		int	pr,
		int	pz,
		int	itheta
	)

charge_interpol_set QQQQ.

Definition at line 1412 of file poisson.c.

Here is the call graph for this function:

void charge_interpol	(	mapper_t *	mapper,
		grid_t *	source,
		grid_t *	target,
		interpol_t *	interpol,
		int	ir,
		int	iz,
		int	itheta
	)

charge_interpol QQQQ.

Definition at line 1432 of file poisson.c.

Here is the call graph for this function:

double pois_phi_at	(	pois_grid_t *	grid,
		double	r,
		double	z,
		double	theta
	)

Finds the best-possible approximation for the potential at a given point (r, z, theta) by interpolating from the finest grid that covers that point.

Definition at line 1452 of file poisson.c.

Here is the call graph for this function: