From 2dc54a6101cd22843342978923636dd7713e7bba Mon Sep 17 00:00:00 2001 From: ih Date: Fri, 25 Feb 2005 20:36:00 +0000 Subject: CarpetAdaptiveRegrid: Code documentation As well as documenting the code a bit, add some fixme's (including noting that the code probably won't work correctly with checkpoint/restart yet!). darcs-hash:20050225203643-0ff1f-ae45e424aa69d52e9381a9a011d2242b58944fe7.gz --- CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc | 220 ++++++++++++++++++++++-------- 1 file changed, 161 insertions(+), 59 deletions(-) (limited to 'CarpetDev') diff --git a/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc b/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc index 207e12c00..157d6f8f1 100644 --- a/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc +++ b/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc @@ -21,19 +21,49 @@ extern "C" { CCTK_FILEVERSION(Carpet_CarpetAdaptiveregrid_regrid_cc); } - +// +// For the moment this is going to live as one large file with one +// large routine. However, it should be broken up into pieces later +// when possible. +// namespace CarpetAdaptiveRegrid { + // + // The following "local_" variables contain copies of the list of + // boxes and the outer boundaries. However, these should be + // PROCESSOR INDEPENDENT; that is, these are never used in calls to + // SplitRegions. All box computation is done with respect to these + // local variables, so should always be independent of the number of + // processors. The results are split across processors and put into + // the "standard" variables before being passed back to Carpet. + // static gh::mexts local_bbsss; static gh::rbnds local_obss; + // + // Keep track of the last iteration on which we were called. This + // means that we can return if we've been called on this timestep + // (because we want to ensure proper nesting we regrid all levels + // finer than the one on which we are called. Therefore if we're + // called on the same iteration, as long as we're called in + // coarse->fine order, we only need to regrid once.) + // + // FIXME: Note that this may well cause problems with checkpoint / + // restart; should this be moved to a grid scalar? + // + static CCTK_INT last_iteration = -1; using namespace std; using namespace Carpet; + // + // The following Fortran helper routines are defined in + // CAR_utils.F90. + // + extern "C" { void CCTK_FCALL CCTK_FNAME(copy_mask) (const CCTK_INT& snx, const CCTK_INT& sny, const CCTK_INT& snz, @@ -51,11 +81,23 @@ namespace CarpetAdaptiveRegrid { CCTK_INT& didit); } + // + // The following helper routines translate between real coordinates + // and integer (Carpet bbox) coordinates on a given refinement + // level. Basically copied from CarpetRegrid, "pos" <-> "position" + // and "int" <-> "integer". Actually defined at the bottom of this + // file. + // + ivect pos2int (const cGH* const cctkGH, const gh& hh, const rvect & rpos, const int rl); rvect int2pos (const cGH* const cctkGH, const gh& hh, const ivect & ipos, const int rl); + // + // The real routine that does all the work + // + CCTK_INT CarpetAdaptiveRegrid_Regrid (CCTK_POINTER_TO_CONST const cctkGH_, CCTK_POINTER const bbsss_, CCTK_POINTER const obss_, @@ -65,7 +107,19 @@ namespace CarpetAdaptiveRegrid { DECLARE_CCTK_PARAMETERS; const cGH * const cctkGH = (const cGH *) cctkGH_; - + + // + // The following variables are the important ones that Carpet uses + // to set up the grids. They contain the bounding boxes (bb), the + // outer boundary specifications (ob) and the processor numbers + // (p). I believe the "s" stands for "set"; so the "obss" is a set + // of sets of outer boundary specifications. However, each "s" is + // implemented as a vector, not a set. For the ob and p the + // indices on the vectors are the refinement levels and the maps; + // for the bounding boxes the order is multigrid levels, + // refinement levels, maps. + // + gh::mexts & bbsss = * (gh::mexts *) bbsss_; gh::rbnds & obss = * (gh::rbnds *) obss_; gh::rprocs & pss = * (gh::rprocs *) pss_; @@ -75,9 +129,15 @@ namespace CarpetAdaptiveRegrid { assert (is_singlemap_mode()); if (local_bbsss.empty()) { // It's the first call - // Is this really the right thing to do on - // multiprocessors? - // local_bbsss = bbsss; + // + // We set up the "local_" variables using just the base grid + // which we get from looking at the base extent; of course, the + // outer boundary set up is trivial. + // + // Note that this will need improving (in fact the whole + // regridding mechanism need minor changes) in order to do mesh + // refinement with multiple maps. + // const ibbox& baseext = vdd.at(Carpet::map)->bases.at(mglevel).at(reflevel).exterior; vector tmp_bbs; @@ -91,6 +151,11 @@ namespace CarpetAdaptiveRegrid { MakeMultigridBoxes(cctkGH, tmp_bbss, tmp_obss, local_bbsss); local_obss = tmp_obss; last_iteration = cctkGH->cctk_iteration; + // + // Having set up the base grid we then set any finer grids + // according to the coordinate parameters, as if we were the + // standard CarpetRegrid. + // CCTK_INT do_recompose = ManualCoordinateList (cctkGH, hh, bbsss, obss, pss, local_bbsss, local_obss); @@ -129,16 +194,14 @@ namespace CarpetAdaptiveRegrid { { return 0; } - - // Return if it's initial data as we can't handle that yet. - // Actually don't as initial data should now be handled - // by the manualcoordinatelist above. - // if (cctkGH->cctk_iteration == 0) { - // return 0; - // } } + // + // Even if force is set, there are times that I'm not going to do + // any regridding, so be told. + // + if (reflevel == maxreflevels - 1) return 0; // Return if we want to regrid regularly, but not at this time @@ -156,9 +219,11 @@ namespace CarpetAdaptiveRegrid { last_iteration = cctkGH->cctk_iteration; } -// cout << "bbsss at start" << endl << bbsss << endl; -// cout << "obss at start" << endl << obss << endl; -// cout << "pss at start" << endl << pss << endl; + // + // Set up a few local variables for later use. In particular we + // need to keep track of the level and map on which we were + // called. + // CCTK_INT do_recompose; do_recompose = 1; @@ -172,8 +237,10 @@ namespace CarpetAdaptiveRegrid { CCTK_INT finest_current_rl = local_bbsss.at(0).size(); finest_current_rl = min(finest_current_rl, maxreflevels - 1); + // // Loop over all levels finer than this one. - + // + leave_singlemap_mode(const_cast (cctkGH)); leave_level_mode(const_cast (cctkGH)); for (CCTK_INT rl = called_on_rl; rl < finest_current_rl; ++rl) { @@ -187,26 +254,6 @@ namespace CarpetAdaptiveRegrid { CCTK_INFO(buf.str().c_str()); } - // So the full algorithm should look something like: - - // Find how big the first bounding box should be on this level - // Do this by finding min lower and max upper bounds of all bboxes - // Allocate box - // Fill errors from local arrays - // If grandchildren exist use their bboxes (expanded) to add to errors - // Reduce errors (MPI sum) - // Set errors to 0/1 - // Define vectors of bboxes final (empty) and todo (contains bbox) - // Define vector of masks (contains error mask) - // Loop over all entries in todo: - // Setup appropriate 1d array memory - // Call fortran routine - // If return is: - // zero: add bbox to final - // one: add new bbox to todo and assoc mask to masklist - // two: add both new bboxs to todo and assoc masks to masklist - // - vector bbs = local_bbsss.at(mglevel).at(reflevel); stack final; @@ -216,6 +263,11 @@ namespace CarpetAdaptiveRegrid { bool did_regrid = false; + // + // For later use get the physical domain specification in real + // coordinates. + // + rvect physical_min, physical_max; rvect interior_min, interior_max; rvect exterior_min, exterior_max; @@ -225,6 +277,10 @@ namespace CarpetAdaptiveRegrid { &interior_min[0], &interior_max[0], &exterior_min[0], &exterior_max[0], &base_spacing[0]); assert (!ierr); + + // + // Loop over all boxes on this level. + // for ( vector::const_iterator bbi = bbs.begin(); bbi != bbs.end(); @@ -253,9 +309,13 @@ namespace CarpetAdaptiveRegrid { << " (points: " << prod(bb.shape()/bb.stride()) << ")"; CCTK_INFO(buf.str().c_str()); } - + + // // Setup the mask. - + // That is, for any errors that we can locally see that exceed + // the threshold, set the mask to 1. + // + const ibbox& baseext = vdd.at(Carpet::map)->bases.at(mglevel).at(reflevel).exterior; ivect imin = (bb.lower() - baseext.lower())/bb.stride(), @@ -328,9 +388,11 @@ namespace CarpetAdaptiveRegrid { } } END_LOCAL_COMPONENT_LOOP; - // Instead check the error on child level, if exists + // + // Check the error on child level, if such a level exists // This should fix the "orphaned grandchild" problem - + // + if (local_bbss.size() > reflevel+1) { CCTK_INT currentml = mglevel; @@ -419,8 +481,13 @@ namespace CarpetAdaptiveRegrid { enter_singlemap_mode(const_cast (cctkGH), currentmap); } - - // Reduce errors (MPI sum) + + // + // Reduce errors (MPI sum). + // This sets up the mask globally, although the mask is now + // non-zero where in error instead of just being 1. It's still + // 0 at points not in error. + // CCTK_INT ierr = CCTK_ReduceLocArrayToArray1D (cctkGH, @@ -438,8 +505,11 @@ namespace CarpetAdaptiveRegrid { CCTK_WARN(0, buf.str().c_str()); } - // Set errors to 0/1 - + // + // Set errors to 0/1; so, where the mask indicates that a + // point is in error, set it precisely to 1. + // + for (CCTK_INT k = 0; k < imax[2] - imin[2] + 1; k++) { for (CCTK_INT j = 0; j < imax[1] - imin[1] + 1; j++) { for (CCTK_INT i = 0; i < imax[0] - imin[0] + 1; i++) { @@ -452,8 +522,10 @@ namespace CarpetAdaptiveRegrid { } } + // // Pad the errors: stage 1 - buffer points marked as 2. - + // + for (CCTK_INT k = 0; k < imax[2] - imin[2] + 1; k++) { for (CCTK_INT j = 0; j < imax[1] - imin[1] + 1; j++) { for (CCTK_INT i = 0; i < imax[0] - imin[0] + 1; i++) { @@ -483,8 +555,10 @@ namespace CarpetAdaptiveRegrid { } } } + // // stage 2: all buffer points marked truly in error. // Also mark if there are any errors. + // bool should_regrid = false; for (CCTK_INT k = 0; k < imax[2] - imin[2] + 1; k++) { for (CCTK_INT j = 0; j < imax[1] - imin[1] + 1; j++) { @@ -511,8 +585,16 @@ namespace CarpetAdaptiveRegrid { CCTK_INFO(buf.str().c_str()); } - // Define vectors of bboxes final (empty) and todo (contains bbox) - + // + // For this box on this level we now have the marked + // points. If there are any errors then we should actually + // create the new boxes. Starting from the actual bbox, placed + // in the "todo" stack, we iterate over the "todo" stack + // creating new boxes which are either accepted (and hence + // placed on the "final" stack) or still need work done to + // them (hence placed on the "todo" stack). + // + if (should_regrid) { stack todo; @@ -524,15 +606,18 @@ namespace CarpetAdaptiveRegrid { stack > masklist; masklist.push(mask); - + + // // Loop over all entries in todo: // Setup appropriate 1d array memory // Call fortran routine // If return is: // zero: add bbox to final // one: add new bbox to todo and assoc mask to masklist - // two: add both new bboxs to todo and assoc masks to masklist - + // two: add both new bboxs to todo and assoc masks to + // masklist + // + while (!todo.empty()) { @@ -565,6 +650,11 @@ namespace CarpetAdaptiveRegrid { } CCTK_INT didit; + + // + // This is the actual Fortran routine that does the work + // of the Berger-Rigoutsos algorithm. + // CCTK_FNAME(check_box)(nx, ny, nz, &mask.front(), @@ -577,7 +667,7 @@ namespace CarpetAdaptiveRegrid { min_width, min_fraction, didit); - if (didit == 0) { + if (didit == 0) { // Box was accepted final.push(bb); @@ -588,8 +678,8 @@ namespace CarpetAdaptiveRegrid { CCTK_INFO(buf.str().c_str()); } } - else if (didit == 1) { - + else if (didit == 1) { // Box was replaced by a new single + // box ibbox newbbox1(ivect::ref(&fbbox1[0][0]), ivect::ref(&fbbox1[1][0]), ivect::ref(&fbbox1[2][0])); @@ -615,7 +705,7 @@ namespace CarpetAdaptiveRegrid { CCTK_INFO(buf.str().c_str()); } } - else if (didit == 2) { + else if (didit == 2) { // Box was replaced with two boxes ibbox newbbox1(ivect::ref(&fbbox1[0][0]), ivect::ref(&fbbox1[1][0]), @@ -669,8 +759,9 @@ namespace CarpetAdaptiveRegrid { } // should regrid. } // Loop over boxes on the parent grid. - if (did_regrid) { - + if (did_regrid) { // If we actually did something, reconvert the + // boxes to correct Carpet style, plus correct + // the boundaries. // Fixup the stride vector newbbs; vector obs; @@ -772,12 +863,15 @@ namespace CarpetAdaptiveRegrid { << "obox is:" << endl << ob; CCTK_INFO(buf.str().c_str()); } - + + // // Convert back to integer coordinates // We have to do this on the fine grid to ensure that // it is correct for an outer boundary with odd numbers - // of ghost zones where the bbox does not align with the parent. - + // of ghost zones where the bbox does not align with the + // parent. + // + ilo = pos2int(cctkGH, hh, newbbcoord.lower(), reflevel+1); ihi = pos2int(cctkGH, hh, newbbcoord.upper(), reflevel+1); ivect istr = bb.stride() / reffact; @@ -834,6 +928,14 @@ namespace CarpetAdaptiveRegrid { bbs = newbbs; // Set local bbss + + // + // This mess ensures that the local bbsss is correctly set for + // multigrid levels, then splits over regions (remember that + // the local_bbsss is processor independent), then goes + // through the whole multigrid procedure with the processor + // dependent bbsss. + // if (bbss.size() < reflevel+2) { if (verbose) { -- cgit v1.2.3