aboutsummaryrefslogtreecommitdiff
path: root/CarpetDev/CarpetAdaptiveRegrid
diff options
context:
space:
mode:
authorIan Hawke <ih@maths.soton.ac.uk>2005-02-08 18:04:00 +0000
committerIan Hawke <ih@maths.soton.ac.uk>2005-02-08 18:04:00 +0000
commit5b4addf1dff82e3d54725385ace11baecacee0de (patch)
treed2fc20448f18392e38260d24df11e268594cc49f /CarpetDev/CarpetAdaptiveRegrid
parent55cf6a666c23377187ac12e26dad69af39882b94 (diff)
Adaptive Regridding might work
Changes to the adaptive regridding routine such that it appears to work for short times on one processor with various changes to, e.g., WaveMoL, before hitting internal Carpet problems. darcs-hash:20050208180418-58c7f-22dc56373df5a208b82d45604267f398bfd3ca0f.gz
Diffstat (limited to 'CarpetDev/CarpetAdaptiveRegrid')
-rw-r--r--CarpetDev/CarpetAdaptiveRegrid/configuration.ccl3
-rw-r--r--CarpetDev/CarpetAdaptiveRegrid/interface.ccl3
-rw-r--r--CarpetDev/CarpetAdaptiveRegrid/par/AMR1.par57
-rw-r--r--CarpetDev/CarpetAdaptiveRegrid/param.ccl21
-rw-r--r--CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc338
-rw-r--r--CarpetDev/CarpetAdaptiveRegrid/src/CAR_Paramcheck.cc6
-rw-r--r--CarpetDev/CarpetAdaptiveRegrid/src/CAR_utils.F9092
-rw-r--r--CarpetDev/CarpetAdaptiveRegrid/src/make.code.defn4
8 files changed, 465 insertions, 59 deletions
diff --git a/CarpetDev/CarpetAdaptiveRegrid/configuration.ccl b/CarpetDev/CarpetAdaptiveRegrid/configuration.ccl
new file mode 100644
index 000000000..d0cc32a17
--- /dev/null
+++ b/CarpetDev/CarpetAdaptiveRegrid/configuration.ccl
@@ -0,0 +1,3 @@
+
+REQUIRES Carpet CarpetLib
+
diff --git a/CarpetDev/CarpetAdaptiveRegrid/interface.ccl b/CarpetDev/CarpetAdaptiveRegrid/interface.ccl
index fe14102e7..f88523f83 100644
--- a/CarpetDev/CarpetAdaptiveRegrid/interface.ccl
+++ b/CarpetDev/CarpetAdaptiveRegrid/interface.ccl
@@ -1,8 +1,7 @@
# Interface definition for thorn CarpetAdaptiveRegrid
# $Header$
-implements: CarpetAdaptiveRegrid
-inherits: Carpet CarpetLib
+implements: CarpetRegrid
uses include header: carpet.hh
diff --git a/CarpetDev/CarpetAdaptiveRegrid/par/AMR1.par b/CarpetDev/CarpetAdaptiveRegrid/par/AMR1.par
new file mode 100644
index 000000000..5d55f2443
--- /dev/null
+++ b/CarpetDev/CarpetAdaptiveRegrid/par/AMR1.par
@@ -0,0 +1,57 @@
+
+ActiveThorns = "coordbase SymBase NaNChecker carpetReduce CartGrid3D carpet carpetlib carpetadaptiveregrid Boundary IOBasic IOUtil carpetIOASCII IDWaveMoL carpetSlab WaveMoL Time MoL"
+
+
+IDWaveMoL::initial_data = "gaussian"
+
+wavemol::bound = "radiation"
+
+mol::initial_data_is_crap = "yes"
+
+##carpet::adaptive_stepsize = "yes"
+
+##carpet::veryverbose = "yes"
+
+grid::domain = "full"
+grid::type = "byspacing"
+grid::avoid_origin = "no"
+driver::global_nx = 51
+driver::global_ny = 51
+driver::global_nz = 51
+grid::dxyz = 0.02
+driver::ghost_size = 1
+time::dtfac = 0.5
+
+carpet::max_refinement_levels = 2
+carpet::init_3_timelevels = yes
+carpetadaptiveregrid::max_error = 1.e-6
+carpetadaptiveregrid::error_var = "mol::errorestimate[0]"
+#carpetadaptiveregrid::error_var = "wavemol::phi"
+carpetadaptiveregrid::regrid_every = 1
+
+
+cactus::terminate = "time"
+cactus::cctk_final_time = 0.5
+
+iobasic::outScalar_every = 1
+iobasic::outScalar_vars = "wavemol::phi mol::errorestimate"
+
+iobasic::outInfo_every = 1
+iobasic::outInfo_vars = "wavemol::phi mol::errorestimate[0]"
+
+ioascii::out1D_every = 1
+ioascii::out1D_vars = "wavemol::scalarevolvemol_scalar mol::errorestimate"
+
+IO::out_dir = $.
+
+#mol::ode_method = "Generic"
+#mol::mol_intermediate_steps = 4
+#mol::mol_num_scratch_levels = 3
+mol::ode_method = "RK45"
+mol::mol_intermediate_steps = 6
+mol::mol_num_scratch_levels = 6
+mol::adaptive_stepsize = "no"
+
+carpet::prolongation_order_space = 3
+carpet::prolongation_order_time = 2
+driver::ghost_size = 2
diff --git a/CarpetDev/CarpetAdaptiveRegrid/param.ccl b/CarpetDev/CarpetAdaptiveRegrid/param.ccl
index d853de204..68690f418 100644
--- a/CarpetDev/CarpetAdaptiveRegrid/param.ccl
+++ b/CarpetDev/CarpetAdaptiveRegrid/param.ccl
@@ -3,28 +3,35 @@
# Refinement criteria for automatic refining
-CCTK_INT minwidth "Minimum width of refined region" STEERABLE=always
+CCTK_INT min_width "Minimum width of refined region" STEERABLE=always
{
1:* :: "must be positive"
} 8
-CCTK_REAL minfraction "Minimum fraction of points in need of refinement in a refined region" STEERABLE=always
+CCTK_REAL min_fraction "Minimum fraction of points in need of refinement in a refined region" STEERABLE=always
{
0:1 :: "must be positive and less than one"
} 0.75
-CCTK_REAL maxerror "Maximum allowed error for non-refined grid points" STEERABLE=always
+CCTK_REAL max_error "Maximum allowed error for non-refined grid points" STEERABLE=always
{
*:* :: "everything goes"
} 1.0
-CCTK_STRING errorvar "Name of grid function that contains the error" STEERABLE=always
+CCTK_STRING error_var "Name of grid function that contains the error" STEERABLE=always
{
".*" :: "must be the name of a grid function"
} ""
-shares: driver
+CCTK_INT regrid_every "Regrid every n time steps" STEERABLE=always
+{
+ -1 :: "regrid never"
+ 0 :: "regrid during initial data calculation only"
+ 1:* :: "regrid every n time steps"
+} 0
-USES INT ghost_size ""
+CCTK_INT pad "Number of cells that the error should be padded" STEERABLE=always
{
-}
+ 0:* :: "Non negative"
+} 2
+
diff --git a/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc b/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc
index e76d20716..08df72d52 100644
--- a/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc
+++ b/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc
@@ -1,7 +1,11 @@
-#include <assert.h>
-
+#include <cassert>
+#include <cmath>
+#include <algorithm>
+#include <iostream>
#include <sstream>
#include <string>
+#include <vector>
+#include <stack>
#include "cctk.h"
#include "cctk_Parameters.h"
@@ -37,7 +41,7 @@ namespace CarpetAdaptiveRegrid {
CCTK_INT* sig_x, CCTK_INT* sig_y, CCTK_INT* sig_z,
const CCTK_INT bbox[3][3],
CCTK_INT newbbox1[3][3], CCTK_INT newbbox2[3][3],
- const CCTK_INT& min_width, const CCTK_REAL& min_density,
+ const CCTK_INT& min_width, const CCTK_REAL& min_fraction,
CCTK_INT& didit);
}
@@ -51,11 +55,11 @@ namespace CarpetAdaptiveRegrid {
const cGH * const cctkGH = (const cGH *) cctkGH_;
- gh<dim>::rexts & bbsss = * (gh<dim>::rexts *) bbsss_;
- gh<dim>::rbnds & obss = * (gh<dim>::rbnds *) obss_;
- gh<dim>::rprocs & pss = * (gh<dim>::rprocs *) pss_;
+ gh::mexts & bbsss = * (gh::mexts *) bbsss_;
+ gh::rbnds & obss = * (gh::rbnds *) obss_;
+ gh::rprocs & pss = * (gh::rprocs *) pss_;
- gh<dim> const & hh = *vhh.at(Carpet::map);
+ gh const & hh = *vhh.at(Carpet::map);
assert (is_singlemap_mode());
@@ -81,9 +85,327 @@ namespace CarpetAdaptiveRegrid {
return 0;
}
+ // Return if it's initial data as we can't handle that yet.
+ if (cctkGH->cctk_iteration ==0) {
+ return 0;
+ }
+
}
+ if (reflevel == maxreflevels - 1) return 0;
+
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<ibbox> bbs = 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());
+
+ vector<CCTK_INT> mask(prod(bb.shape()/bb.stride()), 0);
+
+ // 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));
+
+ 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));
+
+ 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])*(jj + (imax[1] - imin[1]) * kk);
+ mask[mindex] = 1;
+ }
+ }
+ }
+ }
+ } 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])*(j + (imax[1] - imin[1]) * k);
+ if (mask[index] == 1) {
+ for (CCTK_INT kk = max(k-pad, 0);
+ kk < min(k+pad, imax[2] - imin[2]);
+ ++kk)
+ {
+ for (CCTK_INT jj = max(j-pad, 0);
+ jj < min(j+pad, imax[1] - imin[1]);
+ ++jj)
+ {
+ for (CCTK_INT ii = max(i-pad, 0);
+ ii < min(i+pad, imax[0] - imin[0]);
+ ++ii)
+ {
+ CCTK_INT mindex = ii +
+ (imax[0]-imin[0])*(jj + (imax[1] - imin[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])*(j + (imax[1] - imin[1]) * k);
+ if (mask[index] > 1) mask[index] = 1;
+ should_regrid |= (mask[index]);
+ }
+ }
+ }
+
+ // Define vectors of bboxes final (empty) and todo (contains bbox)
+
+ stack<ibbox> final;
+
+ vector<vector<ibbox> > bbss = bbsss.at(0);
+
+
+ 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();
+
+ 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> 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);
+
+ int fbbox[3][3], fbbox1[3][3], fbbox2[3][3];
+
+ 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 (didit == 0) {
+
+ final.push(bb);
+
+ }
+ 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);
+ }
+ 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);
+ }
+ else {
+ CCTK_WARN(0, "The fortran routine must be confused.");
+ }
+
+ } // loop over todo vector (boxes needing to be done).
+
+ // Fixup the stride
+ vector<ibbox> newbbs;
+ vector<bbvect> obs;
+ while (! final.empty()) {
+ ibbox bb = final.top(); final.pop();
+ ivect lo = bb.lower();
+ ivect hi = bb.upper();
+ ivect str = bb.stride() / reffact;
+ newbbs.push_back (ibbox(lo,hi,str));
+
+ bbvect ob(false);
+ obs.push_back(ob);
+ }
+ // FIXME: check if the newbbs is really different from the current bbs
+ // 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) {
+ bbss.resize(reflevel+2);
+ obss.resize(reflevel+2);
+ pss.resize(reflevel+2);
+ }
+
+ bbss.at(reflevel+1) = bbs;
+ obss.at(reflevel+1) = obs;
+ pss.at(reflevel+1) = ps;
+
+ } // should_regrid?
+ else
+ {
+ if (bbss.size() > reflevel+1) {
+ 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.
@@ -105,7 +427,7 @@ namespace CarpetAdaptiveRegrid {
// Split the cbox (find zero in sigma_:, split along that line)
- // if density too low then find zero crossings of signature and
+ // if fraction too low then find zero crossings of signature and
// split there
// if none of these steps are taken the cbox is finished.
diff --git a/CarpetDev/CarpetAdaptiveRegrid/src/CAR_Paramcheck.cc b/CarpetDev/CarpetAdaptiveRegrid/src/CAR_Paramcheck.cc
index 5a6970565..28fe821ed 100644
--- a/CarpetDev/CarpetAdaptiveRegrid/src/CAR_Paramcheck.cc
+++ b/CarpetDev/CarpetAdaptiveRegrid/src/CAR_Paramcheck.cc
@@ -26,9 +26,9 @@ namespace CarpetAdaptiveRegrid {
= (const CCTK_INT *) CCTK_ParameterGet ("domain_from_coordbase", "Carpet", &type);
assert (domain_from_coordbase);
assert (type == PARAMETER_BOOLEAN);
- if (! *domain_from_coordbase) {
- CCTK_PARAMWARN ("CarpetAdaptiveRegrid requires that Carpet::domain_from_coordbase=yes");
- }
+// if (! *domain_from_coordbase) {
+// CCTK_PARAMWARN ("CarpetAdaptiveRegrid requires that Carpet::domain_from_coordbase=yes");
+// }
return 0;
}
diff --git a/CarpetDev/CarpetAdaptiveRegrid/src/CAR_utils.F90 b/CarpetDev/CarpetAdaptiveRegrid/src/CAR_utils.F90
index 64018c999..61bc0bc2f 100644
--- a/CarpetDev/CarpetAdaptiveRegrid/src/CAR_utils.F90
+++ b/CarpetDev/CarpetAdaptiveRegrid/src/CAR_utils.F90
@@ -136,8 +136,8 @@ subroutine prune_box(nx, ny, nz, sum_x, sum_y, sum_z, &
call prune(nx, sum_x, lo, hi)
- newbbox(1,1) = bbox(1,1) + lo - 1
- newbbox(1,2) = bbox(1,2) + hi - nx
+ newbbox(1,1) = bbox(1,1) + (lo - 1) * newbbox(1,3)
+ newbbox(1,2) = bbox(1,2) + (hi - nx) * newbbox(1,3)
!!$ write(*,*) "x dirn",lo,hi,nx,bbox(1,:)
@@ -149,8 +149,8 @@ subroutine prune_box(nx, ny, nz, sum_x, sum_y, sum_z, &
call prune(ny, sum_y, lo, hi)
- newbbox(2,1) = bbox(2,1) + lo - 1
- newbbox(2,2) = bbox(2,2) + hi - ny
+ newbbox(2,1) = bbox(2,1) + (lo - 1) * newbbox(2,3)
+ newbbox(2,2) = bbox(2,2) + (hi - ny) * newbbox(2,3)
!!$ write(*,*) "y dirn",lo,hi,ny,bbox(2,:)
@@ -162,8 +162,8 @@ subroutine prune_box(nx, ny, nz, sum_x, sum_y, sum_z, &
call prune(nz, sum_z, lo, hi)
- newbbox(3,1) = bbox(3,1) + lo - 1
- newbbox(3,2) = bbox(3,2) + hi - nz
+ newbbox(3,1) = bbox(3,1) + (lo - 1) * newbbox(3,3)
+ newbbox(3,2) = bbox(3,2) + (hi - nz) * newbbox(3,3)
!!$ write(*,*) "z dirn",lo,hi,nz,bbox(3,:)
@@ -260,8 +260,8 @@ subroutine split_box_at_hole(nx, ny, nz, sum_x, sum_y, sum_z, &
didit = 1
- newbbox1(3,2) = bbox(3,1) + ihole - 2
- newbbox2(3,1) = newbbox1(3,2) + 1
+ newbbox1(3,2) = bbox(3,1) + (ihole - 1) * bbox(3,3)
+ newbbox2(3,1) = newbbox1(3,2) + bbox(3,3)
return
@@ -277,8 +277,8 @@ subroutine split_box_at_hole(nx, ny, nz, sum_x, sum_y, sum_z, &
didit = 1
- newbbox1(2,2) = bbox(2,1) + ihole - 2
- newbbox2(2,1) = newbbox1(2,2) + 1
+ newbbox1(2,2) = bbox(2,1) + (ihole - 1) * bbox(2,3)
+ newbbox2(2,1) = newbbox1(2,2) + bbox(2,3)
return
@@ -294,8 +294,8 @@ subroutine split_box_at_hole(nx, ny, nz, sum_x, sum_y, sum_z, &
didit = 1
- newbbox1(1,2) = bbox(1,1) + ihole - 2
- newbbox2(1,1) = newbbox1(1,2) + 1
+ newbbox1(1,2) = bbox(1,1) + (ihole - 1) * bbox(1,3)
+ newbbox2(1,1) = newbbox1(1,2) + bbox(1,3)
return
@@ -425,6 +425,8 @@ subroutine split_box_at_sig(nx, ny, nz, sig_x, sig_y, sig_z, &
newbbox1(:,3) = bbox(:,3)
newbbox2(:,3) = bbox(:,3)
+!!$ write(*,*) "Split at sig", nx, ny, nz,bbox
+
!!$ Try splitting in z first
call split(nz, sig_z, min_width, isplit, max_jump)
@@ -435,10 +437,16 @@ subroutine split_box_at_sig(nx, ny, nz, sig_x, sig_y, sig_z, &
didit = 1
- newbbox1(3,2) = bbox(3,1) + isplit - 2
- newbbox2(3,1) = newbbox1(3,2) + 1
+ newbbox1(3,2) = bbox(3,1) + (isplit - 1) * bbox(3,3)
+ newbbox2(3,1) = newbbox1(3,2) + bbox(3,3)
end if
+
+ if (didit > 0) then
+!!$ write(*,*) 'new box 1',newbbox1
+!!$ write(*,*) 'new box 2',newbbox2
+ return
+ end if
!!$ Then try splitting in y next
@@ -452,11 +460,17 @@ subroutine split_box_at_sig(nx, ny, nz, sig_x, sig_y, sig_z, &
newbbox1(3,2) = bbox(3,2)
newbbox2(3,1) = bbox(3,1)
- newbbox1(2,2) = bbox(2,1) + isplit - 2
- newbbox2(2,1) = newbbox1(2,2) + 1
+ newbbox1(2,2) = bbox(2,1) + (isplit - 1) * bbox(2,3)
+ newbbox2(2,1) = newbbox1(2,2) + bbox(2,3)
end if
-
+
+ if (didit > 0) then
+!!$ write(*,*) 'new box 1',newbbox1
+!!$ write(*,*) 'new box 2',newbbox2
+ return
+ end if
+
!!$ Then try splitting in x last
call split(nx, sig_x, min_width, isplit, max_jump)
@@ -471,13 +485,17 @@ subroutine split_box_at_sig(nx, ny, nz, sig_x, sig_y, sig_z, &
newbbox2(3,1) = bbox(3,1)
newbbox1(2,2) = bbox(2,2)
newbbox2(2,1) = bbox(2,1)
- newbbox1(1,2) = bbox(1,1) + isplit - 2
- newbbox2(1,1) = newbbox1(1,2) + 1
+ newbbox1(1,2) = bbox(1,1) + (isplit - 1) * bbox(1,3)
+ newbbox2(1,1) = newbbox1(1,2) + bbox(1,3)
end if
-
- if (didit > 0) return
+ if (didit > 0) then
+!!$ write(*,*) 'new box 1',newbbox1
+!!$ write(*,*) 'new box 2',newbbox2
+ return
+ end if
+
!!$ We should only reach here if we did not find any holes
!!$ So set newbbox2 to dummy values
@@ -509,9 +527,9 @@ subroutine prune_new_box(nx, ny, nz, mask, bbox, newbbox)
CCTK_INT :: tmp_didit
CCTK_INT :: ierr
- tmp_nx = newbbox(1,2) - newbbox(1,1) + 1
- tmp_ny = newbbox(2,2) - newbbox(2,1) + 1
- tmp_nz = newbbox(3,2) - newbbox(3,1) + 1
+ tmp_nx = (newbbox(1,2) - newbbox(1,1))/newbbox(1,3) + 1
+ tmp_ny = (newbbox(2,2) - newbbox(2,1))/newbbox(2,3) + 1
+ tmp_nz = (newbbox(3,2) - newbbox(3,1))/newbbox(3,3) + 1
newbbox(:,3) = bbox(:,3)
@@ -730,10 +748,10 @@ subroutine copy_mask(snx, sny, snz, smask, sbbox, &
CCTK_INT :: snx, sny, snz
CCTK_INT :: smask(snx, sny, snz)
- CCTK_INT :: sbbox(3, 2)
+ CCTK_INT :: sbbox(3, 3)
CCTK_INT :: dnx, dny, dnz
CCTK_INT :: dmask(dnx, dny, dnz)
- CCTK_INT :: dbbox(3, 2)
+ CCTK_INT :: dbbox(3, 3)
CCTK_INT :: i, j, k, si, sj, sk, di, dj, dk
@@ -743,21 +761,21 @@ subroutine copy_mask(snx, sny, snz, smask, sbbox, &
(dbbox(1,2) > sbbox(1,2)) .or. &
(dbbox(2,2) > sbbox(2,2)) .or. &
(dbbox(3,2) > sbbox(3,2)) ) then
- call CCTK_WARN(0, &
- "The destination mask is not contained in the source mask!")
+ call CCTK_WARN(0, \
+ "The destination mask is not contained in the source mask!")
end if
- do k = dbbox(3,1), dbbox(3,2)
- do j = dbbox(2,1), dbbox(2,2)
- do i = dbbox(1,1), dbbox(1,2)
+ do k = dbbox(3,1), dbbox(3,2), dbbox(3,3)
+ do j = dbbox(2,1), dbbox(2,2), dbbox(2, 3)
+ do i = dbbox(1,1), dbbox(1,2), dbbox(1, 3)
- si = 1 + i - sbbox(1,1)
- sj = 1 + j - sbbox(2,1)
- sk = 1 + k - sbbox(3,1)
+ si = 1 + (i - sbbox(1,1)) / dbbox(1, 3)
+ sj = 1 + (j - sbbox(2,1)) / dbbox(2, 3)
+ sk = 1 + (k - sbbox(3,1)) / dbbox(3, 3)
- di = 1 + i - dbbox(1,1)
- dj = 1 + j - dbbox(2,1)
- dk = 1 + k - dbbox(3,1)
+ di = 1 + (i - dbbox(1,1)) / dbbox(1, 3)
+ dj = 1 + (j - dbbox(2,1)) / dbbox(2, 3)
+ dk = 1 + (k - dbbox(3,1)) / dbbox(3, 3)
dmask(di, dj, dk) = smask(si, sj, sk)
diff --git a/CarpetDev/CarpetAdaptiveRegrid/src/make.code.defn b/CarpetDev/CarpetAdaptiveRegrid/src/make.code.defn
index 049bdfb27..8f36e4674 100644
--- a/CarpetDev/CarpetAdaptiveRegrid/src/make.code.defn
+++ b/CarpetDev/CarpetAdaptiveRegrid/src/make.code.defn
@@ -2,8 +2,8 @@
# $Header$
# Source files in this directory
-SRCS = CAR_Paramcheck.cc \\
- CAR.cc \\
+SRCS = CAR_Paramcheck.cc \
+ CAR.cc \
CAR_utils.F90
# Subdirectories containing source files