From 78b1f2df0b1b45d00ba9cbb172790d5d675348bb Mon Sep 17 00:00:00 2001 From: Ian Hawke Date: Tue, 22 Feb 2005 17:44:00 +0000 Subject: CarpetAdaptiveRegrid: Enforce proper nesting Enforce proper nesting both for parents and children by looping over all finer levels regridding them all at once. This completely destroys efficiency... darcs-hash:20050222174437-58c7f-d280bacafa2fb4c634613381b4c7134733348ec5.gz --- CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc | 1083 ++++++++++++++++------------- 1 file changed, 615 insertions(+), 468 deletions(-) (limited to 'CarpetDev') diff --git a/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc b/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc index 6923d9188..f99d44e84 100644 --- a/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc +++ b/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc @@ -105,11 +105,18 @@ namespace CarpetAdaptiveRegrid { if (cctkGH->cctk_iteration ==0) { return 0; } - - } + + } if (reflevel == maxreflevels - 1) return 0; + // Return if we want to regrid regularly, but not at this time + if (regrid_every > 0 && cctkGH->cctk_iteration != 0 + && (cctkGH->cctk_iteration-1) % regrid_every != 0) + { + return 0; + } + // cout << "bbsss at start" << endl << bbsss << endl; // cout << "obss at start" << endl << obss << endl; // cout << "pss at start" << endl << pss << endl; @@ -117,563 +124,703 @@ namespace CarpetAdaptiveRegrid { CCTK_INT do_recompose; do_recompose = 1; - // We assume that we are called in fine->coarse order - // So that on finer levels regridding has already been done - - // 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); + CCTK_INT called_on_ml = mglevel; + CCTK_INT called_on_rl = reflevel; + CCTK_INT called_on_map = carpetGH.map; - stack final; - - vector > bbss = local_bbsss.at(0); - - // ivect low = 100000, upp = -1; - - bool did_regrid = false; + CCTK_INT finest_current_rl = local_bbsss.at(0).size(); + finest_current_rl = min(finest_current_rl, maxreflevels - 1); - rvect physical_min, physical_max; - rvect interior_min, interior_max; - rvect exterior_min, exterior_max; - rvect base_spacing; - int 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); - - for ( vector::const_iterator bbi = bbs.begin(); - bbi != bbs.end(); - ++bbi) - { - // low = min(low, bbi->lower()); - // upp = max(upp, bbi->upper()); - // } + // Loop over all levels finer than this one. - ivect low = bbi->lower(); - ivect upp = bbi->upper(); - - // low and upp now define the starting bbox. + leave_singlemap_mode(const_cast (cctkGH)); + leave_level_mode(const_cast (cctkGH)); + for (CCTK_INT rl = called_on_rl; rl < finest_current_rl; ++rl) { + enter_level_mode(const_cast (cctkGH), rl); + enter_singlemap_mode(const_cast (cctkGH), called_on_map); - ibbox bb(low, upp, bbs.at(0).stride()); - if (verbose) { ostringstream buf; - buf << "Found the local size of the box: " << endl << bb; + buf << "Entering level " << rl << " map " << called_on_map; CCTK_INFO(buf.str().c_str()); } - vector mask(prod(bb.shape()/bb.stride()), 0); + // So the full algorithm should look something like: - if (veryverbose) { - ostringstream buf; - buf << "Allocated mask size: " << bb.shape()/bb.stride() - << " (points: " << prod(bb.shape()/bb.stride()) << ")"; - CCTK_INFO(buf.str().c_str()); - } + // 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 + // - - // Setup the mask. + vector bbs = local_bbsss.at(mglevel).at(reflevel); - 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)); - // FIXME: Why should the following assert be true? - // assert(all(imax < ivect::ref(cctkGH->cctk_lsh))); - assert(all(imin <= imax)); + stack final; - 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); - CCTK_REAL local_error = abs(error_var_ptr[index]); - if (local_error > 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; + vector > bbss = local_bbsss.at(0); - // 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. + bool did_regrid = false; + + rvect physical_min, physical_max; + rvect interior_min, interior_max; + rvect exterior_min, exterior_max; + rvect base_spacing; + int 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); + + for ( vector::const_iterator bbi = bbs.begin(); + bbi != bbs.end(); + ++bbi) + { + + 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")); + - 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++) { - 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] + 1); - ++kk) - { - for (CCTK_INT jj = max(j - pad, 0); - jj < min(j + pad + 1, imax[1] - imin[1] + 1); - ++jj) - { - for (CCTK_INT ii = max(i - pad, 0); - ii < min(i + pad + 1, imax[0] - imin[0] + 1); - ++ii) + // These can actually be negative if the parent shrinks. + // Of course, the final grid should still be properly + // nested... + // 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); + CCTK_REAL local_error = abs(error_var_ptr[index]); + if (local_error > 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); - if (!mask[mindex]) mask[mindex] = 2; + 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); + } } } } } } - } - } - // 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++) { - for (CCTK_INT i = 0; i < imax[0] - imin[0] + 1; 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; + } END_LOCAL_COMPONENT_LOOP; + + // Instead check the error on child level, if exists + // This should fix the "orphaned grandchild" problem + + if (bbss.size() > reflevel+1) { + + CCTK_INT currentml = mglevel; + CCTK_INT currentrl = reflevel; + CCTK_INT currentmap = carpetGH.map; + + leave_singlemap_mode(const_cast (cctkGH)); + leave_level_mode(const_cast (cctkGH)); + + enter_level_mode(const_cast (cctkGH), currentrl + 1); + enter_singlemap_mode(const_cast (cctkGH), currentmap); + + if (verbose) { + ostringstream buf; + buf << "Checking for errors on child level " + << reflevel << " map " << currentmap; + CCTK_INFO(buf.str().c_str()); } - } - } - // Define vectors of bboxes final (empty) and todo (contains bbox) + 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)); + // 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); + CCTK_REAL local_error = abs(error_var_ptr[index]); + if (local_error > 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; + + leave_singlemap_mode(const_cast (cctkGH)); + leave_level_mode(const_cast (cctkGH)); + + enter_level_mode(const_cast (cctkGH), currentrl); + enter_singlemap_mode(const_cast (cctkGH), currentmap); - if (should_regrid) { - - stack todo; + } + + // FIXME - todo.push(bb); + // Reduce errors (MPI sum) + // Set errors to 0/1 - // Define vector of masks (contains error mask) + // Pad the errors: stage 1 - buffer points marked as 2. - stack > masklist; + 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++) { + 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] + 1); + ++kk) + { + for (CCTK_INT jj = max(j - pad, 0); + jj < min(j + pad + 1, imax[1] - imin[1] + 1); + ++jj) + { + for (CCTK_INT ii = max(i - pad, 0); + ii < min(i + pad + 1, imax[0] - imin[0] + 1); + ++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] + 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++) { + 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; + } + } + } - masklist.push(mask); + if (verbose) { + ostringstream buf; + buf << "Finished looking for errors on level " + << reflevel << endl << "should_regrid " << should_regrid + << " did_regrid " << did_regrid; + CCTK_INFO(buf.str().c_str()); + } - // 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 + // Define vectors of bboxes final (empty) and todo (contains bbox) - while (!todo.empty()) - { - - ibbox bb = todo.top(); todo.pop(); - vector mask = masklist.top(); masklist.pop(); - - 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]; - - if (verbose) { - ostringstream buf; - buf << "todo loop. Box: " << endl << bb; - CCTK_INFO(buf.str().c_str()); - } + if (should_regrid) { - 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); + stack todo; - CCTK_INT fbbox[3][3], fbbox1[3][3], fbbox2[3][3]; + todo.push(bb); - 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]; - } + // Define vector of masks (contains error mask) - CCTK_INT didit; + stack > masklist; - 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); + masklist.push(mask); - if (didit == 0) { + // 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 - final.push(bb); - + while (!todo.empty()) + { + + ibbox bb = todo.top(); todo.pop(); + vector mask = masklist.top(); masklist.pop(); + + 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]; + if (verbose) { ostringstream buf; - buf << "todo loop. Box pushed to final: " - << endl << bb; + buf << "todo loop. Box: " << endl << bb; CCTK_INFO(buf.str().c_str()); } - } - else if (didit == 1) { - ibbox newbbox1(ivect::ref(&fbbox1[0][0]), - ivect::ref(&fbbox1[1][0]), - ivect::ref(&fbbox1[2][0])); - todo.push(newbbox1); + 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); - - CCTK_FNAME(copy_mask)(nx, ny, nz, - &mask.front(), fbbox, - dnx, dny, dnz, - &newmask1.front(), fbbox1); - masklist.push(newmask1); + 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_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); - if (verbose) { - ostringstream buf; - buf << "todo loop. New (single) box created: " - << endl << newbbox1; - CCTK_INFO(buf.str().c_str()); + 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) { + + 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."); } - } - 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); + } // 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; + while (! final.empty()) { + ibbox bb = final.top(); final.pop(); + + if (veryverbose) { + ostringstream buf; + buf << "Looping over the final list. Box is:" + << endl << bb; + CCTK_INFO(buf.str().c_str()); + } + + ivect ilo = bb.lower(); + ivect ihi = bb.upper(); + 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); + rbbox newbbcoord(lo, hi, str); + + if (veryverbose) { + ostringstream buf; + buf << "Dealing with boundaries. Coord box is:" + << endl << newbbcoord; + CCTK_INFO(buf.str().c_str()); + } + + // Set the correct ob here. + + bbvect ob(false); + for (int d=0; d - newmask1(prod(newbbox1.shape()/newbbox1.stride()), 0); + rvect const spacing = base_spacing * + ipow((CCTK_REAL)mgfact, basemglevel) / + ipow(reffact, reflevel+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); - CCTK_FNAME(copy_mask)(nx, ny, nz, - &mask.front(), fbbox, - dnx, dny, dnz, - &newmask1.front(), fbbox1); - masklist.push(newmask1); + // If need be clip the domain - dnx = newbbox2.shape()[0]/newbbox2.stride()[0]; - dny = newbbox2.shape()[1]/newbbox2.stride()[1]; - dnz = newbbox2.shape()[2]/newbbox2.stride()[2]; + rvect lo = newbbcoord.lower(); + if (newbbcoord.lower()[d] < physical_min[d]) { + lo[d] = exterior_min[d]; + } + rvect up = newbbcoord.upper(); + if (newbbcoord.upper()[d] > physical_max[d]) { + up[d] = exterior_max[d]; + } + rvect str = newbbcoord.stride(); - vector - newmask2(prod(newbbox2.shape()/newbbox2.stride()), 0); + // Set the ob if outside the physical domain - CCTK_FNAME(copy_mask)(nx, ny, nz, - &mask.front(), fbbox, - dnx, dny, dnz, - &newmask2.front(), fbbox2); - masklist.push(newmask2); + ob[d][0] = + abs(lo[d] - exterior_min[d]) < 1.0e-6 * spacing[d]; + ob[d][1] = + abs(up[d] - exterior_max[d]) < 1.0e-6 * spacing[d]; - if (verbose) { + if (veryverbose) { ostringstream buf; - buf << "todo loop. New (double) box created. Box 1: " - << endl << newbbox1 - << " Box 2: " - << endl << newbbox2; + buf << "Done clipping domain:" + << endl << lo << endl << up << endl << str; CCTK_INFO(buf.str().c_str()); + } + + // Check that the striding is correct. + + rvect remainder = (up - lo) / str - floor( (up - lo) / str ); + + for (int d=0; d < dim; ++d) { + if ( abs(remainder[d]) > 1.e-6 ) { + if (ob[d][0]) { + up[d] += str[d] * (1 - remainder[d]); + } + else if (ob[d][1]) { + lo[d] -= str[d] * remainder[d]; + } + } } + + if (veryverbose) { + ostringstream buf; + buf << "Corrected coords for striding:" + << endl << lo << endl << up << endl << str; + CCTK_INFO(buf.str().c_str()); + } + + newbbcoord = rbbox(lo, up, 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; - while (! final.empty()) { - ibbox bb = final.top(); final.pop(); - - if (veryverbose) { - ostringstream buf; - buf << "Looping over the final list. Box is:" - << endl << bb; - CCTK_INFO(buf.str().c_str()); - } - - ivect ilo = bb.lower(); - ivect ihi = bb.upper(); - 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); - rbbox newbbcoord(lo, hi, str); - - if (veryverbose) { - ostringstream buf; - buf << "Dealing with boundaries. Coord box is:" - << endl << newbbcoord; - CCTK_INFO(buf.str().c_str()); - } - - // FIXME: Set the correct ob here. - - bbvect ob(false); - for (int d=0; d physical_max[d]) { - up[d] = exterior_max[d]; - } - rvect str = newbbcoord.stride(); - - // Set the ob if outside the physical domain - - ob[d][0] = - abs(lo[d] - exterior_min[d]) < 1.0e-6 * spacing[d]; - ob[d][1] = - abs(up[d] - exterior_max[d]) < 1.0e-6 * spacing[d]; - - if (veryverbose) { + if (verbose) { ostringstream buf; - buf << "Done clipping domain:" - << endl << lo << endl << up << endl << str; + buf << "Done dealing with boundaries. Coord box is:" + << endl << newbbcoord << endl + << "obox is:" << endl << ob; CCTK_INFO(buf.str().c_str()); - } + } - // Check that the striding is correct. + // 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. - rvect remainder = (up - lo) / str - floor( (up - lo) / str ); - + ilo = pos2int(cctkGH, hh, newbbcoord.lower(), reflevel+1); + ihi = pos2int(cctkGH, hh, newbbcoord.upper(), reflevel+1); + ivect istr = bb.stride() / reffact; + + // Check that the width is sufficient + // This can only be too small if the domain was clipped for (int d=0; d < dim; ++d) { - if ( abs(remainder[d]) > 1.e-6 ) { + if (ihi[d] - ilo[d] < min_width * istr[d]) { if (ob[d][0]) { - up[d] += str[d] * (1 - remainder[d]); + if (ob[d][1]) { + CCTK_WARN(0, "The domain is too small?!"); + } + ihi[d] = ilo[d] + min_width * istr[d]; } else if (ob[d][1]) { - lo[d] -= str[d] * remainder[d]; + if (ob[d][0]) { + CCTK_WARN(0, "The domain is too small?!"); + } + ilo[d] = ihi[d] - min_width * istr[d]; + } + else { + ostringstream buf; + buf << "The grid is unclipped and too small?" << endl + << ilo << endl << ihi << endl << istr << endl << d; + CCTK_WARN(0, buf.str().c_str()); } } } if (veryverbose) { ostringstream buf; - buf << "Corrected coords for striding:" - << endl << lo << endl << up << endl << str; + buf << "Corrected integer coords for min_width:" + << endl << ilo << endl << ihi << endl << istr; CCTK_INFO(buf.str().c_str()); } - - newbbcoord = rbbox(lo, up, str); - } - if (verbose) { - ostringstream buf; - buf << "Done dealing with boundaries. Coord box is:" - << endl << newbbcoord << endl - << "obox is:" << endl << ob; - CCTK_INFO(buf.str().c_str()); + + ibbox newbb(ilo, ihi, istr); + + if (verbose) { + ostringstream buf; + buf << "After dealing with boundaries. Final box is:" + << endl << newbb; + CCTK_INFO(buf.str().c_str()); + } + + newbbs.push_back (newbb); + obs.push_back(ob); } - // 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. - ilo = pos2int(cctkGH, hh, newbbcoord.lower(), reflevel+1); - ihi = pos2int(cctkGH, hh, newbbcoord.upper(), reflevel+1); - ivect istr = bb.stride() / reffact; + // FIXME: check if the newbbs is really different + // from the current bbs + // if not, set do_recompose = 0 + bbs = newbbs; - // Check that the width is sufficient - // This can only be too small if the domain was clipped - for (int d=0; d < dim; ++d) { - if (ihi[d] - ilo[d] < min_width * istr[d]) { - if (ob[d][0]) { - if (ob[d][1]) { - CCTK_WARN(0, "The domain is too small?!"); - } - ihi[d] = ilo[d] + min_width * istr[d]; - } - else if (ob[d][1]) { - if (ob[d][0]) { - CCTK_WARN(0, "The domain is too small?!"); - } - ilo[d] = ihi[d] - min_width * istr[d]; - } - else { - CCTK_WARN(0, "The grid is unclipped and too small?!"); - } + // Set local bbss + + if (bbss.size() < reflevel+2) { + if (verbose) { + CCTK_INFO("Adding new refinement level"); } + bbss.resize(reflevel+2); + obss.resize(reflevel+2); + pss.resize(reflevel+2); } - - if (veryverbose) { - ostringstream buf; - buf << "Corrected integer coords for min_width:" - << endl << ilo << endl << ihi << endl << istr; - CCTK_INFO(buf.str().c_str()); - } + bbss.at(reflevel+1) = bbs; + obss.at(reflevel+1) = obs; + MakeMultigridBoxes (cctkGH, bbss, obss, local_bbsss); - ibbox newbb(ilo, ihi, istr); + // make multiprocessor aware + gh::cprocs ps; + SplitRegions (cctkGH, bbs, obs, ps); - if (verbose) { - ostringstream buf; - buf << "After dealing with boundaries. Final box is:" - << endl << newbb; - CCTK_INFO(buf.str().c_str()); - } + bbss.at(reflevel+1) = bbs; + obss.at(reflevel+1) = obs; + pss.at(reflevel+1) = ps; - newbbs.push_back (newbb); - obs.push_back(ob); - } - - - // FIXME: check if the newbbs is really different from the current bbs - // if not, set do_recompose = 0 - bbs = newbbs; - - // Set local bbss - - if (bbss.size() < reflevel+2) { - if (verbose) { - CCTK_INFO("Adding new refinement level"); - } - bbss.resize(reflevel+2); - obss.resize(reflevel+2); - pss.resize(reflevel+2); - } - bbss.at(reflevel+1) = bbs; - obss.at(reflevel+1) = obs; - MakeMultigridBoxes (cctkGH, bbss, obss, local_bbsss); - - // make multiprocessor aware - gh::cprocs ps; - SplitRegions (cctkGH, bbs, obs, ps); - - bbss.at(reflevel+1) = bbs; - obss.at(reflevel+1) = obs; - pss.at(reflevel+1) = ps; - - } // did_regrid? - else - { - if (bbss.size() > reflevel+1) { - if (verbose) { - CCTK_INFO("Removing refinement level"); + } // did_regrid? + else + { + if (bbss.size() > reflevel+1) { + if (verbose) { + CCTK_INFO("Removing refinement level"); + } } + bbss.resize(reflevel+1); + obss.resize(reflevel+1); + // Set local bbsss + MakeMultigridBoxes (cctkGH, bbss, obss, local_bbsss); + + pss.resize(reflevel+1); + + do_recompose = 1; } - bbss.resize(reflevel+1); - obss.resize(reflevel+1); - // Set local bbsss - MakeMultigridBoxes (cctkGH, bbss, obss, local_bbsss); + + // make multigrid aware + MakeMultigridBoxes (cctkGH, bbss, obss, bbsss); - pss.resize(reflevel+1); + leave_singlemap_mode(const_cast (cctkGH)); + leave_level_mode(const_cast (cctkGH)); - do_recompose = 1; - } - - // make multigrid aware - MakeMultigridBoxes (cctkGH, bbss, obss, bbsss); + } + + enter_level_mode(const_cast (cctkGH), called_on_rl); + enter_singlemap_mode(const_cast (cctkGH), called_on_map); -// cout << "bbsss done:" << endl << bbsss << endl; -// cout << "obss done:" << endl << obss << endl; -// cout << "pss done:" << endl << pss << endl; + if (verbose) { + ostringstream buf; + buf << "Done with it all. Total list is:" + << endl << local_bbsss; + CCTK_INFO(buf.str().c_str()); + } return do_recompose; } -- cgit v1.2.3