/** @file cdr.c

*   @brief Routines for handling the CDR (convection-diffusion-reaction) equations

*

* Default interpolation methods.\n

* May be changed elsewhere.\n

* \n

* interpol_method_t *cdr_interpol_inside = \&interpol_quadratic;\n

* interpol_method_t *cdr_interpol_bnd = \&interpol_bilin;\n

* \n

*/

#include <stdlib.h>

#include <stdio.h>

#include <math.h>

#include <string.h>

#include <sys/time.h>

#include "cdr.h"

#include "cstream.h"

#include "grid.h"

#include "interpol2.h"

#include "mapper.h"

#include "parameters.h"

#include "poisson.h"

#include "proto.h"

#include "react_table.h"

#include "rz_array.h"

#include "species.h"

#include "tree.h"

/** Declaration of the functions to interpolate the densities. */

decl_mapper_funcs(dens);

mapper_t **dens_mappers, **dens_bnd_mappers;

int cdr_ext_bnd_cond[4];

double z_cutoff = 10000000.0;

double max_err = 0.0;

react_table mu;

react_table diff;

extern pois_problem_t *pois_electrostatic;

static void set_all_bnd (cdr_grid_t *grid,

			 int sign, int ir, int iz, int itheta,

			 int dim0, int inout,

			 int dim1, int dim1_from, int dim1_to,

			 int dim2, int dim2_from, int dim2_to);

static void set_axis_bnd (cdr_grid_t *grid);

static void max_update (double *x, double cmp);

static double f_ad (rz_array_t *dens, double efield, double dc, double mc,

		    int ir, int iz, int itheta, int dim);

static double psi (double theta);

static void prepare_grid (cdr_grid_t *grid);

static double curv_at (cdr_grid_t *grid, rz_array_t *ar,

		       int ir, int iz, int itheta, double (*f) (double));

static int needs_refinement (cdr_grid_t *grid, int ir, int iz, int itheta,

			     int *in_edge);

static int any_needs_refinement (cdr_grid_t *grid, int ir, int iz,

				 int *in_edge);

static int brick_needs_refinement (cdr_grid_t *grid, int r0, int z0,

				   int r1, int z1, int *in_edge);

static void refine_in (cdr_grid_t *grid, int cr0, int cz0, int cr1, int cz1,

		       int contains_edge);

static void restrict_from (cdr_grid_t *parent, cdr_grid_t *child);

static int match_grids (cdr_grid_t *fro, cdr_grid_t *to);

double gauss2_xyz (double x, double y, double z);

int curr_seed = 0;

react_table mu;

react_table diff;

void aux_dump_frames_r (cdr_grid_t *grid, FILE *fp);

/** @brief Some initialization is needed for this module which has to be

 * inside this function.

 */

void

cdr_init (void)

{

  cdr_ext_bnd_cond[BND_BOTTOM] = cdr_bnd_bottom;

  cdr_ext_bnd_cond[BND_TOP] = cdr_bnd_top;

  cdr_ext_bnd_cond[BND_RIGHT] = cdr_bnd_right;

  dens_mappers = cdr_mappers_a (interpol_methods_index[cdr_interp_in]);

  dens_bnd_mappers = cdr_mappers_a (interpol_methods_index[cdr_interp_bnd]);

  react_table_read("input/default.mobility",&mu);

  react_table_read("input/default.diffusion",&diff);

}

/** @brief And the cleaning up has to be here. */

void

cdr_end (void)

{

  cdr_free_mappers (dens_mappers);

  cdr_free_mappers (dens_bnd_mappers);

}

/** @brief Creates a new 3D cdr grid */

cdr_grid_t*

cdr_new_3d_a (int r0, int z0, int r1, int z1, int ntheta)

