diff options
Diffstat (limited to 'Carpet')
-rw-r--r-- | Carpet/CarpetRegrid2/interface.ccl | 2 | ||||
-rw-r--r-- | Carpet/CarpetRegrid2/param.ccl | 27 | ||||
-rw-r--r-- | Carpet/CarpetRegrid2/src/amr.cc | 210 | ||||
-rw-r--r-- | Carpet/CarpetRegrid2/src/amr.hh | 15 | ||||
-rw-r--r-- | Carpet/CarpetRegrid2/src/make.code.defn | 3 |
5 files changed, 256 insertions, 1 deletions
diff --git a/Carpet/CarpetRegrid2/interface.ccl b/Carpet/CarpetRegrid2/interface.ccl index 8202f4072..ee2bb307d 100644 --- a/Carpet/CarpetRegrid2/interface.ccl +++ b/Carpet/CarpetRegrid2/interface.ccl @@ -139,6 +139,8 @@ CCTK_REAL radiixyz[10] TYPE=array DIM=1 SIZE=30 DISTRIB=constant radius_x radius_y radius_z } "Radii of refined regions for each level" +CCTK_REAL level_mask TYPE=gf TAGS='checkpoint="no" prolongation="none"' "Requested refinement level for this grid point" + PRIVATE: CCTK_INT old_active[10] TYPE=scalar "Old whether this centre is active" diff --git a/Carpet/CarpetRegrid2/param.ccl b/Carpet/CarpetRegrid2/param.ccl index c45241de9..70b1be5f1 100644 --- a/Carpet/CarpetRegrid2/param.ccl +++ b/Carpet/CarpetRegrid2/param.ccl @@ -64,6 +64,33 @@ BOOLEAN symmetry_rotating180 "Ensure a 180 degree rotating symmetry about the z +BOOLEAN adaptive_refinement "Use level_mask for adaptive refinement" STEERABLE=always +{ +} "no" + +CCTK_INT adaptive_block_size "Block size for adaptive refinement" STEERABLE=always +{ + 1:* :: "" +} 8 + +CCTK_INT adaptive_block_size_x "Block size in x direction for adaptive refinement" STEERABLE=always +{ + -1 :: "use adaptive_block_size" + 1:* :: "" +} -1 +CCTK_INT adaptive_block_size_y "Block size in y direction for adaptive refinement" STEERABLE=always +{ + -1 :: "use adaptive_block_size" + 1:* :: "" +} -1 +CCTK_INT adaptive_block_size_z "Block size in z direction for adaptive refinement" STEERABLE=always +{ + -1 :: "use adaptive_block_size" + 1:* :: "" +} -1 + + + CCTK_INT num_centres "Number of refinement centres" { 0:10 :: "" diff --git a/Carpet/CarpetRegrid2/src/amr.cc b/Carpet/CarpetRegrid2/src/amr.cc new file mode 100644 index 000000000..572654f6d --- /dev/null +++ b/Carpet/CarpetRegrid2/src/amr.cc @@ -0,0 +1,210 @@ +#include <cctk.h> +#include <cctk_Arguments.h> +#include <cctk_Parameters.h> + +#include <cmath> + +#include <carpet.hh> +#include <mpi_string.hh> + +#include "boundary.hh" + + + +namespace CarpetRegrid2 { + + using namespace std; + using namespace Carpet; + + void evaluate_level_mask (cGH const* restrict const cctkGH, + vector<ibset>& regions, int const rl) + { + DECLARE_CCTK_PARAMETERS; + + assert (is_singlemap_mode()); + int const m = Carpet::map; + assert (rl > 0); + BEGIN_GLOBAL_MODE (cctkGH) { + ENTER_LEVEL_MODE (cctkGH, rl-1) { + ENTER_SINGLEMAP_MODE (cctkGH, m, CCTK_GF) { + if (verbose or veryverbose) { + cout << "Refinement level " << rl << ":\n"; + } + + gh const& hh = *vhh.AT(Carpet::map); + + ivect const gsh = ivect::ref(cctkGH->cctk_gsh); + ivect const nghostzones = ivect::ref(cctkGH->cctk_nghostzones); + if (veryverbose) { + cout << " gsh: " << gsh << "\n" + << " nghostzones: " << nghostzones << "\n"; + } + + // Determine the block size + ivect block_size = adaptive_block_size; + if (adaptive_block_size_x > 0) block_size[0] = adaptive_block_size_x; + if (adaptive_block_size_y > 0) block_size[1] = adaptive_block_size_y; + if (adaptive_block_size_z > 0) block_size[2] = adaptive_block_size_z; + assert (all (block_size > 0)); + + // Determine the block offset, i.e. the grid point index of + // the coordinate origin. We don't start blocks at index + // zero, so that one can change the location of the outer + // boundary without influencing the refinement mechanism. + // Instead, we align the block such that one of the blocks + // begins at the coordinate origin. + ivect origin; + rvect rorigin; + { + DECLARE_CCTK_ARGUMENTS; + for (int d=0; d<dim; ++d) { + // We round up, because the origin is located between + // grid points for cell-centred grids, and we want the + // next upper grid point in this case. We subtract 0.25 + // just for proper rounding. + origin[d] = + ceil(- CCTK_ORIGIN_SPACE(d) / CCTK_DELTA_SPACE(d) - 0.25); + rorigin[d] = + CCTK_ORIGIN_SPACE(d) + origin[d] * CCTK_DELTA_SPACE(d); + } + } + // block_index := ([0...gsh] + block_offset) / block_size + // [0...gsh] = block_index * block_size - block_offset + ivect const block_offset = + (block_size - origin % block_size) % block_size; + ivect const num_blocks = + (gsh + block_offset + block_size - 1) / block_size; + if (veryverbose) { + cout << " origin: " << origin << "\n" + << " coordinates of origin: " << rorigin << "\n" + << " offset: " << block_offset << "\n" + << " num blocks: " << num_blocks << "\n"; + } + assert (all (num_blocks > 0)); + + // Mask containing all refined block indices (in arbitrary + // order, and potentially even duplicated). This is really a + // set, but we need to communicate this data structure, and + // therefore a vector is a better choice. + vector<ivect> mask; + + // For vertex centred grids, the blocks need to overlap by + // one grid point, so that e.g. an 8^3 block is turned into + // a 9^3 block. This will ensure symmetry of refined regions + // in the simulation domain. + int const overlap = hh.refcent == vertex_centered ? 1 : 0; + + BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) { + DECLARE_CCTK_ARGUMENTS; + + if (veryverbose) { + cout << " component " << component << ":\n"; + } + + ivect const lbnd = ivect::ref(cctk_lbnd); + ivect const lssh (CCTK_LSSH(0,0), CCTK_LSSH(0,1), CCTK_LSSH(0,2)); + ivect const bboxlo (cctk_bbox[0], cctk_bbox[2], cctk_bbox[4]); + ivect const bboxhi (cctk_bbox[1], cctk_bbox[3], cctk_bbox[5]); + + ivect const imin = 0 + either(bboxlo, 0, nghostzones); + ivect const imax = lssh - either(bboxhi, 0, nghostzones); + + ivect const bmin = + max (0, + ((lbnd + imin + block_offset + block_size - overlap) / + block_size) - + 1); + ivect const bmax = (lbnd + imax + block_offset) / block_size; + + // Loop over all blocks + int nblocks = 0; + int nrefined = 0; + for (int bk=bmin[2]; bk<bmax[2]; ++bk) { + for (int bj=bmin[1]; bj<bmax[1]; ++bj) { + for (int bi=bmin[0]; bi<bmax[0]; ++bi) { + ivect const bind (bi,bj,bk); + + ivect const bimin = + max(imin, + (bind ) * block_size - block_offset - lbnd ); + ivect const bimax = + min(imax, + (bind+1) * block_size - block_offset - lbnd + overlap); + + bool refine = false; + + // Loop over all points in this block + for (int k=bimin[2]; k<bimax[2]; ++k) { + for (int j=bimin[1]; j<bimax[1]; ++j) { + for (int i=bimin[0]; i<bimax[0]; ++i) { + int const ind3d = CCTK_GFINDEX3D(cctkGH, i,j,k); + refine = refine or level_mask[ind3d] >= rl; + } + } + } + + // Refine this block if any point in this block requires + // refinement + if (refine) { + if (veryverbose) { + cout << " refining block " << bind << "\n"; + } + mask.push_back(bind); + ++ nrefined; + } + ++ nblocks; + + } + } + } + + if (veryverbose) { + cout << " refining " << nrefined << " " + << "of " << nblocks << " blocks on this component\n"; + } + + } END_LOCAL_COMPONENT_LOOP; + + // Combine this mask from all processes, and to all processes + vector<ivect> const fullmask = allgatherv1 (dist::comm(), mask); + + // Convert vector into ibset + const ibbox& cbaseext = hh.baseextent(mglevel, rl-1); + const ibbox& baseext = hh.baseextent(mglevel, rl); + + ibset& region = regions.at(rl); + for (vector<ivect>::const_iterator + it = fullmask.begin(); it != fullmask.end(); ++it) + { + ivect const& bind = *it; + ivect const cstr = cbaseext.stride(); + ivect const clo = + cbaseext.lower() + (bind * block_size - block_offset) * cstr; + ivect const cup = clo + (block_size + overlap - 1) * cstr; + ibbox const cblock (clo, cup, cstr); + ibbox const block = cblock.expanded_for(baseext); + region |= block; + } + + if (verbose or veryverbose) { + if (veryverbose) { + cout << " refined region is " << region << "\n" + << " domain is " << baseext << "\n" + << " next coarser level is " << regions.at(rl-1) << "\n"; + } + ibbox::size_type const nrefined = region.size(); + ibbox::size_type const npoints = baseext.size(); + ibbox::size_type const ncoarse = + regions.at(rl-1).expanded_for(baseext).size(); + cout << " refining " + << (100.0*nrefined/ncoarse) << "% of the next coarser level\n" + << " refining " + << (100.0*nrefined/npoints) << "% of the domain\n"; + } + + } LEAVE_SINGLEMAP_MODE; + } LEAVE_LEVEL_MODE; + } END_GLOBAL_MODE; + } + +} // namespace CarpetRegrid2 diff --git a/Carpet/CarpetRegrid2/src/amr.hh b/Carpet/CarpetRegrid2/src/amr.hh new file mode 100644 index 000000000..27abb6349 --- /dev/null +++ b/Carpet/CarpetRegrid2/src/amr.hh @@ -0,0 +1,15 @@ +#ifndef AMR_HH +#define AMR_HH + +#include <carpet.hh> + + + +namespace CarpetRegrid2 { + + void evaluate_level_mask (cGH const* restrict cctkGH, + vector<ibset>& regions, int rl); + +} // namespace CarpetRegrid2 + +#endif // #ifndef AMR_HH diff --git a/Carpet/CarpetRegrid2/src/make.code.defn b/Carpet/CarpetRegrid2/src/make.code.defn index e8bd87f0b..db068520d 100644 --- a/Carpet/CarpetRegrid2/src/make.code.defn +++ b/Carpet/CarpetRegrid2/src/make.code.defn @@ -1,7 +1,8 @@ # Main make.code.defn file for thorn CarpetRegrid2 # Source files in this directory -SRCS = boundary.cc \ +SRCS = amr.cc \ + boundary.cc \ indexing.cc \ initialise.cc \ paramcheck.cc \ |