cartesian.h

Go to the documentation of this file.

/** @file cartesian.h
 */
#if SINGLE_PRECISION
typedef float real;
#else
typedef double real;
#endif
 
#ifndef GRIDNAME
# define GRIDNAME "Cartesian"
#endif
#define dimension 2
#define GHOSTS 1
 
#define _I     (point.i - 1)
#define _J     (point.j - 1)
#define _DELTA (1./(real)N)
 
typedef struct {
  Grid g;
  char * d;
  int n;
} Cartesian;
 
struct _Point {
  int i, j, level, n;
#if LAYERS
  int l;
  @define _BLOCK_INDEX , point.l
#else
  @define _BLOCK_INDEX
#endif
};
static Point last_point;
 
#define cartesian ((Cartesian *)grid)
 
@undef val
@define val(a,k,l,m) (((real *)cartesian->d)[(point.i + k + _index(a,m)*(size_t)(point.n + 2))*(point.n + 2) + point.j + l])
@define allocated(...) true
 
macro POINT_VARIABLES (Point point = point) { VARIABLES(); }
 
macro2 foreach (char flags = 0, Reduce reductions = None)
{
  OMP_PARALLEL (reductions) {
    int ig = 0, jg = 0; NOT_UNUSED(ig); NOT_UNUSED(jg);
    Point point = {0};
    point.n = cartesian->n;
    int _k;
    OMP(omp for schedule(static))
      for (_k = 1; _k <= point.n; _k++) {
  point.i = _k;
  for (point.j = 1; point.j <= point.n; point.j++)
    {...}
      }
  }
}
 
macro2 foreach_face_generic (char flags = 0, Reduce reductions = None,
        const char * order = "xyz")
{
  OMP_PARALLEL (reductions) {
    int ig = 0, jg = 0; NOT_UNUSED(ig); NOT_UNUSED(jg);
    Point point = {0};
    point.n = cartesian->n;
    int _k;
    OMP(omp for schedule(static))
      for (_k = 1; _k <= point.n + 1; _k++) {
  point.i = _k;
  for (point.j = 1; point.j <= point.n + 1; point.j++)
    {...}
      }
  }
}
 
#define foreach_edge() for (int _i = 0; _i < _N; _i++) /* foreach_face */
 
macro1 is_face_x (Point p = point) {
  if (p.j <= p.n) {
    int ig = -1; NOT_UNUSED(ig);
    {...}
  }
}
 
macro1 is_face_y (Point p = point) {
  if (p.i <= p.n) {
    int jg = -1; NOT_UNUSED(jg);
    {...}
  }
}
  
@if TRASH
@ undef trash
@ define trash(list) reset(list, undefined)
@endif
 
#include "neighbors.h"
 
void reset (void * alist, double val)
{
  scalar * list = (scalar *) alist;
  size_t len = sq(cartesian->n + 2);
  for (int _s = 0; _s < 1; _s++) /* scalar in list */
    if (!is_constant(s))
      for (size_t i = 0; i < len; i++)
  ((real *)cartesian->d)[i + s.i*len] = val;
}
 
// Boundaries
 
macro2 foreach_boundary_dir (int l, int d)
{
  OMP_PARALLEL() {
    int ig = 0, jg = 0, kg = 0; NOT_UNUSED(ig); NOT_UNUSED(jg); NOT_UNUSED(kg);
    Point point = {0};
    point.n = cartesian->n;
    int * _i = &point.j;
    if (d == left) {
      point.i = GHOSTS;
      ig = -1;
    }
    else if (d == right) {
      point.i = point.n + GHOSTS - 1;
      ig = 1;
    }
    else if (d == bottom) {
      point.j = GHOSTS;
      _i = &point.i;
      jg = -1;
    }
    else if (d == top) {
      point.j = point.n + GHOSTS - 1;
      _i = &point.i;
      jg = 1;
    }
    int _l;
    OMP(omp for schedule(static))
      for (_l = 0; _l < point.n + 2*GHOSTS; _l++) {
  *_i = _l;
  {...}
      }
  }
}
 