{

  rz_array_t *ez, *er, *etheta, *eabs, *charge, **dens, **d_dens, *photo;

  cdr_grid_t *grid;

  int i;

  REAL *max_dens;

  debug (2, "cdr_new_3d_a (%d, %d, %d, %d)\n", r0, z0, r1, z1);

  ez = rz_new_3d_a (r0, z0, r1, z1, ntheta);

  er = rz_new_3d_a (r0, z0, r1, z1, ntheta);

  etheta = rz_new_3d_a (r0, z0, r1, z1, ntheta);

  eabs = rz_new_3d_a (r0, z0, r1, z1, ntheta);

  charge = rz_new_3d_a (r0, z0, r1, z1, ntheta);

  photo = has_photoionization? rz_new_3d_a (r0, z0, r1, z1, ntheta): NULL;

  /* Apart from the real species, (electron, ions...) we can create

   * a virtual species that enter into the possible reactions.

   * For example, the absolute value of the electric field, eabs.

   */

  dens = (rz_array_t **) xmalloc (sizeof(rz_array_t*)

				  * (no_species + N_VIRTUAL_SPECIES));

  d_dens = (rz_array_t **) xmalloc (sizeof(rz_array_t*) * no_species);

  max_dens = (REAL *) xmalloc (sizeof(REAL) * no_species);

  for (i = 0; i < no_species; i++) {

    dens[i] = rz_new_3d_a (r0, z0, r1, z1, ntheta);

    d_dens[i] = rz_new_3d_a (r0, z0, r1, z1, ntheta);

  }

  dens[no_species] = eabs;

  grid = (cdr_grid_t *) xmalloc (sizeof(cdr_grid_t));

  grid->r0 = r0;

  grid->r1 = r1;

  grid->z0 = z0;

  grid->z1 = z1;

  grid->charge = charge;

  grid->er = er;

  grid->ez = ez;

  grid->etheta = etheta;

  grid->eabs = eabs;

  grid->photo = photo;

  grid->dens = dens;

  grid->d_dens = d_dens;

  grid->max_dens = max_dens;

  grid->ntheta = ntheta;

  grid->contains_edge = FALSE;

  init_leaf (grid);

  return grid;

}

/** @brief Creates a guest cdr grid:

 *

 *  one that does not have memory of its own but

 *  points to the data of another grid (the @a host grid).  This is useful

 *  to easily set boundary conditions.  Probably, this function and

 *  cdr_new_3d_a should be somehow merged.

 */

cdr_grid_t*

cdr_guest (cdr_grid_t *host, int r0, int z0, int r1, int z1)

{

  rz_array_t *ez, *er, *etheta, *eabs, *charge, **dens, **d_dens, *photo;

  cdr_grid_t *grid;

  int i;

  REAL *max_dens;

  debug (2, "cdr_guest (" grid_printf_str ",%d, %d, %d, %d)\n",

	 grid_printf_args(host), r0, z0, r1, z1);

  ez = rz_guest (host->ez, r0, z0, r1, z1);

  er = rz_guest (host->er, r0, z0, r1, z1);

  etheta = rz_guest (host->etheta, r0, z0, r1, z1);

  eabs = rz_guest (host->eabs, r0, z0, r1, z1);

  charge = rz_guest (host->charge, r0, z0, r1, z1);

  photo = (has_photoionization?

	   rz_guest (host->photo, r0, z0, r1, z1): NULL);

  dens = (rz_array_t **) xmalloc (sizeof(rz_array_t*)

				  * (no_species + N_VIRTUAL_SPECIES));

  d_dens = (rz_array_t **) xmalloc (sizeof(rz_array_t*) * no_species);

  max_dens = (REAL *) xmalloc (sizeof(REAL) * no_species);

  for (i = 0; i < no_species; i++) {

    dens[i] = rz_guest (host->dens[i], r0, z0, r1, z1);

    d_dens[i] = rz_guest (host->d_dens[i], r0, z0, r1, z1);

  }

  dens[no_species] = eabs;

  grid = (cdr_grid_t *) xmalloc (sizeof(cdr_grid_t));

  grid->r0 = r0;

  grid->r1 = r1;

  grid->z0 = z0;

  grid->z1 = z1;

  grid->charge = charge;

  grid->er = er;

  grid->ez = ez;

  grid->etheta = etheta;

  grid->eabs = eabs;

  grid->photo = photo;

  grid->dens = dens;

  grid->d_dens = d_dens;

  grid->max_dens = max_dens;

  grid->ntheta = host->ntheta;

  grid->contains_edge = FALSE;

  set_leaf (grid, host->parent, NULL, NULL, host->level);

  return grid;

}

/** @brief Returns a grid with the same dimensions as grid and at the same

 * level.

 *

 * Realize that the data is _not_ copied into the new grid.

 */

