diff options
Diffstat (limited to 'Carpet/CarpetReduce/src')
-rw-r--r-- | Carpet/CarpetReduce/src/mask_carpet.cc | 104 | ||||
-rw-r--r-- | Carpet/CarpetReduce/src/mask_coords.c | 40 |
2 files changed, 111 insertions, 33 deletions
diff --git a/Carpet/CarpetReduce/src/mask_carpet.cc b/Carpet/CarpetReduce/src/mask_carpet.cc index 81c53e7b7..160471cb4 100644 --- a/Carpet/CarpetReduce/src/mask_carpet.cc +++ b/Carpet/CarpetReduce/src/mask_carpet.cc @@ -16,6 +16,7 @@ +#if 0 #define LOOP_OVER_NEIGHBOURS(dir) \ { \ ivect dir_(-1); \ @@ -33,6 +34,7 @@ } \ } while (not all (dir_ == -1)); \ } +#endif @@ -112,7 +114,15 @@ namespace CarpetMask { ibset test_boxes; ibset test_cfboxes; - LOOP_OVER_NEIGHBOURS (shift) { + for (int neighbour=0; neighbour<ipow(3,dim); ++neighbour) { + ivect shift; + int itmp=neighbour; + for (int d=0; d<dim; ++d) { + shift[d] = itmp % 3 - 1; // [-1 ... +1] + itmp /= 3; + } + assert (itmp==0); + // In this loop, shift [1,1,1] denotes a convex corner of the // region which should be masked out, i.e. a region where only // a small bit (1/8) of the region should be masked out. @@ -122,34 +132,51 @@ namespace CarpetMask { ibset boxes = not_active; ibset fboxes = fine_active; - for (int d=0; d<dim; ++d) { - ivect const dir = ivect::dir(d); - fboxes = fboxes.shift(-dir) & fboxes & fboxes.shift(+dir); - } - for (int d=0; d<dim; ++d) { - // Calculate the boundary in direction d - ivect const dir = ivect::dir(d); - switch (shift[d]) { - case -1: { - // left boundary - boxes = boxes.shift (-dir) - boxes; - fboxes = fboxes.shift(-dir) - fboxes; - break; + switch (hh.refcent) { + case vertex_centered: { + for (int d=0; d<dim; ++d) { + ivect const dir = ivect::dir(d); + fboxes = fboxes.shift(-dir) & fboxes & fboxes.shift(+dir); } - case 0: { - // interior - // do nothing - break; - } - case +1: { - // right boundary - boxes = boxes.shift (+dir) - boxes; - fboxes = fboxes.shift(+dir) - fboxes; - break; + for (int d=0; d<dim; ++d) { + // Calculate the boundary in direction d + ivect const dir = ivect::dir(d); + switch (shift[d]) { + case -1: { + // left boundary + boxes = boxes.shift (-dir) - boxes; + fboxes = fboxes.shift(-dir) - fboxes; + break; + } + case 0: { + // interior + // do nothing + break; + } + case +1: { + // right boundary + boxes = boxes.shift (+dir) - boxes; + fboxes = fboxes.shift(+dir) - fboxes; + break; + } + default: + assert (0); + } } - default: - assert (0); + break; + } + case cell_centered: { + // Assume that all cell boundaries are aligned + if (all(shift == 0)) { + // do nothing + } else { + boxes = ibset(); + fboxes = ibset(); } + break; + } + default: + assert (0); } boxes &= ext; ibset const cfboxes = fboxes.contracted_for(base) & ext; @@ -213,16 +240,29 @@ namespace CarpetMask { } CCTK_ENDLOOP3(CarpetMaskSetup_restriction); } END_LOOP_OVER_BSET; - } END_LOOP_OVER_NEIGHBOURS; + } // for neighbours { - ibset const boxes = not_active.expand(ivect(1), ivect(1)) & ext; + ibset boxes = not_active; ibset fboxes = fine_active; - for (int d=0; d<dim; ++d) { - ivect const dir = ivect::dir(d); - fboxes = fboxes.shift(-dir) & fboxes & fboxes.shift(+dir); + switch (hh.refcent) { + case vertex_centered: { + boxes = boxes.expand(ivect(1), ivect(1)); + for (int d=0; d<dim; ++d) { + ivect const dir = ivect::dir(d); + fboxes = fboxes.shift(-dir) & fboxes & fboxes.shift(+dir); + } + fboxes = fboxes.expand(ivect(1), ivect(1)); + break; + } + case cell_centered: { + // do nothing + break; } - fboxes = fboxes.expand(ivect(1), ivect(1)); + default: + assert (0); + } + boxes &= ext; ibset const cfboxes = fboxes.contracted_for(base) & ext; if (not (test_boxes == boxes ) or not (test_cfboxes == cfboxes)) diff --git a/Carpet/CarpetReduce/src/mask_coords.c b/Carpet/CarpetReduce/src/mask_coords.c index 271e1b957..d861a6d2b 100644 --- a/Carpet/CarpetReduce/src/mask_coords.c +++ b/Carpet/CarpetReduce/src/mask_coords.c @@ -87,6 +87,7 @@ CoordBase_SetupMask (CCTK_ARGUMENTS) /* Loop over all dimensions and faces */ for (int d=0; d<cctk_dim; ++d) { for (int f=0; f<2; ++f) { + /* If this processor has the outer boundary */ if (cctk_bbox[2*d+f]) { @@ -184,7 +185,44 @@ CoordBase_SetupMask (CCTK_ARGUMENTS) } /* if the boundary is not staggered */ - } /* if is outer boundary */ + } else { /* if this is a ghost boundary */ + + /* Calculate the extent of the local part of the domain */ + for (int dd=0; dd<cctk_dim; ++dd) { + imin[dd] = 0; + imax[dd] = cctk_lsh[dd]; + } + + /* Calculate the extent of the boundary */ + for (int dd=0; dd<cctk_dim; ++dd) { + bmin[dd] = imin[dd]; + bmax[dd] = imax[dd]; + } + if (f==0) { + /* lower face */ + bmax[d] = imin[d] + bnd_points[2*d+f]; + } else { + /* upper face */ + bmin[d] = imax[d] - bnd_points[2*d+f]; + } + + /* Loop over the boundary */ + if (verbose) { + CCTK_VInfo (CCTK_THORNSTRING, + "Setting ghost points in direction %d face %d to weight 0 on level %d", d, f, reflevel); + } +#pragma omp parallel + CCTK_LOOP3(CoordBase_SetupMask_ghost, + i,j,k, + bmin[0],bmin[1],bmin[2], bmax[0],bmax[1],bmax[2], + cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) + { + int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k); + iweight[ind] = 0; + } CCTK_ENDLOOP3(CoordBase_SetupMask_ghost); + + } /* if this is a ghost boundary */ + } /* loop over faces */ } /* loop over directions */ |