cdr.c

Go to the documentation of this file.

#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"

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

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);
}

void
cdr_end (void)
{
  cdr_free_mappers (dens_mappers);
  cdr_free_mappers (dens_bnd_mappers);
}

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;
}

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;
}

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;
}

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;
}

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);
}

mk_recursive (cdr_set_periodic, cdr_grid_t)


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);
}

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);
}

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;
     }
   }
      }
    }
  }
}

mk_recursive (cdr_calc_charge, cdr_grid_t)


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);
}

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);
}

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)
  {
#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;
}

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;
}

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;
}

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;
}

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));
    }
  }
}

mk_recursive (cdr_add_ext_field, cdr_grid_t)



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);
  }
}

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

    }
  }
}

mk_recursive (cdr_calc_eabs, cdr_grid_t)

/********************
 * Time integration *
*********************/


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));
         }
       }
     }
}

mk_recursive (cdr_nonegative, cdr_grid_t)


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);
}

mk_recursive (cdr_set_ext_bnd, cdr_grid_t)


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);
  }
}

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

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.
 */

#define pij(p_, stride_)               \
  (((*(p_)) - *((p_) - (stride_))) / (*((p_) + (stride_)) - *(p_)))

static double
f_ad (rz_array_t *dens, double efield, double dc, double mc,
      int ir, int iz, int itheta, int dim)
{
  double sigma, psi_p, sigmasigma;
  int ishift;
  double *p;

  /* a p is the unmodified density-value in the cell (ir,iz,itheta).
   * Normally this should be n_e, as only electrons are mobile.
   */
  p = RZTP (dens, ir, iz, itheta);

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

  /* sigma is just
   * sigma_{i    , j} for efield < 0,
   * sigma_{i + 1, j} for efield > 0
   */
  sigma = *(p + ishift * dens->strides[dim]);

  /* psi_p has to be
   *  psi(  p_{i  , j}) for efield < 0
   * -psi(1/p_{i+1, j}) for efield > 0
   */
  psi_p = pij (p + ishift * dens->strides[dim], dens->strides[dim]);

  /* E > 0            ---> psi_p = -psi( 1 / p_ij )
   * E < 0, p_ij == 0 ---> psi_p = 1
   * E < 0, p_ij != 0 ---> psi_p = psi( 1 / p_ij )
   */
  psi_p = (efield > 0)? -psi (psi_p) : psi_p == 0? 1.0: psi (1.0 / psi_p);

  /* sigmasigma is (sigma_{i+1, j} - sigma_{i, j}) */
  sigmasigma = (*p) - *(p + dens->strides[dim]);

  /* The first term corresponds to the advection term, the second to the
   * diffusion term.
   *
   * See eq. 3.1.5~3.1.8 in the thesis of C. Montijn. Division by cell-size
   * is done in the function calling this one.
   */
  return efield * mc * (sigma + psi_p * sigmasigma) + dc * sigmasigma;
}

void
cdr_advect_diffu (cdr_grid_t *grid)
{
  int ir, iz, itheta, s;
  int ii;
  double gdr_inv, gdz_inv, d_iso0, d_iso;

  double field;
  int dimd, t, sg;
  REAL *d_dens_p;
  double r_inv, er_r;
  rz_array_t *er,*ez,*etheta;

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

  gdr_inv = 1.0 / dr[grid->level];
  gdz_inv = 1.0 / dz[grid->level];
  er     = grid->er;
  ez     = grid->ez;
  etheta = grid->etheta;

  for (s = 0; s < no_species; s++) {
    /* The charge / mass ratio of a species. */
    double cm0, cm;
    int *strides;
    int ir, rmax, dim2;
    rz_array_t *dens, *d_dens, **d_dens_copy;

    if (spec_index[s]->mass <= 0) continue;

    dens = grid->dens[s];
    d_dens = grid->d_dens[s];
    strides = dens->strides;
    rmax = strides[1]; dim2 = strides[2]/strides[1];

    /* Make a copy of d_dens to be able to compare both versions
     * of arcos_cdr_advect
    int nr = grid->r1 - grid->r0 + 4;
    int nz = grid->z1 - grid->z0 + 4;
    d_dens_copy = (rz_array_t **) xmalloc (sizeof(rz_array_t*) * no_species);
    rz_copy(d_dens,grid->r0,grid->z0,*d_dens_copy,grid->r0,grid->z0,nr,nz);

    strides = dens->strides;
    */

    arcos_cdr_advect_vec(spec_index[s]->mass,spec_index[s]->charge,
                     dens->data,d_dens->data,er->data,er->len,
           ez->data,ez->len,diffusion_coeff,
                     dr[grid->level],dz[grid->level],
                     sprite_module,grid->r0,grid->r1,grid->z0,grid->z1);

    /* rz_free(*d_dens_copy); */
  }
}

mk_recursive (cdr_advect_diffu, cdr_grid_t)


static double 
psi (double theta) 
{
  if (theta < 0) return 0;
  if (theta < 0.4) return theta;
  if (theta < 4) return (1 + 0.5 * theta) / 3.0;
  else return 1.0;
}

