#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() foreach_face(y,x)
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 (scalar s 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 (scalar s 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 (scalar s 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 foreach_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 (scalar s 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 (scalar s 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
foreach_dimension()
static void periodic_boundary_level_x (const Boundary * b, scalar * list, int l)
{
scalar * list1 = NULL;
for (scalar s 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 (scalar s 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 (scalar s 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
foreach_dimension() {
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 foreach_vertex (char flags = 0, Reduce reductions = None)
{
foreach_face_generic (flags, reductions) {
int ig = -1, jg = -1; NOT_UNUSED(ig); NOT_UNUSED(jg);
{...}
}
}
