From de1e2516157419f47c65ed2b5a130a6d6111563c Mon Sep 17 00:00:00 2001 From: Ian Hawke Date: Fri, 18 Feb 2005 12:04:00 +0000 Subject: Multilevel AMR Basic support for more than 2 refinement levels with AMR. Basic tests done with shock tubes seem to show that it works. Note that proper nesting is not ensured as yet as we don't check for orphaned grandchildren. darcs-hash:20050218120448-58c7f-5d2e6605ac05f1363596c0e61f6026e8fc057ff8.gz --- CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc | 609 ++++++++++++++++-------------- 1 file changed, 325 insertions(+), 284 deletions(-) (limited to 'CarpetDev') diff --git a/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc b/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc index ce59cb533..eb18a5ca6 100644 --- a/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc +++ b/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc @@ -28,6 +28,8 @@ namespace CarpetAdaptiveRegrid { using namespace std; using namespace Carpet; + static gh::mexts local_bbsss; + extern "C" { void CCTK_FCALL CCTK_FNAME(copy_mask) (const CCTK_INT& snx, const CCTK_INT& sny, const CCTK_INT& snz, @@ -67,6 +69,15 @@ namespace CarpetAdaptiveRegrid { gh const & hh = *vhh.at(Carpet::map); 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; + } + + // FIXME: We should check that the local reflevel "agrees" + // with what is passed in. // In force mode (force == true) we do not check the // CarpetAdaptiveregrid parameters @@ -129,44 +140,15 @@ namespace CarpetAdaptiveRegrid { // two: add both new bboxs to todo and assoc masks to masklist // - vector bbs = bbsss.at(mglevel).at(reflevel); + vector bbs = local_bbsss.at(mglevel).at(reflevel); - ivect low = 100000, upp = -1; - - for ( vector::const_iterator bbi = bbs.begin(); - bbi != bbs.end(); - ++bbi) - { - low = min(low, bbi->lower()); - upp = max(upp, bbi->upper()); - } - - // low and upp now define the starting bbox. - - ibbox bb(low, upp, bbs.at(0).stride()); - - if (verbose) { - ostringstream buf; - buf << "Found the local size of the box: " << endl << bb; - CCTK_INFO(buf.str().c_str()); - } - - vector mask(prod(bb.shape()/bb.stride()), 0); - - if (veryverbose) { - ostringstream buf; - buf << "Allocated mask size: " << bb.shape()/bb.stride() - << " (points: " << prod(bb.shape()/bb.stride()) << ")"; - CCTK_INFO(buf.str().c_str()); - } - + stack final; - // Setup the mask. + vector > bbss = local_bbsss.at(0); - const ibbox& baseext = - vdd.at(Carpet::map)->bases.at(mglevel).at(reflevel).exterior; - ivect imin = (bb.lower() - baseext.lower())/bb.stride(), - imax = (bb.upper() - baseext.lower())/bb.stride(); + // ivect low = 100000, upp = -1; + + bool did_regrid = false; rvect physical_min, physical_max; rvect interior_min, interior_max; @@ -178,268 +160,315 @@ namespace CarpetAdaptiveRegrid { &exterior_min[0], &exterior_max[0], &base_spacing[0]); assert (!ierr); - BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) + for ( vector::const_iterator bbi = bbs.begin(); + bbi != bbs.end(); + ++bbi) { - const CCTK_REAL *error_var_ptr = - static_cast(CCTK_VarDataPtr(cctkGH, 0, error_var)); - const CCTK_REAL *x_var_ptr = - static_cast(CCTK_VarDataPtr(cctkGH, 0, "Grid::x")); - const CCTK_REAL *y_var_ptr = - static_cast(CCTK_VarDataPtr(cctkGH, 0, "Grid::y")); - const CCTK_REAL *z_var_ptr = - static_cast(CCTK_VarDataPtr(cctkGH, 0, "Grid::z")); + // low = min(low, bbi->lower()); + // upp = max(upp, bbi->upper()); + // } + + ivect low = bbi->lower(); + ivect upp = bbi->upper(); + + // low and upp now define the starting bbox. + + ibbox bb(low, upp, bbs.at(0).stride()); + + if (verbose) { + ostringstream buf; + buf << "Found the local size of the box: " << endl << bb; + CCTK_INFO(buf.str().c_str()); + } + + vector mask(prod(bb.shape()/bb.stride()), 0); + + if (veryverbose) { + ostringstream buf; + buf << "Allocated mask size: " << bb.shape()/bb.stride() + << " (points: " << prod(bb.shape()/bb.stride()) << ")"; + CCTK_INFO(buf.str().c_str()); + } + + + // Setup the mask. + + const ibbox& baseext = + vdd.at(Carpet::map)->bases.at(mglevel).at(reflevel).exterior; + ivect imin = (bb.lower() - baseext.lower())/bb.stride(), + imax = (bb.upper() - baseext.lower())/bb.stride(); + + BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) + { + const CCTK_REAL *error_var_ptr = + static_cast(CCTK_VarDataPtr(cctkGH, 0, error_var)); + const CCTK_REAL *x_var_ptr = + static_cast(CCTK_VarDataPtr(cctkGH, 0, "Grid::x")); + const CCTK_REAL *y_var_ptr = + static_cast(CCTK_VarDataPtr(cctkGH, 0, "Grid::y")); + const CCTK_REAL *z_var_ptr = + static_cast(CCTK_VarDataPtr(cctkGH, 0, "Grid::z")); - assert(all(imin >= 0)); - assert(all(imax >= 0)); - assert(all(imin < ivect::ref(cctkGH->cctk_lsh))); - assert(all(imax < ivect::ref(cctkGH->cctk_lsh))); - assert(all(imin <= imax)); + assert(all(imin >= 0)); + assert(all(imax >= 0)); + // FIXME: Why should the following assert be true? + // assert(all(imax < ivect::ref(cctkGH->cctk_lsh))); + assert(all(imin <= imax)); - for (CCTK_INT k = 0; k < cctkGH->cctk_lsh[2]; ++k) { - for (CCTK_INT j = 0; j < cctkGH->cctk_lsh[1]; ++j) { - for (CCTK_INT i = 0; i < cctkGH->cctk_lsh[0]; ++i) { - CCTK_INT index = CCTK_GFINDEX3D(cctkGH, i, j, k); - if (abs(error_var_ptr[index]) > max_error) { - CCTK_INT ii = i + cctkGH->cctk_lbnd[0] - imin[0]; - CCTK_INT jj = j + cctkGH->cctk_lbnd[1] - imin[1]; - CCTK_INT kk = k + cctkGH->cctk_lbnd[2] - imin[2]; - assert (ii >= 0); - assert (jj >= 0); - assert (kk >= 0); - assert (ii <= imax[0] - imin[0]); - assert (jj <= imax[1] - imin[1]); - assert (kk <= imax[2] - imin[2]); - CCTK_INT mindex = ii + - (imax[0] - imin[0] + 1)*(jj + (imax[1] - imin[1] + 1) * kk); - mask[mindex] = 1; - if (veryverbose) { - CCTK_VInfo(CCTK_THORNSTRING, "In error at point" - "\n(%g,%g,%g) [%d,%d,%d] [[%d,%d,%d]]", - x_var_ptr[index], - y_var_ptr[index], - z_var_ptr[index], - ii, jj, kk, i,j,k); + for (CCTK_INT k = 0; k < cctkGH->cctk_lsh[2]; ++k) { + for (CCTK_INT j = 0; j < cctkGH->cctk_lsh[1]; ++j) { + for (CCTK_INT i = 0; i < cctkGH->cctk_lsh[0]; ++i) { + CCTK_INT index = CCTK_GFINDEX3D(cctkGH, i, j, k); + if (abs(error_var_ptr[index]) > max_error) { + CCTK_INT ii = i + cctkGH->cctk_lbnd[0] - imin[0]; + CCTK_INT jj = j + cctkGH->cctk_lbnd[1] - imin[1]; + CCTK_INT kk = k + cctkGH->cctk_lbnd[2] - imin[2]; + // Check that this point actually intersects with + // this box (if this component was actually a + // different grid on the same processor, it need not) + if ( (ii >= 0) and (jj >= 0) and (kk >= 0) and + (ii <= imax[0] - imin[0]) and + (jj <= imax[1] - imin[1]) and + (kk <= imax[2] - imin[2]) ) + { + assert (ii >= 0); + assert (jj >= 0); + assert (kk >= 0); + assert (ii <= imax[0] - imin[0]); + assert (jj <= imax[1] - imin[1]); + assert (kk <= imax[2] - imin[2]); + CCTK_INT mindex = ii + + (imax[0] - imin[0] + 1)*(jj + (imax[1] - imin[1] + 1) * kk); + mask[mindex] = 1; + if (veryverbose) { + CCTK_VInfo(CCTK_THORNSTRING, "In error at point" + "\n(%g,%g,%g) [%d,%d,%d] [[%d,%d,%d]]", + x_var_ptr[index], + y_var_ptr[index], + z_var_ptr[index], + ii, jj, kk, i,j,k); + } + } } } } } - } - } END_LOCAL_COMPONENT_LOOP; - - // FIXME - // If grandchildren exist use their bboxes (expanded) to add to errors - // Reduce errors (MPI sum) - // Set errors to 0/1 - - // Pad the errors: stage 1 - buffer points marked as 2. - - for (CCTK_INT k = 0; k < imax[2] - imin[2]; k++) { - for (CCTK_INT j = 0; j < imax[1] - imin[1]; j++) { - for (CCTK_INT i = 0; i < imax[0] - imin[0]; i++) { - CCTK_INT index = i + - (imax[0] - imin[0] + 1)*(j + (imax[1] - imin[1] + 1) * k); - if (mask[index] == 1) { - for (CCTK_INT kk = max(k-pad, 0); - kk < min(k+pad+1, imax[2] - imin[2]); - ++kk) - { - for (CCTK_INT jj = max(j-pad, 0); - jj < min(j+pad+1, imax[1] - imin[1]); - ++jj) + } END_LOCAL_COMPONENT_LOOP; + + // FIXME + // If grandchildren exist use their bboxes (expanded) to add to errors + // Reduce errors (MPI sum) + // Set errors to 0/1 + + // Pad the errors: stage 1 - buffer points marked as 2. + + for (CCTK_INT k = 0; k < imax[2] - imin[2]; k++) { + for (CCTK_INT j = 0; j < imax[1] - imin[1]; j++) { + for (CCTK_INT i = 0; i < imax[0] - imin[0]; i++) { + CCTK_INT index = i + + (imax[0] - imin[0] + 1)*(j + (imax[1] - imin[1] + 1) * k); + if (mask[index] == 1) { + for (CCTK_INT kk = max(k-pad, 0); + kk < min(k+pad+1, imax[2] - imin[2]); + ++kk) { - for (CCTK_INT ii = max(i-pad, 0); - ii < min(i+pad+1, imax[0] - imin[0]); - ++ii) + for (CCTK_INT jj = max(j-pad, 0); + jj < min(j+pad+1, imax[1] - imin[1]); + ++jj) { - CCTK_INT mindex = ii + - (imax[0] - imin[0] + 1)* - (jj + (imax[1] - imin[1] + 1) * kk); - if (!mask[mindex]) mask[mindex] = 2; + for (CCTK_INT ii = max(i-pad, 0); + ii < min(i+pad+1, imax[0] - imin[0]); + ++ii) + { + CCTK_INT mindex = ii + + (imax[0] - imin[0] + 1)* + (jj + (imax[1] - imin[1] + 1) * kk); + if (!mask[mindex]) mask[mindex] = 2; + } } } } } } } - } - // 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]; k++) { - for (CCTK_INT j = 0; j < imax[1] - imin[1]; j++) { - for (CCTK_INT i = 0; i < imax[0] - imin[0]; i++) { - CCTK_INT index = i + - (imax[0]-imin[0] + 1)*(j + (imax[1] - imin[1] + 1) * k); - if (mask[index] > 1) mask[index] = 1; - if ((veryverbose) and (mask[index])) { - CCTK_VInfo(CCTK_THORNSTRING, "Mask set at point" - "\n[%d,%d,%d]", - i,j,k); - } - should_regrid |= (mask[index]); + // 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]; k++) { + for (CCTK_INT j = 0; j < imax[1] - imin[1]; j++) { + for (CCTK_INT i = 0; i < imax[0] - imin[0]; i++) { + CCTK_INT index = i + + (imax[0]-imin[0] + 1)*(j + (imax[1] - imin[1] + 1) * k); + if (mask[index] > 1) mask[index] = 1; + if ((veryverbose) and (mask[index])) { + CCTK_VInfo(CCTK_THORNSTRING, "Mask set at point" + "\n[%d,%d,%d]", + i,j,k); + } + should_regrid |= (mask[index]); + did_regrid |= should_regrid; + } } - } - } - - // Define vectors of bboxes final (empty) and todo (contains bbox) - - stack final; - - vector > bbss = bbsss.at(0); + } + // Define vectors of bboxes final (empty) and todo (contains bbox) - if (should_regrid) { + if (should_regrid) { - stack todo; - - todo.push(bb); - - // Define vector of masks (contains error mask) - - 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 - - while (!todo.empty()) - { - - ibbox bb = todo.top(); todo.pop(); - vector mask = masklist.top(); masklist.pop(); + stack todo; - CCTK_INT nx = bb.shape()[0]/bb.stride()[0]; - CCTK_INT ny = bb.shape()[1]/bb.stride()[1]; - CCTK_INT nz = bb.shape()[2]/bb.stride()[2]; + todo.push(bb); - if (verbose) { - ostringstream buf; - buf << "todo loop. Box: " << endl << bb; - CCTK_INFO(buf.str().c_str()); - } - - vector sum_x(nx, 0); - vector sig_x(nx, 0); - vector sum_y(ny, 0); - vector sig_y(ny, 0); - vector sum_z(nz, 0); - vector sig_z(nz, 0); + // Define vector of masks (contains error mask) - CCTK_INT fbbox[3][3], fbbox1[3][3], fbbox2[3][3]; + stack > masklist; - for (CCTK_INT d = 0; d < 3; ++d) { - fbbox[0][d] = bb.lower()[d]; - fbbox[1][d] = bb.upper()[d]; - fbbox[2][d] = bb.stride()[d]; - } + masklist.push(mask); - CCTK_INT didit; - - CCTK_FNAME(check_box)(nx, ny, nz, - &mask.front(), - &sum_x.front(), &sum_y.front(), &sum_z.front(), - &sig_x.front(), &sig_y.front(), &sig_z.front(), - fbbox, - fbbox1, fbbox2, - min_width, min_fraction, - didit); + // 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 - if (didit == 0) { - - final.push(bb); - - if (verbose) { - ostringstream buf; - buf << "todo loop. Box pushed to final: " - << endl << bb; - CCTK_INFO(buf.str().c_str()); - } - } - else if (didit == 1) { + while (!todo.empty()) + { - ibbox newbbox1(ivect::ref(&fbbox1[0][0]), - ivect::ref(&fbbox1[1][0]), - ivect::ref(&fbbox1[2][0])); - todo.push(newbbox1); + ibbox bb = todo.top(); todo.pop(); + vector mask = masklist.top(); masklist.pop(); - CCTK_INT dnx = newbbox1.shape()[0]/newbbox1.stride()[0]; - CCTK_INT dny = newbbox1.shape()[1]/newbbox1.stride()[1]; - CCTK_INT dnz = newbbox1.shape()[2]/newbbox1.stride()[2]; + CCTK_INT nx = bb.shape()[0]/bb.stride()[0]; + CCTK_INT ny = bb.shape()[1]/bb.stride()[1]; + CCTK_INT nz = bb.shape()[2]/bb.stride()[2]; - vector - newmask1(prod(newbbox1.shape()/newbbox1.stride()), 0); - - CCTK_FNAME(copy_mask)(nx, ny, nz, - &mask.front(), fbbox, - dnx, dny, dnz, - &newmask1.front(), fbbox1); - masklist.push(newmask1); - if (verbose) { ostringstream buf; - buf << "todo loop. New (single) box created: " - << endl << newbbox1; + buf << "todo loop. Box: " << endl << bb; CCTK_INFO(buf.str().c_str()); } - } - else if (didit == 2) { - ibbox newbbox1(ivect::ref(&fbbox1[0][0]), - ivect::ref(&fbbox1[1][0]), - ivect::ref(&fbbox1[2][0])); - todo.push(newbbox1); - ibbox newbbox2(ivect::ref(&fbbox2[0][0]), - ivect::ref(&fbbox2[1][0]), - ivect::ref(&fbbox2[2][0])); - todo.push(newbbox2); + vector sum_x(nx, 0); + vector sig_x(nx, 0); + vector sum_y(ny, 0); + vector sig_y(ny, 0); + vector sum_z(nz, 0); + vector sig_z(nz, 0); - CCTK_INT dnx = newbbox1.shape()[0]/newbbox1.stride()[0]; - CCTK_INT dny = newbbox1.shape()[1]/newbbox1.stride()[1]; - CCTK_INT dnz = newbbox1.shape()[2]/newbbox1.stride()[2]; + CCTK_INT fbbox[3][3], fbbox1[3][3], fbbox2[3][3]; - vector - newmask1(prod(newbbox1.shape()/newbbox1.stride()), 0); + for (CCTK_INT d = 0; d < 3; ++d) { + fbbox[0][d] = bb.lower()[d]; + fbbox[1][d] = bb.upper()[d]; + fbbox[2][d] = bb.stride()[d]; + } - CCTK_FNAME(copy_mask)(nx, ny, nz, - &mask.front(), fbbox, - dnx, dny, dnz, - &newmask1.front(), fbbox1); - masklist.push(newmask1); + CCTK_INT didit; - dnx = newbbox2.shape()[0]/newbbox2.stride()[0]; - dny = newbbox2.shape()[1]/newbbox2.stride()[1]; - dnz = newbbox2.shape()[2]/newbbox2.stride()[2]; + CCTK_FNAME(check_box)(nx, ny, nz, + &mask.front(), + &sum_x.front(), &sum_y.front(), &sum_z.front(), + &sig_x.front(), &sig_y.front(), &sig_z.front(), + fbbox, + fbbox1, fbbox2, + min_width, min_fraction, + didit); - vector - newmask2(prod(newbbox2.shape()/newbbox2.stride()), 0); + if (didit == 0) { - CCTK_FNAME(copy_mask)(nx, ny, nz, - &mask.front(), fbbox, - dnx, dny, dnz, - &newmask2.front(), fbbox2); - masklist.push(newmask2); + final.push(bb); - if (verbose) { - ostringstream buf; - buf << "todo loop. New (double) box created. Box 1: " - << endl << newbbox1 - << " Box 2: " - << endl << newbbox2; - CCTK_INFO(buf.str().c_str()); + if (verbose) { + ostringstream buf; + buf << "todo loop. Box pushed to final: " + << endl << bb; + CCTK_INFO(buf.str().c_str()); + } } - } - else { - CCTK_WARN(0, "The fortran routine must be confused."); - } - - } // loop over todo vector (boxes needing to be done). - + else if (didit == 1) { + + ibbox newbbox1(ivect::ref(&fbbox1[0][0]), + ivect::ref(&fbbox1[1][0]), + ivect::ref(&fbbox1[2][0])); + todo.push(newbbox1); + + CCTK_INT dnx = newbbox1.shape()[0]/newbbox1.stride()[0]; + CCTK_INT dny = newbbox1.shape()[1]/newbbox1.stride()[1]; + CCTK_INT dnz = newbbox1.shape()[2]/newbbox1.stride()[2]; + + vector + newmask1(prod(newbbox1.shape()/newbbox1.stride()), 0); + + CCTK_FNAME(copy_mask)(nx, ny, nz, + &mask.front(), fbbox, + dnx, dny, dnz, + &newmask1.front(), fbbox1); + masklist.push(newmask1); + + if (verbose) { + ostringstream buf; + buf << "todo loop. New (single) box created: " + << endl << newbbox1; + CCTK_INFO(buf.str().c_str()); + } + } + else if (didit == 2) { + + ibbox newbbox1(ivect::ref(&fbbox1[0][0]), + ivect::ref(&fbbox1[1][0]), + ivect::ref(&fbbox1[2][0])); + todo.push(newbbox1); + ibbox newbbox2(ivect::ref(&fbbox2[0][0]), + ivect::ref(&fbbox2[1][0]), + ivect::ref(&fbbox2[2][0])); + todo.push(newbbox2); + + CCTK_INT dnx = newbbox1.shape()[0]/newbbox1.stride()[0]; + CCTK_INT dny = newbbox1.shape()[1]/newbbox1.stride()[1]; + CCTK_INT dnz = newbbox1.shape()[2]/newbbox1.stride()[2]; + + vector + newmask1(prod(newbbox1.shape()/newbbox1.stride()), 0); + + CCTK_FNAME(copy_mask)(nx, ny, nz, + &mask.front(), fbbox, + dnx, dny, dnz, + &newmask1.front(), fbbox1); + masklist.push(newmask1); + + dnx = newbbox2.shape()[0]/newbbox2.stride()[0]; + dny = newbbox2.shape()[1]/newbbox2.stride()[1]; + dnz = newbbox2.shape()[2]/newbbox2.stride()[2]; + + vector + newmask2(prod(newbbox2.shape()/newbbox2.stride()), 0); + + CCTK_FNAME(copy_mask)(nx, ny, nz, + &mask.front(), fbbox, + dnx, dny, dnz, + &newmask2.front(), fbbox2); + masklist.push(newmask2); + + if (verbose) { + ostringstream buf; + buf << "todo loop. New (double) box created. Box 1: " + << endl << newbbox1 + << " Box 2: " + << endl << newbbox2; + CCTK_INFO(buf.str().c_str()); + } + } + else { + CCTK_WARN(0, "The fortran routine must be confused."); + } + + } // loop over todo vector (boxes needing to be done). + } // should regrid. + } // Loop over boxes on the parent grid. + + if (did_regrid) { + // Fixup the stride vector newbbs; vector obs; @@ -450,7 +479,7 @@ namespace CarpetAdaptiveRegrid { rvect lo = int2pos(cctkGH, hh, ilo, reflevel); rvect hi = int2pos(cctkGH, hh, ihi, reflevel); rvect str = base_spacing * - ipow((CCTK_REAL)mgfact, basemglevel) / ipow(reffact, reflevel); + ipow((CCTK_REAL)mgfact, basemglevel) / ipow(reffact, reflevel); rbbox newbbcoord(lo, hi, str); if (veryverbose) { @@ -459,15 +488,15 @@ namespace CarpetAdaptiveRegrid { << endl << newbbcoord; CCTK_INFO(buf.str().c_str()); } - + // FIXME: Set the correct ob here. - + bbvect ob(false); for (int d=0; d reflevel+1) { @@ -577,8 +615,11 @@ namespace CarpetAdaptiveRegrid { } bbss.resize(reflevel+1); obss.resize(reflevel+1); + // Set local bbsss + MakeMultigridBoxes (cctkGH, bbss, obss, local_bbsss); + pss.resize(reflevel+1); - + do_recompose = 1; } -- cgit v1.2.3