double
cdr_courant (cdr_grid_t *grid)
{
  int s, ir, iz, itheta, depth, is_3d;

  double dr_min, dz_min, rdtheta_min, emax_r, emax_z, emax_theta, dmax, ddmax;
  double tau_a, tau_d, tau_rt, tau_f, max_diff;

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

  is_3d = (grid->ntheta == 1? 0: 1);

  /* Should we better keep track of the level of the finest grid?  */
  depth = grid_max_depth_r ((grid_t *) grid);
  dr_min = dr[depth];
  dz_min = dz[depth];
  rdtheta_min = dtheta * grid_rmin_r ((grid_t *) grid);

  /* This part is not parallel, but if needed it can be parallelized
   * easily by creating vectors xmax[0...ntheta-1], where we would store
   * the maximum for each ntheta.
   */
  emax_r = -1;
  emax_z = -1;
  emax_theta = -1;
  dmax = -1;
  ddmax = 0;

  iter_grid_3d (grid, ir, iz, itheta) {
    double mu;
    if (sprite_module) {
       mu = 1.0 / spr_density_at (z_at (iz, grid->level));
    } else {
       mu = 1.0;
    }
    //mu = sprite_module? 1.0 / spr_density_at (z_at (iz, grid->level)): 1.0;

    max_update (&emax_r, mu * fabs (RZT (grid->er, ir, iz, itheta)));
    max_update (&emax_z, mu * fabs (RZT (grid->ez, ir, iz, itheta)));
    max_update (&emax_theta, mu * fabs (RZT(grid->etheta, ir, iz, itheta)));

    for (s = 0; s < no_species; s++) {
      if (spec_index[s]->mass > 0) {
        max_update (&dmax, mu * RZT(grid->dens[s], ir, iz, itheta));
        max_update (&ddmax, -(RZT(grid->d_dens[s], ir, iz, itheta) / 
              RZT(grid->dens[s], ir, iz, itheta)));
      }
    }
  }

  /* This assert fails when the grid has become completely nan.  In that
   * case, let us not waste cpu cycles any longer.
   */
  assert (emax_r >= 0 && emax_z >= 0 && emax_theta >= 0 && dmax >= 0);

  /* Advection.  This one usually dominates. */
  tau_a = nu_a / (emax_r / dr_min + emax_z / dz_min 
        + is_3d * emax_theta / rdtheta_min);

  /* Diffusion. */
  /* <ISOTROPIC DIFFUSION> */
  if (sprite_module) {
    double dens_z1, dens_z0;
    /* Since we do not know where the density is higher, we check at both
     * extremes of the domain.  We are assuming that the density is monotonic
     * but nothing else (if it is not an exp. it works still).
     */
    dens_z0 = spr_density_at (z_at (grid->z0, grid->level));
    dens_z1 = spr_density_at (z_at (grid->z1, grid->level));
    max_diff = diffusion_coeff / MYMIN(dens_z0, dens_z1);
  } else {
    max_diff = 0.1; /* diffusion_coeff */
  }

  tau_d = nu_d / (max_diff / dr_min / dr_min 
        + max_diff / dz_min / dz_min
        + is_3d * max_diff / rdtheta_min / rdtheta_min);

  /* </ISOTROPIC DIFFUSION> */

  /* Mainly reaction growth / decay time */
  /* tau_f = nu_f / ddmax; */

  /* If sigma is zero somewhere, tau_f turns inf even if nu_f is very large.
   * We correct this here. */
  /*if (tau_f == 0)*/ tau_f = tau_d;  /* Just ignore tau_f please. */

  /* Dielectric relaxation time */
  tau_rt = nu_rt / dmax;

  return MYMIN (tau_a, MYMIN(tau_d, MYMIN(tau_f, tau_rt)));
}

static
void max_update (double *x, double cmp)
{
  if (*x < cmp) *x = cmp;
}

void
cdr_update (cdr_grid_t *orig, cdr_grid_t *dest, double h)
{
  int ir, iz, itheta, s;

  debug (2, "cdr_update (orig = " grid_printf_str ", dest = " 
    grid_printf_str ", %f)\n",
    grid_printf_args(orig), grid_printf_args(dest), h);

  /* We assume that orig and dest have the same shape. */
  assert (orig->r0 == dest->r0 && orig->z0 == dest->z0 && 
     orig->r1 == dest->r1 && orig->z1 == dest->z1 &&
     orig->ntheta == dest->ntheta);

#pragma omp parallel
  {
#pragma omp for private(ir, iz, s)
    iter_grid_3d_n (orig, ir, iz, itheta, 2) {
      for (s = 0; s < no_species; s++) {
   RZT (dest->dens[s], ir, iz, itheta) = 
     RZT (orig->dens[s], ir, iz, itheta) 
     + h * RZT (orig->d_dens[s], ir, iz, itheta);
      }
    }
  }
  cdr_calc_charge (dest);
}

void
cdr_rk2_update (cdr_grid_t *dens_0, cdr_grid_t *d_dens_1, 
      cdr_grid_t *d_dens_2, cdr_grid_t *dest, 
      double h)
{
  int ir, iz, itheta, s;

  debug (2, "cdr_rk2_update (dens_0 = " grid_printf_str 
    ",\n\t\td_dens_1 = " grid_printf_str
    ",\n\t\td_dens_2 = " grid_printf_str
    ",\n\t\tdest = " grid_printf_str ")\n",
    grid_printf_args(dens_0), 
    grid_printf_args(d_dens_1), 
    grid_printf_args(d_dens_2), 
    grid_printf_args(dest));

#pragma omp parallel
  {
#pragma omp for private(ir, iz, s)
    iter_grid_3d_n (dest, ir, iz, itheta, 2) {
      for (s = 0; s < no_species; s++) {
   RZT (dest->dens[s], ir, iz, itheta) = 
     RZT (dens_0->dens[s], ir, iz, itheta) 
     + 0.5 * h * RZT (d_dens_1->d_dens[s], ir, iz, itheta)
     + 0.5 * h * RZT (d_dens_2->d_dens[s], ir, iz, itheta);
      }
    }
  }
  cdr_calc_charge (dest);
}

void
cdr_rk2_update_r (cdr_grid_t *dens_0, cdr_grid_t *d_dens_1, 
        cdr_grid_t *d_dens_2, cdr_grid_t *dest, 
        double h)
{
  cdr_grid_t *c_dens_0, *c_d_dens_1, *c_d_dens_2, *c_dest;

  cdr_rk2_update (dens_0, d_dens_1, d_dens_2, dest, h);

  c_dens_0 = dens_0->first_child;
  c_d_dens_1 = d_dens_1->first_child;
  c_d_dens_2 = d_dens_2->first_child;

  iter_childs (dest, c_dest) {
    cdr_rk2_update_r (c_dens_0, c_d_dens_1, c_d_dens_2, c_dest, h);

    c_dens_0 = c_dens_0->next;
    c_d_dens_1 = c_d_dens_1->next;
    c_d_dens_2 = c_d_dens_2->next;
  }
}

void
cdr_self_update_r (cdr_grid_t *grid, double h)
{
  cdr_grid_t *child;
  cdr_update (grid, grid, h);

  iter_childs (grid, child) {
    cdr_self_update_r (child, h);
  }
}

void 
cdr_like_update_ar (cdr_grid_t *grid, cdr_grid_t *new_grid, double h)
{
  cdr_grid_t *child, *new_child, *prev = NULL;

  cdr_update (grid, new_grid, h);

  iter_childs (grid, child) {
    new_child = cdr_like_a (child);

    /* We can't use add_child here because we read the children of the grid
     * in inverse order as we inserted them.
     */
    if (NULL == prev) {
      new_grid->first_child = new_child;
    } else {
      prev->next = new_child;
    }
    set_leaf (new_child, new_grid, NULL, NULL, new_grid->level + 1);
    grid_inherit_ext_bound ((grid_t *) new_child);

    prev = new_child;

    cdr_like_update_ar (child, new_child, h);
  }
}

