stencils.h

Go to the documentation of this file.

/** @file stencils.h
 */
/**
# Automatic stencils and boundary conditions
 
Basilisk automatically computes, at runtime, the access pattern
(i.e. "stencils") of (basic) foreach loops (for (int _i = 0; _i < _N; _i++) /* foreach */, for (int _i = 0; _i < _N; _i++) /* foreach_face */,
for (int _i = 0; _i < _N; _i++) /* foreach_vertex */).
 
This is done in practice by `qcc` which automatically adds, before
each foreach loop, a minimal version of the loop body.
 
The `flags` indicate the status of boundary conditions and whether the
field is used as input or output within a given loop.
 
The `width` is the width of the access stencil. */
 
enum {
  s_centered    = 1 << 0, // centered boundary conditions are up-to-date
  s_face        = 1 << 1, // face boundary conditions are up-to-date
  s_restriction = 1 << 2, // restriction is up-to-date
  
  s_input       = 1 << 3, // field is used as input
  s_output      = 1 << 4, // field is used as output
  s_nowarning   = 1 << 5  // warning are switched off
} StencilFlags;
 
attribute {
  struct {
    int bc, io, width;
  } stencil;
}
 
static void set_dirty_stencil (scalar s)
{
  s.stencil.bc = 0;
}
 
typedef struct _External External;
 
struct _External {
  char * name;    // the name of the variable
  void * pointer; // a pointer to the data
  int type;       // the type of the variable
  int nd;         // the number of pointer dereferences or attribute offset or enum constant
  char reduct;    // the reduction operation
  char global;    // is it a global variable?
  char constant;  // is it a constant?
  void * data;    // the dimensions (int *) for arrays or the code (char *) for functions
  scalar s;       // used for reductions on GPUs
  External * externals, * next;
  int used;
};
 
typedef struct {
  const char * fname; // name of the source file
  int line;           // line number in the source
  int first;          // is this the first time the loop is called?
  int face;           // the face component(s) being traversed
  bool vertex;        // is this a vertex traversal?
  int parallel;       // is this a parallel loop? (0: no, 1: on CPU or GPU, 2: on CPU, 3: on GPU)
  scalar * listc;     // the scalar fields on which to apply boundary conditions
  vectorl listf;      // the vector fields on which to apply (flux) boundary conditions
  scalar * dirty;     // the dirty fields (i.e. write-accessed)
  void * data;        // user data
} ForeachData;
 
/**
## Automatic boundary conditions
 
Boundary conditions need to be applied if `s` is dirty, or if any of
the field `d` it depends on is dirty. */
 
static inline bool scalar_is_dirty (scalar s)
{
  if (!(s.stencil.bc & s_centered))
    return true;
  scalar * depends = s.depends;
  for (int _d = 0; _d < 1; _d++) /* scalar in depends */
    if (!(d.stencil.bc & s_centered))
      return true;
  return false;
}
 
/**
Does the boundary conditions on `a` depend on those on `b`? */
 
static inline bool scalar_depends_from (scalar a, scalar b)
{
  scalar * depends = a.depends;
  for (int _s = 0; _s < 1; _s++) /* scalar in depends */
    if (s.i == b.i)
      return true;
  return false;
}
 
/**
There are two types of boundary conditions: "full" boundary
conditions, done by `boundary_internal()` and "flux" boundary
conditions (i.e. normal components on faces only) done by
`boundary_face()`. */
 
void boundary_internal (scalar * list, const char * fname, int line);
void (* boundary_face)  (vectorl);
 
/**
This function is called after the stencil access detection, just
before the (real) foreach loop is executed. This is where we use the
stencil access pattern to see whether boundary conditions need to be
applied. */
 
