aboutsummaryrefslogtreecommitdiff
path: root/Carpet/CarpetRegrid/src/regrid.cc
diff options
context:
space:
mode:
Diffstat (limited to 'Carpet/CarpetRegrid/src/regrid.cc')
-rw-r--r--Carpet/CarpetRegrid/src/regrid.cc955
1 files changed, 108 insertions, 847 deletions
diff --git a/Carpet/CarpetRegrid/src/regrid.cc b/Carpet/CarpetRegrid/src/regrid.cc
index 82f6e0954..3c6b8b0a7 100644
--- a/Carpet/CarpetRegrid/src/regrid.cc
+++ b/Carpet/CarpetRegrid/src/regrid.cc
@@ -1,81 +1,70 @@
#include <assert.h>
-#include <math.h>
-#include <stdarg.h>
-#include <stdlib.h>
-#include <string.h>
-#include <algorithm>
-#include <list>
#include <sstream>
#include <string>
-#include <vector>
#include "cctk.h"
#include "cctk_Parameters.h"
-#include "bbox.hh"
-#include "bboxset.hh"
-#include "defs.hh"
#include "gh.hh"
-#include "gf.hh"
#include "vect.hh"
#include "carpet.hh"
#include "regrid.hh"
extern "C" {
- static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/regrid.cc,v 1.33 2004/01/15 09:45:58 cott Exp $";
+ static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/regrid.cc,v 1.34 2004/01/25 14:57:30 schnetter Exp $";
CCTK_FILEVERSION(Carpet_CarpetRegrid_regrid_cc);
}
-#undef AUTOMATIC_BOUNDARIES
-
-
-
namespace CarpetRegrid {
using namespace std;
using namespace Carpet;
- int initial_refinement_levels = 0; // Store the initial number of
- // levels to refine for correct
- // progressive MR.
- // Note that this does not need
- // to be checkpointed as it is
- // reset at the beginning of every
- // iteration.
- int last_regridded_iteration = -1; // The last iteration where any
- // regridding was done.
-
- int CarpetRegridStartup ()
- {
- RegisterRegridRoutine (CarpetRegridRegrid);
- return 0;
- }
-
-
-
- int CarpetRegridRegrid (const cGH * const cctkGH,
- gh<dim>::rexts& bbsss,
- gh<dim>::rbnds& obss,
- gh<dim>::rprocs& pss)
+ CCTK_INT CarpetRegrid_Regrid (CCTK_POINTER_TO_CONST const cctkGH_,
+ CCTK_INT const reflevel,
+ CCTK_INT const map,
+ CCTK_INT const size,
+ CCTK_INT const * const nboundaryzones_,
+ CCTK_INT const * const is_internal_,
+ CCTK_INT const * const is_staggered_,
+ CCTK_INT const * const shiftout_,
+ CCTK_POINTER const bbsss_,
+ CCTK_POINTER const obss_,
+ CCTK_POINTER const pss_)
{
DECLARE_CCTK_PARAMETERS;
+ const cGH * const cctkGH = (const cGH *) cctkGH_;
+
+ assert (reflevel>=0 && reflevel<maxreflevels);
+ assert (map>=0 && map<maps);
+
+ assert (size == 2*dim);
+ jjvect const nboundaryzones (* (jjvect const *) nboundaryzones_);
+ jjvect const is_internal (* (jjvect const *) is_internal_);
+ jjvect const is_staggered (* (jjvect const *) is_staggered_);
+ jjvect const shiftout (* (jjvect const *) shiftout_);
+
+ gh<dim>::rexts & bbsss = * (gh<dim>::rexts *) bbsss_;
+ gh<dim>::rbnds & obss = * (gh<dim>::rbnds *) obss_;
+ gh<dim>::rprocs & pss = * (gh<dim>::rprocs *) pss_;
+
+ gh<dim> const & hh = *vhh.at(map);
+
+ assert (is_meta_mode());
+
+
+
assert (regrid_every == -1 || regrid_every == 0
|| regrid_every % maxmglevelfact == 0);
// Return if no regridding is desired
if (regrid_every == -1) return 0;
-#if 0
- // Return if this is not the main hierarchy
- if (mglevel != 0) return 0;
-#endif
- assert (mglevel == -1);
-
// Return if we want to regrid during initial data only, and this
// is not the time for initial data
if (regrid_every == 0 && cctkGH->cctk_iteration != 0) return 0;
@@ -85,833 +74,105 @@ namespace CarpetRegrid {
&& cctkGH->cctk_iteration % regrid_every != 0) {
return 0;
}
-
- // Find the number of levels that should be active.
- int newnumlevels;
-
- // Is this the first time that we have regridded this iteration?
- if (cctkGH->cctk_iteration > last_regridded_iteration) {
- last_regridded_iteration = cctkGH->cctk_iteration;
- initial_refinement_levels = refinement_levels;
+
+
+
+ // Steer parameters
+ if (activate_newlevels_on_regrid != 0) {
+ if (cctkGH->cctk_iteration >= activate_next) {
+ const int newnumlevels
+ = min(refinement_levels + activate_newlevels_on_regrid,
+ maxreflevels);
+ assert (newnumlevels>0 && newnumlevels<=maxreflevels);
+
+ *const_cast<CCTK_INT*>(&activate_next) = cctkGH->cctk_iteration + 1;
+ ostringstream next;
+ next << activate_next;
+ CCTK_ParameterSet
+ ("activate_next", "CarpetRegrid", next.str().c_str());
+
+ *const_cast<CCTK_INT*>(&refinement_levels) = newnumlevels;
+ ostringstream param;
+ param << refinement_levels;
+ CCTK_ParameterSet
+ ("refinement_levels", "CarpetRegrid", param.str().c_str());
+
+ }
}
- if (cctkGH->cctk_iteration == 0) {
- newnumlevels = refinement_levels;
- } else {
- if (regrid_from_function) {
- if (CCTK_IsFunctionAliased("RegridLevel")) {
- int tempnewnumlevels = 0;
- newnumlevels = 0;
- BEGIN_MGLEVEL_LOOP(cctkGH) {
- tempnewnumlevels =
- RegridLevel(cctkGH, reflevel, refinement_levels,maxreflevels);
- newnumlevels = max(tempnewnumlevels, newnumlevels);
- } END_MGLEVEL_LOOP;
- } else {
- CCTK_WARN(1, "No thorn has provided the function "
- "\"RegridLevel\". Regridding will not be done.");
- }
- } else {
- // The standard progressive MR number of levels.
- newnumlevels = initial_refinement_levels +
- activate_newlevels_on_regrid;
+ if (regrid_from_function) {
+ if (! CCTK_IsFunctionAliased("RegridLevel")) {
+ CCTK_WARN (0, "No thorn has provided the function \"RegridLevel\"");
}
+ const int newnumlevels
+ = RegridLevel (cctkGH, reflevel, refinement_levels, maxreflevels);
+ assert (newnumlevels>0 && newnumlevels<=maxreflevels);
+
+ *const_cast<CCTK_INT*>(&refinement_levels) = newnumlevels;
+ ostringstream param;
+ param << refinement_levels;
+ CCTK_ParameterSet
+ ("refinement_levels", "CarpetRegrid", param.str().c_str());
}
-
- newnumlevels = min(newnumlevels, maxreflevels);
-
- if ( ( newnumlevels >= 1) && (newnumlevels <= maxreflevels )) {
- char numlevelstring[10];
- sprintf(numlevelstring,"%d",newnumlevels);
- CCTK_ParameterSet("refinement_levels","carpetregrid",numlevelstring);
- }
- list<ibbox> bbl;
- list<bvect> obl;
+
+
+ int do_recompose;
if (CCTK_EQUALS(refined_regions, "none")) {
- MakeRegions_BaseLevel (cctkGH, bbl, obl);
-
+ do_recompose = BaseLevel
+ (cctkGH, hh, reflevel, map,
+ size, nboundaryzones, is_internal, is_staggered, shiftout,
+ bbsss, obss, pss);
+
} else if (CCTK_EQUALS(refined_regions, "centre")) {
- MakeRegions_RefineCentre (cctkGH, newnumlevels, bbl, obl);
-
+ do_recompose = Centre
+ (cctkGH, hh, reflevel, map,
+ size, nboundaryzones, is_internal, is_staggered, shiftout,
+ bbsss, obss, pss);
+
} else if (CCTK_EQUALS(refined_regions, "manual-gridpoints")) {
- if (refinement_levels > 4) {
- CCTK_WARN (0, "Cannot currently specify manual refinement regions for more than 4 refinement levels");
- }
- assert (refinement_levels<=4);
- vector<ivect> lower(3), upper(3);
- lower[0] = ivect (l1ixmin, l1iymin, l1izmin);
- upper[0] = ivect (l1ixmax, l1iymax, l1izmax);
- lower[1] = ivect (l2ixmin, l2iymin, l2izmin);
- upper[1] = ivect (l2ixmax, l2iymax, l2izmax);
- lower[2] = ivect (l3ixmin, l3iymin, l3izmin);
- upper[2] = ivect (l3ixmax, l3iymax, l3izmax);
- MakeRegions_AsSpecified (cctkGH, newnumlevels, lower, upper,
- bbl, obl);
-
+ do_recompose = ManualGridpoints
+ (cctkGH, hh, reflevel, map,
+ size, nboundaryzones, is_internal, is_staggered, shiftout,
+ bbsss, obss, pss);
+
} else if (CCTK_EQUALS(refined_regions, "manual-coordinates")) {
- if (refinement_levels > 4) {
- CCTK_WARN (0, "Cannot currently specify manual refinement regions for more than 4 refinement levels");
- }
- assert (refinement_levels<=4);
- vector<rvect> lower(3), upper(3);
- lower[0] = rvect (l1xmin, l1ymin, l1zmin);
- upper[0] = rvect (l1xmax, l1ymax, l1zmax);
- lower[1] = rvect (l2xmin, l2ymin, l2zmin);
- upper[1] = rvect (l2xmax, l2ymax, l2zmax);
- lower[2] = rvect (l3xmin, l3ymin, l3zmin);
- upper[2] = rvect (l3xmax, l3ymax, l3zmax);
- MakeRegions_AsSpecified (cctkGH, newnumlevels, lower, upper,
- bbl, obl);
-
+ do_recompose = ManualCoordinates
+ (cctkGH, hh, reflevel, map,
+ size, nboundaryzones, is_internal, is_staggered, shiftout,
+ bbsss, obss, pss);
+
} else if (CCTK_EQUALS(refined_regions, "manual-gridpoint-list")) {
- vector<vector<ibbox> > bbss;
- if (strcmp(gridpoints, "") != 0) {
- istringstream gp_str(gridpoints);
- gp_str >> bbss;
- }
-
- vector<vector<vect<vect<bool,2>,dim> > > obss;
- if (strcmp(outerbounds, "") !=0 ) {
- istringstream ob_str (outerbounds);
- ob_str >> obss;
- bool good = obss.size() == bbss.size();
- if (good) {
- for (int rl=0; rl<(int)obss.size(); ++rl) {
- good = good && obss[rl].size() == bbss[rl].size();
- }
- }
- if (! good) {
- cout << "gridpoints: " << bbss << endl;
- cout << "outerbounds: " << obss << endl;
- CCTK_WARN (0, "The parameters \"outerbounds\" and \"gridpoints\" must have the same structure");
- }
- } else {
- obss.resize(bbss.size());
- for (int rl=0; rl<(int)obss.size(); ++rl) {
- obss[rl].resize(bbss[rl].size());
- for (int c=0; c<(int)obss[rl].size(); ++c) {
- obss[rl][c] = vect<vect<bool,2>,dim>(vect<bool,2>(false));
- }
- }
- }
-
- MakeRegions_AsSpecified (cctkGH, newnumlevels, bbss, obss,
- bbl, obl);
-
+ do_recompose = ManualGridpointList
+ (cctkGH, hh, reflevel, map,
+ size, nboundaryzones, is_internal, is_staggered, shiftout,
+ bbsss, obss, pss);
+
} else if (CCTK_EQUALS(refined_regions, "manual-coordinate-list")) {
- vector<vector<rbbox> > bbss;
- if (strcmp(coordinates, "") != 0) {
- istringstream co_str(coordinates);
- co_str >> bbss;
- }
-
- vector<vector<vect<vect<bool,2>,dim> > > obss;
- if (strcmp(outerbounds, "") !=0 ) {
- istringstream ob_str (outerbounds);
- ob_str >> obss;
- bool good = obss.size() == bbss.size();
- if (good) {
- for (int rl=0; rl<(int)obss.size(); ++rl) {
- good = good && obss[rl].size() == bbss[rl].size();
- }
- }
- if (! good) {
- cout << "coordinates: " << bbss << endl;
- cout << "outerbounds: " << obss << endl;
- CCTK_WARN (0, "The parameters \"outerbounds\" and \"coordinates\" must have the same structure");
- }
- } else {
- obss.resize(bbss.size());
- for (int rl=0; rl<(int)obss.size(); ++rl) {
- obss[rl].resize(bbss[rl].size());
- for (int c=0; c<(int)obss[rl].size(); ++c) {
- obss[rl][c] = vect<vect<bool,2>,dim>(vect<bool,2>(false));
- }
- }
- }
-
- MakeRegions_AsSpecified (cctkGH, newnumlevels, bbss, obss,
- bbl, obl);
-
+ do_recompose = ManualCoordinateList
+ (cctkGH, hh, reflevel, map,
+ size, nboundaryzones, is_internal, is_staggered, shiftout,
+ bbsss, obss, pss);
+
} else if (CCTK_EQUALS(refined_regions, "automatic")) {
- const int vi = CCTK_VarIndex (errorvar);
- assert (vi>=0 && vi<CCTK_NumVars());
- const int gi = CCTK_GroupIndexFromVarI (vi);
- assert (gi>=0 && gi<CCTK_NumGroups());
- const int v1 = CCTK_FirstVarIndexI(gi);
- assert (v1>=0 && v1<=vi && v1<CCTK_NumVars());
-
- assert (CCTK_GroupTypeI(gi) == CCTK_GF);
- assert (CCTK_VarTypeI(vi) == CCTK_VARIABLE_REAL);
- assert (CCTK_GroupDimI(gi) == dim);
-
- assert (gi < (int)arrdata.size());
- assert (vi-v1 < (int)arrdata[gi].data.size());
- assert (arrdata[gi].data[vi-v1]);
- const gf<CCTK_REAL,dim>& error
- = *dynamic_cast<const gf<CCTK_REAL,dim>*>(arrdata[gi].data[vi-v1]);
-
- MakeRegions_Adaptively (cctkGH, minwidth, minfraction, maxerror, error,
- bbl, obl);
-
+ do_recompose = Automatic
+ (cctkGH, hh, reflevel, map,
+ size, nboundaryzones, is_internal, is_staggered, shiftout,
+ bbsss, obss, pss);
+
} else {
assert (0);
}
-#ifdef AUTOMATIC_BOUNDARIES
- assert (obl.size() == 0);
-#else
- assert (bbl.size() == obl.size());
-#endif
-
- // transform bbox list into bbox vector
- vector<ibbox> bbs;
- bbs.reserve (bbl.size());
- for (list<ibbox>::const_iterator it = bbl.begin();
- it != bbl.end();
- ++it) {
- bbs.push_back (*it);
- }
-#ifdef AUTOMATIC_BOUNDARIES
- vector<bvect> obs;
- obs.resize (bbs.size());
- for (int c=0; c<(int)bbs.size(); ++c) {
- ivect extlo = hh->baseextent.lower();
- ivect extup = hh->baseextent.upper();
- ivect str = bbs[c].stride();
- for (int d=0; d<dim; ++d) {
- // Decide what is an outer boundary.
- // Allow it to be one grid point to the interior.
- assert (bbs[c].lower()[d] >= extlo[d]);
- obs[c][d][0] = bbs[c].lower()[d] <= extlo[d] + str[d];
- assert (bbs[c].upper()[d] <= extup[d]);
- obs[c][d][1] = bbs[c].upper()[d] >= extup[d] - str[d];
- }
- }
-#else
- vector<bvect> obs;
- obs.reserve (obl.size());
- for (list<bvect>::const_iterator it = obl.begin();
- it != obl.end();
- ++it) {
- obs.push_back (*it);
- }
-#endif
-
-
-
- // TODO: ensure nesting properties
-
- // make multiprocessor aware
- SplitRegions (cctkGH, bbs, obs);
-
- // make multigrid aware
- gh<dim>::cexts bbss;
- bbss = hh->make_reflevel_multigrid_boxes(bbs, mglevels);
-
- // insert into grid hierarchy
- bbsss = hh->extents;
- obss = hh->outer_boundaries;
- if (bbss.size() == 0) {
- // remove all finer levels
- // TODO: this might not work
- bbsss.resize(reflevel+1);
- obss.resize(reflevel+1);
- } else {
- assert (reflevel < (int)bbsss.size());
- if (reflevel+1 == (int)bbsss.size()) {
- // add a finer level
- bbsss.push_back (bbss);
- obss.push_back (obs);
- } else {
- // change a finer level
- bbsss[reflevel+1] = bbss;
- obss[reflevel+1] = obs;
- }
- }
-
- MakeProcessors (cctkGH, bbsss, pss);
-
- return 1; // do recompose
- }
-
-
-
- void MakeRegions_BaseLevel (const cGH* cctkGH,
- list<ibbox>& bbl, list<bvect>& obl)
- {
- assert (bbl.empty());
- assert (obl.empty());
- }
-
-
-
- // This is a helpful helper routine. The user can use it to define
- // how the hierarchy should be refined. But the result of this
- // routine is rather arbitrary.
- void MakeRegions_RefineCentre (const cGH* cctkGH, const int reflevels,
- list<ibbox>& bbl, list<bvect>& obl)
- {
- assert (bbl.empty());
- assert (obl.empty());
-
- if (reflevel+1 >= reflevels) return;
-
- ivect rstr = hh->baseextent.stride();
- ivect rlb = hh->baseextent.lower();
- ivect rub = hh->baseextent.upper();
-
- for (int rl=0; rl<reflevel+1; ++rl) {
- // save old values
- const ivect oldrlb = rlb;
- const ivect oldrub = rub;
- // refined boxes have smaller stride
- assert (all(rstr%hh->reffact == 0));
- rstr /= hh->reffact;
- // calculate new extent
- const ivect quarter = (rub - rlb) / 4 / rstr * rstr;
- rlb = oldrlb + quarter;
- rub = oldrub - quarter;
- assert (all(rlb >= oldrlb && rub <= oldrub));
- }
-
- bbl.push_back (ibbox(rlb, rub, rstr));
- obl.push_back (bvect(vect<bool,2>(false)));
- }
-
-
-
- static void
- MakeRegions_AsSpecified_OneLevel (const cGH* cctkGH, const int reflevels,
- const ivect ilower,
- const ivect iupper,
- const bvect obound,
- list<ibbox>& bbl, list<bvect>& obl)
- {
- if (reflevel+1 >= reflevels) return;
-
- const ivect rstr = hh->baseextent.stride();
- const ivect rlb = hh->baseextent.lower();
- const ivect rub = hh->baseextent.upper();
-
- const int rl = reflevel+1;
-
- const int levfac = ipow(hh->reffact, rl);
- assert (all (rstr % levfac == 0));
- const ivect str (rstr / levfac);
- const ivect lb (ilower);
- const ivect ub (iupper);
- if (! all(lb>=rlb && ub<=rub)) {
- ostringstream buf;
- buf << "The refinement region boundaries for refinement level #" << rl << " are not within the main grid. Allowed are the grid point boundaries " << rlb << " - " << rub << "; specified were " << lb << " - " << ub << ends;
- CCTK_WARN (0, buf.str().c_str());
- }
- if (! all(lb<=ub)) {
- CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
- "The refinement region boundaries for refinement level #%d have the upper boundary less than the lower boundary", rl);
- }
- if (! all(lb%str==0 && ub%str==0)) {
- CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
- "The refinement region boundaries for refinement level #%d are not a multiple of the stride for that level", rl);
- }
- assert (all(lb>=rlb && ub<=rub));
- assert (all(lb<=ub));
- assert (all(lb%str==0 && ub%str==0));
-
- bbl.push_back (ibbox(lb, ub, str));
- obl.push_back (obound);
- }
-
-
-
- void MakeRegions_AsSpecified (const cGH* cctkGH, const int reflevels,
- const vector<ivect> lower,
- const vector<ivect> upper,
- list<ibbox>& bbl, list<bvect>& obl)
- {
- assert (lower.size() == upper.size());
-
- if (reflevel+1 >= reflevels) return;
-
- const int rl = reflevel+1;
-
- const ivect ilower = lower[rl-1];
- const ivect iupper = upper[rl-1];
- const bvect obound = bvect(vect<bool,2>(false));
-
- MakeRegions_AsSpecified_OneLevel (cctkGH, reflevels,
- ilower, iupper, obound,
- bbl, obl);
- }
-
-
-
- static ivect delta2int (const cGH* cctkGH, const rvect & rpos, const int rl)
- {
- rvect global_lower, global_upper;
- for (int d=0; d<dim; ++d) {
- const int ierr = CCTK_CoordRange
- (cctkGH, &global_lower[d], &global_upper[d], d+1, 0, "cart3d");
- if (ierr<0) {
- global_lower[d] = 0;
- global_upper[d] = 1;
- }
- }
- const ivect global_extent (hh->baseextent.upper() - hh->baseextent.lower());
-
- CCTK_REAL (* const rfloor) (CCTK_REAL const) = floor;
-
- const rvect scale = rvect(global_extent) / (global_upper - global_lower);
- const int levfac = ipow(hh->reffact, rl);
- assert (all (hh->baseextent.stride() % levfac == 0));
- const ivect istride = hh->baseextent.stride() / levfac;
-
- const ivect ipos = ivect(map(rfloor, rpos * scale / rvect(istride) + 0.5)) * istride;
-
- const rvect apos = rpos * scale;
- assert (all(abs(apos - rvect(ipos)) < rvect(istride)*0.01));
-
- return ipos;
- }
-
-
-
- static ivect pos2int (const cGH* cctkGH, const rvect & rpos, const int rl)
- {
- rvect global_lower, global_upper;
- for (int d=0; d<dim; ++d) {
- const int ierr = CCTK_CoordRange
- (cctkGH, &global_lower[d], &global_upper[d], d+1, 0, "cart3d");
- if (ierr<0) {
- global_lower[d] = 0;
- global_upper[d] = 1;
- }
- }
- const ivect global_extent (hh->baseextent.upper() - hh->baseextent.lower());
-
- CCTK_REAL (* const rfloor) (CCTK_REAL const) = floor;
-
- const rvect scale = rvect(global_extent) / (global_upper - global_lower);
- const int levfac = ipow(hh->reffact, rl);
- assert (all (hh->baseextent.stride() % levfac == 0));
- const ivect istride = hh->baseextent.stride() / levfac;
-
- const ivect ipos = ivect(map(rfloor, (rpos - global_lower) * scale / rvect(istride) + 0.5)) * istride;
-
- return ipos;
- }
-
-
-
-#if 0
- void MakeRegions_AsSpecified (const cGH* cctkGH, const int reflevels,
- const vector<rvect> lower,
- const vector<rvect> upper,
- list<ibbox>& bbl, list<bvect>& obl)
- {
- assert (lower.size() == upper.size());
-
- if (reflevel+1 >= reflevels) return;
-
- rvect global_lower, global_upper;
- for (int d=0; d<dim; ++d) {
- const int ierr = CCTK_CoordRange
- (cctkGH, &global_lower[d], &global_upper[d], d+1, 0, "cart3d");
- if (ierr<0) {
- global_lower[d] = 0;
- global_upper[d] = 1;
- }
- }
- const ivect global_extent (hh->baseextent.upper() - hh->baseextent.lower());
-
- const int rl = reflevel+1;
-
- const rvect scale = rvect(global_extent) / (global_upper - global_lower);
- const rvect rlower = (lower[rl-1] - global_lower) * scale;
- const rvect rupper = (upper[rl-1] - global_lower) * scale;
-// const ivect ilower = ivect(map((CCTK_REAL(*)(CCTK_REAL))floor, rlower + 0.5));
-// const ivect iupper = ivect(map((CCTK_REAL(*)(CCTK_REAL))floor, rupper + 0.5));
- CCTK_REAL (* const rfloor) (CCTK_REAL const) = floor;
- const int levfac = ipow(hh->reffact, rl);
- assert (all (hh->baseextent.stride() % levfac == 0));
- const ivect istride = hh->baseextent.stride() / levfac;
- const ivect ilower = ivect(map(rfloor, rlower / rvect(istride) + 0.5)) * istride;
- const ivect iupper = ivect(map(rfloor, rupper / rvect(istride) + 0.5)) * istride;
- const bvect obound = bvect(vect<bool,2>(false));
-
- MakeRegions_AsSpecified_OneLevel (cctkGH, reflevels,
- ilower, iupper, obound,
- bbl, obl);
- }
-#endif
-
-
-
- void MakeRegions_AsSpecified (const cGH* cctkGH, const int reflevels,
- const vector<rvect> lower,
- const vector<rvect> upper,
- list<ibbox>& bbl, list<bvect>& obl)
- {
- assert (lower.size() == upper.size());
- vector<ivect> ilower, iupper;
- ilower.resize (lower.size());
- iupper.resize (upper.size());
- for (int rl=1; rl<lower.size()+1; ++rl) {
- ilower[rl-1] = pos2int (cctkGH, lower[rl-1], rl);
- iupper[rl-1] = pos2int (cctkGH, upper[rl-1], rl);
- }
- MakeRegions_AsSpecified (cctkGH, reflevels, ilower, iupper, bbl, obl);
- }
-
-
-
- void MakeRegions_AsSpecified (const cGH* cctkGH, const int reflevels,
- const vector<vector<ibbox> > bbss,
- const vector<vector<bvect> > obss,
- list<ibbox>& bbl, list<bvect>& obl)
- {
- if (reflevel+1 >= reflevels) return;
-
- const int rl = reflevel+1;
- assert (rl-1 < (int)bbss.size());
- assert (obss.size() == bbss.size());
-
- for (int c=0; c<(int)bbss[rl-1].size(); ++c) {
- assert (obss[rl-1].size() == bbss[rl-1].size());
-
- const ivect ilower = bbss[rl-1][c].lower();
- const ivect iupper = bbss[rl-1][c].upper();
- const bvect obound = obss[rl-1][c];
-
- MakeRegions_AsSpecified_OneLevel (cctkGH, reflevels,
- ilower, iupper, obound,
- bbl, obl);
-
- }
- }
-
-
-
- void MakeRegions_AsSpecified (const cGH* cctkGH, const int reflevels,
- const vector<vector<rbbox> > bbss,
- const vector<vector<bvect> > obss,
- list<ibbox>& bbl, list<bvect>& obl)
- {
- vector<vector<ibbox> > ibbss;
- ibbss.resize (bbss.size());
- for (int rl=1; rl<bbss.size()+1; ++rl) {
- ibbss[rl-1].resize (bbss[rl-1].size());
- for (int c=0; c<bbss[rl-1].size(); ++c) {
- const ivect ilower = pos2int (cctkGH, bbss[rl-1][c].lower(), rl);
- const ivect iupper = pos2int (cctkGH, bbss[rl-1][c].upper(), rl);
- const ivect istride = delta2int (cctkGH, bbss[rl-1][c].stride(), rl);
- ibbss[rl-1][c] = ibbox (ilower, iupper, istride);
- }
- }
- MakeRegions_AsSpecified (cctkGH, reflevels, ibbss, obss, bbl, obl);
- }
-
-
-
- static void
- MakeRegions_Adaptively_Recombine (list<ibbox>& bbl1,
- list<ibbox>& bbl2,
- list<ibbox>& bbl,
- const ibbox& iface,
- const int dir)
- {
- assert (!iface.empty());
- assert (iface.lower()[dir] == iface.upper()[dir]);
-
- const int oldnumboxes = bbl.size() + bbl1.size() + bbl2.size();
- int numcombinedboxes = 0;
-
- int oldnumpoints = 0;
- for (list<ibbox>::const_iterator ibb = bbl.begin(); ibb != bbl.end(); ++ibb) {
- oldnumpoints += ibb->size();
- }
- for (list<ibbox>::const_iterator ibb1 = bbl1.begin(); ibb1 != bbl1.end(); ++ibb1) {
- oldnumpoints += ibb1->size();
- }
- for (list<ibbox>::const_iterator ibb2 = bbl2.begin(); ibb2 != bbl2.end(); ++ibb2) {
- oldnumpoints += ibb2->size();
- }
-
-#if 0
- // remember old bounding boxes
- bboxset<int,dim> oldboxes;
- for (list<ibbox>::const_iterator ibb1 = bbl1.begin(); ibb1 != bbl1.end(); ++ibb1) {
- oldboxes += *ibb1;
- }
- for (list<ibbox>::const_iterator ibb2 = bbl2.begin(); ibb2 != bbl2.end(); ++ibb2) {
- oldboxes += *ibb2;
- }
-#endif
-#if 0
- cout << endl;
- cout << "MakeRegions_Adaptively_Recombine: initial list:" << endl;
- cout << bbl << endl;
- cout << "MakeRegions_Adaptively_Recombine: initial list 1:" << endl;
- cout << bbl1 << endl;
- cout << "MakeRegions_Adaptively_Recombine: initial list 2:" << endl;
- cout << bbl2 << endl;
-#endif
-
- const ivect lo = iface.lower();
- const ivect up = iface.upper();
- const ivect str = iface.stride();
-
- {
- // prune boxes on the left
- list<ibbox>::iterator ibb1 = bbl1.begin();
- while (ibb1 != bbl1.end()) {
- // is this bbox just to the left of the interface?
-// const ivect lo1 = ibb1->lower();
- const ivect up1 = ibb1->upper();
- const ivect str1 = ibb1->stride();
- assert (up1[dir]+str1[dir] <= lo[dir]);
- assert (all(str1 == str));
- if (up1[dir]+str1[dir] < lo[dir]) {
- // no: forget it
- bbl.push_back (*ibb1);
- ibb1 = bbl1.erase(ibb1);
- continue;
- }
- ++ibb1;
- } // while
- }
-
- {
- // prune boxes on the right
- list<ibbox>::iterator ibb2 = bbl2.begin();
- while (ibb2 != bbl2.end()) {
- // is this bbox just to the right of the interface?
- const ivect lo2 = ibb2->lower();
-// const ivect up2 = ibb2->upper();
- const ivect str2 = ibb2->stride();
- assert (up[dir] <= lo2[dir]);
- assert (all(str2 == str));
- if (up[dir] < lo2[dir]) {
- // no: forget it
- bbl.push_back (*ibb2);
- ibb2 = bbl2.erase(ibb2);
- continue;
- }
- ++ibb2;
- } // while
- }
-
- {
- // walk all boxes on the left
- list<ibbox>::iterator ibb1 = bbl1.begin();
- while (ibb1 != bbl1.end()) {
- ivect lo1 = ibb1->lower();
- ivect up1 = ibb1->upper();
- ivect str1 = ibb1->stride();
- assert (up1[dir]+str1[dir] == lo[dir]);
- lo1[dir] = lo[dir];
- up1[dir] = up[dir];
- const ibbox iface1 (lo1,up1,str1);
-
- {
- // walk all boxes on the right
- list<ibbox>::iterator ibb2 = bbl2.begin();
- while (ibb2 != bbl2.end()) {
- ivect lo2 = ibb2->lower();
- ivect up2 = ibb2->upper();
- ivect str2 = ibb2->stride();
- assert (lo2[dir] == up[dir]);
- lo2[dir] = lo[dir];
- up2[dir] = up[dir];
- const ibbox iface2 (lo2,up2,str2);
-
- // check for a match
- if (iface1 == iface2) {
- const ibbox combined (ibb1->lower(), ibb2->upper(), str);
- bbl.push_back (combined);
- ibb1 = bbl1.erase(ibb1);
- ibb2 = bbl2.erase(ibb2);
- ++numcombinedboxes;
-// cout << "MRA: Combining along " << "xyz"[dir] << " to " << bbl.back() << " size " << bbl.back().size() << endl;
- goto continue_search;
- }
-
- ++ibb2;
- } // while
- }
-
- ++ibb1;
- continue_search:;
- } // while
- }
-
- bbl.splice (bbl.end(), bbl1);
- bbl.splice (bbl.end(), bbl2);
-
- assert (bbl1.empty() && bbl2.empty());
-
- const int newnumboxes = bbl.size();
- assert (newnumboxes + numcombinedboxes == oldnumboxes);
-
- int newnumpoints = 0;
- for (list<ibbox>::const_iterator ibb = bbl.begin(); ibb != bbl.end(); ++ibb) {
- newnumpoints += ibb->size();
- }
- assert (newnumpoints == oldnumpoints);
-
-#if 0
- // find new bounding boxes
- bboxset<int,dim> newboxes;
- for (list<ibbox>::const_iterator ibb = bbl.begin(); ibb != bbl.end(); ++ibb) {
- newboxes += *ibb;
- }
-
- // Check that they are equal
- assert (newboxes.size() <= oldboxes.size());
- assert ((newboxes.size()==0) == (oldboxes.size()==0));
- assert (oldboxes == newboxes);
-#endif
-#if 0
- cout << "MakeRegions_Adaptively_Recombine: final list:" << endl;
- cout << bbl << endl;
- cout << endl;
-#endif
- }
-
- static void
- MakeRegions_Adaptively (const cGH* const cctkGH,
- const int minwidth,
- const CCTK_REAL minfraction,
- const CCTK_REAL maxerror,
- const data<CCTK_REAL,dim>& error,
- list<ibbox>& bbl,
- const ibbox& region)
- {
- // Just to be sure
- assert (! region.empty());
-
- // Count grid points that need to be refined
- // (this doesn't work yet on multiple processors)
- assert (CCTK_nProcs(cctkGH)==1);
- int cnt = 0;
- for (ibbox::iterator it=region.begin(); it!=region.end(); ++it) {
- if (error[*it] > maxerror) ++cnt;
- }
- const CCTK_REAL fraction = (CCTK_REAL)cnt / region.size();
- const int width = maxval(region.shape() / region.stride());
-
- if (cnt == 0) {
- // Don't refine
- } else if (width < 2*minwidth || fraction >= minfraction) {
- // Refine the whole region
- const ivect lo(region.lower());
- const ivect up(region.upper());
- const ivect str(region.stride());
- bbl.push_back (ibbox(lo,up+str-str/reffact,str/reffact));
-// cout << "MRA: Refining to " << bbl.back() << " size " << bbl.back().size() << " fraction " << fraction << endl;
- } else {
- // Split the region and check recursively
- const int dir = maxloc(region.shape());
- const ivect lo(region.lower());
- const ivect up(region.upper());
- const ivect str(region.stride());
- ivect lo1(lo), lo2(lo);
- ivect up1(up), up2(up);
- const int mgstr = ipow(hh->mgfact, mglevels); // honour multigrid factors
- const int step = str[dir]*mgstr;
- lo2[dir] = ((up[dir]+lo[dir])/2/step)*step;
- up1[dir] = lo2[dir]-str[dir];
- const ibbox region1(lo1,up1,str);
- const ibbox region2(lo2,up2,str);
- assert (region1.is_contained_in(region));
- assert (region2.is_contained_in(region));
- assert ((region1 & region2).empty());
- assert (region1 + region2 == region);
- list<ibbox> bbl1, bbl2;
- MakeRegions_Adaptively (cctkGH, minwidth, minfraction, maxerror,
- error, bbl1, region1);
- MakeRegions_Adaptively (cctkGH, minwidth, minfraction, maxerror,
- error, bbl2, region2);
- // Combine regions if possible
- up2 += str-str/reffact;
- up2[dir] = lo2[dir];
- const ibbox iface(lo2,up2,str/reffact);
- MakeRegions_Adaptively_Recombine (bbl1, bbl2, bbl, iface, dir);
- }
-
- }
-
- void
- MakeRegions_Adaptively (const cGH* const cctkGH,
- const int minwidth,
- const CCTK_REAL minfraction,
- const CCTK_REAL maxerror,
- const gf<CCTK_REAL,dim>& error,
- list<ibbox>& bbl, list<bvect>& obl)
- {
- assert (bbl.empty());
-
- if (reflevel+1 >= maxreflevels) return;
-
- assert (component == -1);
-
- const int rl = reflevel;
-
- // Arbitrary
- const int tl = 0;
- const int ml = 0;
-
-// cout << endl << "MRA: Choosing regions to refine in " << hh->components(rl) << " components" << endl;
-
- for (int c=0; c<hh->components(rl); ++c) {
- const ibbox region = hh->extents[rl][c][ml];
- assert (! region.empty());
-
- const data<CCTK_REAL,dim>& errdata = *error(tl,rl,c,ml);
-
- MakeRegions_Adaptively (cctkGH, minwidth, minfraction, maxerror,
- errdata, bbl, region);
- }
-
-// int numpoints = 0;
-// for (list<ibbox>::const_iterator ibb = bbl.begin(); ibb != bbl.end(); ++ibb) {
-// numpoints += ibb->size();
-// }
-// cout << "MRA: Chose " << bbl.size() << " regions with a total size of " << numpoints << " to refine." << endl << endl;
-
- // Remove grid points outside the outer boundary
- vect<bool,dim> obp;
- {
- DECLARE_CCTK_PARAMETERS;
- assert (sizeof outside_boundary_points == dim * sizeof (CCTK_INT));
- obp = outside_boundary_points;
- }
- for (list<ibbox>::iterator it = bbl.begin();
- it != bbl.end();
- ++it) {
- const ivect ub = obp.ifthen (it->upper(),
- min (it->upper(), hh->baseextent.upper()));
- *it = ibbox(it->lower(), ub, it->stride());
- }
-
- // Create obl from bbl
- for (list<ibbox>::const_iterator it = bbl.begin();
- it != bbl.end();
- ++it) {
- obl.push_back (zip ((vect<bool,2> (*) (bool, bool)) &vect<bool,2>::make,
- it->lower() == hh->baseextent.lower(),
- it->upper() == hh->baseextent.upper()));
- }
-
+ return do_recompose;
}
} // namespace CarpetRegrid
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-09 07:23:51 +0000