From 06abca2346dae4d881a7733230579a7efb9a4587 Mon Sep 17 00:00:00 2001 From: Ian Hawke Date: Wed, 9 Feb 2005 16:19:00 +0000 Subject: Adaptive regridding symmetry / indexing partial fix Fixed a bug where the array indexing in C++ was wrong. Now something that should be symmetric almost is. There is still an assertion failure when it tries to derefine - for the moment this appears to be a Carpet internal problem. darcs-hash:20050209161923-58c7f-7d13dd7c72c44d935bff9050e3d9281f5f39ec90.gz --- CarpetDev/CarpetAdaptiveRegrid/par/AMR1.par | 4 + CarpetDev/CarpetAdaptiveRegrid/param.ccl | 7 + CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc | 160 +++++++++++++---------- CarpetDev/CarpetAdaptiveRegrid/src/CAR_utils.F90 | 23 +++- 4 files changed, 124 insertions(+), 70 deletions(-) (limited to 'CarpetDev/CarpetAdaptiveRegrid') diff --git a/CarpetDev/CarpetAdaptiveRegrid/par/AMR1.par b/CarpetDev/CarpetAdaptiveRegrid/par/AMR1.par index 5d55f2443..359c35d88 100644 --- a/CarpetDev/CarpetAdaptiveRegrid/par/AMR1.par +++ b/CarpetDev/CarpetAdaptiveRegrid/par/AMR1.par @@ -25,10 +25,14 @@ time::dtfac = 0.5 carpet::max_refinement_levels = 2 carpet::init_3_timelevels = yes carpetadaptiveregrid::max_error = 1.e-6 +carpetadaptiveregrid::pad = 2 +carpetadaptiveregrid::min_width = 6 carpetadaptiveregrid::error_var = "mol::errorestimate[0]" #carpetadaptiveregrid::error_var = "wavemol::phi" carpetadaptiveregrid::regrid_every = 1 +#carpetadaptiveregrid::verbose = yes +#carpetadaptiveregrid::veryverbose = yes cactus::terminate = "time" cactus::cctk_final_time = 0.5 diff --git a/CarpetDev/CarpetAdaptiveRegrid/param.ccl b/CarpetDev/CarpetAdaptiveRegrid/param.ccl index 68690f418..9660d454c 100644 --- a/CarpetDev/CarpetAdaptiveRegrid/param.ccl +++ b/CarpetDev/CarpetAdaptiveRegrid/param.ccl @@ -35,3 +35,10 @@ CCTK_INT pad "Number of cells that the error should be padded" STEERABLE=always 0:* :: "Non negative" } 2 +CCTK_BOOLEAN verbose "Should we be verbose" +{ +} "no" + +CCTK_BOOLEAN veryverbose "Should we be very verbose" +{ +} "no" diff --git a/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc b/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc index 08df72d52..4af9e3879 100644 --- a/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc +++ b/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc @@ -94,7 +94,11 @@ namespace CarpetAdaptiveRegrid { if (reflevel == maxreflevels - 1) return 0; - int do_recompose; + cout << "bbsss at start" << endl << bbsss << endl; + cout << "obss at start" << endl << obss << endl; + cout << "pss at start" << endl << pss << endl; + + CCTK_INT do_recompose; do_recompose = 1; // We assume that we are called in fine->coarse order @@ -136,8 +140,22 @@ namespace CarpetAdaptiveRegrid { 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 = @@ -149,6 +167,12 @@ namespace CarpetAdaptiveRegrid { { 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)); @@ -171,8 +195,16 @@ namespace CarpetAdaptiveRegrid { assert (jj < imax[1] - imin[1]); assert (kk < imax[2] - imin[2]); CCTK_INT mindex = ii + - (imax[0]-imin[0])*(jj + (imax[1] - imin[1]) * kk); + (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); + } } } } @@ -190,22 +222,23 @@ namespace CarpetAdaptiveRegrid { 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])*(j + (imax[1] - imin[1]) * k); + (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, imax[2] - imin[2]); + kk < min(k+pad+1, imax[2] - imin[2]); ++kk) { for (CCTK_INT jj = max(j-pad, 0); - jj < min(j+pad, imax[1] - imin[1]); + jj < min(j+pad+1, imax[1] - imin[1]); ++jj) { for (CCTK_INT ii = max(i-pad, 0); - ii < min(i+pad, imax[0] - imin[0]); + ii < min(i+pad+1, imax[0] - imin[0]); ++ii) { CCTK_INT mindex = ii + - (imax[0]-imin[0])*(jj + (imax[1] - imin[1]) * kk); + (imax[0] - imin[0] + 1)* + (jj + (imax[1] - imin[1] + 1) * kk); if (!mask[mindex]) mask[mindex] = 2; } } @@ -221,8 +254,13 @@ namespace CarpetAdaptiveRegrid { 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])*(j + (imax[1] - imin[1]) * k); + (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]); } } @@ -265,6 +303,12 @@ namespace CarpetAdaptiveRegrid { 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()); + } + vector sum_x(nx, 0); vector sig_x(nx, 0); vector sum_y(ny, 0); @@ -272,7 +316,7 @@ namespace CarpetAdaptiveRegrid { vector sum_z(nz, 0); vector sig_z(nz, 0); - int fbbox[3][3], fbbox1[3][3], fbbox2[3][3]; + CCTK_INT fbbox[3][3], fbbox1[3][3], fbbox2[3][3]; for (CCTK_INT d = 0; d < 3; ++d) { fbbox[0][d] = bb.lower()[d]; @@ -281,7 +325,7 @@ namespace CarpetAdaptiveRegrid { } CCTK_INT didit; - + CCTK_FNAME(check_box)(nx, ny, nz, &mask.front(), &sum_x.front(), &sum_y.front(), &sum_z.front(), @@ -293,8 +337,14 @@ namespace CarpetAdaptiveRegrid { if (didit == 0) { - final.push(bb); - + 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) { @@ -307,7 +357,7 @@ namespace CarpetAdaptiveRegrid { CCTK_INT dny = newbbox1.shape()[1]/newbbox1.stride()[1]; CCTK_INT dnz = newbbox1.shape()[2]/newbbox1.stride()[2]; - vector + vector newmask1(prod(newbbox1.shape()/newbbox1.stride()), 0); CCTK_FNAME(copy_mask)(nx, ny, nz, @@ -315,6 +365,13 @@ namespace CarpetAdaptiveRegrid { 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) { @@ -331,7 +388,7 @@ namespace CarpetAdaptiveRegrid { CCTK_INT dny = newbbox1.shape()[1]/newbbox1.stride()[1]; CCTK_INT dnz = newbbox1.shape()[2]/newbbox1.stride()[2]; - vector + vector newmask1(prod(newbbox1.shape()/newbbox1.stride()), 0); CCTK_FNAME(copy_mask)(nx, ny, nz, @@ -344,7 +401,7 @@ namespace CarpetAdaptiveRegrid { dny = newbbox2.shape()[1]/newbbox2.stride()[1]; dnz = newbbox2.shape()[2]/newbbox2.stride()[2]; - vector + vector newmask2(prod(newbbox2.shape()/newbbox2.stride()), 0); CCTK_FNAME(copy_mask)(nx, ny, nz, @@ -352,6 +409,15 @@ namespace CarpetAdaptiveRegrid { 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."); @@ -376,14 +442,15 @@ namespace CarpetAdaptiveRegrid { // if not, set do_recompose = 0 bbs = newbbs; - cout << bbs << endl; - // make multiprocessor aware gh::cprocs ps; SplitRegions (cctkGH, bbs, obs, ps); 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); @@ -397,62 +464,21 @@ namespace CarpetAdaptiveRegrid { else { if (bbss.size() > reflevel+1) { - bbss.resize(reflevel+1); - obss.resize(reflevel+1); - pss.resize(reflevel+1); + if (verbose) { + CCTK_INFO("Removing refinement level"); + } } + bbss.resize(reflevel+1); + obss.resize(reflevel+1); + pss.resize(reflevel+1); } // make multigrid aware MakeMultigridBoxes (cctkGH, bbss, obss, bbsss); - - // Make a cboxlist (or set) from the boxes on the level. - - // cboxlist is list of cbox's. - - // cbox is a bbox - // plus a dim-D array of booleans mask(i,j,k) - // plus 3 sum arrays (sigma_i = sum_{j,k} mask(i,:,:) etc) - // plus 3 signature arrays (1D laplacian of sigma arrays) - - // First thing to do; make a cbox containing all active boxes - // on this level. - // cbox list contains this cbox. - // Create the mask for this cbox by looking at the error GF. - - // Then recursively for every cbox in the list until no changes: - - // Prune the cbox (remove all zeros from the ends) - - // Split the cbox (find zero in sigma_:, split along that line) - - // if fraction too low then find zero crossings of signature and - // split there - - // if none of these steps are taken the cbox is finished. - - // Most of the actual work is done by the utility functions in - // CAR_utils.F77. - // What we have to do is - - // loop over every box in the list - // for each box call check_box - // if it returns zero the box is done - // if it returns one then the box needs shrinking - // if it returns two then the box needs splitting - - // so, the method will be: - - // Create a cboxlist with one member - - // loop over the list: - // call check_box - // if the return is two then create a new box from newbbox2 and add it - // to the end of the list - // if the return value is one or two shrink the current box to the - // values given in newbbox1 - // if the return value is one or two go redo this loop again. + cout << "bbsss done:" << endl << bbsss << endl; + cout << "obss done:" << endl << obss << endl; + cout << "pss done:" << endl << pss << endl; return do_recompose; } diff --git a/CarpetDev/CarpetAdaptiveRegrid/src/CAR_utils.F90 b/CarpetDev/CarpetAdaptiveRegrid/src/CAR_utils.F90 index 61bc0bc2f..1b2c9f80f 100644 --- a/CarpetDev/CarpetAdaptiveRegrid/src/CAR_utils.F90 +++ b/CarpetDev/CarpetAdaptiveRegrid/src/CAR_utils.F90 @@ -39,6 +39,10 @@ subroutine count_points(nx, ny, nz, mask, sum_x, sum_y, sum_z) sum_x(i) = sum_x(i) + mask(i, j, k) sum_y(j) = sum_y(j) + mask(i, j, k) sum_z(k) = sum_z(k) + mask(i, j, k) + +!!$ if (mask(i,j,k) > 0) then +!!$ write(*,*) "count",i,j,k,mask(i,j,k) +!!$ end if end do end do @@ -591,16 +595,29 @@ subroutine check_box(nx, ny, nz, mask, sum_x, sum_y, sum_z, & CCTK_REAL :: min_density CCTK_INT :: didit - CCTK_INT :: i + CCTK_INT :: i, j, k CCTK_REAL :: density -!!$ write(*,*) nx, ny, nz - +!!$ write(*,*) "sizes:", nx, ny, nz +!!$ write(*,*) "mask size", shape(mask) +!!$ do k = 1, nz +!!$ do j = 1, ny +!!$ do i = 1, nx +!!$ if (mask(i,j,k) == 1) write(*,*) 'fmask set at',i,j,k +!!$ end do +!!$ end do +!!$ end do + !!$ First set up the sums call count_points(nx, ny, nz, mask, sum_x, sum_y, sum_z) +!!$ write(*,*) "sums:" +!!$ write(*,*) "x",sum_x +!!$ write(*,*) "y",sum_y +!!$ write(*,*) "z",sum_z + !!$ Then prune the box call prune_box(nx, ny, nz, sum_x, sum_y, sum_z, bbox, newbbox1, & -- cgit v1.2.3