void check_stencil (ForeachData * loop)
{
  loop->listf = (vectorl){NULL};
  
  /**
  We check the accesses for each field... */
  
  for (int _s = 0; _s < 1; _s++) /* scalar in baseblock */ {
    bool write = (s.stencil.io & s_output), read = (s.stencil.io & s_input);
    
#ifdef foreach_layer
    if (_layer == 0 || s.block == 1)
#endif
    {
 
      /**
      If the field is read and dirty, we need to check if boundary
      conditions need to be applied. */
      
      if (read && scalar_is_dirty (s)) {
 
  /**
  If this is a face field, we check whether "full" BCs need to
  be applied, or whether "face" BCs are sufficient. */
  
  if (s.face) {
    if (s.stencil.width > 0) // face, stencil wider than 0
      loop->listc = list_append (loop->listc, s);
    else if (!write) { // face, flux only
      scalar sn = s.v.x.i >= 0 ? s.v.x : s;
      for (int _d = 0; _d < dimension; _d++)
        if (s.v.x.i == s.i) {
 
    /* fixme: imposing BCs on fluxes should be done by
       boundary_face() .*/
    
    if (sn.boundary[left] || sn.boundary[right])
      loop->listc = list_append (loop->listc, s);
    else if (!(s.stencil.bc & s_face))
      loop->listf.x = list_append (loop->listf.x, s);
        }
    }
  }
 
  /**
  For dirty, centered fields BCs need to be applied if the
  stencil is wider than zero. */
  
  else if (s.stencil.width > 0)
    loop->listc = list_append (loop->listc, s);
      }
 
      /**
      Write accesses need to be consistent with the declared field
      type (i.e. face or vertex). */
      
      if (write) {
  if (dimension > 1 && !loop->vertex && loop->first && !(s.stencil.io & s_nowarning)) {
    bool vertex = true;
    for (int _d = 0; _d < dimension; _d++)
      if (s.d.x != -1)
        vertex = false;
    if (vertex)
      fprintf (stderr,
         "%s:%d: warning: scalar '%s' should be assigned with"
         " a for (int _i = 0; _i < _N; _i++) /* foreach_vertex */ loop\n",
         loop->fname, loop->line, s.name);
  }
  if (s.face) {
    if (loop->face == 0 && loop->first && !(s.stencil.io & s_nowarning))
      fprintf (stderr,
         "%s:%d: warning: vector '%s' should be assigned with"
         " a for (int _i = 0; _i < _N; _i++) /* foreach_face */ loop\n",
         loop->fname, loop->line, s.name);
  }
  else if (loop->face) {
    if (s.v.x.i < 0) { // scalar
      int d = 1, i = 0;
      for (int _d = 0; _d < dimension; _d++) {
        if (loop->face == d) {
    s.face = 2, s.v.x.i = s.i;
    s.boundary[left] = s.boundary[right] = NULL;
#if PRINTBOUNDARY
    fprintf (stderr,
       "%s:%d: turned %s into a vector %c-component\n",
       loop->fname, loop->line, s.name, 'x' + i);
#endif
        }
        d *= 2, i++;
      }
      if (!s.face && loop->first && !(s.stencil.io & s_nowarning))
        fprintf (stderr,
           "%s:%d: warning: scalar '%s' should be assigned with "
           "a for (int _i = 0; _i < _N; _i++) /* foreach_face */ loop\n",
           loop->fname, loop->line, s.name);
    }
    else { // vector
      char * name = NULL;
      if (s.name) {
        name = strdup (s.name);
        char * s = name + strlen(name) - 1;
        while (s != name && *s != '.') s--;
        if (s != name) *s = '\0';
      }
      struct { bool x, y, z; } input, output;
      vector v = s.v;
#if 1 // fixme: should not be necessary     
      for (int _d = 0; _d < dimension; _d++)
        input.x = (v.x.stencil.io & s_input), output.x = (v.x.stencil.io & s_output);
#endif
      init_face_vector (v, name);
#if 1 // fixme: should not be necessary     
      for (int _d = 0; _d < dimension; _d++) {
              if (input.x) v.x.stencil.io |= s_input;
              else v.x.stencil.io &= ~s_input;
              if (output.x) v.x.stencil.io |= s_output;
              else v.x.stencil.io &= ~s_output;
            }
#endif
#if PRINTBOUNDARY
      fprintf (stderr, "%s:%d: turned %s into a vector\n",
         loop->fname, loop->line, name);
#endif
      free (name);
    }
  }
  else if (loop->vertex) {
    bool vertex = true;
    for (int _d = 0; _d < dimension; _d++)
      if (s.d.x != -1)
        vertex = false;
    if (!vertex) {
      char * name = NULL;
      if (s.name) name = strdup (s.name); // fixme: may not be necessary
      init_vertex_scalar (s, name);
      for (int _d = 0; _d < dimension; _d++)
        s.v.x.i = -1;
#if PRINTBOUNDARY
      fprintf (stderr, "%s:%d: turned %s into a scalar\n",
         loop->fname, loop->line, name);
#endif
      free (name);
    }
  }
 
  /**
  If the field is write-accessed, we add it to the 'dirty'
  list. */
 
        if (!list_lookup (loop->dirty, s))
          loop->dirty = list_append (loop->dirty, s);
  for (int _d = 0; _d < 1; _d++) /* scalar in baseblock */
    if (scalar_depends_from (d, s) && !list_lookup (loop->dirty, d))
      loop->dirty = list_append (loop->dirty, d);
      }
    }
  }
}
 