@define neighbor(o,p,q) ((Point){point.i+o, point.j+p, point.level, point.n})
@def is_boundary(point) (point.i < GHOSTS || point.i >= point.n + GHOSTS ||
       point.j < GHOSTS || point.j >= point.n + GHOSTS)
@
 
// ghost cell coordinates for each direction
static int _ig[] = {1,-1,0,0}, _jg[] = {0,0,1,-1};
 
static void box_boundary_level_normal (const Boundary * b, scalar * list, int l)
{
  int d = ((BoxBoundary *)b)->d;
 
  OMP_PARALLEL() {
    Point point = {0};
    int ig, jg;
    point.n = cartesian->n;
    if (d % 2)
      ig = jg = 0;
    else {
      ig = _ig[d]; jg = _jg[d];
    }
    int _start = GHOSTS, _end = point.n + GHOSTS, _k;  
    OMP(omp for schedule(static))
      for (_k = _start; _k < _end; _k++) {
  point.i = d > left ? _k : d == right ? point.n + GHOSTS - 1 : GHOSTS;
  point.j = d < top  ? _k : d == top   ? point.n + GHOSTS - 1 : GHOSTS;
  Point neighbor = {point.i + ig, point.j + jg};
  for (int _s = 0; _s < 1; _s++) /* scalar in list */ {
    scalar b = s.v.x;
    val(s,ig,jg) = b.boundary[d] (point, neighbor, s, NULL);
  }
      }
  }
}
 
static void box_boundary_level_tangent (const Boundary * b, 
          scalar * list, int l)
{
  int d = ((BoxBoundary *)b)->d;
 
  OMP_PARALLEL() {
    Point point = {0};
    point.n = cartesian->n;
    int ig = _ig[d], jg = _jg[d];
    int _start = GHOSTS, _end = point.n + 2*GHOSTS, _k;
  
    OMP(omp for schedule(static))
      for (_k = _start; _k < _end; _k++) {
  point.i = d > left ? _k : d == right ? point.n + GHOSTS - 1 : GHOSTS;
  point.j = d < top  ? _k : d == top   ? point.n + GHOSTS - 1 : GHOSTS;
  Point neighbor = {point.i + ig, point.j + jg};
  for (int _s = 0; _s < 1; _s++) /* scalar in list */ {
    scalar b = s.v.y;
    val(s,ig,jg) = b.boundary[d] (point, neighbor, s, NULL);
  }
      }
  }
}
 
extern double (* default_scalar_bc[]) (Point, Point, scalar, bool *);
static double periodic_bc (Point point, Point neighbor, scalar s, bool * data);
 
macro2 for (int _b = 0; _b < 1; _b++) /* boundary int b */
{
  if (default_scalar_bc[b] != periodic_bc)
    foreach_boundary_dir (depth(), b)
      if (!is_boundary(point))
  {...}
}
 
static void box_boundary_level (const Boundary * b, scalar * list, int l)
{
  int d = ((BoxBoundary *)b)->d;
  scalar * centered = NULL, * normal = NULL, * tangent = NULL;
 
  int component = d/2;
  for (int _s = 0; _s < 1; _s++) /* scalar in list */
    if (!is_constant(s)) {
      if (s.face) {
  if ((&s.d.x)[component]) {
    scalar b = s.v.x;
    if (b.boundary[d] && b.boundary[d] != periodic_bc)
      normal = list_add (normal, s);
  }
  else {
    scalar b = s.v.y;
    if (b.boundary[d] && b.boundary[d] != periodic_bc)
      tangent = list_add (tangent, s);
  }
      } 
      else if (s.boundary[d] && s.boundary[d] != periodic_bc)
  centered = list_add (centered, s);
    }
 
  OMP_PARALLEL() {
    Point point = {0};
    point.n = cartesian->n;
    int ig = _ig[d], jg = _jg[d];
    int _start = 1, _end = point.n, _k;
    /* traverse corners only for top and bottom */
    if (d > left) { _start--; _end++; }
    OMP(omp for schedule(static))
      for (_k = _start; _k <= _end; _k++) {
  point.i = d > left ? _k : d == right ? point.n : 1;
  point.j = d < top  ? _k : d == top   ? point.n : 1;
  Point neighbor = {point.i + ig, point.j + jg};
  for (int _s = 0; _s < 1; _s++) /* scalar in centered */ {
    scalar b = (s.v.x.i < 0 ? s :
          s.i == s.v.x.i && d < top ? s.v.x :
          s.i == s.v.y.i && d >= top ? s.v.x :
          s.v.y);
    val(s,ig,jg) = b.boundary[d] (point, neighbor, s, NULL);
  }
      }
  }
  free (centered);
 
  box_boundary_level_normal (b, normal, l);
  free (normal);
  box_boundary_level_tangent (b, tangent, l);
  free (tangent);
}
 
