diff options
Diffstat (limited to 'Carpet/CarpetRegrid')
| -rw-r--r-- | Carpet/CarpetRegrid/interface.ccl | 29 | ||||
| -rw-r--r-- | Carpet/CarpetRegrid/param.ccl | 38 | ||||
| -rw-r--r-- | Carpet/CarpetRegrid/src/automatic.cc | 60 | ||||
| -rw-r--r-- | Carpet/CarpetRegrid/src/baselevel.cc | 12 | ||||
| -rw-r--r-- | Carpet/CarpetRegrid/src/centre.cc | 21 | ||||
| -rw-r--r-- | Carpet/CarpetRegrid/src/make.code.defn | 12 | ||||
| -rw-r--r-- | Carpet/CarpetRegrid/src/manualcoordinatelist.cc | 132 | ||||
| -rw-r--r-- | Carpet/CarpetRegrid/src/manualcoordinates.cc | 51 | ||||
| -rw-r--r-- | Carpet/CarpetRegrid/src/manualgridpointlist.cc | 20 | ||||
| -rw-r--r-- | Carpet/CarpetRegrid/src/manualgridpoints.cc | 38 | ||||
| -rw-r--r-- | Carpet/CarpetRegrid/src/regrid.cc | 955 | ||||
| -rw-r--r-- | Carpet/CarpetRegrid/src/regrid.hh | 219 |
12 files changed, 512 insertions, 1075 deletions
diff --git a/Carpet/CarpetRegrid/interface.ccl b/Carpet/CarpetRegrid/interface.ccl index 3c6352fda..efc701a73 100644 --- a/Carpet/CarpetRegrid/interface.ccl +++ b/Carpet/CarpetRegrid/interface.ccl @@ -1,5 +1,5 @@ # Interface definition for thorn CarpetRegrid -# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/interface.ccl,v 1.5 2004/01/15 09:45:58 cott Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/interface.ccl,v 1.6 2004/01/25 14:57:30 schnetter Exp $ implements: CarpetRegrid @@ -14,8 +14,31 @@ uses include header: vect.hh uses include header: gf.hh uses include header: gh.hh -CCTK_INT FUNCTION RegridLevel (CCTK_POINTER_TO_CONST IN cctkGH, \ - CCTK_INT IN current_reflevel, \ + + +# The true prototype of the routine below: +# int Carpet_Regrid (const cGH * cctkGH, +# gh<dim>::rexts * bbsss, +# gh<dim>::rbnds * obss, +# gh<dim>::rprocs * pss); +CCTK_INT FUNCTION Carpet_Regrid (CCTK_POINTER_TO_CONST IN cctkGH, \ + CCTK_INT IN reflevel, \ + CCTK_INT IN map, \ + CCTK_INT IN size, \ + CCTK_INT IN ARRAY nboundaryzones, \ + CCTK_INT IN ARRAY is_internal, \ + CCTK_INT IN ARRAY is_staggered, \ + CCTK_INT IN ARRAY shiftout, \ + CCTK_POINTER IN bsss, \ + CCTK_POINTER IN obss, \ + CCTK_POINTER IN pss) +PROVIDES FUNCTION Carpet_Regrid WITH CarpetRegrid_Regrid LANGUAGE C + + + + +CCTK_INT FUNCTION RegridLevel (CCTK_POINTER_TO_CONST IN cctkGH, \ + CCTK_INT IN current_reflevel, \ CCTK_INT IN current_max_reflevel, \ CCTK_INT IN max_reflevels) USES FUNCTION RegridLevel diff --git a/Carpet/CarpetRegrid/param.ccl b/Carpet/CarpetRegrid/param.ccl index 9f2dc9bf9..1de876c73 100644 --- a/Carpet/CarpetRegrid/param.ccl +++ b/Carpet/CarpetRegrid/param.ccl @@ -1,5 +1,5 @@ # Parameter definitions for thorn CarpetRegrid -# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/param.ccl,v 1.13 2004/01/13 13:50:05 hawke Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/param.ccl,v 1.14 2004/01/25 14:57:30 schnetter Exp $ @@ -17,16 +17,19 @@ CCTK_INT regrid_every "Regrid every n time steps" STEERABLE=always 1:* :: "regrid every n time steps" } 0 -CCTK_INT activate_newlevels_on_regrid "When regridding, activate this many new levels (if possible)" -# Note that this feature will STEER the parameter refinement_levels!!! -# It does this by adding activate_newlevels_on_regrid to refinement levels + + +CCTK_INT activate_newlevels_on_regrid "When regridding, activate this many new levels (if possible). Note that this will steer the parameter refinement_levels." STEERABLE=always { - : :: "Number of new levels to activate. Negative numbers mean to de-activate" + : :: "Number of new levels to activate (negative numbers deactivate)" } 0 +CCTK_INT activate_next "The next iteration at which new levels should be activated" STEERABLE=always +{ + 0: :: "Note that this parameter is steered when new levels are activated" +} 1 - -BOOLEAN outside_boundary_points[3] "On finer grids, where the upper grid boundary is adjacent to the outer boundary, put points outside the outer boundary (needed e.g. for periodicity)" +BOOLEAN regrid_from_function "Regridding criteria from the aliased function?" { } "no" @@ -45,6 +48,20 @@ KEYWORD refined_regions "Regions where the grid is refined" STEERABLE=always +# Region specifications for centre refinement + +BOOLEAN symmetry_x "Refine the lower half in x-direction" +{ +} "no" +BOOLEAN symmetry_y "Refine the lower half in y-direction" +{ +} "no" +BOOLEAN symmetry_z "Refine the lower half in z-direction" +{ +} "no" + + + # Region specifications for manual gridpoint refinement CCTK_INT l1ixmin "Lower boundary of level 1 box in x-direction" STEERABLE=always @@ -268,10 +285,3 @@ CCTK_STRING errorvar "Name of grid function that contains the error" STEERABLE=a { ".*" :: "must be the name of a grid function" } "" - - -# Should we set the regridding criteria from an aliased function? - -BOOLEAN regrid_from_function "Regridding criteria from the aliased function?" -{ -} "no" diff --git a/Carpet/CarpetRegrid/src/automatic.cc b/Carpet/CarpetRegrid/src/automatic.cc index a3b4da055..6add8c770 100644 --- a/Carpet/CarpetRegrid/src/automatic.cc +++ b/Carpet/CarpetRegrid/src/automatic.cc @@ -16,7 +16,7 @@ #include "regrid.hh" extern "C" { - static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/automatic.cc,v 1.5 2004/08/04 16:25:58 schnetter Exp $"; + static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/automatic.cc,v 1.1 2004/01/25 14:57:30 schnetter Exp $"; CCTK_FILEVERSION(Carpet_CarpetRegrid_automatic_cc); } @@ -31,12 +31,22 @@ namespace CarpetRegrid { int Automatic (cGH const * const cctkGH, gh<dim> const & hh, + int const reflevel, + int const map, + int const size, + jjvect const & nboundaryzones, + jjvect const & is_internal, + jjvect const & is_staggered, + jjvect const & shiftout, gh<dim>::rexts & bbsss, gh<dim>::rbnds & obss, gh<dim>::rprocs & pss) { DECLARE_CCTK_PARAMETERS; + assert (reflevel>=0 && reflevel<maxreflevels); + assert (map>=0 && map<maps); + assert (refinement_levels >= 1); assert (bbsss.size() >= 1); @@ -54,19 +64,15 @@ namespace CarpetRegrid { assert (CCTK_VarTypeI(vi) == CCTK_VARIABLE_REAL); assert (CCTK_GroupDimI(gi) == dim); - assert (arrdata.at(gi).at(Carpet::map).data.at(vi-v1)); - const gf<CCTK_REAL,dim>& errorgf + assert (arrdata.at(gi).at(map).data.at(vi-v1)); + const gf<CCTK_REAL,dim>& errorvar = (*dynamic_cast<const gf<CCTK_REAL,dim>*> - (arrdata.at(gi).at(Carpet::map).data.at(vi-v1))); - - assert (! smart_outer_boundaries); + (arrdata.at(gi).at(map).data.at(vi-v1))); vector<ibbox> bbs; gh<dim>::cbnds obs; Automatic_OneLevel - (cctkGH, hh, - reflevel, min(reflevels+1, maxreflevels), - minwidth, minfraction, maxerror, errorgf, + (cctkGH, hh, reflevel, minwidth, minfraction, maxerror, errorvar, bbs, obs); // make multiprocessor aware @@ -75,7 +81,10 @@ namespace CarpetRegrid { // make multigrid aware vector<vector<ibbox> > bbss; - MakeMultigridBoxes (cctkGH, bbs, obs, bbss); + MakeMultigridBoxes + (cctkGH, + size, nboundaryzones, is_internal, is_staggered, shiftout, + bbs, obs, bbss); @@ -106,32 +115,31 @@ namespace CarpetRegrid { void Automatic_OneLevel (const cGH * const cctkGH, const gh<dim> & hh, - const int rl, - const int numrl, + const int reflevel, const int minwidth, const CCTK_REAL minfraction, const CCTK_REAL maxerror, - const gf<CCTK_REAL,dim> & errorgf, + const gf<CCTK_REAL,dim> & errorvar, vector<ibbox> & bbs, vector<bbvect> & obs) { - if (rl+1 >= numrl) return; + if (reflevel+1 >= maxreflevels) return; // Arbitrary const int tl = 0; const int ml = 0; -// cout << endl << "MRA: Choosing regions to refine in " << hh.components(rl) << " components" << endl; +// cout << endl << "MRA: Choosing regions to refine in " << hh.components(reflevel) << " components" << endl; list<ibbox> bbl; - for (int c=0; c<hh.components(rl); ++c) { - const ibbox region = hh.extents.at(rl).at(c).at(ml); + for (int c=0; c<hh.components(reflevel); ++c) { + const ibbox region = hh.extents.at(reflevel).at(c).at(ml); assert (! region.empty()); - const data<CCTK_REAL,dim>& errordata = *errorgf(tl,rl,c,ml); + const data<CCTK_REAL,dim>& errdata = *errorvar(tl,reflevel,c,ml); Automatic_Recursive (cctkGH, hh, minwidth, minfraction, maxerror, - errordata, bbl, region); + errdata, bbl, region); } // int numpoints = 0; @@ -177,7 +185,7 @@ namespace CarpetRegrid { const int minwidth, const CCTK_REAL minfraction, const CCTK_REAL maxerror, - const data<CCTK_REAL,dim> & errordata, + const data<CCTK_REAL,dim> & errorvar, list<ibbox> & bbl, const ibbox & region) { @@ -188,12 +196,8 @@ namespace CarpetRegrid { // (this doesn't work yet on multiple processors) assert (CCTK_nProcs(cctkGH)==1); int cnt = 0; - { - ibbox::iterator it=region.begin(); - do { - if (errordata[*it] > maxerror) ++cnt; - ++it; - } while (it!=region.end()); + for (ibbox::iterator it=region.begin(); it!=region.end(); ++it) { + if (errorvar[*it] > maxerror) ++cnt; } const CCTK_REAL fraction = (CCTK_REAL)cnt / region.size(); const int width = maxval(region.shape() / region.stride()); @@ -227,9 +231,9 @@ namespace CarpetRegrid { assert (region1 + region2 == region); list<ibbox> bbl1, bbl2; Automatic_Recursive (cctkGH, hh, minwidth, minfraction, maxerror, - errordata, bbl1, region1); + errorvar, bbl1, region1); Automatic_Recursive (cctkGH, hh, minwidth, minfraction, maxerror, - errordata, bbl2, region2); + errorvar, bbl2, region2); // Combine regions if possible up2 += str-str/reffact; up2[dir] = lo2[dir]; diff --git a/Carpet/CarpetRegrid/src/baselevel.cc b/Carpet/CarpetRegrid/src/baselevel.cc index c380c6844..18cc8a95c 100644 --- a/Carpet/CarpetRegrid/src/baselevel.cc +++ b/Carpet/CarpetRegrid/src/baselevel.cc @@ -9,7 +9,7 @@ #include "regrid.hh" extern "C" { - static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/baselevel.cc,v 1.2 2004/04/18 13:29:43 schnetter Exp $"; + static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/baselevel.cc,v 1.1 2004/01/25 14:57:30 schnetter Exp $"; CCTK_FILEVERSION(Carpet_CarpetRegrid_baselevel_cc); } @@ -24,12 +24,22 @@ namespace CarpetRegrid { int BaseLevel (cGH const * const cctkGH, gh<dim> const & hh, + int const reflevel, + int const map, + int const size, + jjvect const & nboundaryzones, + jjvect const & is_internal, + jjvect const & is_staggered, + jjvect const & shiftout, gh<dim>::rexts & bbsss, gh<dim>::rbnds & obss, gh<dim>::rprocs & pss) { DECLARE_CCTK_PARAMETERS; + assert (reflevel>=0 && reflevel<maxreflevels); + assert (map>=0 && map<maps); + assert (refinement_levels == 1); assert (bbsss.size() == 1); diff --git a/Carpet/CarpetRegrid/src/centre.cc b/Carpet/CarpetRegrid/src/centre.cc index 7599324ea..46a23f794 100644 --- a/Carpet/CarpetRegrid/src/centre.cc +++ b/Carpet/CarpetRegrid/src/centre.cc @@ -9,7 +9,7 @@ #include "regrid.hh" extern "C" { - static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/centre.cc,v 1.4 2004/04/28 15:45:25 schnetter Exp $"; + static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/centre.cc,v 1.1 2004/01/25 14:57:30 schnetter Exp $"; CCTK_FILEVERSION(Carpet_CarpetRegrid_centre_cc); } @@ -24,16 +24,26 @@ namespace CarpetRegrid { int Centre (cGH const * const cctkGH, gh<dim> const & hh, + int const reflevel, + int const map, + int const size, + jjvect const & nboundaryzones, + jjvect const & is_internal, + jjvect const & is_staggered, + jjvect const & shiftout, gh<dim>::rexts & bbsss, gh<dim>::rbnds & obss, gh<dim>::rprocs & pss) { DECLARE_CCTK_PARAMETERS; + assert (reflevel>=0 && reflevel<maxreflevels); + assert (map>=0 && map<maps); + assert (refinement_levels >= 1); // do nothing if the levels already exist - if (reflevel == refinement_levels) return 0; + if (bbsss.size() == refinement_levels) return 0; assert (bbsss.size() >= 1); @@ -48,8 +58,6 @@ namespace CarpetRegrid { ivect rlb = hh.baseextent.lower(); ivect rub = hh.baseextent.upper(); - assert (! smart_outer_boundaries); - for (size_t rl=1; rl<bbsss.size(); ++rl) { // save old values @@ -81,7 +89,10 @@ namespace CarpetRegrid { // make multigrid aware vector<vector<ibbox> > bbss; - MakeMultigridBoxes (cctkGH, bbs, obs, bbss); + MakeMultigridBoxes + (cctkGH, + size, nboundaryzones, is_internal, is_staggered, shiftout, + bbs, obs, bbss); bbsss.at(rl) = bbss; obss.at(rl) = obs; diff --git a/Carpet/CarpetRegrid/src/make.code.defn b/Carpet/CarpetRegrid/src/make.code.defn index ec45ed348..0b692a316 100644 --- a/Carpet/CarpetRegrid/src/make.code.defn +++ b/Carpet/CarpetRegrid/src/make.code.defn @@ -1,8 +1,16 @@ # Main make.code.defn file for thorn CarpetRegrid -# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/make.code.defn,v 1.2 2002/05/16 23:25:54 schnetter Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/make.code.defn,v 1.3 2004/01/25 14:57:30 schnetter Exp $ # Source files in this directory -SRCS = regrid.cc paramcheck.cc +SRCS = automatic.cc \ + baselevel.cc \ + centre.cc \ + manualcoordinatelist.cc \ + manualcoordinates.cc \ + manualgridpointlist.cc \ + manualgridpoints.cc \ + regrid.cc \ + paramcheck.cc # Subdirectories containing source files SUBDIRS = diff --git a/Carpet/CarpetRegrid/src/manualcoordinatelist.cc b/Carpet/CarpetRegrid/src/manualcoordinatelist.cc index 10e40ecbd..92ab6900a 100644 --- a/Carpet/CarpetRegrid/src/manualcoordinatelist.cc +++ b/Carpet/CarpetRegrid/src/manualcoordinatelist.cc @@ -1,7 +1,6 @@ -#include <cassert> -#include <cmath> -#include <cstring> -#include <sstream> +#include <assert.h> +#include <string.h> + #include <vector> #include "cctk.h" @@ -13,7 +12,7 @@ #include "regrid.hh" extern "C" { - static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/manualcoordinatelist.cc,v 1.12 2004/08/14 07:42:00 schnetter Exp $"; + static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/manualcoordinatelist.cc,v 1.1 2004/01/25 14:57:30 schnetter Exp $"; CCTK_FILEVERSION(Carpet_CarpetRegrid_manualcoordinatelist_cc); } @@ -28,35 +27,29 @@ namespace CarpetRegrid { int ManualCoordinateList (cGH const * const cctkGH, gh<dim> const & hh, + int const reflevel, + int const map, + int const size, + jjvect const & nboundaryzones, + jjvect const & is_internal, + jjvect const & is_staggered, + jjvect const & shiftout, gh<dim>::rexts & bbsss, gh<dim>::rbnds & obss, gh<dim>::rprocs & pss) { DECLARE_CCTK_PARAMETERS; - int ierr; + + assert (reflevel>=0 && reflevel<maxreflevels); + assert (map>=0 && map<maps); assert (refinement_levels >= 1); // do nothing if the levels already exist - if (reflevel == refinement_levels) return 0; + if (bbsss.size() == refinement_levels) return 0; assert (bbsss.size() >= 1); - jjvect nboundaryzones, is_internal, is_staggered, shiftout; - ierr = GetBoundarySpecification - (2*dim, &nboundaryzones[0][0], &is_internal[0][0], - &is_staggered[0][0], &shiftout[0][0]); - assert (!ierr); - rvect physical_min, physical_max; - rvect interior_min, interior_max; - rvect exterior_min, exterior_max; - rvect base_spacing; - ierr = GetDomainSpecification - (dim, &physical_min[0], &physical_max[0], - &interior_min[0], &interior_max[0], - &exterior_min[0], &exterior_max[0], &base_spacing[0]); - assert (!ierr); - bbsss.resize (refinement_levels); obss.resize (refinement_levels); pss.resize (refinement_levels); @@ -70,79 +63,39 @@ namespace CarpetRegrid { CCTK_WARN (0, "Could not parse parameter \"coordinates\""); } } - + vector<vector<bbvect> > newobss; - if (smart_outer_boundaries) { - // TODO: - // assert (domain_from_coordbase); - + if (strcmp(outerbounds, "") !=0 ) { + istringstream ob_str (outerbounds); + try { + ob_str >> newobss; + } catch (input_error) { + CCTK_WARN (0, "Could not parse parameter \"outerbounds\""); + } + bool good = newobss.size() == newbbss.size(); + if (good) { + for (size_t rl=0; rl<newobss.size(); ++rl) { + good = good && newobss.at(rl).size() == newbbss.at(rl).size(); + } + } + if (! good) { + cout << "coordinates: " << newbbss << endl; + cout << "outerbounds: " << newobss << endl; + CCTK_WARN (0, "The parameters \"outerbounds\" and \"coordinates\" must have the same structure"); + } + } else { newobss.resize(newbbss.size()); for (size_t rl=0; rl<newobss.size(); ++rl) { newobss.at(rl).resize(newbbss.at(rl).size()); for (size_t c=0; c<newobss.at(rl).size(); ++c) { - for (int d=0; d<dim; ++d) { - assert (mglevel==0); - rvect const spacing = base_spacing * pow(CCTK_REAL(mgfact), basemglevel) / ipow(reffact, rl+1); - ierr = ConvertFromPhysicalBoundary - (dim, &physical_min[0], &physical_max[0], - &interior_min[0], &interior_max[0], - &exterior_min[0], &exterior_max[0], &spacing[0]); - assert (!ierr); - newobss.at(rl).at(c)[d][0] = abs(newbbss.at(rl).at(c).lower()[d] - physical_min[d]) < 1.0e-6 * spacing[d]; - if (newobss.at(rl).at(c)[d][0]) { - rvect lo = newbbss.at(rl).at(c).lower(); - rvect up = newbbss.at(rl).at(c).upper(); - rvect str = newbbss.at(rl).at(c).stride(); - lo[d] = exterior_min[d]; - newbbss.at(rl).at(c) = rbbox(lo, up, str); - } - newobss.at(rl).at(c)[d][1] = abs(newbbss.at(rl).at(c).upper()[d] - physical_max[d]) < 1.0e-6 * base_spacing[d] / ipow(reffact, rl); - if (newobss.at(rl).at(c)[d][1]) { - rvect lo = newbbss.at(rl).at(c).lower(); - rvect up = newbbss.at(rl).at(c).upper(); - rvect str = newbbss.at(rl).at(c).stride(); - up[d] = exterior_max[d]; - newbbss.at(rl).at(c) = rbbox(lo, up, str); - } - } + newobss.at(rl).at(c) = bbvect(false); } } - - } else { // if ! smart_outer_boundaries - - if (strcmp(outerbounds, "") !=0 ) { - istringstream ob_str (outerbounds); - try { - ob_str >> newobss; - } catch (input_error) { - CCTK_WARN (0, "Could not parse parameter \"outerbounds\""); - } - bool good = newobss.size() == newbbss.size(); - if (good) { - for (size_t rl=0; rl<newobss.size(); ++rl) { - good = good && newobss.at(rl).size() == newbbss.at(rl).size(); - } - } - if (! good) { - cout << "coordinates: " << newbbss << endl; - cout << "outerbounds: " << newobss << endl; - CCTK_WARN (0, "The parameters \"outerbounds\" and \"coordinates\" must have the same structure"); - } - } else { - newobss.resize(newbbss.size()); - for (size_t rl=0; rl<newobss.size(); ++rl) { - newobss.at(rl).resize(newbbss.at(rl).size()); - for (size_t c=0; c<newobss.at(rl).size(); ++c) { - newobss.at(rl).at(c) = bbvect(false); - } - } - } - - } // if ! smart_outer_boundaries + } if (newbbss.size() < refinement_levels-1) { CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, - "The parameter \"coordinates\" must contain at least \"refinement_levels-1\" (here: %d) levels", int(refinement_levels-1)); + "The parameter \"coordinates\" must contain at least \"refinement_levels-1\" (here: %d) levels", (int)refinement_levels-1); } for (size_t rl=1; rl<refinement_levels; ++rl) { @@ -156,10 +109,6 @@ namespace CarpetRegrid { for (size_t c=0; c<newbbss.at(rl-1).size(); ++c) { rbbox const & ext = newbbss.at(rl-1).at(c); bbvect const & ob = newobss.at(rl-1).at(c); - // TODO: - // assert (domain_from_coordbase); - rvect const spacing = base_spacing * pow(CCTK_REAL(mgfact), basemglevel) / ipow(reffact, rl); - assert (all(abs(ext.stride() - spacing) < spacing * 1.0e-10)); ManualCoordinates_OneLevel (cctkGH, hh, rl, refinement_levels, ext.lower(), ext.upper(), ob, bbs, obs); @@ -171,7 +120,10 @@ namespace CarpetRegrid { // make multigrid aware vector<vector<ibbox> > bbss; - MakeMultigridBoxes (cctkGH, bbs, obs, bbss); + MakeMultigridBoxes + (cctkGH, + size, nboundaryzones, is_internal, is_staggered, shiftout, + bbs, obs, bbss); bbsss.at(rl) = bbss; obss.at(rl) = obs; diff --git a/Carpet/CarpetRegrid/src/manualcoordinates.cc b/Carpet/CarpetRegrid/src/manualcoordinates.cc index 91018cf95..1e3719f18 100644 --- a/Carpet/CarpetRegrid/src/manualcoordinates.cc +++ b/Carpet/CarpetRegrid/src/manualcoordinates.cc @@ -11,7 +11,7 @@ #include "regrid.hh" extern "C" { - static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/manualcoordinates.cc,v 1.5 2004/04/28 15:45:25 schnetter Exp $"; + static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/manualcoordinates.cc,v 1.1 2004/01/25 14:57:30 schnetter Exp $"; CCTK_FILEVERSION(Carpet_CarpetRegrid_manualcoordinates_cc); } @@ -26,19 +26,29 @@ namespace CarpetRegrid { int ManualCoordinates (cGH const * const cctkGH, gh<dim> const & hh, + int const reflevel, + int const map, + int const size, + jjvect const & nboundaryzones, + jjvect const & is_internal, + jjvect const & is_staggered, + jjvect const & shiftout, gh<dim>::rexts & bbsss, gh<dim>::rbnds & obss, gh<dim>::rprocs & pss) { DECLARE_CCTK_PARAMETERS; + assert (reflevel>=0 && reflevel<maxreflevels); + assert (map>=0 && map<maps); + if (refinement_levels > 4) { CCTK_WARN (0, "Cannot currently specify manual refinement regions for more than 4 refinement levels"); } assert (refinement_levels >= 1 && refinement_levels <= 4); // do nothing if the levels already exist - if (reflevel == refinement_levels) return 0; + if (bbsss.size() == refinement_levels) return 0; assert (bbsss.size() >= 1); @@ -54,8 +64,6 @@ namespace CarpetRegrid { lower.at(2) = rvect (l3xmin, l3ymin, l3zmin); upper.at(2) = rvect (l3xmax, l3ymax, l3zmax); - assert (! smart_outer_boundaries); - for (size_t rl=1; rl<bbsss.size(); ++rl) { bbvect const ob (false); @@ -73,7 +81,10 @@ namespace CarpetRegrid { // make multigrid aware vector<vector<ibbox> > bbss; - MakeMultigridBoxes (cctkGH, bbs, obs, bbss); + MakeMultigridBoxes + (cctkGH, + size, nboundaryzones, is_internal, is_staggered, shiftout, + bbs, obs, bbss); bbsss.at(rl) = bbss; obss.at(rl) = obs; @@ -88,27 +99,27 @@ namespace CarpetRegrid { void ManualCoordinates_OneLevel (const cGH * const cctkGH, const gh<dim> & hh, - const int rl, - const int numrl, + const int reflevel, + const int reflevels, const rvect lower, const rvect upper, const bbvect obound, vector<ibbox> & bbs, vector<bbvect> & obs) { - if (rl >= numrl) return; + if (reflevel >= reflevels) return; - jvect const ilower = pos2int (cctkGH, hh, lower, rl); - jvect const iupper = pos2int (cctkGH, hh, upper, rl); + jvect const ilower = pos2int (cctkGH, hh, lower, reflevel); + jvect const iupper = pos2int (cctkGH, hh, upper, reflevel); ManualGridpoints_OneLevel - (cctkGH, hh, rl, numrl, ilower, iupper, obound, bbs, obs); + (cctkGH, hh, reflevel, reflevels, ilower, iupper, obound, bbs, obs); } ivect delta2int (const cGH * const cctkGH, const gh<dim>& hh, - const rvect & rpos, const int rl) + const rvect & rpos, const int reflevel) { rvect global_lower, global_upper; for (int d=0; d<dim; ++d) { @@ -121,13 +132,15 @@ namespace CarpetRegrid { } 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); + const int levfac = ipow(hh.reffact, reflevel); assert (all (hh.baseextent.stride() % levfac == 0)); const ivect istride = hh.baseextent.stride() / levfac; const ivect ipos - = ivect(floor(rpos * scale / rvect(istride) + 0.5)) * istride; + = 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)); @@ -138,7 +151,7 @@ namespace CarpetRegrid { ivect pos2int (const cGH* const cctkGH, const gh<dim>& hh, - const rvect & rpos, const int rl) + const rvect & rpos, const int reflevel) { rvect global_lower, global_upper; for (int d=0; d<dim; ++d) { @@ -151,14 +164,16 @@ namespace CarpetRegrid { } 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); + const int levfac = ipow(hh.reffact, reflevel); assert (all (hh.baseextent.stride() % levfac == 0)); const ivect istride = hh.baseextent.stride() / levfac; const ivect ipos - = (ivect(floor((rpos - global_lower) * scale / rvect(istride) + 0.5)) - * istride); + = ivect(::map(rfloor, (rpos - global_lower) * scale / rvect(istride) + + 0.5)) * istride; return ipos; } diff --git a/Carpet/CarpetRegrid/src/manualgridpointlist.cc b/Carpet/CarpetRegrid/src/manualgridpointlist.cc index e53866e6c..6976b982b 100644 --- a/Carpet/CarpetRegrid/src/manualgridpointlist.cc +++ b/Carpet/CarpetRegrid/src/manualgridpointlist.cc @@ -1,7 +1,6 @@ #include <assert.h> #include <string.h> -#include <sstream> #include <vector> #include "cctk.h" @@ -13,7 +12,7 @@ #include "regrid.hh" extern "C" { - static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/manualgridpointlist.cc,v 1.4 2004/07/02 10:14:51 tradke Exp $"; + static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/manualgridpointlist.cc,v 1.1 2004/01/25 14:57:30 schnetter Exp $"; CCTK_FILEVERSION(Carpet_CarpetRegrid_manualgridpointlist_cc); } @@ -28,16 +27,26 @@ namespace CarpetRegrid { int ManualGridpointList (cGH const * const cctkGH, gh<dim> const & hh, + int const reflevel, + int const map, + int const size, + jjvect const & nboundaryzones, + jjvect const & is_internal, + jjvect const & is_staggered, + jjvect const & shiftout, gh<dim>::rexts & bbsss, gh<dim>::rbnds & obss, gh<dim>::rprocs & pss) { DECLARE_CCTK_PARAMETERS; + assert (reflevel>=0 && reflevel<maxreflevels); + assert (map>=0 && map<maps); + assert (refinement_levels >= 1); // do nothing if the levels already exist - if (reflevel == refinement_levels) return 0; + if (bbsss.size() == refinement_levels) return 0; assert (bbsss.size() >= 1); @@ -111,7 +120,10 @@ namespace CarpetRegrid { // make multigrid aware vector<vector<ibbox> > bbss; - MakeMultigridBoxes (cctkGH, bbs, obs, bbss); + MakeMultigridBoxes + (cctkGH, + size, nboundaryzones, is_internal, is_staggered, shiftout, + bbs, obs, bbss); bbsss.at(rl) = bbss; obss.at(rl) = obs; diff --git a/Carpet/CarpetRegrid/src/manualgridpoints.cc b/Carpet/CarpetRegrid/src/manualgridpoints.cc index 2a8bedb5e..e8cd6c17a 100644 --- a/Carpet/CarpetRegrid/src/manualgridpoints.cc +++ b/Carpet/CarpetRegrid/src/manualgridpoints.cc @@ -1,6 +1,5 @@ #include <assert.h> -#include <sstream> #include <vector> #include "cctk.h" @@ -12,7 +11,7 @@ #include "regrid.hh" extern "C" { - static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/manualgridpoints.cc,v 1.5 2004/07/02 10:14:51 tradke Exp $"; + static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/manualgridpoints.cc,v 1.1 2004/01/25 14:57:30 schnetter Exp $"; CCTK_FILEVERSION(Carpet_CarpetRegrid_manualgridpoints_cc); } @@ -27,19 +26,29 @@ namespace CarpetRegrid { int ManualGridpoints (cGH const * const cctkGH, gh<dim> const & hh, + int const reflevel, + int const map, + int const size, + jjvect const & nboundaryzones, + jjvect const & is_internal, + jjvect const & is_staggered, + jjvect const & shiftout, gh<dim>::rexts & bbsss, gh<dim>::rbnds & obss, gh<dim>::rprocs & pss) { DECLARE_CCTK_PARAMETERS; + assert (reflevel>=0 && reflevel<maxreflevels); + assert (map>=0 && map<maps); + if (refinement_levels > 4) { CCTK_WARN (0, "Cannot currently specify manual refinement regions for more than 4 refinement levels"); } assert (refinement_levels >= 1 && refinement_levels <= 4); // do nothing if the levels already exist - if (reflevel == refinement_levels) return 0; + if (bbsss.size() == refinement_levels) return 0; assert (bbsss.size() >= 1); @@ -55,8 +64,6 @@ namespace CarpetRegrid { ilower.at(2) = ivect (l3ixmin, l3iymin, l3izmin); iupper.at(2) = ivect (l3ixmax, l3iymax, l3izmax); - assert (! smart_outer_boundaries); - for (size_t rl=1; rl<bbsss.size(); ++rl) { bbvect const ob (false); @@ -65,7 +72,7 @@ namespace CarpetRegrid { gh<dim>::cbnds obs; ManualGridpoints_OneLevel - (cctkGH, hh, rl,refinement_levels, + (cctkGH, hh, rl, refinement_levels, ilower.at(rl-1), iupper.at(rl-1), ob, bbs, obs); // make multiprocessor aware @@ -74,7 +81,10 @@ namespace CarpetRegrid { // make multigrid aware vector<vector<ibbox> > bbss; - MakeMultigridBoxes (cctkGH, bbs, obs, bbss); + MakeMultigridBoxes + (cctkGH, + size, nboundaryzones, is_internal, is_staggered, shiftout, + bbs, obs, bbss); bbsss.at(rl) = bbss; obss.at(rl) = obs; @@ -89,36 +99,38 @@ namespace CarpetRegrid { void ManualGridpoints_OneLevel (const cGH * const cctkGH, const gh<dim> & hh, - const int rl, - const int numrl, + const int reflevel, + const int reflevels, const ivect ilower, const ivect iupper, const bbvect obound, vector<ibbox> & bbs, vector<bbvect> & obs) { + if (reflevel >= reflevels) return; + const ivect rstr = hh.baseextent.stride(); const ivect rlb = hh.baseextent.lower(); const ivect rub = hh.baseextent.upper(); - const int levfac = ipow(hh.reffact, rl); + const int levfac = ipow(hh.reffact, reflevel); 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; + buf << "The refinement region boundaries for refinement level #" << reflevel << " 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)) { ostringstream buf; - buf << "The refinement region boundaries for refinement level #" << rl << " have the upper boundary (" << ub << ") less than the lower boundary (" << lb << ")" << ends; + buf << "The refinement region boundaries for refinement level #" << reflevel << " have the upper boundary (" << ub << ") less than the lower boundary (" << lb << ")" << ends; CCTK_WARN (0, buf.str().c_str()); } 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); + "The refinement region boundaries for refinement level #%d are not a multiple of the stride for that level", reflevel); } assert (all(lb>=rlb && ub<=rub)); assert (all(lb<=ub)); 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 diff --git a/Carpet/CarpetRegrid/src/regrid.hh b/Carpet/CarpetRegrid/src/regrid.hh index 45a390a97..27edc99c9 100644 --- a/Carpet/CarpetRegrid/src/regrid.hh +++ b/Carpet/CarpetRegrid/src/regrid.hh @@ -1,4 +1,4 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/regrid.hh,v 1.10 2003/11/13 16:04:37 schnetter Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/regrid.hh,v 1.11 2004/01/25 14:57:30 schnetter Exp $ #ifndef CARPETREGRID_HH #define CARPETREGRID_HH @@ -6,14 +6,15 @@ #include <list> #include "cctk.h" +#include "cctk_Arguments.h" +#include "bbox.hh" #include "gf.hh" #include "gh.hh" +#include "vect.hh" #include "carpet.hh" -#include "regrid.h" - namespace CarpetRegrid { @@ -23,53 +24,171 @@ namespace CarpetRegrid { - typedef vect<int,dim> ivect; - typedef bbox<int,dim> ibbox; - - typedef vect<vect<bool,2>,dim> bvect; - - typedef vect<CCTK_REAL,dim> rvect; - typedef bbox<CCTK_REAL,dim> rbbox; - - - - int CarpetRegridRegrid (const cGH * const cctkGH, - gh<dim>::rexts& bbsss, - gh<dim>::rbnds& obss, - gh<dim>::rprocs& pss); - - - - void MakeRegions_BaseLevel (const cGH* cctkGH, - list<ibbox>& bbl, list<bvect>& obl); - - void MakeRegions_RefineCentre (const cGH* cctkGH, const int reflevels, - list<ibbox>& bbl, list<bvect>& obl); - - void MakeRegions_AsSpecified (const cGH* cctkGH, const int reflevels, - const vector<ivect> lower, - const vector<ivect> upper, - list<ibbox>& bbl, list<bvect>& obl); - void MakeRegions_AsSpecified (const cGH* cctkGH, const int reflevels, - const vector<rvect> lower, - const vector<rvect> upper, - list<ibbox>& bbl, list<bvect>& 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); - 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); - - void MakeRegions_Adaptively (const cGH* cctkGH, - const int minwidth, - const CCTK_REAL minfraction, - const CCTK_REAL maxerror, - const gf<CCTK_REAL,dim>& error, - list<ibbox>& bbl, list<bvect>& obl); + extern "C" { + + /* Scheduled functions */ + int CarpetRegridParamcheck (CCTK_ARGUMENTS); + + /* Aliased functions */ +// CCTK_INT CarpetRegrid_Regrid (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_); + } + + + + int BaseLevel (cGH const * const cctkGH, + gh<dim> const & hh, + int const reflevel, + int const map, + int const size, + jjvect const & nboundaryzones, + jjvect const & is_internal, + jjvect const & is_staggered, + jjvect const & shiftout, + gh<dim>::rexts & bbsss, + gh<dim>::rbnds & obss, + gh<dim>::rprocs & pss); + + int Centre (cGH const * const cctkGH, + gh<dim> const & hh, + int const reflevel, + int const map, + int const size, + jjvect const & nboundaryzones, + jjvect const & is_internal, + jjvect const & is_staggered, + jjvect const & shiftout, + gh<dim>::rexts & bbsss, + gh<dim>::rbnds & obss, + gh<dim>::rprocs & pss); + + int ManualGridpoints (cGH const * const cctkGH, + gh<dim> const & hh, + int const reflevel, + int const map, + int const size, + jjvect const & nboundaryzones, + jjvect const & is_internal, + jjvect const & is_staggered, + jjvect const & shiftout, + gh<dim>::rexts & bbsss, + gh<dim>::rbnds & obss, + gh<dim>::rprocs & pss); + + void ManualGridpoints_OneLevel (const cGH * const cctkGH, + const gh<dim> & hh, + const int reflevel, + const int reflevels, + const ivect ilower, + const ivect iupper, + const bbvect obound, + vector<ibbox> & bbs, + vector<bbvect> & obs); + + int ManualCoordinates (cGH const * const cctkGH, + gh<dim> const & hh, + int const reflevel, + int const map, + int const size, + jjvect const & nboundaryzones, + jjvect const & is_internal, + jjvect const & is_staggered, + jjvect const & shiftout, + gh<dim>::rexts & bbsss, + gh<dim>::rbnds & obss, + gh<dim>::rprocs & pss); + + void ManualCoordinates_OneLevel (const cGH * const cctkGH, + const gh<dim> & hh, + const int reflevel, + const int reflevels, + const rvect ilower, + const rvect iupper, + const bbvect obound, + vector<ibbox> & bbs, + vector<bbvect> & obs); + + ivect delta2int (const cGH * const cctkGH, const gh<dim>& hh, + const rvect & rpos, const int reflevel); + ivect pos2int (const cGH* const cctkGH, const gh<dim>& hh, + const rvect & rpos, const int reflevel); + + int ManualGridpointList (cGH const * const cctkGH, + gh<dim> const & hh, + int const reflevel, + int const map, + int const size, + jjvect const & nboundaryzones, + jjvect const & is_internal, + jjvect const & is_staggered, + jjvect const & shiftout, + gh<dim>::rexts & bbsss, + gh<dim>::rbnds & obss, + gh<dim>::rprocs & pss); + + int ManualCoordinateList (cGH const * const cctkGH, + gh<dim> const & hh, + int const reflevel, + int const map, + int const size, + jjvect const & nboundaryzones, + jjvect const & is_internal, + jjvect const & is_staggered, + jjvect const & shiftout, + gh<dim>::rexts & bbsss, + gh<dim>::rbnds & obss, + gh<dim>::rprocs & pss); + + int Automatic (cGH const * const cctkGH, + gh<dim> const & hh, + int const reflevel, + int const map, + int const size, + jjvect const & nboundaryzones, + jjvect const & is_internal, + jjvect const & is_staggered, + jjvect const & shiftout, + gh<dim>::rexts & bbsss, + gh<dim>::rbnds & obss, + gh<dim>::rprocs & pss); + + void Automatic_OneLevel (const cGH * const cctkGH, + const gh<dim> & hh, + const int reflevel, + const int minwidth, + const CCTK_REAL minfraction, + const CCTK_REAL maxerror, + const gf<CCTK_REAL,dim> & errorvar, + vector<ibbox> & bbs, + vector<bbvect> & obs); + + void Automatic_Recursive (const cGH * const cctkGH, + const gh<dim> & hh, + const int minwidth, + const CCTK_REAL minfraction, + const CCTK_REAL maxerror, + const data<CCTK_REAL,dim> & errorvar, + list<ibbox> & bbl, + const ibbox & region); + + void Automatic_Recombine (list<ibbox> & bbl1, + list<ibbox> & bbl2, + list<ibbox> & bbl, + const ibbox & iface, + const int dir); } // namespace CarpetRegrid |