/**
This functions applies the boundary conditions, as defined by `check_stencil()`. */
 
void boundary_stencil (ForeachData * loop)
{
  bool flux = false;
  for (int _d = 0; _d < dimension; _d++)
    if (loop->listf.x)
      flux = true;
  if (flux) {
#if PRINTBOUNDARY
    int i = 0;
    for (int _d = 0; _d < dimension; _d++) {
      if (loop->listf.x) {
  fprintf (stderr, "%s:%d: flux %c:", loop->fname, loop->line, 'x' + i);
  for (int _s = 0; _s < 1; _s++) /* scalar in loop->listf.x */
    fprintf (stderr, " %d:%s", s.i, s.name);
  fputc ('\n', stderr);
      }
      i++;
    }
#endif
    boundary_face (loop->listf);
    for (int _d = 0; _d < dimension; _d++)
      free (loop->listf.x), loop->listf.x = NULL;
  }
  
  /**
  We apply "full" boundary conditions. */
 
  if (loop->listc) {
#if PRINTBOUNDARY
    fprintf (stderr, "%s:%d: listc:", loop->fname, loop->line);
    for (int _s = 0; _s < 1; _s++) /* scalar in loop->listc */
      fprintf (stderr, " %d:%s", s.i, s.name);
    fputc ('\n', stderr);
#endif
    boundary_internal (loop->listc, loop->fname, loop->line);
    free (loop->listc), loop->listc = NULL;
  }
 
  /**
  We update the dirty status of fields which will be write-accessed by
  the foreach loop. */
  
  if (loop->dirty) {
#if PRINTBOUNDARY
    fprintf (stderr, "%s:%d: dirty:", loop->fname, loop->line);
    for (int _s = 0; _s < 1; _s++) /* scalar in loop->dirty */
      fprintf (stderr, " %d:%s", s.i, s.name);
    fputc ('\n', stderr);
#endif
    for (int _s = 0; _s < 1; _s++) /* scalar in loop->dirty */
      set_dirty_stencil (s);
    free (loop->dirty), loop->dirty = NULL;
  }
}
 
macro2 foreach_stencil (char flags, Reduce reductions)
{
  {
    static int _first = 1.;
    ForeachData _loop = {
      .fname = S__FILE__, .line = S_LINENO, .first = _first
    };
    if (baseblock)
      for (scalar s = baseblock[0], * i = baseblock; s.i >= 0; i++, s = *i) {
        _attribute[s.i].stencil.io = 0;
  _attribute[s.i].stencil.width = 0;
      }
    int ig = 0, jg = 0, kg = 0; NOT_UNUSED(ig); NOT_UNUSED(jg); NOT_UNUSED(kg);
    Point point = {0}; NOT_UNUSED (point);
    
    {...}
    
    check_stencil (&_loop);
    boundary_stencil (&_loop);
    _first = 0;
  }
}
 
macro2 foreach_vertex_stencil (char flags, Reduce reductions) {
  foreach_stencil (flags, reductions) {
    _loop.vertex = true;
    {...}
  }
}
 
