aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Carpet/Carpet/param.ccl7
-rw-r--r--Carpet/Carpet/src/Recompose.cc19
-rw-r--r--Carpet/Carpet/src/SetupGH.cc10
-rw-r--r--Carpet/Carpet/src/carpet_public.hh8
-rw-r--r--Carpet/Carpet/src/helpers.cc5
-rw-r--r--Carpet/Carpet/src/variables.cc10
-rw-r--r--Carpet/CarpetLib/src/defs.hh18
-rw-r--r--Carpet/CarpetLib/src/dh.cc3
-rw-r--r--Carpet/CarpetLib/src/gh.cc4
-rw-r--r--Carpet/CarpetRegrid/src/regrid.cc150
10 files changed, 201 insertions, 33 deletions
diff --git a/Carpet/Carpet/param.ccl b/Carpet/Carpet/param.ccl
index a688d6c49..6d9369396 100644
--- a/Carpet/Carpet/param.ccl
+++ b/Carpet/Carpet/param.ccl
@@ -1,5 +1,5 @@
# Parameter definitions for thorn Carpet
-# $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/param.ccl,v 1.15 2001/12/14 16:39:04 schnetter Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/param.ccl,v 1.16 2002/01/09 13:56:24 schnetter Exp $
shares: Cactus
@@ -67,6 +67,11 @@ CCTK_INT refinement_factor "Refinement factor"
1:* :: "must be positive"
} 2
+CCTK_INT multigrid_levels "Number of multigrid levels (including the base level)"
+{
+ 1:* :: "must be positive"
+} 1
+
CCTK_INT multigrid_factor "Multigrid factor"
{
1:* :: "must be positive"
diff --git a/Carpet/Carpet/src/Recompose.cc b/Carpet/Carpet/src/Recompose.cc
index c71c3f112..3792ee082 100644
--- a/Carpet/Carpet/src/Recompose.cc
+++ b/Carpet/Carpet/src/Recompose.cc
@@ -1,5 +1,4 @@
#include <assert.h>
-#include <math.h>
#include <stdlib.h>
#include <list>
@@ -10,12 +9,13 @@
#include "Carpet/CarpetLib/src/bbox.hh"
#include "Carpet/CarpetLib/src/bboxset.hh"
+#include "Carpet/CarpetLib/src/defs.hh"
#include "Carpet/CarpetLib/src/gh.hh"
#include "Carpet/CarpetLib/src/vect.hh"
#include "carpet.hh"
-static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Recompose.cc,v 1.14 2002/01/01 16:48:29 schnetter Exp $";
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Recompose.cc,v 1.15 2002/01/09 13:56:25 schnetter Exp $";
@@ -57,11 +57,17 @@ namespace Carpet {
// No empty levels
assert (bbsss[rl].size() > 0);
for (int c=0; c<(int)bbsss[rl].size(); ++c) {
- // Just one multigrid level
- assert (bbsss[rl][c].size() == 1);
+ // At least one multigrid level
+ assert (bbsss[rl][c].size() > 0);
for (int ml=0; ml<(int)bbsss[rl][c].size(); ++ml) {
- assert (all(bbsss[rl][c][ml].stride()
- == floor(pow((double)reffact,maxreflevels-rl-1)+0.5)));
+ // Check sizes
+ assert (all(bbsss[rl][c][ml].lower() <= bbsss[rl][c][ml].upper()));
+ // Check strides
+ const int str = ipow(reffact, maxreflevels-rl-1) * ipow(mgfact, ml);
+ assert (all(bbsss[rl][c][ml].stride() == str));
+ // Check alignments
+ assert (all(bbsss[rl][c][ml].lower() % str == 0));
+ assert (all(bbsss[rl][c][ml].upper() % str == 0));
}
}
}
@@ -133,7 +139,6 @@ namespace Carpet {
static void Adapt (const cGH* cgh, const int reflevels, gh<dim>* hh)
{
const int nprocs = CCTK_nProcs(cgh);
- const int mglevels = 1; // for now
vector<vector<bbox<int,dim> > > bbss(reflevels);
// note: what this routine calls "ub" is "ub+str" elsewhere
vect<int,dim> rstr = hh->baseextent.stride();
diff --git a/Carpet/Carpet/src/SetupGH.cc b/Carpet/Carpet/src/SetupGH.cc
index 431514da7..1cb27c766 100644
--- a/Carpet/Carpet/src/SetupGH.cc
+++ b/Carpet/Carpet/src/SetupGH.cc
@@ -5,13 +5,14 @@
#include "cctk.h"
#include "cctk_Parameters.h"
+#include "Carpet/CarpetLib/src/defs.hh"
#include "Carpet/CarpetLib/src/dist.hh"
#include "Carpet/CarpetLib/src/ggf.hh"
#include "Carpet/CarpetLib/src/vect.hh"
#include "carpet.hh"
-static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/SetupGH.cc,v 1.17 2001/12/17 13:34:01 schnetter Exp $";
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/SetupGH.cc,v 1.18 2002/01/09 13:56:25 schnetter Exp $";
@@ -40,7 +41,11 @@ namespace Carpet {
// Refinement information
maxreflevels = max_refinement_levels;
reffact = refinement_factor;
- maxreflevelfact = floor(pow((double)reffact, maxreflevels-1) + 0.5);
+ maxreflevelfact = ipow(reffact, maxreflevels-1);
+
+ // Multigrid information
+ mglevels = multigrid_levels;
+ mgfact = multigrid_factor;
// Ghost zones
vect<int,dim> lghosts, ughosts;
@@ -250,7 +255,6 @@ namespace Carpet {
bbss[0] = bbs;
gh<dim>::rexts bbsss;
- const int mglevels = 1;
bbsss = hh->make_multigrid_boxes(bbss, mglevels);
gh<dim>::rprocs pss;
diff --git a/Carpet/Carpet/src/carpet_public.hh b/Carpet/Carpet/src/carpet_public.hh
index 8c967501c..5f097ae33 100644
--- a/Carpet/Carpet/src/carpet_public.hh
+++ b/Carpet/Carpet/src/carpet_public.hh
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/carpet_public.hh,v 1.12 2002/01/02 17:14:08 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/carpet_public.hh,v 1.13 2002/01/09 13:56:25 schnetter Exp $
// It is assumed that the number of components of all arrays is equal
// to the number of components of the grid functions, and that their
@@ -41,6 +41,12 @@ namespace Carpet {
// Refinement factor on finest grid
extern int maxreflevelfact;
+ // Multigrid levels
+ extern int mglevels;
+
+ // Multigrid factor
+ extern int mgfact;
+
// Current iteration per refinement level
extern vector<int> iteration;
diff --git a/Carpet/Carpet/src/helpers.cc b/Carpet/Carpet/src/helpers.cc
index fcc294cb6..6209d1256 100644
--- a/Carpet/Carpet/src/helpers.cc
+++ b/Carpet/Carpet/src/helpers.cc
@@ -7,11 +7,12 @@
#include "cctk.h"
#include "cctk_Parameters.h"
+#include "Carpet/CarpetLib/src/defs.hh"
#include "Carpet/CarpetLib/src/dist.hh"
#include "carpet.hh"
-static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/helpers.cc,v 1.14 2001/12/17 13:34:01 schnetter Exp $";
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/helpers.cc,v 1.15 2002/01/09 13:56:26 schnetter Exp $";
@@ -184,7 +185,7 @@ namespace Carpet {
// Change
reflevel = rl;
const bbox<int,dim>& base = hh->baseextent;
- reflevelfact = floor(pow((double)hh->reffact, reflevel) + 0.5);
+ reflevelfact = ipow(hh->reffact, reflevel);
cgh->cctk_delta_time = base_delta_time / reflevelfact;
for (int d=0; d<dim; ++d) {
cgh->cctk_gsh[d]
diff --git a/Carpet/Carpet/src/variables.cc b/Carpet/Carpet/src/variables.cc
index b1d4e5a0d..a47073da2 100644
--- a/Carpet/Carpet/src/variables.cc
+++ b/Carpet/Carpet/src/variables.cc
@@ -5,7 +5,7 @@
#include "carpet.hh"
-static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/variables.cc,v 1.3 2001/12/17 13:34:02 schnetter Exp $";
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/variables.cc,v 1.4 2002/01/09 13:56:26 schnetter Exp $";
@@ -27,6 +27,12 @@ namespace Carpet {
// Refinement factor on finest grid
int maxreflevelfact;
+ // Multigrid levels
+ int mglevels;
+
+ // Multigrid factor
+ int mgfact;
+
// Current iteration per refinement level
vector<int> iteration;
@@ -35,7 +41,7 @@ namespace Carpet {
int reflevel;
int component;
- // refinement factor of current level: pow(refinement_factor, reflevel)
+ // refinement factor of current level: ipow(refinement_factor, reflevel)
int reflevelfact;
// Time step on base grid
diff --git a/Carpet/CarpetLib/src/defs.hh b/Carpet/CarpetLib/src/defs.hh
index f311b9bc8..a85b164c1 100644
--- a/Carpet/CarpetLib/src/defs.hh
+++ b/Carpet/CarpetLib/src/defs.hh
@@ -5,7 +5,7 @@
copyright : (C) 2000 by Erik Schnetter
email : schnetter@astro.psu.edu
- $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/defs.hh,v 1.5 2001/03/27 22:26:31 eschnett Exp $
+ $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/defs.hh,v 1.6 2002/01/09 13:56:26 schnetter Exp $
***************************************************************************/
@@ -25,6 +25,8 @@
#include <config.h>
#endif
+#include <assert.h>
+
#include <algorithm>
#include <iostream>
#include <list>
@@ -46,6 +48,20 @@ enum centering { vertex_centered, cell_centered };
template<class T>
inline T square (const T& x) { return x*x; }
+// Another useful helper
+template<class T>
+inline T ipow (const T& x, const int y) {
+ if (y<0) {
+ return T(1)/ipow(x,-y);
+ } else if (y==0) {
+ return T(1);
+ } else if (y%2) {
+ return x * ipow(x*x,y/2);
+ } else {
+ return ipow(x*x,y/2);
+ }
+}
+
// Container output
template<class T> ostream& output (ostream& os, const list<T>& l);
template<class T> ostream& output (ostream& os, const set<T>& s);
diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc
index 2f3863dac..d89ed9bc4 100644
--- a/Carpet/CarpetLib/src/dh.cc
+++ b/Carpet/CarpetLib/src/dh.cc
@@ -6,7 +6,7 @@
copyright : (C) 2000 by Erik Schnetter
email : schnetter@astro.psu.edu
- $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.cc,v 1.17 2001/12/09 16:43:09 schnetter Exp $
+ $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.cc,v 1.18 2002/01/09 13:56:27 schnetter Exp $
***************************************************************************/
@@ -280,6 +280,7 @@ void dh<D>::recompose () {
}
#ifdef DEBUG_OUTPUT
+ cout << endl << h << endl;
for (int rl=0; rl<h.reflevels(); ++rl) {
for (int c=0; c<h.components(rl); ++c) {
for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
diff --git a/Carpet/CarpetLib/src/gh.cc b/Carpet/CarpetLib/src/gh.cc
index 7c307860b..95bf16350 100644
--- a/Carpet/CarpetLib/src/gh.cc
+++ b/Carpet/CarpetLib/src/gh.cc
@@ -7,7 +7,7 @@
copyright : (C) 2000 by Erik Schnetter
email : schnetter@astro.psu.edu
- $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.cc,v 1.10 2001/12/14 16:39:42 schnetter Exp $
+ $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.cc,v 1.11 2002/01/09 13:56:27 schnetter Exp $
***************************************************************************/
@@ -189,7 +189,7 @@ gh<D>::cexts gh<D>::make_reflevel_multigrid_boxes (const vector<ibbox>& exts,
up = up - (up - lo) % str;
ext = ibbox(lo,up,str);
- } // for ml
+ } // for ml
} // for c
return mexts;
diff --git a/Carpet/CarpetRegrid/src/regrid.cc b/Carpet/CarpetRegrid/src/regrid.cc
index 5e7e28588..1f3b96b16 100644
--- a/Carpet/CarpetRegrid/src/regrid.cc
+++ b/Carpet/CarpetRegrid/src/regrid.cc
@@ -1,5 +1,4 @@
#include <assert.h>
-#include <math.h>
#include <stdarg.h>
#include <stdlib.h>
@@ -12,6 +11,7 @@
#include "Carpet/CarpetLib/src/bbox.hh"
#include "Carpet/CarpetLib/src/bboxset.hh"
+#include "Carpet/CarpetLib/src/defs.hh"
#include "Carpet/CarpetLib/src/gh.hh"
#include "Carpet/CarpetLib/src/gf.hh"
#include "Carpet/CarpetLib/src/vect.hh"
@@ -19,7 +19,7 @@
#include "carpet.hh"
#include "regrid.hh"
-static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/regrid.cc,v 1.2 2001/12/17 13:32:43 schnetter Exp $";
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetRegrid/src/regrid.cc,v 1.3 2002/01/09 13:56:28 schnetter Exp $";
@@ -112,7 +112,9 @@ namespace CarpetRegrid {
vector<bbox<int,dim> > bbs;
bbs.reserve (bbl.size());
- for (list<bbox<int,dim> >::const_iterator it = bbl.begin(); it != bbl.end(); ++it) {
+ for (list<bbox<int,dim> >::const_iterator it = bbl.begin();
+ it != bbl.end();
+ ++it) {
bbs.push_back (*it);
}
@@ -121,7 +123,6 @@ namespace CarpetRegrid {
SplitRegions (cctkGH, bbs);
gh<dim>::cexts bbss;
- const int mglevels = 1; // arbitrary value
bbss = hh->make_reflevel_multigrid_boxes(bbs, mglevels);
gh<dim>::rexts bbsss = hh->extents;
@@ -184,8 +185,13 @@ namespace CarpetRegrid {
rstr /= hh->reffact;
// refine (arbitrarily) around the center only
rlb = rcentre - (newrextent/2 / rstr) * rstr;
- // // refine (arbitrarily) around the lower boundary only
- // rlb = rlb;
+#if 0
+ // refine (arbitrarily) around the lower boundary only
+ rlb = rlb;
+#endif
+ // honour multigrid factors
+ const int mgstr = ipow(hh->mgfact, mglevels);
+ rlb = (rlb / mgstr) * mgstr;
rub = rlb + newrextent;
// require rub<oldrub because we really want rub-rstr<=oldrub-oldstr
assert (all(rlb >= oldrlb && rub < oldrub));
@@ -211,7 +217,7 @@ namespace CarpetRegrid {
const int rl = reflevel+1;
- const int levfac = floor(pow((double)hh->reffact, rl) + 0.5);
+ const int levfac = ipow(hh->reffact, rl);
assert (all (rstr % levfac == 0));
const vect<int,dim> str (rstr / levfac);
const vect<int,dim> lb (lower[rl-1]);
@@ -238,22 +244,135 @@ namespace CarpetRegrid {
static void
+ MakeRegions_Adaptively_Recombine (list<bbox<int,dim> >& bbl1,
+ list<bbox<int,dim> >& bbl2,
+ list<bbox<int,dim> >& bbl,
+ const bbox<int,dim>& iface)
+ {
+#if 0
+ // remember old bounding boxes
+ bboxset<int,dim> oldboxes;
+ for (list<bbox<int,dim> >::const_iterator ibb1 = bbl1.begin(); ibb1 != bbl1.end(); ++ibb1) {
+ oldboxes += *ibb1;
+ }
+ for (list<bbox<int,dim> >::const_iterator ibb2 = bbl2.begin(); ibb2 != bbl2.end(); ++ibb2) {
+ oldboxes += *ibb2;
+ }
+#endif
+#if 0
+ cout << endl;
+ cout << "MakeRegions_Adaptively_Recombine: initial list:" << endl;
+ cout << bbl << endl;
+ cout << "MakeRegions_Adaptively_Recombine: initial list 1:" << endl;
+ cout << bbl1 << endl;
+ cout << "MakeRegions_Adaptively_Recombine: initial list 2:" << endl;
+ cout << bbl2 << endl;
+#endif
+
+ {
+ // prune boxes on the left
+ list<bbox<int,dim> >::iterator ibb1 = bbl1.begin();
+ while (ibb1 != bbl1.end()) {
+ // is this bbox just to the left of the interface?
+ const bbox<int,dim> iface1 (ibb1->upper()+ibb1->stride(), ibb1->upper()+ibb1->stride(), ibb1->stride());
+ assert (!iface1.empty());
+ if (!iface1.is_contained_in(iface)) {
+ // no: forget it
+ bbl.push_back (*ibb1);
+ ibb1 = bbl1.erase(ibb1);
+ continue;
+ }
+ ++ibb1;
+ } // while
+ }
+
+ {
+ // prune boxes on the right
+ list<bbox<int,dim> >::iterator ibb2 = bbl2.begin();
+ while (ibb2 != bbl2.end()) {
+ // is this bbox just to the right of the interface?
+ const bbox<int,dim> iface2 (ibb2->lower(), ibb2->lower(), ibb2->stride());
+ assert (!iface2.empty());
+ if (! iface2.is_contained_in(iface)) {
+ // no: forget it
+ bbl.push_back (*ibb2);
+ ibb2 = bbl2.erase(ibb2);
+ continue;
+ }
+ ++ibb2;
+ } // while
+ }
+
+ {
+ // walk all boxes on the left
+ list<bbox<int,dim> >::iterator ibb1 = bbl1.begin();
+ while (ibb1 != bbl1.end()) {
+ const bbox<int,dim> iface1 (ibb1->upper()+ibb1->stride(), ibb1->upper()+ibb1->stride(), ibb1->stride());
+ assert (iface1.is_contained_in(iface));
+
+ {
+ // walk all boxes on the right
+ list<bbox<int,dim> >::iterator ibb2 = bbl2.begin();
+ while (ibb2 != bbl2.end()) {
+ const bbox<int,dim> iface2 (ibb2->lower(), ibb2->lower(), ibb2->stride());
+ assert (iface2.is_contained_in(iface));
+
+ // check for a match
+ if (iface1 == iface2) {
+ assert (all(ibb1->stride() == ibb2->stride()));
+ const bbox<int,dim> combined (ibb1->lower(), ibb2->upper(), ibb1->stride());
+ bbl.push_back (combined);
+ ibb1 = bbl1.erase(ibb1);
+ ibb2 = bbl2.erase(ibb2);
+ goto continue_search;
+ }
+
+ ++ibb2;
+ } // while
+ }
+
+ ++ibb1;
+ continue_search:;
+ } // while
+ }
+
+ bbl.splice (bbl.end(), bbl1);
+ bbl.splice (bbl.end(), bbl2);
+
+#if 0
+ // find new bounding boxes
+ bboxset<int,dim> newboxes;
+ for (list<bbox<int,dim> >::const_iterator ibb = bbl.begin(); ibb != bbl.end(); ++ibb) {
+ newboxes += *ibb;
+ }
+
+ // Check that they are equal
+ assert (newboxes.size() <= oldboxes.size());
+ assert ((newboxes.size()==0) == (oldboxes.size()==0));
+ assert (oldboxes == newboxes);
+#endif
+#if 0
+ cout << "MakeRegions_Adaptively_Recombine: final list:" << endl;
+ cout << bbl << endl;
+ cout << endl;
+#endif
+ }
+
+ static void
MakeRegions_Adaptively (const cGH* const cctkGH,
const int minwidth, const double minfraction,
const CCTK_REAL maxerror,
const data<CCTK_REAL,dim>& error,
list<bbox<int,dim> >& bbl,
- const bbox<int,dim> region)
+ const bbox<int,dim>& region)
{
// Just to be sure
assert (! region.empty());
// Count grid points that need to be refined
int cnt = 0;
- CCTK_REAL me=-999;
for (bbox<int,dim>::iterator it=region.begin(); it!=region.end(); ++it) {
if (error[*it] > maxerror) ++cnt;
- me=max(me,error[*it]);
}
const CCTK_REAL fraction = (CCTK_REAL)cnt / region.num_points();
const int width = minval(region.shape() / region.stride());
@@ -274,7 +393,8 @@ namespace CarpetRegrid {
const vect<int,dim> str(region.stride());
vect<int,dim> lo1(lo), lo2(lo);
vect<int,dim> up1(up), up2(up);
- up1[d] = (up[d]+lo[d])/2/str[d]*str[d];
+ const int mgstr = ipow(hh->mgfact, mglevels); // honour multigrid factors
+ up1[d] = ((up[d]+lo[d])/2/str[d]/mgstr)*str[d]*mgstr;
lo2[d] = up1[d]+str[d];
const bbox<int,dim> region1(lo1,up1,str);
const bbox<int,dim> region2(lo2,up2,str);
@@ -282,10 +402,14 @@ namespace CarpetRegrid {
assert (region2.is_contained_in(region));
assert ((region1 & region2).empty());
assert (region1 + region2 == region);
+ list<bbox<int,dim> > bbl1, bbl2;
MakeRegions_Adaptively (cctkGH, minwidth, minfraction, maxerror,
- error, bbl, region1);
+ error, bbl1, region1);
MakeRegions_Adaptively (cctkGH, minwidth, minfraction, maxerror,
- error, bbl, region2);
+ error, bbl2, region2);
+ // Combine regions if possible
+ const bbox<int,dim> iface(lo2,lo2,str);
+ MakeRegions_Adaptively_Recombine (bbl1, bbl2, bbl, iface);
}
}
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-30 10:35:33 +0000