//double
//cdr_grid_error (cdr_grid_t *grid, cdr_grid_t *inter)
//{
//  int ir, iz, itheta;
//  int spec_nr;
//  double vol, err, val;
//
//  spec_nr = find_species_by_name("electrons");
//
//  assert(spec_nr > -1);
//
//  vol = dz[grid->level] * dr[grid->level];
//
//  iter_grid_3d(grid, ir, iz, itheta)
//    {
//
//      val = RZT(grid->dens[spec_nr], ir, iz, itheta);
//      err = val > 1e-12 ? (val - RZT(inter->dens[spec_nr], ir, iz, itheta)) / val : 0;
//     err = err > 0 ? err : -err;
//      if (max_err < err) 
// {
//   max_err = err;
// }
//    }
//    return err;
//}
//
//void
//cdr_grid_error_r (cdr_grid_t *grid, cdr_grid_t *inter)
//{
//  cdr_grid_t *g_child, *i_child;
//
//  cdr_grid_error(grid, inter);
//
//  i_child = inter->first_child;
//
//  for (g_child = grid->first_child; g_child; g_child = g_child->next)
//    {
//      cdr_grid_error_r(g_child, i_child);
//
//      i_child = i_child->next;
//    }
//}

double
cdr_rk2 (cdr_grid_t *grid, double h, double t)
{
  /* intermediate step. */
  cdr_grid_t *inter, *child;
  double courant_h;

  /* After this number of warnings about small timesteps, the code stops. */
  static int count_min_timestep = 5;

  debug (2, "cdr_rk2 (" grid_printf_str ", %f)\n",
    grid_printf_args(grid), h);

  cdr_nonegative_r (grid);

  /* We want to be sure that the boundaries are properly set. */
  prepare_grid (grid);

  /* Experimental: Moving z_max boundary -- DISABLED */

  cdr_calc_field_r (grid, /*return_pois = */ FALSE);

  /* If we are using the sprite module, we set here the proper magnitudes
     corresponding to the head altitude here. */
  if (sprite_module) spr_hook (grid);

  cdr_calc_d_dens (grid);

  /* The actual timestep is limited by the Courant stability criterium.
   *
   * Note that cdr_courant has to be called AFTER the computation of the fields.
   */
  courant_h = cdr_courant (grid);
  h = MYMIN (h, courant_h);

  if (h < warn_min_timestep) {
    warning ("Time step [h = %g] is rather small.  If this is ok, reduce "
        "warn_time_step to remove this warning\n", h);

    if (0 == count_min_timestep--) {
      fatal ("Too many warnings about small timesteps.\n");
    }
  }

  inter = cdr_like_a (grid);

  /* We create a new tree and put there y + dy/dt * h at each level. */
  cdr_like_update_ar (grid, inter, h);

  /* Set the proper boundary conditions also in the intermediate step. */
  prepare_grid (inter);

  cdr_calc_field_r (inter, /* return_pois = */ FALSE);
  cdr_calc_d_dens (inter);

  /* This is something like grid = grid + h / 2 * F(grid) + h / 2 F(inter) */
  cdr_rk2_update_r (grid, grid, inter, grid, h);

  /* Also, for the refinement functions, the boundary conditions have to be
   * appropriately set.
   */
  prepare_grid (grid);

  cdr_free_r (inter);

  return h;
}

static void
prepare_grid (cdr_grid_t *grid)
{
  /* Restrict the densities from finer to coarser grids. */
  cdr_restrict_r (grid);

  /* Sets the external boundaries. */
  cdr_set_ext_bnd_r (grid);

  /* Sets the boundaries of children grids interpolating from their parents. */
  cdr_set_bnd_r (grid);

  /* If two grids share a boundary, use the values of one to set the 
     boundaries of the other. */
  cdr_match_r (grid, NULL);

  /* Sets periodic boundary conditions in theta. */
  cdr_set_periodic_r (grid);
}  

/*****************************************************************
 * Refinement and mapping. 
 *****************************************************************/

void
cdr_update_refined (cdr_grid_t **ptree)
{
  cdr_grid_t *old, *new;

  old = *ptree;

  new = cdr_clone_a (old);
  cdr_calc_charge (new);
  cdr_refine_r (new, old);

  *ptree = new;

  cdr_free_r (old);
}

void
cdr_calc_maxs (cdr_grid_t *grid)
{
  int start = TRUE;
  int s, ir, iz, itheta;

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

  /* If we have a parent, we inherit the maxs from him. */
  if (grid->parent) {
    grid->max_charge = grid->parent->max_charge;
    for (s = 0; s < no_species; s++) {
      grid->max_dens[s] = grid->parent->max_dens[s];
    }
    return;  
  }

  iter_grid_3d (grid, ir, iz, itheta) {
    if (start) {
      start = FALSE;
      grid->max_charge = fabs(RZT (grid->charge, ir, iz, itheta));
      for (s = 0; s < no_species; s++) {
   grid->max_dens[s] = RZT (grid->dens[s], ir, iz, itheta);
      }
    } else {
      max_update (&grid->max_charge, fabs(RZT (grid->charge, ir, iz, itheta)));
      for (s = 0; s < no_species; s++) {
   max_update (&grid->max_dens[s], RZT (grid->dens[s], ir, iz, itheta));
      }
    }
  }

  debug (2, "->max_charge = %g, ->max_dens[0] = %g, ->max_dens[1] = %g\n",
    grid->max_charge, grid->max_dens[0], grid->max_dens[1]);
}

static double
curv_at (cdr_grid_t *grid, rz_array_t *ar, int ir, int iz, int itheta,
    double (*f) (double))
{
  int i, level;
  REAL x[3];
  double dur, duz;

  for (i = 0; i < 3; i++) {
    x[i] = RZT (ar, ir - 1 + i, iz, itheta);
    x[i] = f? f(x[i]): x[i];
  }

  level = grid->level;
  /* <CYLINDRICAL> */
  dur = (cyl_er_r_at (ir, level) * (x[2] - x[1]) 
    - cyl_er_r_at (ir - 1, level) * (x[1] - x[0])) / cyl_r_at (ir, level);
  /* </CYLINDRICAL> */

  for (i = 0; i < 3; i++) {
    x[i] = RZT (ar, ir, iz - 1 + i, itheta);
    x[i] = f? f(x[i]): x[i];
  }

  duz = (x[2] - 2 * x[1] + x[0]);

  return fabs(dur) + fabs(duz);
}