cdr_grid_t*

cdr_like_a (cdr_grid_t *grid)

{

  cdr_grid_t *n;

  debug (2, "cdr_like_a(" grid_printf_str ")\n",

	 grid_printf_args(grid));

  n = cdr_new_3d_a (grid->r0, grid->z0, grid->r1, grid->z1, grid->ntheta);

  n->level = grid->level;

  n->ext_bound = grid->ext_bound;

  n->contains_edge = grid->contains_edge;

  return n;

}

/** @brief Returns a grid exactly like te one received, with the densities

 * copied.

 *

 * But the descendants are ignored.

 */

cdr_grid_t*

cdr_clone_a (cdr_grid_t *grid)

{

  cdr_grid_t *n;

  int itheta, s, nr, nz, nspecs;

  n = cdr_like_a (grid);

  nr = grid->r1 - grid->r0 + 4;

  nz = grid->z1 - grid->z0 + 4;

  nspecs = no_species;

  /* When we use the efield as a refinement criterium, we also have to copy

   * eabs (which is a @i virtual density).

   */

  if (ref_level_eabs >= 0 && ref_threshold_eabs >= 0.0) nspecs++;

#pragma omp parallel private(s)

  {

#pragma omp for

    iter_grid_theta_n (grid, itheta, 2) {

      for (s = 0; s < nspecs; s++) {

	rz_copy_modes (grid->dens[s], grid->r0 - 2, grid->z0 - 2,

		       n->dens[s], n->r0 - 2, n->z0 - 2,

		       nr, nz, itheta, itheta);

      }

    }

  }

  return n;

}

/** @brief Impose periodic boundary conditions in \f$\theta\f$ on a grid. */

void

cdr_set_periodic (cdr_grid_t *grid)

{

  int i;

  rz_set_periodic (grid->etheta);

  for (i = 0; i < no_species; i++) {

    /* For negative r, we have the values

     * of the densities on an opposite itheta.

     */

    rz_set_periodic (grid->dens[i]);

    rz_set_periodic (grid->d_dens[i]);

  }

  /* I am not completely sure that this is actually neccesary */

  rz_set_periodic (grid->er);

  rz_set_periodic (grid->ez);

}

/** @brief Recursive version of cdr_set_periodic. */

mk_recursive (cdr_set_periodic, cdr_grid_t)

/** @brief Frees a cdr grid. */

void

cdr_free (cdr_grid_t *grid)

{

  int i;

  debug (2, "cdr_free\n");

  rz_free (grid->charge);

  rz_free (grid->er);

  rz_free (grid->ez);

  rz_free (grid->etheta);

  rz_free (grid->eabs);

  if (has_photoionization) rz_free (grid->photo);

  for (i = 0; i < no_species; i++) {

    rz_free (grid->dens[i]);

    rz_free (grid->d_dens[i]);

  }

  free (grid->dens);

  free (grid->d_dens);

  free (grid->max_dens);

  free (grid);

}

/** @brief Recursive version of cdr_free. */

void

cdr_free_r (cdr_grid_t *grid)

{

  cdr_grid_t *leaf;

  debug (2, "cdr_free_r (" grid_printf_str ")\n", grid_printf_args (grid));

  free_childs (grid, leaf, cdr_free_r);

  cdr_free (grid);

}

/** @brief Calculates the charge in a CDR grid.

 *

 *  Important Note:  The charge is here divided by ntheta because

 *  our FFT functions are unnormalized and hence when they are applied

 *  forth and back they return the initial value multiplied

 *  by ntheta.  I believe that this is the best point to normalize the

 *  data.

 */

void

cdr_calc_charge (cdr_grid_t *grid)

{

  int s, r, z, t;

  double s_charge;

  debug (3, "cdr_calc_charge\n");

  rz_set_zero (grid->charge);

  for (s = 0; s < no_species; s++) {

    s_charge = spec_index[s]->charge;

#pragma omp parallel private(r, z)

    {

#pragma omp for

      iter_grid_theta (grid, t) {

	for (r = grid->r0 - 2; r < grid->r1 + 2; r++) {

	  for (z = grid->z0 - 2; z < grid->z1 + 2; z++) {

	    RZT (grid->charge, r, z, t) += (RZT (grid->dens[s], r, z, t)

					    * s_charge) / grid->ntheta;

	  }

	}

      }

    }

  }

}

