photo.c

Go to the documentation of this file.

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

#include "cdr.h"
#include "interpol2.h"
#include "mapper.h"
#include "parameters.h"
#include "photo.h"
#include "poisson.h"
#include "proto.h"
#include "rz_array.h"
#include "species.h"

void photo_copy (mapper_t *mapper, grid_t *source, grid_t *target, 
       int ir, int iz, int itheta);
void photo_coarsen (mapper_t *mapper, grid_t *source, grid_t *target, 
          int ir, int iz, int itheta);
int photo_interpol_set (mapper_t *mapper, grid_t *source, 
          interpol_t *interpol,
          int pr, int pz, int itheta);
void photo_interpol (mapper_t *mapper, grid_t *source, grid_t *target, 
           interpol_t *interpol, 
           int ir, int iz, int itheta);

photo_term_t *photo_terms = NULL;

mapper_t photo_mapper = {&interpol_bilin, photo_coarsen, photo_copy,
          photo_interpol_set, photo_interpol};

mapper_t *photo_mappers[] = {&photo_mapper, NULL};

pois_problem_t *pois_photo_1;
pois_problem_t *pois_photo_2;
extern pois_problem_t *pois_electrostatic;

void
photo_init ()
{
  /* Since it makes no sense, we use photo_extra_levels < 0 to say that
     the conditions for the photoionization are the same as for
     the electrostatic problem.  This is needed for backwards compatibility. */
  if (extra_photo_levels < 0) 
  {
    pois_photo_1 = pois_electrostatic;
    pois_photo_2 = pois_electrostatic;
    return;
  }

  /* Two different sets of refinement criteria for the two photoionization
     terms. If extra_photo_levels_2 < 0, then pois_photo_2 = pois_photo_1 */
  pois_photo_1 = xmalloc (sizeof (pois_problem_t));

  pois_photo_1->max_level = photo_max_level;
  pois_photo_1->extra_levels = extra_photo_levels;
  pois_photo_1->max_error = photo_max_error;

  pois_photo_1->bnd_right = photo_bnd_right;
  pois_photo_1->bnd_top = photo_bnd_top;
  pois_photo_1->bnd_bottom = photo_bnd_bottom;

  if (extra_photo_levels_2 < 0)
  {
     pois_photo_2 = pois_photo_1;
     return;
  }

  pois_photo_2 = xmalloc (sizeof (pois_problem_t));

  pois_photo_2->max_level = photo_max_level_2;
  pois_photo_2->extra_levels = extra_photo_levels_2;
  pois_photo_2->max_error = photo_max_error_2;

  pois_photo_2->bnd_right = photo_bnd_right_2;
  pois_photo_2->bnd_top = photo_bnd_top_2;
  pois_photo_2->bnd_bottom = photo_bnd_bottom_2;
}


void
photo_copy (mapper_t *mapper, grid_t *source, grid_t *target, 
       int ir, int iz, int itheta)
{
  cdr_grid_t *cdr;
  pois_grid_t *pois;

  cdr = (cdr_grid_t*) target;
  pois = (pois_grid_t*) source;

  RZT (cdr->photo, ir, iz, itheta) = RZ (pois->phi, ir, iz);
}

void 
photo_coarsen (mapper_t *mapper, grid_t *source, grid_t *target, 
          int ir, int iz, int itheta)
{
  cdr_grid_t *cdr;
  pois_grid_t *pois;
  int level_diff, z, r;

  cdr = (cdr_grid_t*) target;
  pois = (pois_grid_t*) source;

  level_diff = pois->level - cdr->level;

  z = (iz << level_diff) + (1 << (level_diff - 1));
  r = (ir << level_diff) + (1 << (level_diff - 1));

  if (grid_contains (source, r, z, GRID_INSIDE) && 
      grid_contains (source, r - 1, z, GRID_INSIDE) &&
      grid_contains (source, r, z - 1, GRID_INSIDE) &&
      grid_contains (source, r - 1, z - 1, GRID_INSIDE)) {
    
    RZT (cdr->photo, ir, iz, itheta)  = 0.25 *
      (RZ (pois->phi, r, z) 
       + RZ (pois->phi, r - 1, z) 
       + RZ (pois->phi, r, z - 1) 
       + RZ (pois->phi, r - 1, z - 1));
  }

}

int 
photo_interpol_set (mapper_t *mapper, grid_t *source, interpol_t *interpol,
            int pr, int pz, int itheta)
{
  pois_grid_t *pois;

  pois = (pois_grid_t*) source;

  interpol_set_stencil (interpol,
         r_at (pr, pois->level),
         z_at (pz, pois->level),
         RZ (pois->phi, pr, pz), 
         RZ (pois->phi, pr, pz + 1),
         RZ (pois->phi, pr + 1, pz), 
         RZ (pois->phi, pr + 1, pz + 1));

  return TRUE;
}

void
photo_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->photo, ir, iz, itheta) = interpol_apply (interpol, r, z);
}


void
photo_register (double A, double lambda)
{
  photo_term_t *t;

  t = (photo_term_t *) xmalloc (sizeof (photo_term_t));

  t->A = A;
  t->lambda = lambda;
  t->next = photo_terms;

  photo_terms = t;
}