static int
needs_refinement (cdr_grid_t *grid, int ir, int iz, int itheta, int *in_edge)
{
  double curv;
  int s;
  /* For sprites, the eabs refinement criterium may depend on the location. */
  int s_ref_level_eabs;
  double s_ref_threshold_eabs;

  /* If the maximum charge is below this number, we ignore the refinement
   * criterium.
   */
  const double epsilon_charge = 1e-8;

  if (!sprite_module || ref_level_eabs < 0) {
    s_ref_level_eabs = ref_level_eabs;
    s_ref_threshold_eabs = ref_threshold_eabs;
  } else {
    double z, back_dens;
    z = z_at (iz, grid->level);
    back_dens = spr_density_at (z);

    s_ref_threshold_eabs = ref_threshold_eabs * back_dens;

    /* log2 appears because whenever we increase the background density by a
     * factor 2, we should refine up to one more level.
     */
    s_ref_level_eabs = ref_level_eabs + (int) floor (log2 (back_dens));
  }

  /* Electric field refinement. If eabs > ref_threshold_eabs, refine up to
   * some predefined level.
   */

  if (grid->level < s_ref_level_eabs && 
      RZT (grid->eabs, ir, iz, itheta) > s_ref_threshold_eabs) {
    debug (4, "Refine grid " grid_printf_str 
      " at ir = %d, iz = %d, itheta = %d\n"
      "\t because too high EABS [eabs = %g]\n",
      grid_printf_args(grid), ir, iz, itheta, 
      RZT (grid->eabs, ir, iz, itheta));

    *in_edge = FALSE;
    return TRUE;
  }

  curv = curv_at (grid, grid->charge, ir, iz, itheta, NULL);

  if (grid->max_charge > epsilon_charge 
      && curv / grid->max_charge > ref_threshold_charge) {
    debug (4, "Refine grid " grid_printf_str 
      " at ir = %d, iz = %d, itheta = %d\n"
      "\t because too high CHARGE curvature [curv / max_charge = %g]\n",
      grid_printf_args(grid), ir, iz, itheta, curv / grid->max_charge);

    *in_edge = FALSE;
    return TRUE;
  }

  for (s = 0; s < no_species; s++) {
    /* We ignore immobile species: */
    if (spec_index[s]->mass <= 0) continue;

    curv = curv_at (grid, grid->dens[s], ir, iz, itheta, NULL);
    if (curv / grid->max_dens[s] > ref_threshold_dens) {
      debug (4, "Refine grid " grid_printf_str 
        " at ir = %d, iz = %d, itheta = %d\n"
        "\t because too high DENS[%s] curvature [curv / max_dens = %g]\n",
        grid_printf_args(grid), ir, iz, itheta, spec_index[s]->name, 
        curv / grid->max_dens[s]);

      *in_edge = FALSE;
      return TRUE;
    }

    /* If the grid contains the leading edge, the criterium of
     * refinement also takes into account the density.
     */
    if (*in_edge && RZT (grid->dens[s], ir, iz, itheta) > ref_threshold_edge) {
      debug (4, "Refine grid " grid_printf_str 
        " at ir = %d, iz = %d, itheta = %d\n"
        "\t because too high DENS[%s] in leading edge [dens = %g]\n",
        grid_printf_args(grid), ir, iz, itheta, spec_index[s]->name, 
        RZT (grid->dens[s], ir, iz, itheta));
      return TRUE;
    }
  }    

  /* If we are close to the border, we also check the border itself. 
   *
   * NOTE: Does this have any sense?  We are anyway checking for
   * the boundaries, aren't we?
   */
  if (ir == grid->r0 + 1 && needs_refinement (grid, ir - 1, iz, itheta, 
                     in_edge)) 
    return TRUE;

  if (iz == grid->z0 + 1 && needs_refinement (grid, ir, iz - 1, itheta,
                     in_edge)) 
    return TRUE;

  if (ir == grid->r1 - 2 && needs_refinement (grid, ir + 1, iz, itheta,
                     in_edge)) 
    return TRUE;

  if (iz == grid->z1 - 2 && needs_refinement (grid, ir, iz + 1, itheta,
                     in_edge)) 
    return TRUE;

  return FALSE;
} 

static int
any_needs_refinement (cdr_grid_t *grid, int ir, int iz, int *in_edge)
{
  int itheta;
  int ret = FALSE;
  iter_grid_theta (grid, itheta) {
    if (needs_refinement (grid, ir, iz, itheta, in_edge)) {
      debug (6, grid_printf_str 
        " needs refinement at ir = %d, iz = %d itheta = %d\n",
        grid_printf_args(grid), ir, iz, itheta);
      ret = TRUE;
    }
    if (ret && ! *in_edge) return TRUE;
  }
  return ret;
}

static int
brick_needs_refinement (cdr_grid_t *grid, int r0, int z0, int r1, int z1,
         int *in_edge)
{
  int ir, iz;
  int ret = FALSE;

  for (ir = r0; ir < r1; ir ++) {
    for (iz = z1 - 1; iz >= z0; iz--) {
      ret = (any_needs_refinement (grid, ir, iz, in_edge) || ret);
      if (ret && ! *in_edge) return TRUE;
    }
  }
  return ret;
}

void
cdr_refine (cdr_grid_t *grid)
{
  int ir;

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

  assert (grid->level <= cdr_max_level);

  /* Checks if the maximum refinement level has been reached. */
  if (grid->level == cdr_max_level) return;

  cdr_calc_maxs (grid);

  for (ir = grid->r0; ir < grid->r1; ir += cdr_brick_dr) {
    int irmax = MYMIN (grid->r1, ir + cdr_brick_dr);
    int tainted, building, z0 = -1, z1 = -1;
    int iz;
    int in_edge;
    int first_edge;

    first_edge = in_edge = grid->contains_edge;

    /* The leading edge extends downwards. */
    for (iz = grid->z0, building = FALSE; 
    iz < grid->z1; iz += cdr_brick_dz) {
      int izmax = MYMIN (grid->z1, iz + cdr_brick_dz);

      tainted = brick_needs_refinement (grid, ir, iz, irmax, izmax, 
               &first_edge);
      if (tainted) {
   z1 = izmax;
   if (!building) {
     z0 = iz;
     building = TRUE;
   }
      } else /* ! tainted */ {
   if (building) {
     assert (z0 >= 0 && z1 > 0);
     /* If a grid does not satisfy the curvature criterium anywhere
      * and is here only because of the density threshold (edge) 
      * criterium, he does not deserve to be refined.
      */
     building = FALSE;
     if (!first_edge) {
       refine_in (grid, ir, z0, irmax, z1, in_edge);
       in_edge = FALSE;
     }
   }
      } 
    }
    if (building) {
      assert (z0 >= 0 && z1 > 0);
      if (!first_edge) refine_in (grid, ir, z0, irmax, z1, in_edge);
    }
  }
}