/** @brief Recursive version of cdr_calc_charge. */

mk_recursive (cdr_calc_charge, cdr_grid_t)

/** @brief Fourier-transforms the charge in grid and all its descendants. */

void

cdr_dft_charge_r (cdr_grid_t *grid, int sign)

{

  cdr_grid_t *leaf;

  debug (2, "cdr_dft_charge_r(" grid_printf_str ")\n",

	 grid_printf_args(grid));

  iter_childs (grid, leaf) {

    cdr_dft_charge_r (leaf, sign);

  }

  dft_transform (grid->charge, grid->charge, sign);

}

/** @brief Fourier-transforms the electric field in grid and all its descendants. */

void

cdr_dft_field_r (cdr_grid_t *grid, int sign)

{

  cdr_grid_t *leaf;

  debug (2, "cdr_dft_field_r(" grid_printf_str ")\n",

	 grid_printf_args(grid));

  iter_childs (grid, leaf) {

    cdr_dft_field_r (leaf, sign);

  }

  dft_transform (grid->er, grid->er, sign);

  dft_transform (grid->ez, grid->ez, sign);

  dft_diff ((grid_t *) grid, grid->etheta);

  dft_transform (grid->etheta, grid->etheta, sign);

}

/** @brief Creates a grid one level coarser than GRID that covers the same area

 *    or somewhat more.

 *

 *  Note:  The new grid contains information only about the charge,

 *  not the other functions.

 *  That's why we don't use cdr_restrict.

 */

cdr_grid_t*

cdr_create_coarser_a (cdr_grid_t *grid)

{

  int r0_half, r1_half, z0_half, z1_half, itheta;

  int r, z;

  cdr_grid_t *new_grid;

  debug (2, "cdr_create_coarser_a\n");

  r0_half = grid->r0 >> 1;

  r1_half = (grid->r1 + 1) >> 1;

  z0_half = grid->z0 >> 1;

  z1_half = (grid->z1 + 1) >> 1;

  new_grid = cdr_new_3d_a (r0_half, z0_half, r1_half, z1_half, grid->ntheta);

#pragma omp parallel private (r, z)

  {

    /** Note that here we may be going one cell out of the aparent boundaries

       of the input grids but we allocated enough memory for that. */

#pragma omp for

    iter_grid_3d_n (new_grid, r, z, itheta, 1) {

      /* Start cylindrical */

      /* Here we are interpolating the masses, but in many other places

       * we interpolate the densities. Why?

       */

      RZT (new_grid->charge, r, z, itheta) =

	0.25

	* (cyl_q (r + 0.25) * (RZT (grid->charge, 2 * r, 2 * z, itheta)

			      + RZT (grid->charge, 2 * r, 2 * z + 1, itheta))

	   + cyl_q (r + 0.75) * (RZT (grid->charge, 2 * r + 1, 2 * z, itheta)

				+ RZT (grid->charge, 2 * r + 1, 2 * z + 1,

				       itheta)))

	/ cyl_q (r + 0.5);

      /* End cylindrical */

    }

  }

  return new_grid;

}

/** @brief Takes a tree of input nodes and adds n coarser grids, returning a

 *  pointer to the root leaf.

 */

cdr_grid_t*

cdr_add_coarser_grids_a (cdr_grid_t *prev_root, int n)

{

  int i;

  cdr_grid_t *grid = NULL;

  debug (2, "cdr_add_coarser_grids_a\n");

  assert (n > 0);

  for (i = 0; i < n; i++) {

    grid = cdr_create_coarser_a (prev_root);

    set_leaf (grid, NULL, NULL, prev_root, prev_root->level - 1);

    prev_root->parent = grid;

    prev_root = grid;

  }

  return grid;

}

/** @brief Frees the memory allocated by cdr_add_coarser_grids_a.

 *

 * Note: be careful to pass here only a tree returned by

 * cdr_add_coarser_grids_a, since some grids may lose their references and

 * stay allocated forever otherwise.

 * (see assert below).

 */

void

cdr_free_coarser_grids (cdr_grid_t *prev_root, int n)

