aboutsummaryrefslogtreecommitdiff
path: root/CarpetDev/CarpetAdaptiveRegrid
diff options
context:
space:
mode:
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
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-10 09:41:36 +0000