rz_array.c

Go to the documentation of this file.

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

#include "parameters.h"
#include "proto.h"
#include "rz_array.h"
#include "species.h"

long int used_disk_space = 0;
static void check_disk_space(void);

rz_array_t *
rz_new_3d_a (int r0, int z0, int r1, int z1, int ntheta)
{
  int rmax, zmax;
  rz_array_t *array;

  debug (3, "rz_new_3d_a (%d, %d, %d, %d, %d)\n", r0, z0, r1, z1, ntheta);

  /* We allow a margin of two cells to set the boundary conditions (see below)*/
  r0 -= 2; z0 -= 2; r1 += 2; z1 += 2;

  rmax = r1 - r0;
  zmax = z1 - z0;

  assert (rmax > 0 && zmax > 0 && ntheta > 0);

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

  array = (rz_array_t *) xmalloc (sizeof(rz_array_t));

  if (ntheta != 1) {
    array->len = rmax * zmax * (ntheta + 4);
    array->theta0 = -2;
  } else {
    array->len = rmax * zmax;
    array->theta0 = 0;
  }

  /*  Initially all arrays are zeroed. */
  array->data = (REAL *) xcalloc (array->len, sizeof(REAL));


  /* These strides are to make the array compatible with array(ir, iz) 
     in FORTRAN. */
  array->strides[R_INDX] = 1;
  array->strides[Z_INDX] = rmax;
  array->strides[THETA_INDX] = rmax * zmax;

  array->dim = ntheta == 1? 2: 3;

  debug (3, "rmax = %d\n", rmax);
  debug (3, "strides = {%d, %d}\n", array->strides[0], array->strides[1]);

  array->r0 = r0;
  array->z0 = z0;
  /* Note that these are NOT the r0, z0 we received.  The idea here is that
     RZ(array, r0, z0) will produce the correct result, where the allocated
     array looks like (in the picture, r0, z0 ARE the received ones).
     Note: in the picture there is only one buffer cell, while we actually
     allocate two.

       +--------------+--------------+...+--------------+--------------+
       |r0 - 1, z0 - 1|              |   |              |              |
       +--------------+--------------+...+--------------+--------------+
       |              |   r0, z0     |   |              |              |
       +--------------+--------------+...+--------------+--------------+
                     ...             |...|             ...
       +--------------+--------------+...+--------------+--------------+
       |              |              |   |r1 - 1, z1 - 1|              |
       +--------------+--------------+...+--------------+--------------+
       |              |              |   |              |    r1, z1    |
       +--------------+--------------+...+--------------+--------------+

     but the rows z0 +/- 1 and the columns r0 +/- 1 do not belong to the
     physical space and are used only to specify boundary conditions.
  */

  array->nr = r1 - r0 - 4;
  array->nz = z1 - z0 - 4;
  array->ntheta = ntheta;

  array->host = NULL;

  return array;
}

rz_array_t *
rz_guest (rz_array_t *host, int r0, int z0, int r1, int z1)
{
  int ntheta, rmax, zmax;
  rz_array_t *array;

  ntheta = host->ntheta;
  array = (rz_array_t *) xmalloc (sizeof(rz_array_t));

  rmax = r1 - r0;
  zmax = z1 - z0;

  if (ntheta != 1) {
    array->len = rmax * zmax * (ntheta + 4);
  } else {
    array->len = rmax * zmax;
  }
  
  array->r0 = r0;
  array->z0 = z0;

  array->theta0 = host->theta0;

  array->strides[R_INDX] = host->strides[R_INDX];
  array->strides[Z_INDX] = host->strides[Z_INDX];
  array->strides[THETA_INDX] = host->strides[THETA_INDX];

  array->dim = host->dim;

  array->nr = r1 - r0 - 4;
  array->nz = z1 - z0 - 4;
  array->ntheta = host->ntheta;

  array->data = RZTP (host, r0, z0, array->theta0);
  array->host = host;

  return array;
}