{

  int i;

  cdr_grid_t *grid;

  debug (2, "cdr_free_coarser_grids_a\n");

  for (i = 0; i < n; i++) {

    assert (NULL == prev_root->next);

    grid = prev_root->first_child;

    cdr_free (prev_root);

    prev_root = grid;

  }

  prev_root->parent = NULL;

}

/** @brief The interface between the cdr and the poisson parts

 *

 *  If return_pois is true, returns the Poisson trees and they have

 *  to be de-allocated by the calling program.  If not, takes care himself

 *  of deallocation and returns NULL.

 */

pois_grid_t**

cdr_calc_field_r (cdr_grid_t *cdr, int return_pois)

{

  pois_grid_t **pois_modes;

  debug (2, "cdr_calc_field_r (" grid_printf_str ", %d)\n",

	 grid_printf_args(cdr), return_pois);

  if (cdr->ntheta != 1)

    cdr_dft_charge_r (cdr, 1);

  pois_modes = pois_solve_a (cdr, pois_electrostatic);

  if (cdr->ntheta != 1)

    cdr_dft_field_r (cdr, -1);

  cdr_add_ext_field_r (cdr);

  /* Formerly, I added the inhomogeneous field here.  But this is problematic:

   *  see pois_add_inhom_phi_r in poisson.c for details.

   *

   *  To return back to adding electric fields and not potentials, uncomment

   *  these lines:

   *  if (pois_inhom) {

   *     q_factor = pois_inhom_q_factor (pois_modes[0]);

   *     debug (1, "q_factor = %g\n", q_factor);

   *

   *     cdr_add_inhom_field_r (cdr, q_factor);

   *  }

   */

  /* Calculates the absolute value of the field. */

  cdr_calc_eabs_r (cdr);

  /* Impose periodic boundary conditions in theta. */

  if (cdr->ntheta != 1)

    cdr_set_periodic_r (cdr);

  if (!return_pois) {

    int i;

    for (i = 0; i < max_ntheta; i++) {

      pois_free_r (pois_modes[i]);

    }

    free (pois_modes);

    pois_modes = NULL;

  }

  return pois_modes;

}

/** @brief Adds an external electric field to a cdr @a grid.

 *

 * Really, it would be more efficient to include that field in some other

 * routine and we will avoid one more grid sweep. However, for the sake

 * of code clarity, extendability and maintenability I prefer to add it here.

 *

 * ext_e_r, ext_e_z and ext_e_theta are functions of three real numbers

 * (r, z, \f$\theta\f$) that return each of the components of the electric field.

 *

 *  Note that it is the responsibility of e_theta to behave correctly

 *  (i.e. periodically) when it receives thetas outside \f$[0, 2 \pi]\f$.

 */

void

cdr_add_ext_field (cdr_grid_t *grid)

{

  int ir, iz, itheta;

  debug (2, "cdr_add_ext_field (" grid_printf_str ")\n",

	 grid_printf_args(grid));

#pragma omp parallel private(ir, iz)

  {

    #pragma omp for

    iter_grid_3d_n (grid, ir, iz, itheta, 2) {

      RZT(grid->er, ir, iz, itheta) +=

	ext_e_r (er_r_at (ir, grid->level),

		 er_z_at (iz, grid->level),

		 theta_at (itheta));

      RZT(grid->ez, ir, iz, itheta) +=

	ext_e_z (ez_r_at (ir, grid->level),

		 ez_z_at (iz, grid->level),

		 theta_at (itheta));

      RZT(grid->etheta, ir, iz, itheta) +=

	ext_e_theta (r_at (ir, grid->level),

		     z_at (iz, grid->level),

		     etheta_theta_at (itheta));

    }

  }

}

/** @brief Recursive version of cdr_add_ext_field. */

mk_recursive (cdr_add_ext_field, cdr_grid_t)

/** @brief Adds to a cdr grid the inhomogeneous field created by a charge

 * q located at (r = 0, z = pois_inhom_z).

 *

 * THIS FUNCTION IS KEPT FOR REFERENCE ONLY.  It was superseded by

 * pois_add_inhom_phi_r (see poisson.c)

 */

void

cdr_add_inhom_field_r (cdr_grid_t *grid, double q)