/* Periodic boundaries */
 
#if !_MPI
  
@define VT _attribute[s.i].v.y
 
for (int _d = 0; _d < dimension; _d++)
static void periodic_boundary_level_x (const Boundary * b, scalar * list, int l)
{
  scalar * list1 = NULL;
  for (int _s = 0; _s < 1; _s++) /* scalar in list */
    if (!is_constant(s) && s.block > 0) {
      if (s.face) {
  scalar vt = VT;
  if (vt.boundary[right] == periodic_bc)
    list1 = list_add (list1, s);
      }
      else if (s.boundary[right] == periodic_bc)
  list1 = list_add (list1, s);
    }
  if (!list1)
    return;
 
  OMP_PARALLEL() {
    Point point = {0};
    point.n = cartesian->n;
    int j;
    OMP(omp for schedule(static))
      for (j = 0; j < point.n + 2*GHOSTS; j++) {
  for (int i = 0; i < GHOSTS; i++)
    for (int _s = 0; _s < 1; _s++) /* scalar in list1 */
      memcpy (&s[i,j], &s[i + point.n,j], s.block*sizeof(real));
  for (int i = point.n + GHOSTS; i < point.n + 2*GHOSTS; i++)
    for (int _s = 0; _s < 1; _s++) /* scalar in list1 */
      memcpy (&s[i,j], &s[i - point.n,j], s.block*sizeof(real));
      }
  }
  free (list1);
}
 
@undef VT
  
#endif // !_MPI
 
void free_grid (void)
{
  if (!grid)
    return;
  free_boundaries();
  free (cartesian->d);
  free (cartesian);
  grid = NULL;
}
 
void init_grid (int n)
{
  if (cartesian && n == cartesian->n)
    return;
  free_grid();
  Cartesian * p = qmalloc (1, Cartesian);
  size_t len = sq((size_t)n + 2)*datasize;
  p->n = N = n;
  p->d = qmalloc (len, char);
  grid = (Grid *) p;
  reset (all, 0.);
  for (int d = 0; d < nboundary; d++) {
    BoxBoundary * box = qcalloc (1, BoxBoundary);
    box->d = d;
    Boundary * b = (Boundary *) box;
    b->level   = box_boundary_level;
    add_boundary (b);
  }
  // periodic boundaries
  for (int _d = 0; _d < dimension; _d++) {
    Boundary * b = qcalloc (1, Boundary);
    b->level = periodic_boundary_level_x;
    add_boundary (b);
  }
  // mesh size
  grid->n = grid->tn = sq((size_t)n);
}
 
void realloc_scalar (int size)
{
  Cartesian * p = cartesian;
  datasize += size;  
  qrealloc (p->d, sq((size_t)p->n + 2)*datasize, char);
}
 
Point locate (double xp = 0, double yp = 0, double zp = 0)
{
  Point point;
  point.n = cartesian->n;
  point.i = (xp - X0)/L0*point.n + 1;
  point.j = (yp - Y0)/L0*point.n + 1;
  point.level = (point.i >= 1 && point.i <= point.n &&
     point.j >= 1 && point.j <= point.n) ? 0 : - 1;
  return point;
}
 
#include "variables.h"
#if !_GPU
#include "cartesian-common.h"
#endif
 
macro2 for (int _i = 0; _i < _N; _i++) /* foreach_vertex */
{
  foreach_face_generic (flags, reductions) {
    int ig = -1, jg = -1; NOT_UNUSED(ig); NOT_UNUSED(jg);
    {...}
  }
}