void
rz_set_zero (rz_array_t *array)
{
  debug (3, "rz_set_zero\n");
  memset (array->data, 0, array->len * sizeof(REAL));
}

void
rz_set_periodic (rz_array_t *array)
{
  /* If the array is 2D, we do nothing. */
  if (array->ntheta == 1) return;

  rz_copy_modes (array, array->r0, array->z0, 
       array, array->r0, array->z0, 
       array->nr + 4, array->nz + 4, array->ntheta - 1, -1);

  rz_copy_modes (array, array->r0, array->z0, 
       array, array->r0, array->z0, 
       array->nr + 4, array->nz + 4, array->ntheta - 2, -2);

  rz_copy_modes (array, array->r0, array->z0, 
       array, array->r0, array->z0, 
       array->nr + 4, array->nz + 4, 0, array->ntheta);

  rz_copy_modes (array, array->r0, array->z0, 
       array, array->r0, array->z0, 
       array->nr + 4, array->nz + 4, 1, array->ntheta + 1);

}


void
rz_copy_bnd (rz_array_t *from, rz_array_t *to, 
        int sign, REAL *start_from, REAL *start_to, 
        int dim0, int inout_from, int inout_to,
        int dim1, int dim1_0, int dim1_1, 
        int dim2, int dim2_0, int dim2_1)
{
  int i, j;
  REAL *pfrom, *pto;

  for (i = dim1_0; i < dim1_1; i++) {
    for (j = dim2_0; j < dim2_1; j++) {
      pfrom = start_from + i * from->strides[dim1] + j * from->strides[dim2];
      pto = start_to + i * to->strides[dim1] + j * to->strides[dim2];

      *(pto + inout_to * to->strides[dim0]) = sign * (*pfrom);

      *(pto + 2 * inout_to * to->strides[dim0]) =  
   sign * (*(pfrom + inout_from * from->strides[dim0]));      
    }
  }
}


void
rz_set_bnd (rz_array_t *array, int sign, REAL *astart, int dim0, int inout,
       int dim1, int dim1_0, int dim1_1, 
       int dim2, int dim2_0, int dim2_1)
{
  rz_copy_bnd (array, array, sign, astart, astart, dim0, -inout, inout, 
          dim1, dim1_0, dim1_1, dim2, dim2_0, dim2_1); 
}


void
rz_free (rz_array_t *array)
{
  debug (3, "rz_free\n");

  if (NULL == array->host) free (array->data);
  free (array);
}

void 
rz_copy (rz_array_t *fro, int rfro, int zfro, 
        rz_array_t *to, int rto, int zto,
        int rn, int zn)
{
  int i;
  REAL *pfro, *pto;
  
  debug (3, "rz_copy\n");

  pfro = RZP (fro, rfro, zfro);
  pto = RZP (to, rto, zto);
  
  for (i = 0; i < zn; i++) {
    memcpy (pto, pfro, sizeof(REAL) * rn);

    pfro += fro->strides[Z_INDX];
    pto += to->strides[Z_INDX];
  }
}

void 
rz_copy_modes (rz_array_t *fro, int rfro, int zfro, 
          rz_array_t *to, int rto, int zto,
          int rn, int zn, int nmode_fro, int nmode_to)
{
  int i, j;
  REAL *pfro, *pto;
  
  debug (3, "rz_copy\n");

  pfro = RZTP (fro, rfro, zfro, nmode_fro);
  pto = RZTP (to, rto, zto, nmode_to);
  
  for (i = 0; i < zn; i++) {
    /* memcpy (pto, pfro, sizeof(REAL) * rn); Maybe not thread safe? */
    for (j = 0; j < rn; j++ ) *(pto + j) = *(pfro + j);

    pfro += fro->strides[Z_INDX];
    pto += to->strides[Z_INDX];
  }
}