{

  int ir, iz, itheta;

  cdr_grid_t *child;

  debug (2, "cdr_add_inhom_field_r (" grid_printf_str ", q = %f)\n",

	 grid_printf_args(grid), q);

  assert (0 != pois_inhom);

#pragma omp parallel private(ir, iz)

  {

    #pragma omp for

    iter_grid_3d_n (grid, ir, iz, itheta, 2) {

      RZT(grid->er, ir, iz, itheta) +=

	q * pois_inhom_er (er_r_at (ir, grid->level),

			   er_z_at (iz, grid->level));

      RZT(grid->ez, ir, iz, itheta) +=

	q * pois_inhom_ez (ez_r_at (ir, grid->level),

			   ez_z_at (iz, grid->level));

    }

  }

  iter_childs (grid, child) {

    cdr_add_inhom_field_r (child, q);

  }

}

/** @brief Calculates the absolute value of the electric field by linear

 *  interpolation of its components.

 *

 *  Of course, it requires that all the components are already present. */

void

cdr_calc_eabs (cdr_grid_t *grid)

{

  int ir, iz, itheta;

  double er, ez, etheta;

#pragma omp parallel

  {

#pragma omp for private(ir, iz, er, ez, etheta)

    iter_grid_3d (grid, ir, iz, itheta) {

      er = 0.5 * (RZT (grid->er, ir, iz, itheta)

		  + RZT (grid->er, ir - 1, iz, itheta));

      ez = 0.5 * (RZT (grid->ez, ir, iz, itheta)

		  + RZT (grid->ez, ir, iz - 1, itheta));

      if (grid->ntheta != 1) {

	etheta = 0.5 * (RZT (grid->etheta, ir, iz, itheta)

			+ RZT (grid->etheta, ir, iz, itheta - 1));

      } else {

	etheta = 0;

      }

      RZT (grid->eabs, ir, iz, itheta) = sqrt (er * er + ez * ez

					     + etheta * etheta);

    }

  }

}

/** @brief Recursive version of cdr_calc_eabs. */

mk_recursive (cdr_calc_eabs, cdr_grid_t)

/********************

 * Time integration *

*********************/

/** @brief Makes sure that all species of one grid are non-negative. */

void

cdr_nonegative (cdr_grid_t *grid)

{

  int s, ir, iz, itheta;

#pragma omp parallel

  {

#pragma omp for private(ir, iz, s)

    iter_grid_theta_n (grid, itheta, 2) {

      iter_grid_n(grid, ir, iz, 2) {

	   for (s = 0; s < no_species; s++)

	        *RZTP (grid->dens[s], ir, iz, itheta) =

	        MYMAX(0, *RZTP (grid->dens[s], ir, iz, itheta));

         }

       }

     }

}

/** @brief Recursive version of cdr_nonegative. */

mk_recursive (cdr_nonegative, cdr_grid_t)

/** @brief Sets the boundary values for the root grid and those grids that

   have a boundary coinciding with any of the boundaries of the root grid. */

void

cdr_set_ext_bnd (cdr_grid_t *grid)

{

  int r0, z0, r1, z1, ntheta;

  debug (2, "cdr_set_ext_bnd (" grid_printf_str " [grid->ext_bound = 0x%x])\n",

	 grid_printf_args(grid), grid->ext_bound);

  r0 = grid->r0;

  z0 = grid->z0;

  r1 = grid->r1;

  z1 = grid->z1;

  ntheta = grid->ntheta;

  /* Matching conditions (reduced to Hom. Neumann for ntheta = 1 at

     r = r0 */

  if (grid->ext_bound & BND_MASK (BND_LEFT)) {

    set_axis_bnd (grid);

  }

  /* At r = r1. */

  if (grid->ext_bound & BND_MASK (BND_RIGHT))

    set_all_bnd (grid, cdr_ext_bnd_cond[BND_RIGHT], r1 - 1, z0 - 2, 0,

		 R_INDX, 1,

		 Z_INDX, 0, z1 - z0 + 4,

		 THETA_INDX, 0, ntheta);

  /* At z = z0. */

  if (grid->ext_bound & BND_MASK (BND_BOTTOM))

    set_all_bnd (grid, cdr_ext_bnd_cond[BND_BOTTOM], r0 - 2, z0, 0,

		 Z_INDX, -1,

		 R_INDX, 0, r1 - r0 + 4,

		 THETA_INDX, 0, ntheta);

  /* At z = z1. */

  if (grid->ext_bound & BND_MASK (BND_TOP))

    set_all_bnd (grid, cdr_ext_bnd_cond[BND_TOP], r0 - 2, z1 - 1, 0,

		 Z_INDX, 1,

		 R_INDX, 0, r1 - r0 + 4,

		 THETA_INDX, 0, ntheta);

}