void
photo_copy_list (photo_term_t *src, photo_term_t **dest)
{
  photo_term_t *ptr_term;
  photo_term_t *p = NULL, *p1;

  *dest = NULL;

  while (src) {
    p1 = (photo_term_t *) xmalloc (sizeof (photo_term_t));
    if (p) {
      p->next = p1;
    } else {
      *dest = p1;
    }

    p = p1;
    p->A = src->A;
    p->lambda = src->lambda;
    src = src->next;
  }

  if (p) p->next = NULL;
}


void
photo_unregister_all (void)
{
  photo_term_t *t, *next;

  for (t = photo_terms; t; t = next) {
    next = t->next;
    free (t);
  }

  photo_terms = NULL;
}

void
photo_dft_r (cdr_grid_t *grid, int sign)
{
  cdr_grid_t *leaf;

  debug (3, "photo_dft_r(" grid_printf_str ", %d)\n", 
    grid_printf_args(grid), sign);

  iter_childs (grid, leaf) {
    photo_dft_r (leaf, sign);
  }

  dft_transform (grid->photo, grid->photo, sign);
}

void
photo_copy_source (cdr_grid_t *grid)
{
  int ir, iz, itheta;

   debug(3, "photo_copy_source (" grid_printf_str ")\n",
    grid_printf_args(grid));

#pragma omp parallel
  {
#pragma omp for private (ir, iz)
    iter_grid_3d_n (grid, ir, iz, itheta, 2) {
      RZT (grid->charge, ir, iz, itheta) = 
   RZT (grid->d_dens[ions], ir, iz, itheta) / grid->ntheta;
    }
  }
}

mk_recursive (photo_copy_source, cdr_grid_t)


void
photo_add_term (photo_term_t *term, cdr_grid_t *cdr)
{
  int s, ir, iz, itheta;
  int updated[2] = {electrons, photo_ions};
  
  debug (3, "photo_add_term (" photo_printf_str ", " grid_printf_str ")\n",
    photo_printf_args(term), grid_printf_args(cdr));
    
#pragma omp parallel
    {
#pragma omp for private (ir, iz, s)
      iter_grid_3d_n (cdr, ir, iz, itheta, 2) {
   for (s = 0; s < 2; s++) {
     RZT (cdr->d_dens[updated[s]], ir, iz, itheta) += 
       term->A * RZT (cdr->photo, ir, iz, itheta);
   }
      }
    }  
}

void
photo_add_term_r (photo_term_t *term, cdr_grid_t *cdr)
{
  cdr_grid_t *child;

  photo_add_term (term, cdr);

  iter_childs (cdr, child) {
    photo_add_term_r (term, child);
  }
}

pois_grid_t **
photo_calc_term (photo_term_t *term, cdr_grid_t *cdr, int i)
{
  /* Call to the generic Poisson/Helmholtz solver */
  if (i == 1) {
   return pois_gen_solve_a (cdr, pois_photo_2, photo_mappers, term->lambda);
  } else {
   return pois_gen_solve_a (cdr, pois_photo_1, photo_mappers, term->lambda);
  }
}

void
photo_calc (photo_term_t *terms, cdr_grid_t *cdr)
{
  pois_grid_t **pois_modes;
  photo_term_t *term;
  // photo_term_t *ptr_term;

  photo_copy_source_r (cdr);
  
  if (cdr->ntheta != 1)
    cdr_dft_charge_r (cdr, 1);
  
  // for (ptr_term = terms; ptr_term != NULL; ptr_term = ptr_term->next) {
  //   printf("terms: A=%g lambda=%g ptr=%d\n", ptr_term->A, ptr_term->lambda, (int)ptr_term->next);
  // }

  int j = 0;
  for (term = terms; term; term = term->next ) {
    int i;
    pois_modes = photo_calc_term (term, cdr, j);
    j++;
    if (cdr->ntheta != 1)
      photo_dft_r (cdr, -1);

    photo_add_term_r (term, cdr);
    debug (3, "photo_calc (" photo_printf_str ", " grid_printf_str ")\n",
          photo_printf_args(term), grid_printf_args(cdr));

    /* Free the allocated memory. */
    for (i = 0; i < max_ntheta; i++) {      
      pois_free_r (pois_modes[i]);
    }
    free (pois_modes);
  }
}

void
photo_load_file (char *fname)
{
  FILE *fp;
  double A, lambda;
  int c;
  int i;

  printf("\n");
  printf ("Loading photoionization data from `%s'...\n", fname);
  printf("\n");
  fp = fopen (fname, "r");

  if (NULL == fp) {
    warning ("Unable to open photoionization file `%s'\n", fname);
    exit (-1);
    return;
  }

  
  i=0;
  do {
    c = fscanf (fp, "%lf %lf", &A, &lambda);
    if (c != 2) {
      break;
    }

    printf ("Registring photoionization term A = %.5g, lambda = %.4g\n",
       A, lambda);

    photo_register (A, lambda);
    i++;
  } while (TRUE);
  printf("\n");

  fclose (fp);
}