static void
refine_in (cdr_grid_t *grid, int cr0, int cz0, int cr1, int cz1, 
      int contains_edge)
{
  cdr_grid_t *child;
  int nr0, nr1, nz0, nz1;

  nr0 = cr0 << 1;
  nz0 = cz0 << 1;
  nr1 = cr1 << 1;
  nz1 = cz1 << 1;

  child = cdr_new_3d_a (nr0, nz0, nr1, nz1, grid->ntheta);

  add_child (grid, child);
  grid_inherit_ext_bound ((grid_t*) child);

  child->contains_edge = contains_edge;

  debug (3, "new child created {r0 = %d, z0 = %d, r1 = %d, z1 = %d, "
    "level = %d}\n", child->r0, child->z0, child->r1, child->z1, 
    child->level);
}

void
cdr_refine_r (cdr_grid_t *grid, cdr_grid_t *source)
{
  cdr_grid_t *child;
  int itheta;

  cdr_calc_charge (grid);

  cdr_refine (grid);

  iter_childs (grid, child) {
    if (NULL != source) {
#pragma omp parallel
      {    
#pragma omp for
   iter_grid_theta (child, itheta) {
     map_grid_r (dens_mappers, (grid_t *) source, (grid_t *) child, 
            itheta, 
            /* copy = */ TRUE, 
            /* interpol = */ TRUE, 
            /* coarsen = */ FALSE, 
            /* s_buff = */ 0, /* t_buff = */ 2);
   }
      }
    } else {
      cdr_init_dens (child);
    }
    cdr_refine_r (child, source);
  }
}

/* Matching of the boundaries:
 *
 * The boundary conditions of a grid that has a common border with another
 * one at the same level have to be matched. This is the purpose of this set
 * of functions.
 */

void
cdr_match_r (cdr_grid_t *grid1, cdr_grid_t *grid2)
{
  cdr_grid_t *child1, *child2;
  int match = TRUE;

  if (grid2) {
    debug (2, "cdr_match_r (grid1 = " grid_printf_str 
      ", grid2 = " grid_printf_str ")\n",
      grid_printf_args (grid1), grid_printf_args (grid2));
  } else {
    debug (2, "cdr_match_r (grid1 = " grid_printf_str 
      ", grid2 = NULL)\n", grid_printf_args (grid1));
  }

  if (NULL != grid2) {
    match = match_grids (grid1, grid2);
  }

  if (!match) return;

  iter_childs (grid1, child1) {
    debug (4, "\tchild1 = " grid_printf_str "\n",
      grid_printf_args (child1));
    iter_childs (grid1, child2) {
      debug (4, "\t\tchild2 = " grid_printf_str "\n",
        grid_printf_args (child2));

      cdr_match_r (child1, (child1 != child2)? child2: NULL);
    }

    if (NULL == grid2)
      continue;

    iter_childs (grid2, child2) {
      cdr_match_r (child1, child2);
    }
  }
}

static int
match_grids (cdr_grid_t *fro, cdr_grid_t *to)
{
  int s, x0, x1;

  debug (2, "match_grids (to = " grid_printf_str 
    ", fro = " grid_printf_str ")\n",
    grid_printf_args(to), grid_printf_args(fro));

  assert (to->level == fro->level);

  if (to->r1 == fro->r0) {
    x0 = MYMAX (to->z0 - 2, fro->z0 - 2);
    x1 = MYMIN (to->z1 + 2, fro->z1 + 2);
    if (x1 <= x0) return FALSE;

    for (s = 0; s < no_species; s++) {
      rz_copy_bnd (fro->dens[s], to->dens[s], 1,
         RZTP (fro->dens[s], fro->r0, x0, 0),
         RZTP (to->dens[s], to->r1 - 1, x0, 0),
         R_INDX, BND_OUTWARD, BND_OUTWARD,
         Z_INDX, 0, x1 - x0,
         THETA_INDX, 0, fro->ntheta);
    }
    return TRUE;
  }

  if (to->r0 == fro->r1) {
    x0 = MYMAX (to->z0 - 2, fro->z0 - 2);
    x1 = MYMIN (to->z1 + 2, fro->z1 + 2);
    if (x1 <= x0) return FALSE;

    for (s = 0; s < no_species; s++) {
      rz_copy_bnd (fro->dens[s], to->dens[s], 1,
         RZTP (fro->dens[s], fro->r1 - 1, x0, 0),
         RZTP (to->dens[s], to->r0, x0, 0),
         R_INDX, BND_INWARD, BND_INWARD,
         Z_INDX, 0, x1 - x0,
         THETA_INDX, 0, fro->ntheta);
    }
    return TRUE;
  }

  if (to->z1 == fro->z0) {
    x0 = MYMAX (to->r0 - 2, fro->r0 - 2);
    x1 = MYMIN (to->r1 + 2, fro->r1 + 2);
    if (x1 <= x0) return FALSE;

    for (s = 0; s < no_species; s++) {
      rz_copy_bnd (fro->dens[s], to->dens[s], 1,
         RZTP (fro->dens[s], x0, fro->z0, 0),
         RZTP (to->dens[s], x0, to->z1 - 1, 0),
         Z_INDX, BND_OUTWARD, BND_OUTWARD,
         R_INDX, 0, x1 - x0,
         THETA_INDX, 0, fro->ntheta);
    }
    return TRUE;
  }

  if (to->z0 == fro->z1) {
    x0 = MYMAX (to->r0 - 2, fro->r0 - 2);
    x1 = MYMIN (to->r1 + 2, fro->r1 + 2);
    if (x1 <= x0) return FALSE;

    for (s = 0; s < no_species; s++) {
      rz_copy_bnd (fro->dens[s], to->dens[s], 1,
         RZTP (fro->dens[s], x0, fro->z1 - 1, 0),
         RZTP (to->dens[s], x0, to->z0, 0),
         Z_INDX, BND_INWARD, BND_INWARD,
         R_INDX, 0, x1 - x0,
         THETA_INDX, 0, fro->ntheta);
    }
    return TRUE;
  }
  return FALSE;
}