/** @brief Recursive version of cdr_set_ext_bnd. */

mk_recursive (cdr_set_ext_bnd, cdr_grid_t)

/** @brief Sets some boundary conditions for all the relevant variables in a

 * grid, to understand the parameters, see rz_array.c:rz_set_bnd(..)

 */

static void

set_all_bnd (cdr_grid_t *grid,

	     int sign, int ir, int iz, int itheta, int dim0, int inout,

	     int dim1, int dim1_from, int dim1_to,

	     int dim2, int dim2_from, int dim2_to)

{

  int s;

  for (s = 0; s < no_species; s++) {

    rz_set_bnd (grid->dens[s], sign, RZTP (grid->dens[s], ir, iz, itheta),

		dim0, inout,

		dim1, dim1_from, dim1_to,

		dim2, dim2_from, dim2_to);

  }

}

/** @brief Sets the proper boundary conditions for a grid that has its

 * left boundary at r = 0 (e.g. the root grid).

 *

 * We don't set Hom. Neumann at r = r0 since r = r0 is

 * not a real boundary: the proper condition here is matching the

 * opposite grid.  But note that this condition reduces to Hom.

 * Neumann in the case ntheta = 1 (because itheta2 = itheta = 0).

 */

static void

set_axis_bnd (cdr_grid_t *grid)

{

  int s, itheta, itheta2, ntheta;

  debug (2, "set_axis_bnd (" grid_printf_str ")\n",

	 grid_printf_args(grid));

  assert (0 == grid->r0);

  ntheta = grid->ntheta;

  for (s = 0; s < no_species; s++) {

    iter_grid_theta (grid, itheta) {

      itheta2 = (itheta + ntheta / 2) % ntheta;

      rz_copy_bnd (grid->dens[s], grid->dens[s],

		   /* sign  = */ 1,

		   /* start_from =  */ RZTP (grid->dens[s],

					     grid->r0, grid->z0 - 2, itheta2),

		   /* start_to = */    RZTP (grid->dens[s],

					     grid->r0, grid->z0 - 2, itheta),

		   R_INDX, BND_OUTWARD, BND_INWARD,

		   Z_INDX, 0, grid->z1 - grid->z0 + 4,

		   THETA_INDX, 0, 1 /* <- relative to start_xxx */);

    }

  }

}

/** @brief Calculates the derivatives of the densities for each point of the

 *  grid.

 */

void

cdr_calc_d_dens (cdr_grid_t *grid)

{

  react_apply_all (grid);

  cdr_advect_diffu_r (grid);

}

/* See "An Adaptive Grid Refinement Strategy for the Simulation

 * of Negative Streamers", C. Montijn, W. Hundsdorfer and U. Ebert

 * for the computations and notation used in this function.

 */

/** @brief This thing looks really ugly, but what does it mean?

 * For example, for stride = dens[electrons]->strides[R_INDX],

 *  it is the \f$p_{ij}\f$ of the aforementioned paper:

 *  \f$p_{i,j} = (\sigma_{i,j} - \sigma_{i-1,j}) / (\sigma_{i+1,j} - \sigma_{i,j})\f$

 */

#define pij(p_, stride_)					\

  (((*(p_)) - *((p_) - (stride_))) / (*((p_) + (stride_)) - *(p_)))

/** @brief This is the F^a + F^d in the paper (times \f$\delta\f${r, z, \f$\theta\f$}). */

static double

f_ad (rz_array_t *dens, double efield, double dc, double mc,

      int ir, int iz, int itheta, int dim)

{

  double sigma, psi_p, sigmasigma;

  int ishift;

  double *p;

  /* a p is the unmodified density-value in the cell (ir,iz,itheta).

   * Normally this should be n_e, as only electrons are mobile.

   */

  p = RZTP (dens, ir, iz, itheta);

  ishift = (efield > 0? 0: 1);