#define QUADTREE 1 // OBSOLETE: for backward-compatibility
#define TREE 1
#include "multigrid-common.h"
// scalar attributes
attribute {
void (* refine) (Point, scalar);
}
// Tree methods
#define periodic_clamp(a, level) do { \
if (a < GHOSTS) a += 1 << level; \
else if (a >= GHOSTS + (1 << level)) a -= 1 << level; } while(0)
/*
A cache of refined cells is maintained (if not NULL).
*/
int refine_cell (Point point, scalar * list, int flag, Cache * refined)
{
int nr = 0;
#if TWO_ONE
/* refine neighborhood if required */
if (level > 0)
for (int k = 0; k != 2*child.x; k += child.x)
#if dimension > 1
for (int l = 0; l != 2*child.y; l += child.y)
#endif
#if dimension > 2
for (int m = 0; m != 2*child.z; m += child.z)
#endif
if (aparent(k,l,m).pid >= 0 && is_leaf(aparent(k,l,m))) {
Point p = point;
/* fixme: this should be made
independent from the tree implementation */
p.level = point.level - 1;
p.i = (point.i + GHOSTS)/2 + k;
periodic_clamp (p.i, p.level);
#if dimension > 1
p.j = (point.j + GHOSTS)/2 + l;
periodic_clamp (p.j, p.level);
#endif
#if dimension > 2
p.k = (point.k + GHOSTS)/2 + m;
periodic_clamp (p.k, p.level);
#endif
nr += refine_cell (p, list, flag, refined);
aparent(k,l,m).flags |= flag;
}
#endif
/* update neighborhood */
increment_neighbors (point);
int cflag = is_active(cell) ? (active|leaf) : leaf;
foreach_child()
cell.flags |= cflag;
/* initialise scalars */
for (scalar s in list)
if (is_local(cell) || s.face)
s.refine (point, s);
/* refine */
cell.flags &= ~leaf;
@if _MPI
if (is_border(cell)) {
foreach_child() {
bool bord = false;
foreach_neighbor() {
if (!is_local(cell) || (level > 0 && !is_local(aparent(0)))) {
bord = true; break;
}
if (is_refined_check())
foreach_child()
if (!is_local(cell)) {
bord = true; break;
}
if (bord)
break;
}
if (bord)
cell.flags |= border;
}
if (refined)
cache_append (refined, point, cell.flags);
nr++;
}
@endif
return nr;
}
attribute {
void (* coarsen) (Point, scalar);
}
bool coarsen_cell (Point point, scalar * list)
{
#if TWO_ONE
/* check that neighboring cells are not too fine.
check that children are not different boundaries */
int pid = cell.pid;
foreach_child()
if (cell.neighbors || (cell.pid < 0 && cell.pid != pid))
return false; // cannot coarsen
#endif
/* restriction/coarsening */
for (scalar s in list) {
s.restriction (point, s);
if (s.coarsen)
s.coarsen (point, s);
}
/* coarsen */
cell.flags |= leaf;
/* update neighborhood */
decrement_neighbors (point);
@if _MPI
if (!is_local(cell)) {
cell.flags &= ~(active|border);
foreach_neighbor(1)
if (cell.neighbors)
foreach_child()
if (allocated(0) && is_local(cell) && !is_border(cell))
cell.flags |= border;
}
@endif
return true;
}
void coarsen_cell_recursive (Point point, scalar * list)
{
#if TWO_ONE
/* recursively coarsen children cells */
foreach_child()
if (cell.neighbors)
foreach_neighbor(1)
if (is_refined (cell))
coarsen_cell_recursive (point, list);
#endif
assert (coarsen_cell (point, list));
}
void mpi_boundary_refine (scalar *);
void mpi_boundary_coarsen (int, int);
void mpi_boundary_update (scalar *);
static
scalar * list_add_depend (scalar * list, scalar s)
{
if (is_constant(s) || s.restriction == no_restriction)
return list;
for (scalar t in list)
if (t.i == s.i)
return list;
for (scalar d in s.depends)
list = list_add_depend (list, d);
return list_append (list, s);
}
typedef struct {
int nc, nf;
} astats;
trace
astats adapt_wavelet (scalar * slist, // list of scalars
double * max, // tolerance for each scalar
int maxlevel, // maximum level of refinement
int minlevel = 1, // minimum level of refinement
scalar * list = all) // list of fields to update
{
scalar * ilist = list;
if (is_constant(cm)) {
if (list == NULL || list == all)
list = list_copy (all);
boundary (list);
restriction (slist);
}
else {
if (list == NULL || list == all) {
list = list_copy ({cm, fm});
for (scalar s in all)
list = list_add (list, s);
}
boundary (list);
scalar * listr = list_concat (slist, {cm});
restriction (listr);
free (listr);
}
astats st = {0, 0};
scalar * listc = NULL;
for (scalar s in list)
listc = list_add_depend (listc, s);
// refinement
if (minlevel < 1)
minlevel = 1;
tree->refined.n = 0;
static const int refined = 1 << user, too_fine = 1 << (user + 1);
foreach_cell() {
if (is_active(cell)) {
static const int too_coarse = 1 << (user + 2);
if (is_leaf (cell)) {
if (cell.flags & too_coarse) {
cell.flags &= ~too_coarse;
refine_cell (point, listc, refined, &tree->refined);
st.nf++;
}
continue;
}
else { // !is_leaf (cell)
if (cell.flags & refined) {
// cell has already been refined, skip its children
cell.flags &= ~too_coarse;
continue;
}
// check whether the cell or any of its children is local
bool local = is_local(cell);
if (!local)
foreach_child()
if (is_local(cell)) {
local = true; break;
}
if (local) {
int i = 0;
static const int just_fine = 1 << (user + 3);
for (scalar s in slist) {
double emax = max[i++], sc[1 << dimension];
int c = 0;
foreach_child()
sc[c++] = s[];
s.prolongation (point, s);
c = 0;
foreach_child() {
double e = fabs(sc[c] - s[]);
if (e > emax && level < maxlevel) {
cell.flags &= ~too_fine;
cell.flags |= too_coarse;
}
else if ((e <= emax/1.5 || level > maxlevel) &&
!(cell.flags & (too_coarse|just_fine))) {
if (level >= minlevel)
cell.flags |= too_fine;
}
else if (!(cell.flags & too_coarse)) {
cell.flags &= ~too_fine;
cell.flags |= just_fine;
}
s[] = sc[c++];
}
}
foreach_child() {
cell.flags &= ~just_fine;
if (!is_leaf(cell)) {
cell.flags &= ~too_coarse;
if (level >= maxlevel)
cell.flags |= too_fine;
}
else if (!is_active(cell))
cell.flags &= ~too_coarse;
}
}
}
}
else // inactive cell
continue;
}
mpi_boundary_refine (listc);
// coarsening
// the loop below is only necessary to ensure symmetry of 2:1 constraint
for (int l = depth(); l >= 0; l--) {
foreach_cell()
if (!is_boundary(cell)) {
if (level == l) {
if (!is_leaf(cell)) {
if (cell.flags & refined)
// cell was refined previously, unset the flag
cell.flags &= ~(refined|too_fine);
else if (cell.flags & too_fine) {
if (is_local(cell) && coarsen_cell (point, listc))
st.nc++;
cell.flags &= ~too_fine; // do not coarsen parent
}
}
if (cell.flags & too_fine)
cell.flags &= ~too_fine;
else if (level > 0 && (aparent(0).flags & too_fine))
aparent(0).flags &= ~too_fine;
continue;
}
else if (is_leaf(cell))
continue;
}
mpi_boundary_coarsen (l, too_fine);
}
free (listc);
mpi_all_reduce (st.nf, MPI_INT, MPI_SUM);
mpi_all_reduce (st.nc, MPI_INT, MPI_SUM);
if (st.nc || st.nf)
mpi_boundary_update (list);
if (list != ilist)
free (list);
return st;
}
#define refine(cond) do { \
int refined; \
do { \
boundary (all); \
refined = 0; \
tree->refined.n = 0; \
foreach_leaf() \
if (cond) { \
refine_cell (point, all, 0, &tree->refined); \
refined++; \
continue; \
} \
mpi_all_reduce (refined, MPI_INT, MPI_SUM); \
if (refined) { \
mpi_boundary_refine (all); \
mpi_boundary_update (all); \
} \
} while (refined); \
} while(0)
static void refine_level (int depth)
{
int refined;
do {
refined = 0;
tree->refined.n = 0;
foreach_leaf()
if (level < depth) {
refine_cell (point, NULL, 0, &tree->refined);
refined++;
continue;
}
mpi_all_reduce (refined, MPI_INT, MPI_SUM);
if (refined) {
mpi_boundary_refine (NULL);
mpi_boundary_update (NULL);
}
} while (refined);
}
#define unrefine(cond) do { \
static const int too_fine = 1 << user; \
foreach_cell() { \
if (is_leaf(cell)) \
continue; \
if (is_local(cell) && (cond)) \
cell.flags |= too_fine; \
} \
for (int _l = depth(); _l >= 0; _l--) { \
foreach_cell() { \
if (is_leaf(cell)) \
continue; \
if (level == _l) { \
if (is_local(cell) && (cell.flags & too_fine)) { \
coarsen_cell (point, all); \
cell.flags &= ~too_fine; \
} \
continue; \
} \
} \
mpi_boundary_coarsen (_l, too_fine); \
} \
mpi_boundary_update (all); \
} while (0)
trace
static void halo_face (vectorl vl)
{
foreach_dimension()
for (scalar s in vl.x)
s.dirty = 2;
for (int l = depth() - 1; l >= 0; l--)
foreach_halo (prolongation, l)
foreach_dimension()
if (vl.x) {
#if dimension == 1
if (is_refined (neighbor(-1)))
for (scalar s in vl.x)
s[] = fine(s,0);
if (is_refined (neighbor(1)))
for (scalar s in vl.x)
s[1] = fine(s,2);
#elif dimension == 2
if (is_refined (neighbor(-1)))
for (scalar s in vl.x)
s[] = (fine(s,0,0) + fine(s,0,1))/2.;
if (is_refined (neighbor(1)))
for (scalar s in vl.x)
s[1] = (fine(s,2,0) + fine(s,2,1))/2.;
#else // dimension == 3
if (is_refined (neighbor(-1)))
for (scalar s in vl.x)
s[] = (fine(s,0,0,0) + fine(s,0,1,0) +
fine(s,0,0,1) + fine(s,0,1,1))/4.;
if (is_refined (neighbor(1)))
for (scalar s in vl.x)
s[1] = (fine(s,2,0,0) + fine(s,2,1,0) +
fine(s,2,0,1) + fine(s,2,1,1))/4.;
#endif
}
}
// Cartesian methods
static scalar tree_init_scalar (scalar s, const char * name)
{
s = multigrid_init_scalar (s, name);
s.refine = s.prolongation;
return s;
}
static void prolongation_vertex (Point point, scalar s)
{
#if dimension == 2
fine(s,1,1) = (s[] + s[1] + s[0,1] + s[1,1])/4.;
#else // dimension == 3
fine(s,1,1,1) = (s[] + s[1] + s[0,1] + s[1,1] +
s[0,0,1] + s[1,0,1] + s[0,1,1] + s[1,1,1])/8.;
#endif
for (int i = 0; i <= 1; i++) {
for (int j = 0; j <= 1; j++)
#if dimension == 3
for (int k = 0; k <= 1; k++)
if (allocated_child(2*i,2*j,2*k))
fine(s,2*i,2*j,2*k) = s[i,j,k];
#else // dimension != 3
if (allocated_child(2*i,2*j))
fine(s,2*i,2*j) = s[i,j];
#endif // dimension != 3
foreach_dimension()
if (neighbor(i).neighbors) {
#if dimension == 2
fine(s,2*i,1) = (s[i,0] + s[i,1])/2.;
#elif dimension == 3
fine(s,2*i,1,1) = (s[i,0,0] + s[i,1,0] + s[i,0,1] + s[i,1,1])/4.;
fine(s,2*i,1,0) = (s[i,0,0] + s[i,1,0])/2.;
fine(s,2*i,0,1) = (s[i,0,0] + s[i,0,1])/2.;
if (allocated_child(2*i,1,2))
fine(s,2*i,1,2) = (s[i,0,1] + s[i,1,1])/2.;
if (allocated_child(2*i,2,1))
fine(s,2*i,2,1) = (s[i,1,0] + s[i,1,1])/2.;
#endif // dimension == 3
}
}
}
static scalar tree_init_vertex_scalar (scalar s, const char * name)
{
s = multigrid_init_vertex_scalar (s, name);
s.refine = s.prolongation = prolongation_vertex;
return s;
}
static void tree_setup_vector (vector v)
{
foreach_dimension()
v.x.refine = v.x.prolongation;
}
static vector tree_init_vector (vector v, const char * name)
{
v = multigrid_init_vector (v, name);
tree_setup_vector (v);
return v;
}
foreach_dimension()
static void refine_face_x (Point point, scalar s)
{
vector v = s.v;
#if dimension <= 2
if (!is_refined(neighbor(-1)) &&
(is_local(cell) || is_local(neighbor(-1)))) {
double g1 = (v.x[0,+1] - v.x[0,-1])/8.;
for (int j = 0; j <= 1; j++)
fine(v.x,0,j) = v.x[] + (2*j - 1)*g1;
}
if (!is_refined(neighbor(1)) && neighbor(1).neighbors &&
(is_local(cell) || is_local(neighbor(1)))) {
double g1 = (v.x[1,+1] - v.x[1,-1])/8.;
for (int j = 0; j <= 1; j++)
fine(v.x,2,j) = v.x[1] + (2*j - 1)*g1;
}
if (is_local(cell)) {
double g1 = (v.x[0,+1] - v.x[0,-1] + v.x[1,+1] - v.x[1,-1])/16.;
for (int j = 0; j <= 1; j++)
fine(v.x,1,j) = (v.x[] + v.x[1])/2. + (2*j - 1)*g1;
}
#else // dimension > 2
if (!is_refined(neighbor(-1)) &&
(is_local(cell) || is_local(neighbor(-1)))) {
double g1 = (v.x[0,+1] - v.x[0,-1])/8.;
double g2 = (v.x[0,0,+1] - v.x[0,0,-1])/8.;
for (int j = 0; j <= 1; j++)
for (int k = 0; k <= 1; k++)
fine(v.x,0,j,k) = v.x[] + (2*j - 1)*g1 + (2*k - 1)*g2;
}
if (!is_refined(neighbor(1)) && neighbor(1).neighbors &&
(is_local(cell) || is_local(neighbor(1)))) {
double g1 = (v.x[1,+1] - v.x[1,-1])/8.;
double g2 = (v.x[1,0,+1] - v.x[1,0,-1])/8.;
for (int j = 0; j <= 1; j++)
for (int k = 0; k <= 1; k++)
fine(v.x,2,j,k) = v.x[1] + (2*j - 1)*g1 + (2*k - 1)*g2;
}
if (is_local(cell)) {
double g1 = (v.x[0,+1] - v.x[0,-1] + v.x[1,+1] - v.x[1,-1])/16.;
double g2 = (v.x[0,0,+1] - v.x[0,0,-1] + v.x[1,0,+1] - v.x[1,0,-1])/16.;
for (int j = 0; j <= 1; j++)
for (int k = 0; k <= 1; k++)
fine(v.x,1,j,k) = (v.x[] + v.x[1])/2. + (2*j - 1)*g1 + (2*k - 1)*g2;
}
#endif // dimension > 2
}
void refine_face (Point point, scalar s)
{
vector v = s.v;
foreach_dimension()
v.x.prolongation (point, v.x);
}
void refine_face_solenoidal (Point point, scalar s)
{
refine_face (point, s);
#if dimension > 1
if (is_local(cell)) {
// local projection, see section 3.3 of Popinet, JCP, 2009
vector v = s.v;
double d[1 << dimension], p[1 << dimension];
int i = 0;
foreach_child() {
d[i] = 0.;
foreach_dimension()
d[i] += v.x[1] - v.x[];
i++;
}
#if dimension == 2
p[0] = 0.;
p[1] = (3.*d[3] + d[0])/4. + d[2]/2.;
p[2] = (d[3] + 3.*d[0])/4. + d[2]/2.;
p[3] = (d[3] + d[0])/2. + d[2];
fine(v.x,1,1) += p[1] - p[0];
fine(v.x,1,0) += p[3] - p[2];
fine(v.y,0,1) += p[0] - p[2];
fine(v.y,1,1) += p[1] - p[3];
#elif dimension == 3
static double m[7][7] = {
{7./12,5./24,3./8,5./24,3./8,1./4,1./3},
{5./24,7./12,3./8,5./24,1./4,3./8,1./3},
{3./8,3./8,3./4,1./4,3./8,3./8,1./2},
{5./24,5./24,1./4,7./12,3./8,3./8,1./3},
{3./8,1./4,3./8,3./8,3./4,3./8,1./2},
{1./4,3./8,3./8,3./8,3./8,3./4,1./2},
{1./3,1./3,1./2,1./3,1./2,1./2,5./6}
};
p[0] = 0.;
for (int i = 0; i < 7; i++) {
p[i + 1] = 0.;
for (int j = 0; j < 7; j++)
p[i + 1] += m[i][j]*d[j+1];
}
for (int k = 0; k <= 1; k++) {
fine(v.x,1,0,k) += p[4+k] - p[0+k];
fine(v.x,1,1,k) += p[6+k] - p[2+k];
fine(v.y,0,1,k) += p[2+k] - p[0+k];
fine(v.y,1,1,k) += p[6+k] - p[4+k];
}
fine(v.z,0,0,1) += p[1] - p[0];
fine(v.z,0,1,1) += p[3] - p[2];
fine(v.z,1,0,1) += p[5] - p[4];
fine(v.z,1,1,1) += p[7] - p[6];
#endif // dimension == 3
}
#endif // dimension > 1
}
static vector tree_init_face_vector (vector v, const char * name)
{
v = cartesian_init_face_vector (v, name);
foreach_dimension()
v.x.restriction = v.x.refine = no_restriction;
v.x.restriction = restriction_face;
v.x.refine = refine_face;
foreach_dimension()
v.x.prolongation = refine_face_x;
return v;
}
static tensor tree_init_tensor (tensor t, const char * name)
{
t = multigrid_init_tensor (t, name);
foreach_dimension()
tree_setup_vector (t.x);
return t;
}
trace
static void tree_boundary_level (scalar * list, int l)
{
int depth = l < 0 ? depth() : l;
if (tree_is_full()) {
boundary_iterate (level, list, depth);
return;
}
scalar * listdef = NULL, * listc = NULL, * list2 = NULL, * vlist = NULL;
for (scalar s in list)
if (!is_constant (s)) {
if (s.restriction == restriction_average) {
listdef = list_add (listdef, s);
list2 = list_add (list2, s);
}
else if (s.restriction != no_restriction) {
listc = list_add (listc, s);
if (s.face)
foreach_dimension()
list2 = list_add (list2, s.v.x);
else {
list2 = list_add (list2, s);
if (s.restriction == restriction_vertex)
vlist = list_add (vlist, s);
}
}
}
if (vlist) // vertex scalars
#if dimension == 1
foreach_vertex (noauto)
if (is_refined(cell) || is_refined(neighbor(-1)))
for (scalar s in vlist)
s[] = is_vertex (child(0)) ? fine(s) : nodata;
#elif dimension == 2
foreach_vertex (noauto) {
if (is_refined(cell) || is_refined(neighbor(-1)) ||
is_refined(neighbor(0,-1)) || is_refined(neighbor(-1,-1))) {
// corner
for (scalar s in vlist)
s[] = is_vertex (child(0)) ? fine(s) : nodata;
}
else
foreach_dimension()
if (child.y == 1 &&
(is_prolongation(cell) || is_prolongation(neighbor(-1)))) {
// center of refined edge
for (scalar s in vlist)
s[] = is_vertex(neighbor(0,-1)) && is_vertex(neighbor(0,1)) ?
(s[0,-1] + s[0,1])/2. : nodata;
}
}
#else // dimension == 3
foreach_vertex (noauto) {
if (is_refined(cell) || is_refined(neighbor(-1)) ||
is_refined(neighbor(0,-1)) || is_refined(neighbor(-1,-1)) ||
is_refined(neighbor(0,0,-1)) || is_refined(neighbor(-1,0,-1)) ||
is_refined(neighbor(0,-1,-1)) || is_refined(neighbor(-1,-1,-1))) {
// corner
for (scalar s in vlist)
s[] = is_vertex (child(0)) ? fine(s) : nodata;
}
else
foreach_dimension() {
if (child.y == 1 && child.z == 1 &&
(is_prolongation(cell) || is_prolongation(neighbor(-1)))) {
// center of refined face
for (scalar s in vlist)
s[] = is_vertex(neighbor(0,-1,-1)) && is_vertex(neighbor(0,1,-1))
&& is_vertex(neighbor(0,-1,1)) && is_vertex(neighbor(0,1,1)) ?
(s[0,-1,-1] + s[0,1,-1] + s[0,-1,1] + s[0,1,1])/4. : nodata;
}
else if (child.x == -1 && child.z == -1 && child.y == 1 &&
(is_prolongation(cell) || is_prolongation(neighbor(-1)) ||
is_prolongation(neighbor(0,0,-1)) ||
is_prolongation(neighbor(-1,0,-1)))) {
// center of refined edge
for (scalar s in vlist)
s[] = is_vertex(neighbor(0,-1)) && is_vertex(neighbor(0,1)) ?
(s[0,-1] + s[0,1])/2. : nodata;
}
}
}
#endif // dimension == 3
free (vlist);
if (listdef || listc) {
boundary_iterate (restriction, list2, depth);
for (int l = depth - 1; l >= 0; l--) {
foreach_coarse_level(l) {
for (scalar s in listdef)
restriction_average (point, s);
for (scalar s in listc)
s.restriction (point, s);
}
boundary_iterate (restriction, list2, l);
}
free (listdef);
free (listc);
free (list2);
}
scalar * listr = NULL;
vector * listf = NULL;
for (scalar s in list)
if (!is_constant (s) && s.refine != no_restriction) {
if (s.face)
listf = vectors_add (listf, s.v);
else
listr = list_add (listr, s);
}
if (listr || listf) {
boundary_iterate (level, list, 0);
for (int i = 0; i < depth; i++) {
foreach_halo (prolongation, i) {
for (scalar s in listr)
s.prolongation (point, s);
for (vector v in listf)
foreach_dimension()
v.x.prolongation (point, v.x);
}
boundary_iterate (level, list, i + 1);
}
free (listr);
free (listf);
}
}
double treex (Point point) {
if (level == 0)
return 0;
#if dimension == 2
double i = 2*child.x - child.y;
if (i <= 1 && i >= -1) i = -i;
#else
assert (false); // not implemented
double i = 0;
#endif
return treex(parent) + i/(1 << 2*(level - 1));
}
double treey (Point point) {
if (level == 0)
return 0;
return treey(parent) + 4./(1 << 2*(level - 1));
}
void output_tree (FILE * fp)
{
foreach_cell()
if (cell.neighbors)
foreach_child()
if (is_local(cell))
fprintf (fp, "%g %g\n%g %g\n\n",
treex(parent), treey(parent), treex(point), treey(point));
}
trace
void tree_check()
{
// checks the consistency of the tree
long nleaves = 0, nactive = 0;
foreach_cell_all() {
if (is_leaf(cell)) {
assert (cell.pid >= 0); // boundaries cannot be leaves
nleaves++;
}
if (is_local(cell))
assert (is_active(cell) || is_prolongation(cell));
if (is_active(cell))
nactive++;
// check number of refined neighbors
int neighbors = 0;
foreach_neighbor(1)
if (allocated(0) && is_refined(cell))
neighbors++;
assert (cell.neighbors == neighbors);
// checks that prolongation cells do not have children
if (!cell.neighbors)
assert (!allocated_child(0));
}
// checks that all active cells are reachable
long reachable = 0;
foreach_cell() {
if (is_active(cell))
reachable++;
else
continue;
}
assert (nactive == reachable);
// checks that all leaf cells are reachable
reachable = 0;
foreach_cell()
if (is_leaf(cell)) {
reachable++;
continue;
}
assert (nleaves == reachable);
}
trace
static void tree_restriction (scalar * list) {
boundary (list);
if (tree_is_full())
multigrid_restriction (list);
}
void tree_methods()
{
multigrid_methods();
init_scalar = tree_init_scalar;
init_vertex_scalar = tree_init_vertex_scalar;
init_vector = tree_init_vector;
init_face_vector = tree_init_face_vector;
init_tensor = tree_init_tensor;
boundary_level = tree_boundary_level;
boundary_face = halo_face;
restriction = tree_restriction;
}