void
cdr_set_bnd (cdr_grid_t *grid)
{
  cdr_grid_t *top, *bottom, *left, *right, *parent;
  int itheta;

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

  parent = grid->parent;
  if (NULL == parent) {
    assert (0 == grid->level);
    return;
  }

  /* Since our mapping functions are easier to implement when we map
   * between rectangular regions (grids) we map the boundary conditions 
   * by creating from grid four "guest" grids that share memory with him
   * and that contains each of the four boundaries.
   */
  top = cdr_guest (grid, grid->r0 - 2, grid->z1, 
         grid->r1 + 2, grid->z1 + 2);
  bottom = cdr_guest (grid, grid->r0 - 2, grid->z0 - 2, 
            grid->r1 + 2, grid->z0);
  left = cdr_guest (grid, grid->r0 - 2, grid->z0, 
          grid->r0, grid->z1);
  right = cdr_guest (grid, grid->r1, grid->z0, 
             grid->r1 + 2, grid->z1);

#pragma omp parallel
    {    
#pragma omp for
      iter_grid_theta (grid, itheta) {
   if (0 == (grid->ext_bound & BND_MASK (BND_TOP)))
     map_grid (dens_bnd_mappers, (grid_t *) parent, (grid_t *) top, itheta, 
          FALSE, TRUE, FALSE, 1, 0);
   
   if (0 == (grid->ext_bound & BND_MASK (BND_BOTTOM)))
     map_grid (dens_bnd_mappers, (grid_t *) parent, (grid_t *) bottom, itheta,
          FALSE, TRUE, FALSE, 1, 0);

   if (0 == (grid->ext_bound & BND_MASK (BND_LEFT))) {
     assert (grid->r0 != 0);
     map_grid (dens_bnd_mappers, (grid_t *) parent, (grid_t *) left, itheta, 
        FALSE, TRUE, FALSE, 1, 0);
     
   }

   if (0 == (grid->ext_bound & BND_MASK (BND_RIGHT)))
     map_grid (dens_bnd_mappers, (grid_t *) parent, (grid_t *) right, itheta, 
          FALSE, TRUE, FALSE, 1, 0);
      }
    }

    cdr_free (top);
    cdr_free (bottom);
    cdr_free (left);
    cdr_free (right);
}

mk_recursive (cdr_set_bnd, cdr_grid_t)


mapper_t**
cdr_mappers_a (interpol_method_t *interp_method)
{
  mapper_t *mapper, **mappers;
  int s, nmappers;

  nmappers = no_species;

  /* mapper->extra represent the species whose density we will map.
   * There are no_species species, but no_species + 1 represents
   * the abs value of the electric field, which is mapped only to use
   * as an electric field-based refinement criterium.
   */
  if (ref_level_eabs >= 0 && ref_threshold_eabs >= 0.0) nmappers++;

  mappers = (mapper_t **) xmalloc (sizeof(mapper_t*) * (nmappers + 1));

  for (s = 0; s < nmappers; s++) {
    mapper = (mapper_t *) xmalloc (sizeof(mapper_t));
    mapper->extra = s;
    mapper->interpol_method = interp_method;
    mapper->copy = dens_copy;
    mapper->interpol_set = dens_interpol_set;
    mapper->interpol = dens_interpol;
    mapper->coarsen = NULL;
    mapper->shift_r = 0;
    mapper->shift_z = 0;
    mappers[s] = mapper;
  }

  mappers[nmappers] = NULL;

  return mappers;
}

void
cdr_free_mappers (mapper_t **mappers)
{
  int s;

  for (s = 0; s < no_species; s++) {
    free (mappers[s]);
  }
  free (mappers);
}  

void
dens_copy (mapper_t *mapper, grid_t *source, grid_t *target, 
      int ir, int iz, int itheta)
{
  cdr_grid_t *cdr_s, *cdr_d;

  cdr_s = (cdr_grid_t*) source;
  cdr_d = (cdr_grid_t*) target;

  RZT (cdr_d->dens[mapper->extra], ir, iz, itheta) = 
    RZT (cdr_s->dens[mapper->extra], ir, iz, itheta);
}

int
dens_interpol_set (mapper_t *mapper, grid_t *source, interpol_t *interpol,
         int pr, int pz, int itheta)
{
  cdr_grid_t *cdr;

  cdr = (cdr_grid_t*) source;

  interpol_set_stencil_at (source, interpol, 
            r_at (pr, cdr->level),
            z_at (pz, cdr->level),
            cdr->dens[mapper->extra], pr, pz, itheta);
  return TRUE;
}

void
dens_interpol (mapper_t *mapper, grid_t *source, grid_t *target, 
          interpol_t *interpol, 
          int ir, int iz, int itheta)
{
  double r, z;
  cdr_grid_t *cdr;

  cdr = (cdr_grid_t *) target;

  r = r_at (ir, cdr->level);
  z = z_at (iz, cdr->level);

  RZT (cdr->dens[mapper->extra], ir, iz, itheta) = 
    interpol_apply (interpol, r, z);
}

static void
restrict_from (cdr_grid_t *parent, cdr_grid_t *child)
{
  int ir, iz, cr, cz, itheta, s;
  REAL sum;

  debug (2, "restrict_from (parent =" grid_printf_str 
    ", child = " grid_printf_str ")\n",
    grid_printf_args(parent), grid_printf_args(child));

  assert (parent->level == child->level - 1);

#pragma omp parallel private(sum, ir, iz, cr, cz, s)
  {
#pragma omp for
    iter_grid_theta (parent, itheta) {
      iter_grid_parent (child, ir, iz) {
   for (s = 0; s < no_species; s++) {
     sum = 0;
     for (cr = (ir << 1); cr <= (ir << 1) + 1; cr++)
       for (cz = (iz << 1); cz <= (iz << 1) + 1; cz++)
         sum += cyl_r_at (cr, child->level) 
      * RZT (child->dens[s], cr, cz, itheta);
     RZT (parent->dens[s], ir, iz, itheta) = 
       0.25 * sum / cyl_r_at (ir, parent->level);
   }
      }
    }
  }
}

void
cdr_restrict (cdr_grid_t *grid)
{
  cdr_grid_t *child;

  iter_childs(grid, child) {
    restrict_from (grid, child);
  }
}

