mapper.c

Go to the documentation of this file.

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#include "cstream.h"
#include "grid.h"
#include "interpol2.h"
#include "mapper.h"
#include "proto.h"
#include "species.h"

static void map_same (mapper_t **mappers, 
            grid_t *source, grid_t *target, int ntheta,
            int r0, int z0, int r1, int z1);
static void map_interpol (mapper_t **mappers, grid_t *source, grid_t *target, 
           int ntheta, int r0, int z0, int r1, int z1, 
           int level_diff);
static void interpol_inner (mapper_t **mappers, grid_t *source, grid_t *target,
             int pr, int pz, int ntheta,  
             interpol_t **interpols, 
             int *extend_r, int *extend_z);
static void map_coarsen (mapper_t **mappers, grid_t *source, grid_t *target, 
          int ntheta, int r0, int z0, int r1, int z1, 
          int level_diff);


void
map_grid (mapper_t **mappers, grid_t *source, grid_t *target, int ntheta, 
     int copy, int interpol, int coarsen, int s_buff, int t_buff)
{
  int r0, z0, r1, z1, level_diff;
  int br0, bz0, br1, bz1;
  mapper_t **m;
  int overlap;

  debug (3, "map_grid (source = " grid_printf_str 
    ",\n\t dest = " grid_printf_str 
    ",\n\t s_buff = %d, t_buff = %d)\n",
    grid_printf_args (source), grid_printf_args (target),
    s_buff, t_buff);

  overlap = grid_overlap (source, target, s_buff, t_buff, 
           &r0, &z0, &r1, &z1, &level_diff);

  if (!overlap) {
    debug (3, " -> no overlap\n");
    return;
  }
  
  if (level_diff == 0 && copy) {
    map_same (mappers, source, target, ntheta, r0, z0, r1, z1);
  } else if (level_diff < 0 && interpol) {
    /* When we interpolate electric fields, the stencil used in the coarse
       grid does not map exactly to cells in the finer grid by the standard
       rules.  Hence may have to extend the range of fine-grid cells that
       can be affected by a coarse grid. 
       But later we must be careful, because not all fine-grid cells can
       be interpolated.
    */
    for (m = mappers; *m; m++) {
      int max_shift;
     
      max_shift =  MYMAX(abs((*m)->shift_r) << (-level_diff - 1),
          abs((*m)->shift_z) << (-level_diff - 1));

      debug (3, "max_shift = %d, shift_r = %d, shift_z = %d\n", 
        max_shift, (*m)->shift_r, (*m)->shift_z);

      grid_overlap_with_shifts (source, target, 
            s_buff, 
            MYMAX (t_buff, max_shift),
            &br0, &bz0, &br1, &bz1, &level_diff,
            (*m)->shift_r, (*m)->shift_z);
      if (br0 < r0) r0 = br0;
      if (bz0 < z0) z0 = bz0;
      if (br1 > r1) r1 = br1;
      if (bz1 > z1) z1 = bz1;
    }
    map_interpol (mappers, source, target, ntheta, r0, z0, r1, z1, 
        -level_diff);
  } else if (level_diff > 0 && coarsen) {
    map_coarsen (mappers, source, target, ntheta, r0, z0, r1, z1, level_diff);
  }
}


void
map_grid_r (mapper_t **mappers, grid_t *source, grid_t *target, int ntheta,
       int copy, int interpol, int coarsen, int s_buff, int t_buff)
{
  grid_t *child;

  debug (3, "map_grid_r (...)\n");

  map_grid (mappers, source, target, ntheta, copy, interpol, coarsen,
       s_buff, t_buff);

  iter_childs (source, child) {
    map_grid_r (mappers, child, target, ntheta, copy, interpol, coarsen,
      s_buff, t_buff);
  }
}

void
map_trees_r (mapper_t **mappers, grid_t *source, grid_t *target, int ntheta,
        int copy, int interpol, int coarsen, int s_buff, int t_buff)
{
  grid_t *child;

  debug (3, "map_trees_r (...)\n");

  iter_childs (target, child) {
    map_trees_r (mappers, source, child, ntheta, copy, interpol, coarsen,
      s_buff, t_buff);
  }

  map_grid_r (mappers, source, target, ntheta, copy, interpol, coarsen,
         s_buff, t_buff);
}

static void
map_same (mapper_t **mappers, 
     grid_t *source, grid_t *target, int ntheta,
     int r0, int z0, int r1, int z1)
{
  int ir, iz;
  mapper_t **m;

  debug (3, "map_same (source = " grid_printf_str 
    ",\n\t dest = " grid_printf_str 
    ",\n\t r0 = %d, z0 = %d, r1 = %d, z1 = %d)\n",
    grid_printf_args (source), grid_printf_args (target),
    r0, z0, r1, z1);

  for (ir = r0; ir < r1; ir++) {
    for (iz = z0; iz < z1; iz++) {
      for (m = mappers; *m; m++) {
   (*m)->copy (*m, source, target, ir, iz, ntheta);
      }
    }
  }
  debug (3, " < map_same (...)\n");
}

#define MAX_MAPPERS  32

