aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Carpet/Carpet/src/Recompose.cc14
-rw-r--r--Carpet/Carpet/src/SetupGH.cc4
-rw-r--r--Carpet/CarpetIOASCII/src/ioascii.cc111
-rw-r--r--Carpet/CarpetIOASCII/src/ioascii.hh6
-rw-r--r--Carpet/CarpetLib/src/gh.cc53
-rw-r--r--Carpet/CarpetRegrid2/interface.ccl36
-rw-r--r--Carpet/CarpetRegrid2/src/regrid.cc159
7 files changed, 266 insertions, 117 deletions
diff --git a/Carpet/Carpet/src/Recompose.cc b/Carpet/Carpet/src/Recompose.cc
index 3027290de..4da7afb73 100644
--- a/Carpet/Carpet/src/Recompose.cc
+++ b/Carpet/Carpet/src/Recompose.cc
@@ -91,13 +91,17 @@ namespace Carpet {
}
assert ((int)regsss.at(0).size() <= maxreflevels);
for (int ml=0; ml<(int)regsss.size(); ++ml) {
+ int num_regions = 0;
for (int rl=0; rl<(int)regsss.at(0).size(); ++rl) {
// No empty levels
- assert (regsss.at(ml).at(rl).size() > 0);
+ // (but allow some empty maps)
+ // assert (regsss.at(ml).at(rl).size() > 0);
+ num_regions += regsss.at(ml).at(rl).size();
for (int c=0; c<(int)regsss.at(ml).at(rl).size(); ++c) {
// Check sizes
- // Do allow processors with zero grid points
-// assert (all(regsss.at(rl).at(c).at(ml).extent.lower() <= regsss.at(rl).at(c).at(ml).extent.upper()));
+ // (but allow processors with zero grid points)
+ // assert (all(regsss.at(rl).at(c).at(ml).extent.lower() <=
+ // regsss.at(rl).at(c).at(ml).extent.upper()));
// Check strides
const ivect str
= (maxspacereflevelfact / spacereffacts.at(rl) * ipow(mgfact, ml));
@@ -107,6 +111,8 @@ namespace Carpet {
assert (all(regsss.at(ml).at(rl).at(c).extent.upper() % str == 0));
}
}
+ // No empty levels
+ assert (num_regions > 0);
}
}
@@ -782,7 +788,7 @@ namespace Carpet {
}
for (int c=0; c<hh->components(rl); ++c) {
++ num_comps;
- dh::dboxes const & b = dd->boxes.AT(m).AT(rl).AT(c);
+ dh::dboxes const & b = dd->boxes.AT(ml).AT(rl).AT(c);
num_active_mem_points += num_gfs * b.active.size();
num_owned_mem_points += num_gfs * b.owned.size();
num_total_mem_points += num_gfs * b.exterior.size();
diff --git a/Carpet/Carpet/src/SetupGH.cc b/Carpet/Carpet/src/SetupGH.cc
index 109c3bcab..5e34fcd30 100644
--- a/Carpet/Carpet/src/SetupGH.cc
+++ b/Carpet/Carpet/src/SetupGH.cc
@@ -1221,6 +1221,8 @@ namespace Carpet {
} else if (domain_from_coordbase) {
+ assert (not CCTK_IsFunctionAliased ("MultiPatch_GetDomainSpecification"));
+
// Ensure that CartGrid3D::type = "coordbase"
ensure_CartGrid3D_type ();
@@ -1235,6 +1237,8 @@ namespace Carpet {
} else {
// Legacy code
+ assert (not CCTK_IsFunctionAliased ("MultiPatch_GetDomainSpecification"));
+
if (max_refinement_levels > 1) {
// Ensure that CartGrid3D::avoid_origin = no
ensure_CartGrid3D_avoid_origin ();
diff --git a/Carpet/CarpetIOASCII/src/ioascii.cc b/Carpet/CarpetIOASCII/src/ioascii.cc
index ade2299dd..f90770cde 100644
--- a/Carpet/CarpetIOASCII/src/ioascii.cc
+++ b/Carpet/CarpetIOASCII/src/ioascii.cc
@@ -358,7 +358,6 @@ namespace CarpetIOASCII {
::OutputVarAs (const cGH* const cctkGH,
const char* const varname, const char* const alias)
{
- DECLARE_CCTK_ARGUMENTS;
DECLARE_CCTK_PARAMETERS;
int vindex = -1;
@@ -434,7 +433,7 @@ namespace CarpetIOASCII {
for (int ml = 0; ml < mglevels; ++ml) {
last_outputs[ml].resize (maxreflevels);
for (int rl = 0; rl < maxreflevels; ++rl) {
- last_outputs[ml][rl].resize (numelems, cctk_iteration - 1);
+ last_outputs[ml][rl].resize (numelems, cctkGH->cctk_iteration - 1);
}
}
// TODO: this makes a copy of last_outputs, which is expensive
@@ -454,7 +453,7 @@ namespace CarpetIOASCII {
elem = vindex;
}
int& last_output = thisfile->second.at(mglevel).at(reflevel).at(elem);
- if (last_output == cctk_iteration) {
+ if (last_output == cctkGH->cctk_iteration) {
// Has already been output during this iteration
char * const fullname = CCTK_FullName (vindex);
CCTK_VWarn (5, __LINE__, __FILE__, CCTK_THORNSTRING,
@@ -465,8 +464,8 @@ namespace CarpetIOASCII {
free (fullname);
return 0;
}
- assert (last_output < cctk_iteration);
- last_output = cctk_iteration;
+ assert (last_output < cctkGH->cctk_iteration);
+ last_output = cctkGH->cctk_iteration;
// Check for storage
if (not CCTK_QueryGroupStorageI (cctkGH, group)) {
@@ -547,7 +546,9 @@ namespace CarpetIOASCII {
if (desired) {
// Traverse all maps on this refinement and multigrid level
- BEGIN_MAP_LOOP(cctkGH, grouptype) {
+ int const m_min = 0;
+ int const m_max = grouptype == CCTK_GF ? Carpet::maps : 1;
+ for (int m = m_min; m < m_max; ++ m) {
// Find the output offset.
ivect offset(0);
@@ -555,17 +556,17 @@ namespace CarpetIOASCII {
switch (outdim) {
case 0:
offset[0] = GetGridOffset
- (cctkGH, 1,
+ (cctkGH, m, 1,
"out%dD_point_xi", /*"out_point_xi"*/ NULL,
"out%dD_point_x", /*"out_point_x"*/ NULL,
/*out_point_x*/ 0.0);
offset[1] = GetGridOffset
- (cctkGH, 2,
+ (cctkGH, m, 2,
"out%dD_point_yi", /*"out_point_yi"*/ NULL,
"out%dD_point_y", /*"out_point_y"*/ NULL,
/*out_point_y*/ 0.0);
offset[2] = GetGridOffset
- (cctkGH, 3,
+ (cctkGH, m, 3,
"out%dD_point_zi", /*"out_point_zi"*/ NULL,
"out%dD_point_z", /*"out_point_z"*/ NULL,
/*out_point_z*/ 0.0);
@@ -573,31 +574,31 @@ namespace CarpetIOASCII {
case 1:
switch (dirs[0]) {
case 0:
- offset[1] = GetGridOffset (cctkGH, 2,
+ offset[1] = GetGridOffset (cctkGH, m, 2,
"out%dD_xline_yi", "out_xline_yi",
"out%dD_xline_y", "out_xline_y",
out_xline_y);
- offset[2] = GetGridOffset (cctkGH, 3,
+ offset[2] = GetGridOffset (cctkGH, m, 3,
"out%dD_xline_zi", "out_xline_zi",
"out%dD_xline_z", "out_xline_z",
out_xline_z);
break;
case 1:
- offset[0] = GetGridOffset (cctkGH, 1,
+ offset[0] = GetGridOffset (cctkGH, m, 1,
"out%dD_yline_xi", "out_yline_xi",
"out%dD_yline_x", "out_yline_x",
out_yline_x);
- offset[2] = GetGridOffset (cctkGH, 3,
+ offset[2] = GetGridOffset (cctkGH, m, 3,
"out%dD_yline_zi", "out_yline_zi",
"out%dD_yline_z", "out_yline_z",
out_yline_z);
break;
case 2:
- offset[0] = GetGridOffset (cctkGH, 1,
+ offset[0] = GetGridOffset (cctkGH, m, 1,
"out%dD_zline_xi", "out_zline_xi",
"out%dD_zline_x", "out_zline_x",
out_zline_x);
- offset[1] = GetGridOffset (cctkGH, 2,
+ offset[1] = GetGridOffset (cctkGH, m, 2,
"out%dD_zline_yi", "out_zline_yi",
"out%dD_zline_y", "out_zline_y",
out_zline_y);
@@ -615,19 +616,19 @@ namespace CarpetIOASCII {
case 2:
if (dirs[0]==0 and dirs[1]==1) {
offset[2] = GetGridOffset
- (cctkGH, 3,
+ (cctkGH, m, 3,
"out%dD_xyplane_zi", "out_xyplane_zi",
"out%dD_xyplane_z", "out_xyplane_z",
out_xyplane_z);
} else if (dirs[0]==0 and dirs[1]==2) {
offset[1] = GetGridOffset
- (cctkGH, 2,
+ (cctkGH, m, 2,
"out%dD_xzplane_yi", "out_xzplane_yi",
"out%dD_xzplane_y", "out_xzplane_y",
out_xzplane_y);
} else if (dirs[0]==1 and dirs[1]==2) {
offset[0] = GetGridOffset
- (cctkGH, 1,
+ (cctkGH, m, 1,
"out%dD_yzplane_xi", "out_yzplane_xi",
"out%dD_yzplane_x", "out_yzplane_x",
out_yzplane_x);
@@ -653,7 +654,7 @@ namespace CarpetIOASCII {
ostringstream filenamebuf;
filenamebuf << my_out_dir << "/" << alias;
if (maps > 1 and grouptype == CCTK_GF) {
- filenamebuf << "." << Carpet::map;
+ filenamebuf << "." << m;
}
filenamebuf << ".";
if (new_filename_scheme) {
@@ -766,7 +767,12 @@ namespace CarpetIOASCII {
// Traverse and components on this multigrid and
// refinement level and map
- BEGIN_COMPONENT_LOOP(cctkGH, grouptype) {
+ int const c_min = 0;
+ int const c_max =
+ grouptype == CCTK_GF ?
+ vhh.at(m)->components(reflevel) :
+ CCTK_nProcs(cctkGH);
+ for (int c = c_min; c < c_max; ++ c) {
cGroup groupdata;
{
@@ -781,18 +787,14 @@ namespace CarpetIOASCII {
assert (not ierr);
}
- if (groupdata.disttype != CCTK_DISTRIB_CONSTANT
- or component == 0)
- {
+ if (groupdata.disttype != CCTK_DISTRIB_CONSTANT or c == 0) {
- const ggf* const ff
- = arrdata.at(group).at(Carpet::map).data.at(var);
+ const ggf* const ff = arrdata.at(group).at(m).data.at(var);
const int maxtl = output_all_timelevels ? num_tl : 1;
for (int tl=0; tl<maxtl; ++tl) {
- const gdata* const data
- = (*ff) (tl, my_reflevel, component, mglevel);
+ const gdata* const data = (*ff) (tl, my_reflevel, c, mglevel);
ibbox ext = data->extent();
ivect lo = ext.lower();
@@ -855,13 +857,17 @@ namespace CarpetIOASCII {
rvect global_lower;
rvect coord_delta;
if (grouptype == CCTK_GF) {
+ rvect const cctk_origin_space =
+ origin_space.at(m).at(mglevel);
+ rvect const cctk_delta_space =
+ delta_space.at(m) * rvect (mglevelfact);
for (int d=0; d<dim; ++d) {
// lower boundary of Carpet's integer indexing
- global_lower[d] = cctkGH->cctk_origin_space[d];
+ global_lower[d] = cctk_origin_space[d];
// grid spacing of Carpet's integer indexing
coord_delta[d] =
- cctkGH->cctk_delta_space[d] /
- vhh.at(Carpet::map)->baseextents.at(0).at(0).stride()[d];
+ cctk_delta_space[d] /
+ vhh.at(m)->baseextents.at(0).at(0).stride()[d];
}
} else {
for (int d=0; d<dim; ++d) {
@@ -877,7 +883,7 @@ namespace CarpetIOASCII {
ivect offset1;
if (grouptype == CCTK_GF) {
const ibbox& baseext
- = vhh.at(Carpet::map)->baseextents.at(mglevel).at(reflevel);
+ = vhh.at(m)->baseextents.at(mglevel).at(reflevel);
offset1 = baseext.lower() + offset * ext.stride();
} else {
offset1 = offset * ext.stride();
@@ -893,9 +899,8 @@ namespace CarpetIOASCII {
int const numvars = CCTK_NumVarsInGroupI(group);
datas.resize (numvars);
for (int n=0; n<numvars; ++n) {
- const ggf* const ff1
- = arrdata.at(group).at(Carpet::map).data.at(n);
- datas.at(n) = (*ff1) (tl, my_reflevel, component, mglevel);
+ const ggf* const ff1 = arrdata.at(group).at(m).data.at(n);
+ datas.at(n) = (*ff1) (tl, my_reflevel, c, mglevel);
}
} else {
datas.resize (1);
@@ -903,7 +908,7 @@ namespace CarpetIOASCII {
}
WriteASCII (file, datas, ext, vindex, cctkGH->cctk_iteration,
offset1, dirs,
- my_reflevel, mglevel, Carpet::map, component, tl,
+ my_reflevel, mglevel, m, c, tl,
coord_time, coord_lower, coord_upper);
// Append EOL after every component
@@ -917,9 +922,9 @@ namespace CarpetIOASCII {
} // for tl
- } // if distrib!=CONSTANT or component==0
+ } // if distrib!=CONSTANT or c==0
- } END_COMPONENT_LOOP;
+ } // for c
// Append EOL after every complete set of components
if (dist::rank()==ioproc) {
@@ -937,7 +942,7 @@ namespace CarpetIOASCII {
CCTK_REAL const io_bytes = io_bytes_end - io_bytes_begin;
EndTimingIO (cctkGH, io_files, io_bytes, false);
- } END_MAP_LOOP;
+ } // for m
} // if (desired)
@@ -993,7 +998,7 @@ namespace CarpetIOASCII {
template<int outdim>
int IOASCII<outdim>
- ::GetGridOffset (const cGH* const cctkGH, const int dir,
+ ::GetGridOffset (const cGH* const cctkGH, const int m, const int dir,
const char* const itempl, const char* const iglobal,
const char* const ctempl, const char* const cglobal,
const CCTK_REAL cfallback)
@@ -1011,7 +1016,7 @@ namespace CarpetIOASCII {
assert (pcoord);
const CCTK_REAL coord = *pcoord;
assert (ptype == PARAMETER_REAL);
- return CoordToOffset (cctkGH, dir, coord, 0);
+ return CoordToOffset (cctkGH, m, dir, coord, 0);
}
// Second choice: explicit index
@@ -1042,7 +1047,7 @@ namespace CarpetIOASCII {
assert (pcoord);
const CCTK_REAL coord = *pcoord;
assert (ptype == PARAMETER_REAL);
- return CoordToOffset (cctkGH, dir, coord, 0);
+ return CoordToOffset (cctkGH, m, dir, coord, 0);
}
}
@@ -1062,22 +1067,30 @@ namespace CarpetIOASCII {
}
// Fifth choice: default coordinate
- return CoordToOffset (cctkGH, dir, cfallback, 0);
+ return CoordToOffset (cctkGH, m, dir, cfallback, 0);
}
template<int outdim>
int IOASCII<outdim>
- ::CoordToOffset (const cGH* cctkGH, const int dir, const CCTK_REAL coord,
- const int ifallback)
+ ::CoordToOffset (const cGH* cctkGH, const int m, const int dir,
+ const CCTK_REAL coord, const int ifallback)
{
- assert (dir>=1 and dir<=dim);
+ assert (m>=0 and m<Carpet::maps and dir>=1 and dir<=dim);
+
+ assert (mglevel!=-1 and reflevel!=-1 and Carpet::map==-1);
- assert (mglevel!=-1 and reflevel!=-1 and Carpet::map!=-1);
+ rvect const cctk_origin_space = origin_space.at(m).at(mglevel);
+ rvect const cctk_delta_space = delta_space.at(m) * rvect (mglevelfact);
+ ivect const cctk_levfac = spacereffacts.at (reflevel);
+ ibbox const & coarseext = vhh.at(m)->baseextents.at(mglevel).at(0 );
+ ibbox const & baseext = vhh.at(m)->baseextents.at(mglevel).at(reflevel);
+ ivect const cctk_levoff = baseext.lower() - coarseext.lower();
+ ivect const cctk_levoffdenom = baseext.stride();
- const CCTK_REAL delta = cctkGH->cctk_delta_space[dir-1] / cctkGH->cctk_levfac[dir-1];
- const CCTK_REAL lower = cctkGH->cctk_origin_space[dir-1] + cctkGH->cctk_delta_space[dir-1] / cctkGH->cctk_levfac[dir-1] * cctkGH->cctk_levoff[dir-1] / cctkGH->cctk_levoffdenom[dir-1];
+ const CCTK_REAL delta = cctk_delta_space[dir-1] / cctk_levfac[dir-1];
+ const CCTK_REAL lower = cctk_origin_space[dir-1] + cctk_delta_space[dir-1] / cctk_levfac[dir-1] * cctk_levoff[dir-1] / cctk_levoffdenom[dir-1];
const CCTK_REAL rindex = (coord - lower) / delta;
int cindex = (int)floor(rindex + 0.75);
@@ -1317,7 +1330,7 @@ namespace CarpetIOASCII {
const vect<int,3> str = gfext.stride();
const bbox<int,3> ext(lo,up,str);
- gh const & hh = *vhh.at(Carpet::map);
+ gh const & hh = *vhh.at(m);
ibbox const & base = hh.baseextents.at(mglevel).at(reflevel);
assert (base.stride()[0] == base.stride()[1]
diff --git a/Carpet/CarpetIOASCII/src/ioascii.hh b/Carpet/CarpetIOASCII/src/ioascii.hh
index 921490af9..6f0bb74dc 100644
--- a/Carpet/CarpetIOASCII/src/ioascii.hh
+++ b/Carpet/CarpetIOASCII/src/ioascii.hh
@@ -46,12 +46,12 @@ namespace CarpetIOASCII {
static int TimeToOutput (const cGH* cctkGH, int vindex);
static int TriggerOutput (const cGH* cctkGH, int vindex);
- static int GetGridOffset (const cGH* cctkGH, int dir,
+ static int GetGridOffset (const cGH* cctkGH, int m, int dir,
const char* itempl, const char* iglobal,
const char* ctempl, const char* cglobal,
CCTK_REAL cfallback);
- static int CoordToOffset (const cGH* cctkGH, int dir, CCTK_REAL coord,
- int ifallback);
+ static int CoordToOffset (const cGH* cctkGH, int m, int dir,
+ CCTK_REAL coord, int ifallback);
static void CheckSteerableParameters (const cGH* cctkGH);
static const char* GetStringParameter (const char* parametertemplate);
diff --git a/Carpet/CarpetLib/src/gh.cc b/Carpet/CarpetLib/src/gh.cc
index 6a7738048..fa9ff2125 100644
--- a/Carpet/CarpetLib/src/gh.cc
+++ b/Carpet/CarpetLib/src/gh.cc
@@ -108,7 +108,7 @@ regrid (mregs const & regs)
// Check component consistency
for (int ml=0; ml<mglevels(); ++ml) {
for (int rl=0; rl<reflevels(); ++rl) {
- assert (components(rl)>0);
+ assert (components(rl)>=0);
for (int c=0; c<components(rl); ++c) {
ibbox const & b = extent(ml,rl,c);
ibbox const & b0 = extent(ml,rl,0);
@@ -136,31 +136,34 @@ regrid (mregs const & regs)
bool have_error = false;
for (int ml=0; ml<mglevels(); ++ml) {
for (int rl=1; rl<reflevels(); ++rl) {
- assert (all (extent(ml,rl,0).stride() * reffacts.AT(rl) ==
- extent(ml,rl-1,0).stride() * reffacts.AT(rl-1)));
- // Check contained-ness:
- // first take all coarse grids
- ibset coarse;
- for (int c=0; c<components(rl-1); ++c) {
- coarse += extent(ml,rl-1,c);
- }
- coarse.normalize();
- // then check all fine grids
- for (int c=0; c<components(rl); ++c) {
- ibbox const & fine =
- extent(ml,rl,c).contracted_for(extent(ml,rl-1,0));
- if (not (fine <= coarse)) {
- if (not have_error) {
- cout << "The following components are not properly nested, i.e.," << endl
- << "they are not contained within the next coarser level's components:" << endl;
- have_error = true;
- }
- cout << " ml " << ml << " rl " << rl << " c " << c << ": "
- << fine << endl;
+ if (components(rl)>0) {
+ assert (components(rl-1)>0);
+ assert (all (extent(ml,rl,0).stride() * reffacts.AT(rl) ==
+ extent(ml,rl-1,0).stride() * reffacts.AT(rl-1)));
+ // Check contained-ness:
+ // first take all coarse grids
+ ibset coarse;
+ for (int c=0; c<components(rl-1); ++c) {
+ coarse += extent(ml,rl-1,c);
}
- } // for c
- } // for rl
- } // for ml
+ coarse.normalize();
+ // then check all fine grids
+ for (int c=0; c<components(rl); ++c) {
+ ibbox const & fine =
+ extent(ml,rl,c).contracted_for(extent(ml,rl-1,0));
+ if (not (fine <= coarse)) {
+ if (not have_error) {
+ cout << "The following components are not properly nested, i.e.," << endl
+ << "they are not contained within the next coarser level's components:" << endl;
+ have_error = true;
+ }
+ cout << " ml " << ml << " rl " << rl << " c " << c << ": "
+ << fine << endl;
+ }
+ } // for c
+ } // if c
+ } // for rl
+ } // for ml
if (have_error) {
cout << "The grid hierarchy is:" << endl;
for (int ml=0; ml<mglevels(); ++ml) {
diff --git a/Carpet/CarpetRegrid2/interface.ccl b/Carpet/CarpetRegrid2/interface.ccl
index 7b55ed523..aa6f799da 100644
--- a/Carpet/CarpetRegrid2/interface.ccl
+++ b/Carpet/CarpetRegrid2/interface.ccl
@@ -47,6 +47,42 @@ CCTK_INT FUNCTION ConvertFromPhysicalBoundary \
CCTK_REAL IN ARRAY spacing)
USES FUNCTION ConvertFromPhysicalBoundary
+# The location of the boundary points
+CCTK_INT FUNCTION MultiPatch_GetBoundarySpecification \
+ (CCTK_INT IN map, \
+ CCTK_INT IN size, \
+ CCTK_INT OUT ARRAY nboundaryzones, \
+ CCTK_INT OUT ARRAY is_internal, \
+ CCTK_INT OUT ARRAY is_staggered, \
+ CCTK_INT OUT ARRAY shiftout)
+USES FUNCTION MultiPatch_GetBoundarySpecification
+
+# The overall size of the domain
+CCTK_INT FUNCTION MultiPatch_GetDomainSpecification \
+ (CCTK_INT IN map, \
+ CCTK_INT IN size, \
+ CCTK_REAL OUT ARRAY physical_min, \
+ CCTK_REAL OUT ARRAY physical_max, \
+ CCTK_REAL OUT ARRAY interior_min, \
+ CCTK_REAL OUT ARRAY interior_max, \
+ CCTK_REAL OUT ARRAY exterior_min, \
+ CCTK_REAL OUT ARRAY exterior_max, \
+ CCTK_REAL OUT ARRAY spacing)
+USES FUNCTION MultiPatch_GetDomainSpecification
+
+# Conversion between boundary types
+CCTK_INT FUNCTION MultiPatch_ConvertFromPhysicalBoundary \
+ (CCTK_INT IN map, \
+ CCTK_INT IN size, \
+ CCTK_REAL IN ARRAY physical_min, \
+ CCTK_REAL IN ARRAY physical_max, \
+ CCTK_REAL OUT ARRAY interior_min, \
+ CCTK_REAL OUT ARRAY interior_max, \
+ CCTK_REAL OUT ARRAY exterior_min, \
+ CCTK_REAL OUT ARRAY exterior_max, \
+ CCTK_REAL IN ARRAY spacing)
+USES FUNCTION MultiPatch_ConvertFromPhysicalBoundary
+
# The true prototype of the routine below:
diff --git a/Carpet/CarpetRegrid2/src/regrid.cc b/Carpet/CarpetRegrid2/src/regrid.cc
index 983f7acba..c08705d29 100644
--- a/Carpet/CarpetRegrid2/src/regrid.cc
+++ b/Carpet/CarpetRegrid2/src/regrid.cc
@@ -174,6 +174,68 @@ namespace CarpetRegrid2 {
+ void
+ get_physical_boundary (rvect & physical_lower,
+ rvect & physical_upper,
+ rvect & spacing)
+ {
+ rvect interior_lower, interior_upper;
+ rvect exterior_lower, exterior_upper;
+ if (CCTK_IsFunctionAliased ("MultiPatch_GetDomainSpecification")) {
+ assert (Carpet::map >= 0);
+ CCTK_INT const ierr = MultiPatch_GetDomainSpecification
+ (Carpet::map, dim,
+ & physical_lower[0], & physical_upper[0],
+ & interior_lower[0], & interior_upper[0],
+ & exterior_lower[0], & exterior_upper[0],
+ & spacing[0]);
+ assert (not ierr);
+ } else if (CCTK_IsFunctionAliased ("GetDomainSpecification")) {
+ CCTK_INT const ierr = GetDomainSpecification
+ (dim,
+ & physical_lower[0], & physical_upper[0],
+ & interior_lower[0], & interior_upper[0],
+ & exterior_lower[0], & exterior_upper[0],
+ & spacing[0]);
+ assert (not ierr);
+ } else {
+ assert (0);
+ }
+ }
+
+ void
+ calculate_exterior_boundary (rvect const & physical_lower,
+ rvect const & physical_upper,
+ rvect & exterior_lower,
+ rvect & exterior_upper,
+ rvect const & spacing)
+ {
+ rvect interior_lower, interior_upper;
+ if (CCTK_IsFunctionAliased ("MultiPatch_ConvertFromPhysicalBoundary")) {
+ assert (Carpet::map >= 0);
+ CCTK_INT const ierr = MultiPatch_ConvertFromPhysicalBoundary
+ (Carpet::map, dim,
+ & physical_lower[0], & physical_upper[0],
+ & interior_lower[0], & interior_upper[0],
+ & exterior_lower[0], & exterior_upper[0],
+ & spacing[0]);
+ assert (not ierr);
+ } else if (CCTK_IsFunctionAliased ("ConvertFromPhysicalBoundary")) {
+ CCTK_INT const ierr =
+ ConvertFromPhysicalBoundary
+ (dim,
+ & physical_lower[0], & physical_upper[0],
+ & interior_lower[0], & interior_upper[0],
+ & exterior_lower[0], & exterior_upper[0],
+ & spacing[0]);
+ assert (not ierr);
+ } else {
+ assert (0);
+ }
+ }
+
+
+
extern "C" {
CCTK_INT
CarpetRegrid2_Regrid (CCTK_POINTER_TO_CONST const cctkGH_,
@@ -208,6 +270,7 @@ namespace CarpetRegrid2 {
#warning "TODO: check this (for Carpet, and maybe also for CartGrid3D)"
#warning "TODO: (the check that these two are consistent should be in Carpet)"
+#if 0
typedef vect<vect<CCTK_INT,2>,3> jjvect;
jjvect nboundaryzones;
jjvect is_internal;
@@ -222,7 +285,9 @@ namespace CarpetRegrid2 {
& shiftout[0][0]);
assert (not ierr);
}
+#endif
+#if 0
rvect physical_lower, physical_upper;
rvect interior_lower, interior_upper;
rvect exterior_lower, exterior_upper;
@@ -236,10 +301,16 @@ namespace CarpetRegrid2 {
& spacing[0]);
assert (not ierr);
}
+#endif
+ rvect physical_lower, physical_upper;
+ rvect spacing;
+ get_physical_boundary (physical_lower, physical_upper,
+ spacing);
// Adapt spacing for convergence level
spacing *= ipow ((CCTK_REAL) Carpet::mgfact, Carpet::basemglevel);
+#if 0
{
CCTK_INT const ierr =
ConvertFromPhysicalBoundary
@@ -250,6 +321,11 @@ namespace CarpetRegrid2 {
& spacing[0]);
assert (not ierr);
}
+#endif
+ rvect exterior_lower, exterior_upper;
+ calculate_exterior_boundary (physical_lower, physical_upper,
+ exterior_lower, exterior_upper,
+ spacing);
//
// Calculate the union of the bounding boxes for all levels
@@ -274,38 +350,43 @@ namespace CarpetRegrid2 {
}
regions.at(0).normalize();
- // Loop over all centres
- for (int n = 0; n < num_centres; ++ n) {
- centre_description centre (cctkGH, n);
- if (centre.active) {
-
- // Loop over all levels for this centre
- for (int rl = 1; rl < centre.num_levels; ++ rl) {
-
- // Calculate a bbox for this region
- rvect const rmin = centre.position - centre.radius.at(rl);
- rvect const rmax = centre.position + centre.radius.at(rl);
-
- // Convert to an integer bbox
- ivect const istride = hh.baseextents.at(0).at(rl).stride();
- ivect const imin =
- rpos2ipos (rmin, origin, scale, hh, rl)
- - boundary_shiftout * istride;
- ivect const imax =
- rpos2ipos1 (rmax, origin, scale, hh, rl)
- + boundary_shiftout * istride;
-
- ibbox const region (imin, imax, istride);
+ // Refine only patch 0
+ if (Carpet::map == 0) {
+ // Loop over all centres
+ for (int n = 0; n < num_centres; ++ n) {
+ centre_description centre (cctkGH, n);
+ if (centre.active) {
- // Add this region to the list of regions
- if (static_cast <int> (regions.size()) < rl+1) regions.resize (rl+1);
- regions.at(rl) |= region;
- regions.at(rl).normalize();
+ // Loop over all levels for this centre
+ for (int rl = 1; rl < centre.num_levels; ++ rl) {
+
+ // Calculate a bbox for this region
+ rvect const rmin = centre.position - centre.radius.at(rl);
+ rvect const rmax = centre.position + centre.radius.at(rl);
+
+ // Convert to an integer bbox
+ ivect const istride = hh.baseextents.at(0).at(rl).stride();
+ ivect const imin =
+ rpos2ipos (rmin, origin, scale, hh, rl)
+ - boundary_shiftout * istride;
+ ivect const imax =
+ rpos2ipos1 (rmax, origin, scale, hh, rl)
+ + boundary_shiftout * istride;
+
+ ibbox const region (imin, imax, istride);
+
+ // Add this region to the list of regions
+ if (static_cast <int> (regions.size()) < rl+1) {
+ regions.resize (rl+1);
+ }
+ regions.at(rl) |= region;
+ regions.at(rl).normalize();
+
+ } // for rl
- } // for rl
-
- } // if centre is active
- } // for n
+ } // if centre is active
+ } // for n
+ } // if map==0
@@ -448,6 +529,7 @@ namespace CarpetRegrid2 {
rvect const level_physical_lower = physical_lower;
rvect const level_physical_upper = physical_upper;
rvect const level_spacing = spacing / rvect (hh.reffacts.at(rl));
+#if 0
rvect level_interior_lower, level_interior_upper;
rvect level_exterior_lower, level_exterior_upper;
{
@@ -460,6 +542,11 @@ namespace CarpetRegrid2 {
& level_spacing[0]);
assert (not ierr);
}
+#endif
+ rvect level_exterior_lower, level_exterior_upper;
+ calculate_exterior_boundary (level_physical_lower, level_physical_upper,
+ level_exterior_lower, level_exterior_upper,
+ level_spacing);
ivect const level_exterior_ilower =
rpos2ipos (level_exterior_lower, origin, scale, hh, rl);
@@ -975,20 +1062,20 @@ namespace CarpetRegrid2 {
for (int m = 0; m < maps; ++ m) {
maxrl = max (maxrl, rls.at(m));
}
+ // All maps must have the same number of levels
+ for (int m = 0; m < maps; ++ m) {
+ regsss.at(m).resize (maxrl);
+ }
// Make multiprocessor aware
for (int rl = 0; rl < maxrl; ++ rl) {
vector <vector <region_t> > regss (maps);
for (int m = 0; m < maps; ++ m) {
- if (rl < rls.at(m)) {
- regss.at(m) = regsss.at(m).at(rl);
- }
+ regss.at(m) = regsss.at(m).at(rl);
}
Carpet::SplitRegionsMaps (cctkGH, regss);
for (int m = 0; m < maps; ++ m) {
- if (rl < rls.at(m)) {
- regsss.at(m).at(rl) = regss.at(m);
- }
+ regsss.at(m).at(rl) = regss.at(m);
}
} // for rl
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-25 05:19:16 +0000