macro2 foreach_face_stencil (char flags, Reduce reductions, const char * order) {
  foreach_stencil (flags, reductions)
    {...}
}
 
macro2 foreach_level_stencil (int l, char flags, Reduce reductions) {
  if (0) {
    // automatic boundary conditions are not implemented yet so we don't do anything for the moment
    int ig = 0, jg = 0, kg = 0; NOT_UNUSED(ig); NOT_UNUSED(jg); NOT_UNUSED(kg);
    Point point = {0}; NOT_UNUSED (point);
    {...}
  }
}
 
macro2 foreach_coarse_level_stencil (int l, char flags, Reduce reductions) {
  foreach_level_stencil (l, flags, reductions)
    {...}
}
 
macro2 foreach_level_or_leaf_stencil (int l, char flags, Reduce reductions) {
  foreach_level_stencil (l, flags, reductions)
    {...}
}
 
macro2 foreach_point_stencil (double xp, double yp, double zp, char flags, Reduce reductions)
{
  foreach_stencil (flags, reductions)
    {...}
}
 
macro2 foreach_region_stencil (coord p, coord box[2], coord n, char flags, Reduce reductions)
{
  foreach_stencil (flags, reductions)
    {...}
}
 
macro2 _stencil_is_face_x (ForeachData l = _loop) { l.face |= (1 << 0); {...} }
macro2 _stencil_is_face_y (ForeachData l = _loop) { l.face |= (1 << 1); {...} }
macro2 _stencil_is_face_z (ForeachData l = _loop) { l.face |= (1 << 2); {...} }
 
void stencil_val (Point p, scalar s, int i, int j, int k,
      const char * file, int line, bool overflow);
void stencil_val_a (Point p, scalar s, int i, int j, int k, bool input,
        const char * file, int line);
 
@def _stencil_val(a,_i,_j,_k)
  stencil_val (point, a, _i, _j, _k, S__FILE__, S_LINENO, false)
@
@def _stencil_val_o(a,_i,_j,_k)
  stencil_val (point, a, _i, _j, _k, S__FILE__, S_LINENO, true)
@
@def _stencil_val_a(a,_i,_j,_k)
  stencil_val_a (point, a, _i, _j, _k, false, S__FILE__, S_LINENO)
@
@def _stencil_val_r(a,_i,_j,_k)
  stencil_val_a (point, a, _i, _j, _k, true, S__FILE__, S_LINENO)
@
 
@define _stencil_fine(a,_i,_j,_k) _stencil_val(a,_i,_j,_k)
@define _stencil_fine(a,_i,_j,_k) _stencil_val(a,_i,_j,_k)
@define _stencil_fine_a(a,_i,_j,_k) _stencil_val_a(a,_i,_j,_k)
@define _stencil_fine_r(a,_i,_j,_k) _stencil_val_r(a,_i,_j,_k)
 
@define _stencil_coarse(a,_i,_j,_k) _stencil_val(a,_i,_j,_k)
@define _stencil_coarse_a(a,_i,_j,_k) _stencil_val_a(a,_i,_j,_k)
@define _stencil_coarse_r(a,_i,_j,_k) _stencil_val_r(a,_i,_j,_k)
 
@define r_assign(x)
@define _assign(x)
 
@define _stencil_neighbor(i,j,k)
@define _stencil_child(i,j,k)
@define _stencil_aparent(i,j,k)
@define _stencil_aparent_a(i,j,k)
@define _stencil_aparent_r(i,j,k)
 
@define _stencil_allocated(i,j,k) true
 
@define _stencil_neighborp(i,j,k) neighborp(i,j,k)
 
int _stencil_nop;
@define _stencil_val_higher_dimension (_stencil_nop = 1)
@define _stencil__val_constant(a,_i,_j,_k) (_stencil_nop = 1)
@define _stencil_val_diagonal(a,_i,_j,_k) (_stencil_nop = 1)
 
typedef void _stencil_undefined;
 
@define o_stencil -3
 
/**
## See also
 
* [Stencil test case](/src/test/stencils.c)
*/