mk_tail_recursive (cdr_restrict, cdr_grid_t)



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

  q =  invpi32 * 1 / (seed_index[curr_seed]->sigma_x * seed_index[curr_seed]->sigma_y *
            seed_index[curr_seed]->sigma_z)
    * exp (- SQ(x - seed_index[curr_seed]->x0) / SQ(seed_index[curr_seed]->sigma_x)
      - SQ(y - seed_index[curr_seed]->y0) / SQ(seed_index[curr_seed]->sigma_y)
      - SQ(z - seed_index[curr_seed]->z0) / SQ(seed_index[curr_seed]->sigma_z));
  return q;
}

void 
cdr_set_dens (cdr_grid_t *cdr, int species, int mode, double factor,
          double (*f) (double, double, double))
{
  int ir, iz, itheta;
  double ctheta, stheta, x, y;
  debug (2, "cdr_set_dens (" grid_printf_str ", %s, ...)\n",
    grid_printf_args(cdr), spec_index[species]->name);

#pragma omp parallel private(ir, iz, ctheta, stheta, x, y)
  {
#pragma omp for
    iter_grid_theta_n (cdr, itheta, 2) {
      //sincos (theta_at(itheta), &stheta, &ctheta);
      stheta=sin(theta_at(itheta));
      ctheta=cos(theta_at(itheta));
      iter_grid_r_n(cdr, ir, 2) {
   x = r_at(ir, cdr->level) * ctheta;
   y = r_at(ir, cdr->level) * stheta;
   iter_grid_z_n(cdr, iz, 2) {
     switch (mode) {
     case SET_DENS_OVERWRITE:
       *RZTP(cdr->dens[species], ir, iz, itheta) = 
         factor * f (x, y, z_at(iz, cdr->level));
       break;
     case SET_DENS_ADD:
       *RZTP(cdr->dens[species], ir, iz, itheta) += 
         factor * f (x, y, z_at(iz, cdr->level));
       break;
     case SET_DENS_SUB:
       *RZTP(cdr->dens[species], ir, iz, itheta) -= 
         factor * f (x, y, z_at(iz, cdr->level));
       break;
     }
   }
      }
    }
  }
}

double f_one(double x, double y, double z)
{
  return 1.0;
}

void
cdr_init_dens (cdr_grid_t *cdr)
{
  int cnt;

  debug (2, "cdr_init_dens (" grid_printf_str ")\n",
    grid_printf_args(cdr));

  for (cnt = 0; cnt < no_species; cnt++)
  {
    cdr_set_dens(cdr, cnt, SET_DENS_OVERWRITE, 0.0, f_one);
  }

  for (cnt = 0; cnt < no_seed; cnt++)
  {
    curr_seed = cnt;
    if (seed_index[cnt]->type == 0)
      cdr_set_dens(cdr, seed_index[cnt]->species, SET_DENS_ADD, seed_index[cnt]->value, gauss2_xyz);
    if (seed_index[cnt]->type == 1)
      cdr_set_dens(cdr, seed_index[cnt]->species, SET_DENS_ADD, seed_index[cnt]->value, f_one);
  }
}

cdr_grid_t*
cdr_scratch_init (void)
{
  cdr_grid_t *cdr;

  cdr = cdr_new_3d_a (0, 0, gridpoints_r, gridpoints_z, max_ntheta);

  cdr->level = 0;
  cdr->ext_bound = BND_MASK_ALL;
  cdr->contains_edge = TRUE;

  cdr_init_dens (cdr);
  cdr_refine_r (cdr, NULL);

  if (perturb_epsilon > 0.0 && max_ntheta > 1) {
    dft_dens_perturb_r (cdr, electrons, NULL);

    /* Maybe the perturbation has messed up the boundary conditions, so we
     * have to repair them.
     */
    cdr_set_ext_bnd_r (cdr);
    cdr_set_periodic_r (cdr);
  }
  return cdr;
}

void
cdr_dump (cdr_grid_t *grid, const char *prefix, const char *name)
{
  char *fname;
  int s, nt;

  int m = cdr_output_margin;

  /* To make it easier to produce plots in theta, we save one extra
   *  theta that has to give the same data as the first one.  Except if
   *  we are working in 2D, where anyway we have only one possible theta (0)
   */
  nt = grid->ntheta == 1? 1: grid->ntheta + 1;


#pragma omp parallel sections private(fname)
  {
#pragma omp section
    {
      asprintf (&fname, "%s/r.%s.tsv", prefix, name);
      rz_axis_dump (fname, grid->r0 - m, grid->r1 + m, dr[grid->level]);
      free (fname);
    }
#pragma omp section
    {
      asprintf (&fname, "%s/z.%s.tsv", prefix, name);
      rz_axis_dump (fname, grid->z0 - m , grid->z1 + m, dz[grid->level]);
      free (fname);
    }
#pragma omp section
    {
      /* In all this function we use ntheta + 1 to make matlab's life easier.
       * Note that some of these variables have not been made periodic
       * and hence the value for ntheta may be undefined. 
       */
      asprintf (&fname, "%s/theta.%s.tsv", prefix, name);
      rz_axis_dump (fname, 0, nt, dtheta);
      free (fname);
    }
#pragma omp section
    {
      asprintf (&fname, "%s/charge.%s.tsv", prefix, name);
      rz_dump_3d (grid->charge, fname, "w", grid->r0 - m, grid->z0 - m, 
        grid->r1 + m, grid->z1 + m, nt);
      free (fname);
    }
#pragma omp section
    {
      asprintf (&fname, "%s/er.%s.tsv", prefix, name);
      rz_dump_3d (grid->er, fname, "w", grid->r0 - m, grid->z0 - m, 
        grid->r1 + m, grid->z1 + m, nt);
      free (fname);
    }
#pragma omp section
    {
      asprintf (&fname, "%s/ez.%s.tsv", prefix, name);
      rz_dump_3d (grid->ez, fname, "w", grid->r0 - m, grid->z0 - m, 
        grid->r1 + m, grid->z1 + m, nt);
      free (fname);
    }
#pragma omp section
    {
      asprintf (&fname, "%s/etheta.%s.tsv", prefix, name);
      rz_dump_3d (grid->etheta, fname, "w", grid->r0 - m, grid->z0 - m, 
        grid->r1 + m, grid->z1 + m, nt);
      free (fname);
    }

#pragma omp section
    {
      asprintf (&fname, "%s/eabs.%s.tsv", prefix, name);
      rz_dump_3d (grid->eabs, fname, "w", grid->r0 - m, grid->z0 - m, 
        grid->r1 + m, grid->z1 + m, nt);
      free (fname);
    }

#pragma omp section
    {
      if (has_photoionization) {
   asprintf (&fname, "%s/photo.%s.tsv", prefix, name);
   rz_dump_3d (grid->photo, fname, "w", grid->r0 - m, grid->z0 - m, 
          grid->r1 + m, grid->z1 + m, nt);
   free (fname);
      }
    }

#pragma omp section
    {
      for (s = 0; s < no_species; s++) {
   /* The densities */
   asprintf (&fname, "%s/%s.%s.tsv", prefix, spec_index[s]->name, name);
   rz_dump_3d (grid->dens[s], fname, "w", grid->r0 - m, grid->z0 - m, 
          grid->r1 + m, grid->z1 + m, nt);
   free (fname);

   /* ...and their derivatives. */
   asprintf (&fname, "%s/d_%s.%s.tsv", prefix, spec_index[s]->name, name);
   rz_dump_3d (grid->d_dens[s], fname, "w", grid->r0 - m, grid->z0 - m, 
          grid->r1 + m, grid->z1 + m, nt);
   free (fname);
      }
    }
  }
}

