aboutsummaryrefslogtreecommitdiff
path: root/CarpetDev
diff options
context:
space:
mode:
Diffstat (limited to 'CarpetDev')
-rw-r--r--CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc609
1 files changed, 325 insertions, 284 deletions
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<ibbox> bbs = bbsss.at(mglevel).at(reflevel);
+ vector<ibbox> bbs = local_bbsss.at(mglevel).at(reflevel);
- ivect low = 100000, upp = -1;
-
- for ( vector<ibbox>::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<CCTK_INT> 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<ibbox> final;
- // Setup the mask.
+ vector<vector<ibbox> > 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<ibbox>::const_iterator bbi = bbs.begin();
+ bbi != bbs.end();
+ ++bbi)
{
- const CCTK_REAL *error_var_ptr =
- static_cast<const CCTK_REAL*>(CCTK_VarDataPtr(cctkGH, 0, error_var));
- const CCTK_REAL *x_var_ptr =
- static_cast<const CCTK_REAL*>(CCTK_VarDataPtr(cctkGH, 0, "Grid::x"));
- const CCTK_REAL *y_var_ptr =
- static_cast<const CCTK_REAL*>(CCTK_VarDataPtr(cctkGH, 0, "Grid::y"));
- const CCTK_REAL *z_var_ptr =
- static_cast<const CCTK_REAL*>(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<CCTK_INT> 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<const CCTK_REAL*>(CCTK_VarDataPtr(cctkGH, 0, error_var));
+ const CCTK_REAL *x_var_ptr =
+ static_cast<const CCTK_REAL*>(CCTK_VarDataPtr(cctkGH, 0, "Grid::x"));
+ const CCTK_REAL *y_var_ptr =
+ static_cast<const CCTK_REAL*>(CCTK_VarDataPtr(cctkGH, 0, "Grid::y"));
+ const CCTK_REAL *z_var_ptr =
+ static_cast<const CCTK_REAL*>(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<ibbox> final;
-
- vector<vector<ibbox> > bbss = bbsss.at(0);
+ }
+ // Define vectors of bboxes final (empty) and todo (contains bbox)
- if (should_regrid) {
+ if (should_regrid) {
- stack<ibbox> todo;
-
- todo.push(bb);
-
- // Define vector of masks (contains error mask)
-
- stack<vector<CCTK_INT> > 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<CCTK_INT> mask = masklist.top(); masklist.pop();
+ stack<ibbox> 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<CCTK_INT> sum_x(nx, 0);
- vector<CCTK_INT> sig_x(nx, 0);
- vector<CCTK_INT> sum_y(ny, 0);
- vector<CCTK_INT> sig_y(ny, 0);
- vector<CCTK_INT> sum_z(nz, 0);
- vector<CCTK_INT> sig_z(nz, 0);
+ // Define vector of masks (contains error mask)
- CCTK_INT fbbox[3][3], fbbox1[3][3], fbbox2[3][3];
+ stack<vector<CCTK_INT> > 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<CCTK_INT> 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<CCTK_INT>
- 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<CCTK_INT> sum_x(nx, 0);
+ vector<CCTK_INT> sig_x(nx, 0);
+ vector<CCTK_INT> sum_y(ny, 0);
+ vector<CCTK_INT> sig_y(ny, 0);
+ vector<CCTK_INT> sum_z(nz, 0);
+ vector<CCTK_INT> 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<CCTK_INT>
- 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<CCTK_INT>
- 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<CCTK_INT>
+ 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<CCTK_INT>
+ 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<CCTK_INT>
+ 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<ibbox> newbbs;
vector<bbvect> 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<dim; ++d) {
assert (mglevel==0);
-
+
// Find the size of the physical domain
-
+
rvect const spacing = base_spacing *
ipow((CCTK_REAL)mgfact, basemglevel) / ipow(reffact, reflevel+1);
ierr = ConvertFromPhysicalBoundary
@@ -475,9 +504,9 @@ namespace CarpetAdaptiveRegrid {
&interior_min[0], &interior_max[0],
&exterior_min[0], &exterior_max[0], &spacing[0]);
assert (!ierr);
-
+
// If need be clip the domain
-
+
rvect lo = newbbcoord.lower();
if (newbbcoord.lower()[d] < physical_min[d]) {
lo[d] = exterior_min[d];
@@ -487,30 +516,33 @@ namespace CarpetAdaptiveRegrid {
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];
-
+
newbbcoord = rbbox(lo, up, str);
}
- if (veryverbose) {
+ 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());
}
-
+
// Convert back to integer coordinates
-
- ilo = pos2int(cctkGH, hh, newbbcoord.lower(), reflevel);
- ihi = pos2int(cctkGH, hh, newbbcoord.upper(), reflevel);
+ // 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;
-
+
// Check that the width is sufficient
// This can only be too small if the domain was clipped
for (int d=0; d < dim; ++d) {
@@ -532,28 +564,27 @@ namespace CarpetAdaptiveRegrid {
}
}
}
-
+
ibbox newbb(ilo, ihi, istr);
-
- if (veryverbose) {
+
+ 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);
}
+
+
// FIXME: check if the newbbs is really different from the current bbs
// if not, set do_recompose = 0
bbs = newbbs;
-
- // make multiprocessor aware
- gh::cprocs ps;
- SplitRegions (cctkGH, bbs, obs, ps);
-
-
+
+ // Set local bbss
+
if (bbss.size() < reflevel+2) {
if (verbose) {
CCTK_INFO("Adding new refinement level");
@@ -562,12 +593,19 @@ namespace CarpetAdaptiveRegrid {
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;
- } // should_regrid?
+ } // did_regrid?
else
{
if (bbss.size() > 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;
}
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-01 23:31:01 +0000