void
rz_dump (rz_array_t *rz_array, const char *fname, const char *mode,
    int r0, int z0, int r1, int z1)
{
  FILE *fp;
  int ir, iz;
  int writing;
  double f;

  debug (3, "rz_dump(\"%s\", \"%s\", %d, %d, %d, %d)\n", fname, mode, 
    r0, z0, r1, z1);
 
  if (0 == strcmp (mode, "w")) {
    writing = TRUE;
  } else if (0 == strcmp (mode, "r")) {
    writing = FALSE;
  } else {
    fatal ("Unknown mode for rz_dump\n");
  }

  fp = fopen (fname, mode);

  if (fp == NULL) {
    warning("Unable to open %s for %s\n", fname, 
       writing? "writing": "reading");
    return;
  }

  for (ir = r0; ir < r1; ir++)
    for (iz = z0; iz < z1; iz++) {
      if (writing)
   fprintf (fp, "%15.5e\n", RZ(rz_array, ir, iz));
      else {
   if (fscanf (fp, "%lf", &f) != 1) {
     warning ("Error reading file %s, at line %d, (ir, iz) = (%d, %d)\n",
         fname, iz + ir * (z1 - z0), ir, iz);
     f = 0.0;
   }
   *RZP (rz_array, ir, iz) = f;
      }
    }
  
  used_disk_space += ftell (fp);
  fclose(fp);

  check_disk_space ();
}

void
rz_dump_3d (rz_array_t *rz_array, const char *fname, const char *mode,
       int r0, int z0, int r1, int z1, int ntheta)
{
  FILE *fp;
  int ir, iz, itheta;
  int writing;
  double f;

  debug (3, "rz_dump_3d(\"%s\", \"%s\", %d, %d, %d, %d, %d)\n", 
    fname, mode, r0, z0, r1, z1, ntheta);
 
  if (0 == strcmp (mode, "w")) {
    writing = TRUE;
  } else if (0 == strcmp (mode, "r")) {
    writing = FALSE;
  } else {
    fatal ("Unknown mode for rz_dump_3d\n");
  }

  fp = fopen (fname, mode);

  if (fp == NULL) {
    warning ("Unable to open %s for %s\n", fname, 
        writing? "writing": "reading");
    return;
  }

  for (itheta = 0; itheta < ntheta; itheta++) {
    for (ir = r0; ir < r1; ir++) {
      for (iz = z0; iz < z1; iz++) {
   if (writing)
     fprintf (fp, "%15.5e\n", RZT (rz_array, ir, iz, itheta));
   else {
     if (fscanf (fp, "%lf", &f) != 1) {
       warning ("Error reading file %s, at line %d, "
           "(ir, iz, itheta) = (%d, %d, %d)\n",
           fname, iz + (ir + itheta * (r1 - r0)) * (z1 - z0), 
           ir, iz, itheta);
       f = 0.0;
     }
     *RZTP (rz_array, ir, iz, itheta) = f;
   }
      }
    }
  }
  
  /* Add the disk space used to the total disk space that we have used so far
   */
  
  used_disk_space += ftell (fp);
  fclose(fp);

  check_disk_space ();
}
 
void
 rz_axis_dump (const char *fname, int x0, int x1, double delta)
{
  FILE *fp;
  int i;

  debug (3, "rz_axis_dump(\"%s\", %d, %d, %f)\n", fname, x0, x1, delta);
 
  fp = fopen (fname, "w");

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

  for (i = x0; i < x1; i++)
    fprintf (fp, "%15.5e\n", ((double) i + 0.5) * delta);


  used_disk_space += ftell (fp);
  fclose(fp);

  check_disk_space ();
}

static void
check_disk_space (void)
{
  if (used_disk_space > (((long int) max_disk_space_mb) << 20)) {
    fatal ("The disk space limit has been surpassed. "
      "Increase max_disk_space_mb if you really need more space\n");
  }
}