void
cdr_dump_r (cdr_grid_t *grid, const char *prefix, const char *name, 
       FILE *infp, double sim_time)
{
  cdr_grid_t *child;
  char *cname;
  int i, nchilds;
  char codes[] = "abcdefghijklmnopqrstuvwxyz";
  FILE *fp;

  if (NULL == infp) {
    /* Root call */
    asprintf (&cname, "%s/tree.%s.dat", prefix, name);
    fp = fopen (cname, "w");
    if (NULL == fp) {
      fatal ("Could not open file %s/tree.%s.dat to write\n", 
        prefix, name);
    }
    free (cname);
    fprintf (fp, "time %g\n", sim_time);
  } else {
    fp = infp;
  }

  /* We want to be sure that the charge that we plot is actually the
   * charge and not its Fourier transform.
   */
  cdr_calc_charge_r (grid);

  fprintf (fp, "open %s %d %d %d %d %d %d %d\n", name, grid->r0, grid->z0, 
      grid->r1, grid->z1, grid->level, grid->ext_bound,
      cdr_output_margin);

  cdr_dump (grid, prefix, name);

  nchilds = grid_howmany_children ((grid_t *) grid);

  for (i = 0; i < nchilds; i++) {
    child = (cdr_grid_t *) grid_get_child ((grid_t*) grid, nchilds - i - 1);

    assert (NULL != child);

    if (i > sizeof(codes) - 2)
      asprintf (&cname, "%s{%d}", name, i);
    else
      asprintf (&cname, "%s%c", name, codes[i]);

    cdr_dump_r (child, prefix, cname, fp, sim_time);
    free (cname);
  }

  fprintf (fp, "close %s\n", name);

  if (NULL == infp) {
    /* Root call */
    fclose (fp);
  }
}

cdr_grid_t *
cdr_load_tree_r (const char *prefix, const char *name, FILE *infp)
{
  char command[16], gridname[32], *fname;
  int r0, z0, r1, z1, level, ext_bound, margin, nt, s;
  int open_close;
  cdr_grid_t *grid = NULL, *leaf;
  FILE *fp;

  debug (2, "cdr_load_tree_r (prefix = \"%s\", name = \"%s\", ...)\n",
    prefix, name);

  nt = max_ntheta == 1? 1: max_ntheta + 1;

  if (NULL == infp) {
    asprintf (&fname, "%s/tree.%s.dat", prefix, name);
    fp = fopen (fname, "r");
    if (NULL == fp) {
      fatal ("Could not open file %s/tree.%s.dat to read\n", 
        prefix, name);
    }
    free (fname);
  } else {
    fp = infp;
  }

  do {
    fscanf (fp, "%15s %15s", command, gridname);
    open_close = TRUE;

    /* We analyze some commands that do not create/close grids, now
     * restricted to set time.
     */
    if (0 == strcmp (command, "time")) {
      double new_time;
      sscanf (gridname, "%lf\n", &new_time);
      warning ("Starting time [%f] was read from %s/tree.%s.dat\n",
          new_time, prefix, name);

      start_t = new_time;
      open_close = FALSE;
    }
  } while (!open_close);

  if (0 == strcmp (command, "open")) {
    debug (3, "opening %s\n", gridname);

    fscanf (fp, "%d %d %d %d %d %d %d\n", 
       &r0, &z0, &r1, &z1, &level, &ext_bound, &margin);

    grid = cdr_new_3d_a (r0, z0, r1, z1, max_ntheta);
    grid->ext_bound = ext_bound;
    grid->level = level;

    for (s = 0; s < no_species; s++) {
      /* The densities */
      asprintf (&fname, "%s/%s.%s.tsv", prefix, spec_index[s]->name, gridname);
      debug (3, "Loading %s\n", fname);
      rz_dump_3d (grid->dens[s], fname, "r", 
        grid->r0 - margin, grid->z0 - margin, 
        grid->r1 + margin, grid->z1 + margin, nt);
      free (fname);
    }

    do {
      leaf = cdr_load_tree_r (prefix, name, fp);
      if (NULL != leaf) {
   assert (leaf->level == grid->level + 1);
   add_child (grid, leaf);
      }

    } while (NULL != leaf);
  } else if (0 == strcmp (command, "close")) {
    debug (3, "closing %s\n", gridname);
    grid = NULL;
  }

  if (NULL == infp) {
    fclose (fp);
  }
  return grid;
}

void
cdr_dump_frames (cdr_grid_t *grid, const char *prefix, const char *name)
{
  FILE *fp;
  char *fname;

  asprintf (&fname, "%s/frames.%s.tsv", prefix, name);
  fp = fopen (fname, "w");
  free (fname);

  if (fp == NULL) {
    warning ("Unable to open %s\n", fname);
    return;
  }
  aux_dump_frames_r (grid, fp);

  fclose (fp);
}

void
aux_dump_frames_r (cdr_grid_t *grid, FILE *fp)
{
  cdr_grid_t *child;
  int level;

  level = grid->level;

  fprintf (fp, "%g %g %g %g %d\n", 
      er_r_at (grid->r0 - 1, level), ez_z_at (grid->z0 - 1, level),
      er_r_at (grid->r1 - 1, level), ez_z_at (grid->z1 - 1, level), 
      level);

  iter_childs(grid, child) {
    aux_dump_frames_r (child, fp);
  }
}