static void
map_interpol (mapper_t **mappers, grid_t *source, grid_t *target, int ntheta,
          int r0, int z0, int r1, int z1, int level_diff)
{
  int pr, pz;
  int i, extend_r, extend_z;
  interpol_t *interpols[MAX_MAPPERS];
  mapper_t **m;

  debug (3, "map_interpol (source = " grid_printf_str 
    ",\n\t dest = " grid_printf_str 
    ",\n\t r0 = %d, z0 = %d, r1 = %d, z1 = %d, level_diff = %d)\n",
    grid_printf_args (source), grid_printf_args (target),
    r0, z0, r1, z1, level_diff);

  for (m = mappers, i = 0; *m; m++, i++) 
    interpols[i] = interpol_new_a (dr[source->level], dz[source->level], 
               (*m)->interpol_method);

  /* The +- 1 are there because even when we are including the
     boundaries to calculate the overlap, we have to let one buffer cell
     for the 9-point interpolation. */
  r0 = MYMAX (r0 >> level_diff, source->r0 - 1);

  /* Explanation of the formula used: r1 is the index of the first cell not
     in the overlap 
     -> r1 - 1 is the last cell in the overlap 
     -> (r1 - 1) >> level_diff indexes the last cell in the coarser grid
        that contains part of the overlap.
     -> since the loop has to include the former index, we add one.
  */
  r1 = MYMIN (((r1 - 1) >> level_diff) + 1, source->r1 + 1);
  z0 = MYMAX (z0 >> level_diff, source->z0 - 1);
  z1 = MYMIN (((z1 - 1) >> level_diff) + 1, source->z1 + 1);


  debug (4, "r0 = %d, z0 = %d, r1 = %d, z1 = %d\n",
    r0, z0, r1, z1);

  /* assert (r0 < r1 && z0 < z1); */

  for (pr = r0, extend_r = TRUE; pr < r1; pr++) {
    /* The 0 initialization is to shut up the compiler. */
    int extend_r_in = 0;
    for (pz = z0, extend_z = TRUE; pz < z1; pz++) {
      debug (6, "pr = %d, pz = %d, extend_r = %d, extend_z = %d\n", 
        pr, pz, extend_r, extend_z);
      extend_r_in = extend_r;

      interpol_inner (mappers, source, target, pr, pz, ntheta,
            interpols, &extend_r_in, &extend_z);
    }     
    extend_r = extend_r_in;
  }

  for (m = mappers, i= 0; *m; m++, i++) 
    interpol_free (interpols[i]);

  debug (3, "< map_interpol (...)\n");
}

static void
interpol_inner (mapper_t **mappers, grid_t *source, grid_t *target, 
      int pr, int pz, int ntheta,
      interpol_t **interpols, int *extend_r, int *extend_z) 
{
  int level_diff;
  int i, ir, iz;
  int rmin, rmax, zmin, zmax;
  int inside;

  mapper_t **m;

  level_diff = target->level - source->level;

  debug (6, "interpol_inner (source = " grid_printf_str ",\n\t\t target = "
    grid_printf_str ",\n\t\t pr = %d, pz = %d, ntheta = %d)\n",
    grid_printf_args (source), grid_printf_args (target), pr, pz, ntheta);

  for (m = mappers, i = 0; *m; m++, i++) {
    int fine_shift_r, fine_shift_z;

    inside = (*m)->interpol_set (*m, source, interpols[i], pr, pz, ntheta);
  
    if (!inside) continue;

    fine_shift_r = (*m)->shift_r << (level_diff - 1);
    fine_shift_z = (*m)->shift_z << (level_diff - 1);

    debug (6, "shift_r = %d, shift_z = %d\n", (*m)->shift_r, (*m)->shift_z);

    rmin = (pr << level_diff) + fine_shift_r;
    /* (*m)->shift_z tells us if the grid is staggered.  If it is not,
       it doesn't make sense to extend it. */
    if (*extend_r && (*m)->shift_z != 0) {
      *extend_r = FALSE;
      rmin--;
    }

    zmin = (pz << level_diff) + fine_shift_z;
    if (*extend_z && (*m)->shift_r != 0) {
      *extend_z = FALSE;
      zmin--;
    }

    /* The 2 are there because there is no problem in writing in the
       buffer margin of the target grid. */
    rmin = MYMAX (target->r0 - 2, rmin);
    rmax = MYMIN (target->r1 + 2, ((pr + 1) << level_diff) + fine_shift_r);
    zmin = MYMAX (target->z0 - 2, zmin);
    zmax = MYMIN (target->z1 + 2, ((pz + 1) << level_diff) + fine_shift_z);
    

    debug(6, "\trmin = %d, rmax = %d, zmin = %d, zmax = %d\n",
     rmin, rmax, zmin, zmax);

    for (ir = rmin; ir < rmax; ir++) {
      for (iz = zmin; iz < zmax; iz++) {
   (*m)->interpol (*m, source, target, interpols[i], 
         ir, iz, ntheta);
      }   
    }
  }
}

static void
map_coarsen (mapper_t **mappers, grid_t *source, grid_t *target, int ntheta,
        int r0, int z0, int r1, int z1, int level_diff)
{
  int i, ir, iz;
  int rmin, rmax, zmin, zmax;
  mapper_t **m;

  debug (3, "map_coarsen (source = " grid_printf_str 
    ",\n\t dest = " grid_printf_str 
    ",\n\t r0 = %d, z0 = %d, r1 = %d, z1 = %d, level_diff = %d)\n",
    grid_printf_args (source), grid_printf_args (target),
    r0, z0, r1, z1, level_diff);

  rmin = MYMAX (r0 >> level_diff, target->r0 - 2);

  /* See above for an explanation for this formula. */
  rmax = MYMIN (((r1 - 1) >> level_diff) + 1, target->r1 + 2);
  zmin = MYMAX (z0 >> level_diff, target->z0 - 2);
  zmax = MYMIN (((z1 - 1) >> level_diff) + 1, target->z1 + 2);
  
  for (iz = zmin; iz < zmax; iz++) {
    for (ir = rmin; ir < rmax; ir++) {
      for (m = mappers, i = 0; *m; m++, i++) {
   (*m)->coarsen (*m, source, target, ir, iz, ntheta);
      }
    } 
  }
  debug (3, "< map_coarsen (...)\n");
}