aboutsummaryrefslogtreecommitdiff
path: root/Carpet
diff options
context:
space:
mode:
Diffstat (limited to 'Carpet')
-rw-r--r--Carpet/Carpet/src/CallFunction.cc17
-rw-r--r--Carpet/Carpet/src/Checksum.cc201
-rw-r--r--Carpet/Carpet/src/Comm.cc62
-rw-r--r--Carpet/Carpet/src/Evolve.cc34
-rw-r--r--Carpet/Carpet/src/Initialise.cc98
-rw-r--r--Carpet/Carpet/src/Poison.cc34
-rw-r--r--Carpet/Carpet/src/Restrict.cc6
-rw-r--r--Carpet/Carpet/src/SetupGH.cc17
-rw-r--r--Carpet/Carpet/src/Shutdown.cc10
-rw-r--r--Carpet/Carpet/src/Storage.cc28
-rw-r--r--Carpet/Carpet/src/carpet_public.hh146
-rw-r--r--Carpet/Carpet/src/helpers.cc69
-rw-r--r--Carpet/Carpet/src/variables.cc12
-rw-r--r--Carpet/CarpetIOASCII/param.ccl8
-rw-r--r--Carpet/CarpetIOASCII/src/ioascii.cc160
-rw-r--r--Carpet/CarpetInterp/src/interp.cc68
-rw-r--r--Carpet/CarpetLib/src/ggf.cc4
-rw-r--r--Carpet/CarpetReduce/src/reduce.cc90
-rw-r--r--Carpet/CarpetSlab/src/slab.cc178
19 files changed, 607 insertions, 635 deletions
diff --git a/Carpet/Carpet/src/CallFunction.cc b/Carpet/Carpet/src/CallFunction.cc
index b1cd011d5..d16675e53 100644
--- a/Carpet/Carpet/src/CallFunction.cc
+++ b/Carpet/Carpet/src/CallFunction.cc
@@ -1,6 +1,8 @@
#include <assert.h>
#include <stdlib.h>
+#include <algorithm>
+
#include "cctk.h"
#include "cctki_GHExtensions.h"
@@ -9,7 +11,7 @@
#include "carpet.hh"
extern "C" {
- static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/CallFunction.cc,v 1.8 2003/05/08 15:35:49 schnetter Exp $";
+ static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/CallFunction.cc,v 1.9 2003/06/18 18:24:27 schnetter Exp $";
CCTK_FILEVERSION(Carpet_Carpet_CallFunction_cc);
}
@@ -32,10 +34,17 @@ namespace Carpet {
// Global operation: call once
if (do_global_mode) {
- Waypoint ("%*sGlobal mode call at %s to %s::%s", 2*reflevel, "",
+ assert (component == -1);
+ const int saved_mglevel = mglevel;
+ if (mglevel!=-1) set_mglevel (cgh, -1);
+ const int saved_reflevel = reflevel;
+ if (reflevel!=-1) set_reflevel (cgh, -1);
+ Waypoint ("Global mode call at %s to %s::%s",
attribute->where, attribute->thorn, attribute->routine);
const int res = CCTK_CallFunction (function, attribute, data);
assert (res==0);
+ if (reflevel!=saved_reflevel) set_reflevel (cgh, saved_reflevel);
+ if (mglevel!=saved_mglevel) set_mglevel (cgh, saved_mglevel);
}
} else if (attribute->level) {
@@ -49,7 +58,7 @@ namespace Carpet {
} else {
// Local operation: call once per component
- BEGIN_LOCAL_COMPONENT_LOOP(cgh) {
+ BEGIN_LOCAL_COMPONENT_LOOP(cgh, CCTK_GF) {
Waypoint ("%*sLocal mode call on component %d at %s to %s::%s",
2*reflevel, "", component,
@@ -57,7 +66,7 @@ namespace Carpet {
const int res = CCTK_CallFunction (function, attribute, data);
assert (res==0);
- } END_LOCAL_COMPONENT_LOOP(cgh);
+ } END_LOCAL_COMPONENT_LOOP;
}
diff --git a/Carpet/Carpet/src/Checksum.cc b/Carpet/Carpet/src/Checksum.cc
index 78b59cde6..ad3972a35 100644
--- a/Carpet/Carpet/src/Checksum.cc
+++ b/Carpet/Carpet/src/Checksum.cc
@@ -9,7 +9,7 @@
#include "carpet.hh"
extern "C" {
- static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Checksum.cc,v 1.8 2002/10/24 10:39:38 schnetter Exp $";
+ static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Checksum.cc,v 1.9 2003/06/18 18:24:27 schnetter Exp $";
CCTK_FILEVERSION(Carpet_Carpet_Checksum_cc);
}
@@ -31,55 +31,61 @@ namespace Carpet {
checksums.resize(CCTK_NumVars());
for (int group=0; group<CCTK_NumGroups(); ++group) {
- for (int var=0; var<CCTK_NumVarsInGroupI(group); ++var) {
-
- const int n = CCTK_FirstVarIndexI(group) + var;
- const int num_tl = CCTK_NumTimeLevelsFromVarI(n);
- assert (num_tl>0);
- const int min_tl = mintl(where, num_tl);
- const int max_tl = maxtl(where, num_tl);
-
- checksums[n].resize(maxreflevels);
- checksums[n][reflevel].resize(num_tl);
- for (int tl=min_tl; tl<=max_tl; ++tl) {
- checksums[n][reflevel][tl].resize(hh->components(reflevel));
- BEGIN_COMPONENT_LOOP(cgh) {
- checksums[n][reflevel][tl][component].valid = false;
- } END_COMPONENT_LOOP(cgh);
- }
+ const int grouptype = CCTK_GroupTypeI(group);
+ if (grouptype == CCTK_GF || reflevel==0) {
+ for (int var=0; var<CCTK_NumVarsInGroupI(group); ++var) {
+
+ const int n = CCTK_FirstVarIndexI(group) + var;
+ const int num_tl = CCTK_NumTimeLevelsFromVarI(n);
+ assert (num_tl>0);
+ const int min_tl = mintl(where, num_tl);
+ const int max_tl = maxtl(where, num_tl);
+
+ checksums[n].resize(maxreflevels);
+ checksums[n][reflevel].resize(num_tl);
+ for (int tl=min_tl; tl<=max_tl; ++tl) {
+ checksums[n][reflevel][tl].resize(hh->components(reflevel));
+ BEGIN_COMPONENT_LOOP(cgh, grouptype) {
+ checksums[n][reflevel][tl][component].valid = false;
+ } END_COMPONENT_LOOP;
+ }
+ }
}
}
for (int group=0; group<CCTK_NumGroups(); ++group) {
- if (CCTK_QueryGroupStorageI(cgh, group)) {
- for (int var=0; var<CCTK_NumVarsInGroupI(group); ++var) {
-
- const int n = CCTK_FirstVarIndexI(group) + var;
- const int sz = CCTK_VarTypeSize(CCTK_VarTypeI(n));
- assert (sz>0);
- const int num_tl = CCTK_NumTimeLevelsFromVarI(n);
- assert (num_tl>0);
- const int min_tl = mintl(where, num_tl);
- const int max_tl = maxtl(where, num_tl);
-
- for (int tl=min_tl; tl<=max_tl; ++tl) {
- BEGIN_LOCAL_COMPONENT_LOOP(cgh) {
- const int gpdim = arrdata[group].info.dim;
- int np = 1;
- for (int d=0; d<gpdim; ++d) {
- np *= *CCTK_ArrayGroupSizeI(cgh, d, group);
- }
- const void* data = cgh->data[n][tl];
- int chk = 0;
- for (int i=0; i<np*sz/(int)sizeof chk; ++i) {
- chk += ((const int*)data)[i];
- }
- checksums[n][reflevel][tl][component].sum = chk;
- checksums[n][reflevel][tl][component].valid = true;
- } END_LOCAL_COMPONENT_LOOP(cgh);
- } // for tl
- } // for var
- } // if has storage
+ const int grouptype = CCTK_GroupTypeI(group);
+ if (grouptype == CCTK_GF || reflevel==0) {
+ if (CCTK_QueryGroupStorageI(cgh, group)) {
+ for (int var=0; var<CCTK_NumVarsInGroupI(group); ++var) {
+
+ const int n = CCTK_FirstVarIndexI(group) + var;
+ const int sz = CCTK_VarTypeSize(CCTK_VarTypeI(n));
+ assert (sz>0);
+ const int num_tl = CCTK_NumTimeLevelsFromVarI(n);
+ assert (num_tl>0);
+ const int min_tl = mintl(where, num_tl);
+ const int max_tl = maxtl(where, num_tl);
+
+ for (int tl=min_tl; tl<=max_tl; ++tl) {
+ BEGIN_LOCAL_COMPONENT_LOOP(cgh, grouptype) {
+ const int gpdim = arrdata[group].info.dim;
+ int np = 1;
+ for (int d=0; d<gpdim; ++d) {
+ np *= *CCTK_ArrayGroupSizeI(cgh, d, group);
+ }
+ const void* data = cgh->data[n][tl];
+ int chk = 0;
+ for (int i=0; i<np*sz/(int)sizeof chk; ++i) {
+ chk += ((const int*)data)[i];
+ }
+ checksums[n][reflevel][tl][component].sum = chk;
+ checksums[n][reflevel][tl][component].valid = true;
+ } END_LOCAL_COMPONENT_LOOP;
+ } // for tl
+ } // for var
+ } // if has storage
+ } // if reflevel
} // for group
}
@@ -95,56 +101,59 @@ namespace Carpet {
assert ((int)checksums.size()==CCTK_NumVars());
for (int group=0; group<CCTK_NumGroups(); ++group) {
- if (CCTK_QueryGroupStorageI(cgh, group)) {
- for (int var=0; var<CCTK_NumVarsInGroupI(group); ++var) {
-
- const int n = CCTK_FirstVarIndexI(group) + var;
- const int sz = CCTK_VarTypeSize(CCTK_VarTypeI(n));
- assert (sz>0);
- const int num_tl = CCTK_NumTimeLevelsFromVarI(n);
- assert (num_tl>0);
- const int min_tl = mintl(where, num_tl);
- const int max_tl = maxtl(where, num_tl);
-
- assert ((int)checksums[n].size()==maxreflevels);
- assert ((int)checksums[n][reflevel].size()==num_tl);
-
- for (int tl=min_tl; tl<=max_tl; ++tl) {
-
- bool unexpected_change = false;
-
- assert ((int)checksums[n][reflevel][tl].size()
- == hh->components(reflevel));
- BEGIN_LOCAL_COMPONENT_LOOP(cgh) {
- if (checksums[n][reflevel][tl][component].valid) {
- const int gpdim = arrdata[group].info.dim;
- int np = 1;
- for (int d=0; d<gpdim; ++d) {
- np *= *CCTK_ArrayGroupSizeI(cgh, d, group);
- }
- const void* data = cgh->data[n][tl];
- int chk = 0;
- for (int i=0; i<np*sz/(int)sizeof chk; ++i) {
- chk += ((const int*)data)[i];
- }
- unexpected_change
- = (unexpected_change
- || chk != checksums[n][reflevel][tl][component].sum);
- }
- } END_LOCAL_COMPONENT_LOOP(cgh);
-
- if (unexpected_change) {
- char* fullname = CCTK_FullName(n);
- CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
- "Timelevel %d of the variable \"%s\" has changed unexpectedly.",
- tl, fullname);
- free (fullname);
- }
-
- } // for tl
-
- } // for var
- } // if has storage
+ const int grouptype = CCTK_GroupTypeI(group);
+ if (grouptype == CCTK_GF || reflevel==0) {
+ if (CCTK_QueryGroupStorageI(cgh, group)) {
+ for (int var=0; var<CCTK_NumVarsInGroupI(group); ++var) {
+
+ const int n = CCTK_FirstVarIndexI(group) + var;
+ const int sz = CCTK_VarTypeSize(CCTK_VarTypeI(n));
+ assert (sz>0);
+ const int num_tl = CCTK_NumTimeLevelsFromVarI(n);
+ assert (num_tl>0);
+ const int min_tl = mintl(where, num_tl);
+ const int max_tl = maxtl(where, num_tl);
+
+ assert ((int)checksums[n].size()==maxreflevels);
+ assert ((int)checksums[n][reflevel].size()==num_tl);
+
+ for (int tl=min_tl; tl<=max_tl; ++tl) {
+
+ bool unexpected_change = false;
+
+ assert ((int)checksums[n][reflevel][tl].size()
+ == hh->components(reflevel));
+ BEGIN_LOCAL_COMPONENT_LOOP(cgh, grouptype) {
+ if (checksums[n][reflevel][tl][component].valid) {
+ const int gpdim = arrdata[group].info.dim;
+ int np = 1;
+ for (int d=0; d<gpdim; ++d) {
+ np *= *CCTK_ArrayGroupSizeI(cgh, d, group);
+ }
+ const void* data = cgh->data[n][tl];
+ int chk = 0;
+ for (int i=0; i<np*sz/(int)sizeof chk; ++i) {
+ chk += ((const int*)data)[i];
+ }
+ unexpected_change
+ = (unexpected_change
+ || chk != checksums[n][reflevel][tl][component].sum);
+ }
+ } END_LOCAL_COMPONENT_LOOP;
+
+ if (unexpected_change) {
+ char* fullname = CCTK_FullName(n);
+ CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "Timelevel %d of the variable \"%s\" has changed unexpectedly.",
+ tl, fullname);
+ free (fullname);
+ }
+
+ } // for tl
+
+ } // for var
+ } // if has storage
+ } // if reflevel
} // for group
}
diff --git a/Carpet/Carpet/src/Comm.cc b/Carpet/Carpet/src/Comm.cc
index c3f59803b..f893b3c71 100644
--- a/Carpet/Carpet/src/Comm.cc
+++ b/Carpet/Carpet/src/Comm.cc
@@ -10,7 +10,7 @@
#include "carpet.hh"
extern "C" {
- static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Comm.cc,v 1.16 2003/05/21 14:30:24 schnetter Exp $";
+ static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Comm.cc,v 1.17 2003/06/18 18:24:27 schnetter Exp $";
CCTK_FILEVERSION(Carpet_Carpet_Comm_cc);
}
@@ -24,19 +24,31 @@ namespace Carpet {
int SyncGroup (cGH* cgh, const char* groupname)
{
- if (hh->local_components(reflevel) != 1 && component != -1) {
- CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
- "Cannot synchronise in local mode "
- "(Tried to synchronise group \"%s\")",
- groupname);
- }
- if (hh->local_components(reflevel) != 1) assert (component == -1);
-
- Checkpoint ("%*sSyncGroup %s", 2*reflevel, "", groupname);
+ DECLARE_CCTK_PARAMETERS;
const int group = CCTK_GroupIndex(groupname);
assert (group>=0 && group<CCTK_NumGroups());
+ Checkpoint ("%*sSyncGroup %s", 2*reflevel, "", groupname);
+
+ const int grouptype = CCTK_GroupTypeI(group);
+
+ if (grouptype == CCTK_GF) {
+ if (reflevel == -1) {
+ CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "Cannot synchronise grid functions in global mode "
+ "(Tried to synchronise group \"%s\")",
+ groupname);
+ }
+ if (hh->local_components(reflevel) != 1 && component != -1) {
+ CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "Cannot synchronise grid functions in local mode "
+ "(Tried to synchronise group \"%s\")",
+ groupname);
+ }
+ if (hh->local_components(reflevel) != 1) assert (component == -1);
+ }
+
if (! CCTK_QueryGroupStorageI(cgh, group)) {
CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
"Cannot synchronise group \"%s\" because it has no storage",
@@ -44,6 +56,8 @@ namespace Carpet {
return -1;
}
+ if (CCTK_NumVarsInGroupI(group)==0) return 0;
+
const int n0 = CCTK_FirstVarIndexI(group);
assert (n0>=0);
const int num_tl = CCTK_NumTimeLevelsFromVarI(n0);
@@ -52,19 +66,19 @@ namespace Carpet {
assert (group<(int)arrdata.size());
for (int var=0; var<(int)arrdata[group].data.size(); ++var) {
- if (CCTK_GroupTypeI(group) == CCTK_GF) {
+ if (grouptype == CCTK_GF) {
if (arrdata[group].do_transfer) {
if (reflevel>0) {
+ // use the current time here (which may be modified by the
+ // user)
+ const CCTK_REAL time = (cgh->cctk_time - cctk_initial_time) / delta_time;
+ if (false) {
+ const CCTK_REAL time1 = tt->time (tl, reflevel, mglevel);
+ const CCTK_REAL time2 = (cgh->cctk_time - cctk_initial_time) / delta_time;
+ assert (fabs((time1 - time2) / (fabs(time1) + fabs(time2) + fabs(cgh->cctk_delta_time))) < 1e-12);
+ }
+
for (int c=0; c<arrdata[group].hh->components(reflevel); ++c) {
- // use the current time here (which may be modified by the
- // user)
- const CCTK_REAL time = cgh->cctk_time / base_delta_time;
- if (false) {
- const CCTK_REAL time1 = tt->time (tl, reflevel, mglevel);
- const CCTK_REAL time2 = cgh->cctk_time / base_delta_time;
- assert (fabs((time1 - time2) / (fabs(time1) + fabs(time2) + fabs(base_delta_time))) < 1e-12);
- }
-
arrdata[group].data[var]->ref_bnd_prolongate
(tl, reflevel, c, mglevel, time);
}
@@ -73,9 +87,11 @@ namespace Carpet {
Checkpoint ("%*s(no prolongating for group %s)",
2*reflevel, "", groupname);
}
- }
- for (int c=0; c<arrdata[group].hh->components(reflevel); ++c) {
- arrdata[group].data[var]->sync (tl, reflevel, c, mglevel);
+ for (int c=0; c<arrdata[group].hh->components(reflevel); ++c) {
+ arrdata[group].data[var]->sync (tl, reflevel, c, mglevel);
+ }
+ } else { // grouptype != CCTK_GF
+ arrdata[group].data[var]->sync (0, 0, 0, 0);
}
}
diff --git a/Carpet/Carpet/src/Evolve.cc b/Carpet/Carpet/src/Evolve.cc
index b20b7d1cd..ebf6f4a66 100644
--- a/Carpet/Carpet/src/Evolve.cc
+++ b/Carpet/Carpet/src/Evolve.cc
@@ -18,7 +18,7 @@
#include "carpet.hh"
extern "C" {
- static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Evolve.cc,v 1.21 2003/05/08 15:35:49 schnetter Exp $";
+ static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Evolve.cc,v 1.22 2003/06/18 18:24:27 schnetter Exp $";
CCTK_FILEVERSION(Carpet_Carpet_Evolve_cc);
}
@@ -139,13 +139,11 @@ namespace Carpet {
// Advance level times
tt->advance_time (reflevel, mglevel);
-
- // TODO: set cctk_time in set_mglevel
- cgh->cctk_time
- = tt->time (0, reflevel, mglevel) * base_delta_time;
+ cgh->cctk_time = cctk_initial_time + tt->time (0, reflevel, mglevel) * cgh->cctk_delta_time;
- Waypoint ("%*sCurrent time is %g%s", 2*reflevel, "",
+ Waypoint ("%*sCurrent time is %g, delta is %g%s", 2*reflevel, "",
cgh->cctk_time,
+ cgh->cctk_delta_time / cgh->cctk_timefac,
do_global_mode ? " (global time)" : "");
// Cycle time levels
@@ -170,14 +168,14 @@ namespace Carpet {
CalculateChecksums (cgh, currenttime);
}
- } END_MGLEVEL_LOOP(cgh);
+ } END_MGLEVEL_LOOP;
// Regrid
Waypoint ("%*sRegrid", 2*reflevel, "");
Regrid (cgh);
}
- } END_REFLEVEL_LOOP(cgh);
+ } END_REFLEVEL_LOOP;
@@ -191,13 +189,11 @@ namespace Carpet {
do_global_mode = cgh->cctk_iteration >= next_global_mode_iter_loop2;
next_global_mode_iter_loop2 = cgh->cctk_iteration + 1;
-
- cgh->cctk_time
- = tt->time (0, reflevel, mglevel) * base_delta_time;
// Restrict
- Waypoint ("%*sCurrent time is %g%s", 2*reflevel, "",
+ Waypoint ("%*sCurrent time is %g, delta is %g%s", 2*reflevel, "",
cgh->cctk_time,
+ cgh->cctk_delta_time / cgh->cctk_timefac,
do_global_mode ? " (global time)" : "");
Restrict (cgh);
@@ -205,10 +201,10 @@ namespace Carpet {
CCTK_ScheduleTraverse ("POSTRESTRICT", cgh, CallFunction);
}
- } END_MGLEVEL_LOOP(cgh);
+ } END_MGLEVEL_LOOP;
}
- } END_REVERSE_REFLEVEL_LOOP(cgh);
+ } END_REVERSE_REFLEVEL_LOOP;
@@ -223,11 +219,9 @@ namespace Carpet {
do_global_mode = cgh->cctk_iteration >= next_global_mode_iter_loop3;
next_global_mode_iter_loop3 = cgh->cctk_iteration + 1;
- cgh->cctk_time
- = tt->time (0, reflevel, mglevel) * base_delta_time;
-
- Waypoint ("%*sCurrent time is %g%s", 2*reflevel, "",
+ Waypoint ("%*sCurrent time is %g, delta is %g%s", 2*reflevel, "",
cgh->cctk_time,
+ cgh->cctk_delta_time / cgh->cctk_timefac,
do_global_mode ? " (global time)" : "");
// Checkpoint
@@ -246,10 +240,10 @@ namespace Carpet {
CheckChecksums (cgh, alltimes);
}
- } END_MGLEVEL_LOOP(cgh);
+ } END_MGLEVEL_LOOP;
}
- } END_REFLEVEL_LOOP(cgh);
+ } END_REFLEVEL_LOOP;
} // main loop
diff --git a/Carpet/Carpet/src/Initialise.cc b/Carpet/Carpet/src/Initialise.cc
index 91b8599a3..49b1b8da1 100644
--- a/Carpet/Carpet/src/Initialise.cc
+++ b/Carpet/Carpet/src/Initialise.cc
@@ -12,7 +12,7 @@
#include "carpet.hh"
extern "C" {
- static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Initialise.cc,v 1.29 2003/05/27 12:01:11 schnetter Exp $";
+ static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Initialise.cc,v 1.30 2003/06/18 18:24:27 schnetter Exp $";
CCTK_FILEVERSION(Carpet_Carpet_Initialise_cc);
}
@@ -42,7 +42,6 @@ namespace Carpet {
// Initialise stuff
cgh->cctk_iteration = 0;
- cgh->cctk_time = cctk_initial_time;
do_global_mode = true;
// Enable storage and communtication
@@ -65,10 +64,12 @@ namespace Carpet {
BEGIN_MGLEVEL_LOOP(cgh) {
+ cgh->cctk_time = cctk_initial_time + tt->time (0, reflevel, mglevel) * cgh->cctk_delta_time;
do_global_mode = reflevel == 0;
- Waypoint ("%*sCurrent time is %g%s", 2*reflevel, "",
+ Waypoint ("%*sCurrent time is %g, delta is %g%s", 2*reflevel, "",
cgh->cctk_time,
+ cgh->cctk_delta_time / cgh->cctk_timefac,
do_global_mode ? " (global time)" : "");
// Checking
@@ -78,20 +79,6 @@ namespace Carpet {
Waypoint ("%*sScheduling BASEGRID", 2*reflevel, "");
CCTK_ScheduleTraverse ("CCTK_BASEGRID", cgh, CallFunction);
- // Allow the time step to be changed
- if (reflevel==0) {
- // Initialise time and time step on coarse grid
- base_delta_time = cgh->cctk_delta_time;
- for (int d=0; d<dim; ++d) {
- base_origin_space[d] = cgh->cctk_origin_space[d];
- }
- } else {
-// assert (abs(cgh->cctk_delta_time - base_delta_time / reflevelfactor)
-// < 1e-6 * base_delta_time);
- // This circumvents a bug in CactusBase/Time
- cgh->cctk_delta_time = base_delta_time / reflevelfact * mglevelfact;
- }
-
if (! init_each_timelevel) {
// Set up the initial data
@@ -103,7 +90,7 @@ namespace Carpet {
} else {
// init_each_timelevel
- bool const orig_do_global_mode = do_global_mode;
+ bool const saved_do_global_mode = do_global_mode;
tt->set_delta
(reflevel, mglevel, - tt->get_delta (reflevel, mglevel));
@@ -115,18 +102,18 @@ namespace Carpet {
for (int tl=-2; tl<=0; ++tl) {
- do_global_mode = orig_do_global_mode && tl==0;
+ do_global_mode = saved_do_global_mode && tl==0;
// Advance level times
tt->advance_time (reflevel, mglevel);
- cgh->cctk_time
- = tt->time (0, reflevel, mglevel) * base_delta_time;
+ cgh->cctk_time = cctk_initial_time + tt->time (0, reflevel, mglevel) * cgh->cctk_delta_time;
// Cycle time levels
CycleTimeLevels (cgh);
- Waypoint ("%*sCurrent time is %g%s", 2*reflevel, "",
+ Waypoint ("%*sCurrent time is %g, delta is %g%s", 2*reflevel, "",
cgh->cctk_time,
+ cgh->cctk_delta_time / cgh->cctk_timefac,
do_global_mode ? " (global time)" : "");
// Set up the initial data
@@ -137,7 +124,7 @@ namespace Carpet {
} // for tl
- do_global_mode = orig_do_global_mode;
+ do_global_mode = saved_do_global_mode;
} // init_each_timelevel
@@ -148,13 +135,13 @@ namespace Carpet {
// Checking
PoisonCheck (cgh, alltimes);
- } END_MGLEVEL_LOOP(cgh);
+ } END_MGLEVEL_LOOP;
// Regrid
Waypoint ("%*sRegrid", 2*reflevel, "");
Regrid (cgh, reflevel);
- } END_REFLEVEL_LOOP(cgh);
+ } END_REFLEVEL_LOOP;
@@ -164,18 +151,20 @@ namespace Carpet {
do_global_mode = reflevel == 0;
- // Restrict
- Waypoint ("%*sCurrent time is %g%s", 2*reflevel, "",
+ Waypoint ("%*sCurrent time is %g, delta is %g%s", 2*reflevel, "",
cgh->cctk_time,
+ cgh->cctk_delta_time / cgh->cctk_timefac,
do_global_mode ? " (global time)" : "");
+
+ // Restrict
Restrict (cgh);
Waypoint ("%*sScheduling POSTRESTRICT", 2*reflevel, "");
CCTK_ScheduleTraverse ("POSTRESTRICT", cgh, CallFunction);
- } END_MGLEVEL_LOOP(cgh);
+ } END_MGLEVEL_LOOP;
- } END_REVERSE_REFLEVEL_LOOP(cgh);
+ } END_REVERSE_REFLEVEL_LOOP;
@@ -192,8 +181,9 @@ namespace Carpet {
do_global_mode = reflevel == 0;
- Waypoint ("%*sCurrent time is %g%s", 2*reflevel, "",
+ Waypoint ("%*sCurrent time is %g, delta is %g%s", 2*reflevel, "",
cgh->cctk_time,
+ cgh->cctk_delta_time / cgh->cctk_timefac,
do_global_mode ? " (global time)" : "");
// Checking
@@ -216,9 +206,9 @@ namespace Carpet {
// Checking
CheckChecksums (cgh, allbutcurrenttime);
- } END_MGLEVEL_LOOP(cgh);
+ } END_MGLEVEL_LOOP;
- } END_REFLEVEL_LOOP(cgh);
+ } END_REFLEVEL_LOOP;
Waypoint ("done with Initialise.");
@@ -239,7 +229,6 @@ namespace Carpet {
BEGIN_REFLEVEL_LOOP(cgh) {
BEGIN_MGLEVEL_LOOP(cgh) {
- cgh->cctk_time = cctk_initial_time;
// Cycle time levels (ignore arrays)
cout << "3TL rl=" << reflevel << " cycling" << endl;
@@ -247,10 +236,10 @@ namespace Carpet {
// Advance level times
tt->advance_time (reflevel, mglevel);
- cgh->cctk_time += cgh->cctk_delta_time;
+ cgh->cctk_time = cctk_initial_time + tt->time (0, reflevel, mglevel) * cgh->cctk_delta_time;
cout << "3TL rl=" << reflevel << " ml=" << mglevel
<< " time=" << tt->get_time (reflevel, mglevel)
- << " time=" << cgh->cctk_time / base_delta_time << endl;
+ << " time=" << cgh->cctk_time / cgh->cctk_delta_time << endl;
// Evolve forward
Waypoint ("%*sScheduling PRESTEP", 2*reflevel, "");
@@ -271,11 +260,11 @@ namespace Carpet {
tt->advance_time (rl, mglevel);
}
cgh->cctk_delta_time *= -1;
- cgh->cctk_time += cgh->cctk_delta_time;
- cgh->cctk_time += cgh->cctk_delta_time;
+ delta_time *= -1;
+ cgh->cctk_time = cctk_initial_time - tt->time (0, reflevel, mglevel) * cgh->cctk_delta_time;
cout << "3TL rl=" << reflevel << " ml=" << mglevel
<< " time=" << tt->get_time (reflevel, mglevel)
- << " time=" << cgh->cctk_time / base_delta_time << endl;
+ << " time=" << cgh->cctk_time / cgh->cctk_delta_time << endl;
// Evolve backward
Waypoint ("%*sScheduling PRESTEP", 2*reflevel, "");
@@ -296,14 +285,14 @@ namespace Carpet {
tt->advance_time (rl, mglevel);
}
cgh->cctk_delta_time *= -1;
- cgh->cctk_time += cgh->cctk_delta_time;
- cgh->cctk_time += cgh->cctk_delta_time;
+ delta_time *= -1;
+ cgh->cctk_time = cctk_initial_time + tt->time (0, reflevel, mglevel) * cgh->cctk_delta_time;
cout << "3TL rl=" << reflevel << " ml=" << mglevel
<< " time=" << tt->get_time (reflevel, mglevel)
- << " time=" << cgh->cctk_time / base_delta_time << endl;
+ << " time=" << cgh->cctk_time / cgh->cctk_delta_time << endl;
- } END_MGLEVEL_LOOP(cgh);
- } END_REFLEVEL_LOOP(cgh);
+ } END_MGLEVEL_LOOP;
+ } END_REFLEVEL_LOOP;
cout << "Hourglass structure in place" << endl;
@@ -315,7 +304,7 @@ namespace Carpet {
BEGIN_REVERSE_REFLEVEL_LOOP(cgh) {
BEGIN_MGLEVEL_LOOP(cgh) {
- cgh->cctk_time = cctk_initial_time + cgh->cctk_delta_time;
+ cgh->cctk_time = cctk_initial_time + cgh->cctk_delta_time / cgh->cctk_timefac;
// Restrict
cout << "3TL rl=" << reflevel << " restricting" << endl;
@@ -335,11 +324,11 @@ namespace Carpet {
tt->advance_time (rl, mglevel);
}
cgh->cctk_delta_time *= -1;
- cgh->cctk_time += cgh->cctk_delta_time;
- cgh->cctk_time += cgh->cctk_delta_time;
+ delta_time *= -1;
+ cgh->cctk_time = cctk_initial_time - tt->time (0, reflevel, mglevel) * cgh->cctk_delta_time;
cout << "3TL rl=" << reflevel << " ml=" << mglevel
<< " time=" << tt->get_time (reflevel, mglevel)
- << " time=" << cgh->cctk_time / base_delta_time << endl;
+ << " time=" << cgh->cctk_time / cgh->cctk_delta_time << endl;
// Cycle time levels
cout << "3TL rl=" << reflevel << " cycling" << endl;
@@ -347,10 +336,10 @@ namespace Carpet {
// Advance level times
tt->advance_time (reflevel, mglevel);
- cgh->cctk_time += cgh->cctk_delta_time;
+ cgh->cctk_time = cctk_initial_time - tt->time (0, reflevel, mglevel) * cgh->cctk_delta_time;
cout << "3TL rl=" << reflevel << " ml=" << mglevel
<< " time=" << tt->get_time (reflevel, mglevel)
- << " time=" << cgh->cctk_time / base_delta_time << endl;
+ << " time=" << cgh->cctk_time / cgh->cctk_delta_time << endl;
// Evolve backward
Waypoint ("%*sScheduling PRESTEP", 2*reflevel, "");
@@ -371,18 +360,17 @@ namespace Carpet {
tt->advance_time (rl, mglevel);
}
cgh->cctk_delta_time *= -1;
- cgh->cctk_time += cgh->cctk_delta_time;
- cgh->cctk_time += cgh->cctk_delta_time;
+ delta_time *= -1;
+ cgh->cctk_time = cctk_initial_time + tt->time (0, reflevel, mglevel) * cgh->cctk_delta_time;
cout << "3TL rl=" << reflevel << " ml=" << mglevel
<< " time=" << tt->get_time (reflevel, mglevel)
- << " time=" << cgh->cctk_time / base_delta_time << endl;
+ << " time=" << cgh->cctk_time / cgh->cctk_delta_time << endl;
- } END_MGLEVEL_LOOP(cgh);
- } END_REVERSE_REFLEVEL_LOOP(cgh);
+ } END_MGLEVEL_LOOP;
+ } END_REVERSE_REFLEVEL_LOOP;
// Reset stuff
cgh->cctk_iteration = 0;
- cgh->cctk_time = cctk_initial_time;
CCTK_INFO ("Finished initialising three timelevels");
}
diff --git a/Carpet/Carpet/src/Poison.cc b/Carpet/Carpet/src/Poison.cc
index a698e1574..8cf7df051 100644
--- a/Carpet/Carpet/src/Poison.cc
+++ b/Carpet/Carpet/src/Poison.cc
@@ -8,7 +8,7 @@
#include "carpet.hh"
extern "C" {
- static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Poison.cc,v 1.11 2003/05/23 23:51:17 schnetter Exp $";
+ static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Poison.cc,v 1.12 2003/06/18 18:24:27 schnetter Exp $";
CCTK_FILEVERSION(Carpet_Carpet_Poison_cc);
}
@@ -38,6 +38,7 @@ namespace Carpet {
void PoisonGroup (cGH* cgh, const int group, const checktimes where)
{
DECLARE_CCTK_PARAMETERS;
+ int ierr;
assert (group>=0 && group<CCTK_NumGroups());
@@ -63,23 +64,20 @@ namespace Carpet {
const int min_tl = mintl(where, num_tl);
const int max_tl = maxtl(where, num_tl);
- int np = 1;
const int gpdim = CCTK_GroupDimI(group);
- for (int d=0; d<gpdim; ++d) {
- np *= *CCTK_ArrayGroupSizeI(cgh, d, group);
- }
- if (CCTK_GroupTypeI(group) == CCTK_GF) {
+ const int grouptype = CCTK_GroupTypeI(group);
+
+ BEGIN_LOCAL_COMPONENT_LOOP(cgh, grouptype) {
- BEGIN_LOCAL_COMPONENT_LOOP(cgh) {
- for (int n=var0; n<var0+nvar; ++n) {
- for (int tl=min_tl; tl<=max_tl; ++tl) {
- memset (cgh->data[n][tl], poison_value, np*sz);
- } // for tl
- } // for var
- } END_LOCAL_COMPONENT_LOOP(cgh);
+ vector<int> lsh(gpdim);
+ ierr = CCTK_GrouplshGI (cgh, gpdim, &lsh[0], group);
+ assert (!ierr);
- } else {
+ int np = 1;
+ for (int d=0; d<gpdim; ++d) {
+ np *= lsh[d];
+ }
for (int n=var0; n<var0+nvar; ++n) {
for (int tl=min_tl; tl<=max_tl; ++tl) {
@@ -87,7 +85,7 @@ namespace Carpet {
} // for tl
} // for var
- } // group type != gf
+ } END_LOCAL_COMPONENT_LOOP;
}
@@ -104,6 +102,8 @@ namespace Carpet {
if (CCTK_QueryGroupStorageI(cgh, group)) {
for (int var=0; var<CCTK_NumVarsInGroupI(group); ++var) {
+ const int grouptype = CCTK_GroupTypeI(group);
+
const int n = CCTK_FirstVarIndexI(group) + var;
const int num_tl = CCTK_NumTimeLevelsFromVarI(n);
@@ -113,7 +113,7 @@ namespace Carpet {
for (int tl=min_tl; tl<=max_tl; ++tl) {
- BEGIN_LOCAL_COMPONENT_LOOP(cgh) {
+ BEGIN_LOCAL_COMPONENT_LOOP(cgh, grouptype) {
vect<int,dim> size(1);
const int gpdim = arrdata[group].info.dim;
for (int d=0; d<gpdim; ++d) {
@@ -161,7 +161,7 @@ namespace Carpet {
fullname, numpoison, tl);
free (fullname);
}
- } END_LOCAL_COMPONENT_LOOP(cgh);
+ } END_LOCAL_COMPONENT_LOOP;
} // for tl
diff --git a/Carpet/Carpet/src/Restrict.cc b/Carpet/Carpet/src/Restrict.cc
index c96ec56a7..810174f18 100644
--- a/Carpet/Carpet/src/Restrict.cc
+++ b/Carpet/Carpet/src/Restrict.cc
@@ -10,7 +10,7 @@
#include "carpet.hh"
extern "C" {
- static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Restrict.cc,v 1.15 2003/05/14 08:33:37 schnetter Exp $";
+ static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Restrict.cc,v 1.16 2003/06/18 18:24:27 schnetter Exp $";
CCTK_FILEVERSION(Carpet_Carpet_Restrict_cc);
}
@@ -42,8 +42,8 @@ namespace Carpet {
const CCTK_REAL time = tt->time (tl, reflevel, mglevel);
if (tl==0) {
const CCTK_REAL time1 = tt->time (tl, reflevel, mglevel);
- const CCTK_REAL time2 = cgh->cctk_time / base_delta_time;
- assert (fabs((time1 - time2) / (fabs(time1) + fabs(time2) + fabs(base_delta_time))) < 1e-12);
+ const CCTK_REAL time2 = cgh->cctk_time / cgh->cctk_delta_time;
+ assert (fabs((time1 - time2) / (fabs(time1) + fabs(time2) + fabs(cgh->cctk_delta_time))) < 1e-12);
}
if (mglevel > 0) {
diff --git a/Carpet/Carpet/src/SetupGH.cc b/Carpet/Carpet/src/SetupGH.cc
index 786320839..809a1e43a 100644
--- a/Carpet/Carpet/src/SetupGH.cc
+++ b/Carpet/Carpet/src/SetupGH.cc
@@ -19,7 +19,7 @@
#include "carpet.hh"
extern "C" {
- static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/SetupGH.cc,v 1.45 2003/05/23 23:51:17 schnetter Exp $";
+ static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/SetupGH.cc,v 1.46 2003/06/18 18:24:27 schnetter Exp $";
CCTK_FILEVERSION(Carpet_Carpet_SetupGH_cc);
}
@@ -461,8 +461,15 @@ namespace Carpet {
// Initialise time step on coarse grid
- base_delta_time = 1.0;
- base_origin_space = vect<CCTK_REAL,dim>((CCTK_REAL)0);
+ cgh->cctk_timefac = 1;
+ for (int d=0; d<dim; ++d) {
+ cgh->cctk_levoff[d] = 0;
+ cgh->cctk_levoffdenom[d] = 1;
+ }
+
+ refleveltimes.resize (maxreflevels);
+ cgh->cctk_time = 0xdeadbeef;
+ cgh->cctk_delta_time = 0xdeadbeef;
@@ -473,8 +480,8 @@ namespace Carpet {
for (int group=0; group<CCTK_NumGroups(); ++group) {
EnableGroupStorage (cgh, CCTK_GroupName(group));
}
- } END_MGLEVEL_LOOP(cgh);
- } END_REFLEVEL_LOOP(cgh);
+ } END_MGLEVEL_LOOP;
+ } END_REFLEVEL_LOOP;
}
Waypoint ("done with SetupGH.");
diff --git a/Carpet/Carpet/src/Shutdown.cc b/Carpet/Carpet/src/Shutdown.cc
index 67759e342..24e41a912 100644
--- a/Carpet/Carpet/src/Shutdown.cc
+++ b/Carpet/Carpet/src/Shutdown.cc
@@ -10,7 +10,7 @@
#include "carpet.hh"
extern "C" {
- static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Shutdown.cc,v 1.9 2003/05/08 15:35:49 schnetter Exp $";
+ static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Shutdown.cc,v 1.10 2003/06/18 18:24:27 schnetter Exp $";
CCTK_FILEVERSION(Carpet_Carpet_Shutdown_cc);
}
@@ -34,7 +34,6 @@ namespace Carpet {
Waypoint ("Current time is %g", cgh->cctk_time);
BEGIN_REFLEVEL_LOOP(cgh) {
-
BEGIN_MGLEVEL_LOOP(cgh) {
do_global_mode = reflevel == 0;
@@ -47,14 +46,13 @@ namespace Carpet {
Waypoint ("%*sScheduling TERMINATE", 2*reflevel, "");
CCTK_ScheduleTraverse ("CCTK_TERMINATE", cgh, CallFunction);
- } END_MGLEVEL_LOOP(cgh);
-
- } END_REFLEVEL_LOOP(cgh);
+ } END_MGLEVEL_LOOP;
+ } END_REFLEVEL_LOOP;
do_global_mode = true;
// Shutdown
- Waypoint ("%*sScheduling SHUTDOWN", 2*reflevel, "");
+ Waypoint ("Scheduling SHUTDOWN");
CCTK_ScheduleTraverse ("CCTK_SHUTDOWN", cgh, CallFunction);
CCTK_PRINTSEPARATOR;
diff --git a/Carpet/Carpet/src/Storage.cc b/Carpet/Carpet/src/Storage.cc
index 86ddd3817..ef2477333 100644
--- a/Carpet/Carpet/src/Storage.cc
+++ b/Carpet/Carpet/src/Storage.cc
@@ -10,7 +10,7 @@
#include "carpet.hh"
extern "C" {
- static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Storage.cc,v 1.20 2003/05/27 12:01:11 schnetter Exp $";
+ static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Storage.cc,v 1.21 2003/06/18 18:24:27 schnetter Exp $";
CCTK_FILEVERSION(Carpet_Carpet_Storage_cc);
}
@@ -44,7 +44,10 @@ namespace Carpet {
const int grouptype = CCTK_GroupTypeI(group);
// No storage change in local mode
- assert (! (component!=-1 && grouptype==CCTK_GF));
+ if (grouptype == CCTK_GF) {
+ assert (reflevel == -1
+ || hh->local_components(reflevel) == 1 || component == -1);
+ }
if (CCTK_QueryGroupStorageI(cgh, group)) {
// storage was enabled previously
@@ -94,11 +97,12 @@ namespace Carpet {
for (int var=0; var<(int)arrdata[group].data.size(); ++var) {
const int n = n0 + var;
switch (CCTK_VarTypeI(n)) {
-#define TYPECASE(N,T) \
- case N: \
- arrdata[group].data[var] = new gf<T,dim> \
- (CCTK_VarName(n), *arrdata[group].tt, *arrdata[group].dd, \
- tmin, tmax, my_prolongation_order_time); \
+#define TYPECASE(N,T) \
+ case N: \
+ assert (! arrdata[group].data[var]); \
+ arrdata[group].data[var] = new gf<T,dim> \
+ (CCTK_VarName(n), *arrdata[group].tt, *arrdata[group].dd, \
+ tmin, tmax, my_prolongation_order_time); \
break;
#include "typecase"
#undef TYPECASE
@@ -115,7 +119,7 @@ namespace Carpet {
} // for
- PoisonGroup (cgh, group, alltimes);
+// PoisonGroup (cgh, group, alltimes);
// storage was not enabled previously
return 0;
@@ -152,9 +156,11 @@ namespace Carpet {
for (int var=0; var<(int)arrdata[group].data.size(); ++var) {
const int n = n0 + var;
switch (CCTK_VarTypeI(n)) {
-#define TYPECASE(N,T) \
- case N: \
- delete (gf<T,dim>*)arrdata[group].data[var]; \
+#define TYPECASE(N,T) \
+ case N: \
+ assert (arrdata[group].data[var]); \
+ delete (gf<T,dim>*)arrdata[group].data[var]; \
+ arrdata[group].data[var] = 0; \
break;
#include "typecase"
#undef TYPECASE
diff --git a/Carpet/Carpet/src/carpet_public.hh b/Carpet/Carpet/src/carpet_public.hh
index f74274b8f..11aa9390b 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.28 2003/05/23 23:51:17 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/carpet_public.hh,v 1.29 2003/06/18 18:24:28 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
@@ -73,13 +73,9 @@ namespace Carpet {
// Is this the time for a global mode call?
extern bool do_global_mode;
-
-
- // Time step on base grid
- extern CCTK_REAL base_delta_time;
-
- // Spatial origin on base grid
- extern vect<CCTK_REAL,dim> base_origin_space;
+ // Current times on the refinement levels
+ extern vector<CCTK_REAL> refleveltimes;
+ extern CCTK_REAL delta_time;
@@ -163,99 +159,117 @@ namespace Carpet {
// Refinement level iterator
-#define BEGIN_REFLEVEL_LOOP(cgh) \
- do { \
- int _rll; \
- assert (reflevel==-1); \
- for (int _rl=0; _rl<hh->reflevels(); ++_rl) { \
- set_reflevel ((cGH*)(cgh), _rl); \
+#define BEGIN_REFLEVEL_LOOP(cgh) \
+ do { \
+ int _rll; \
+ cGH * const _cgh = const_cast<cGH*>(cgh); \
+ assert (reflevel==-1); \
+ for (int _rl=0; _rl<hh->reflevels(); ++_rl) { \
+ set_reflevel (_cgh, _rl); \
{
-#define END_REFLEVEL_LOOP(cgh) \
- } \
- } \
- set_reflevel ((cGH*)(cgh), -1); \
- assert (reflevel==-1); \
- _rll = 0; \
+#define END_REFLEVEL_LOOP \
+ } \
+ } \
+ set_reflevel (_cgh, -1); \
+ _rll = 0; \
} while (0)
// Reverse refinement level iterator
-#define BEGIN_REVERSE_REFLEVEL_LOOP(cgh) \
- do { \
- int _rrll; \
- assert (reflevel==-1); \
- for (int _rl=hh->reflevels()-1; _rl>=0; --_rl) { \
- set_reflevel ((cGH*)(cgh), _rl); \
+#define BEGIN_REVERSE_REFLEVEL_LOOP(cgh) \
+ do { \
+ int _rrll; \
+ cGH * const _cgh = const_cast<cGH*>(cgh); \
+ assert (reflevel==-1); \
+ for (int _rl=hh->reflevels()-1; _rl>=0; --_rl) { \
+ set_reflevel (_cgh, _rl); \
{
-#define END_REVERSE_REFLEVEL_LOOP(cgh) \
- } \
- } \
- set_reflevel ((cGH*)(cgh), -1); \
- assert (reflevel==-1); \
- _rrll = 0; \
+#define END_REVERSE_REFLEVEL_LOOP \
+ } \
+ } \
+ set_reflevel (_cgh, -1); \
+ _rrll = 0; \
} while (0)
// Multigrid level iterator
-#define BEGIN_MGLEVEL_LOOP(cgh) \
- do { \
- int _mgl; \
- assert (reflevel>=0 && reflevel<hh->reflevels()); \
- assert (mglevel==-1); \
- for (int _ml=mglevels-1; _ml>=0; --_ml) { \
- set_mglevel ((cGH*)(cgh), _ml); \
+#define BEGIN_MGLEVEL_LOOP(cgh) \
+ do { \
+ int _mgl; \
+ cGH * const _cgh = const_cast<cGH*>(cgh); \
+ assert (reflevel>=0 && reflevel<hh->reflevels()); \
+ assert (mglevel==-1); \
+ for (int _ml=mglevels-1; _ml>=0; --_ml) { \
+ set_mglevel (_cgh, _ml); \
{
-#define END_MGLEVEL_LOOP(cgh) \
- } \
- } \
- set_mglevel ((cGH*)(cgh), -1); \
- assert (mglevel==-1); \
- _mgl = 0; \
+#define END_MGLEVEL_LOOP \
+ } \
+ } \
+ set_mglevel (_cgh, -1); \
+ _mgl = 0; \
} while (0)
// Component iterator
-#define BEGIN_COMPONENT_LOOP(cgh) \
+#define BEGIN_COMPONENT_LOOP(cgh, grouptype) \
do { \
int _cl; \
- assert (reflevel>=0 && reflevel<hh->reflevels()); \
- assert (mglevel>=0 && mglevel<mglevels); \
- assert (hh->local_components(reflevel)==1 || component==-1); \
- int const _saved_component = component; \
- for (int _c=0; _c<hh->components(reflevel); ++_c) { \
- set_component ((cGH*)(cgh), _c); \
+ cGH * const _cgh = const_cast<cGH*>(cgh); \
+ int const _grouptype = (grouptype); \
+ int _mincl, _maxcl; \
+ if (_grouptype == CCTK_GF) { \
+ assert (reflevel>=0 && reflevel<hh->reflevels()); \
+ assert (mglevel>=0 && mglevel<mglevels); \
+ assert (hh->local_components(reflevel)==1 || component==-1); \
+ _mincl=0; \
+ _maxcl=hh->components(reflevel); \
+ } else { \
+ _mincl=component; \
+ _maxcl=component; \
+ } \
+ for (int _c=_mincl; _c<_maxcl; ++_c) { \
+ if (_grouptype==CCTK_GF) set_component (_cgh, _c); \
{
-#define END_COMPONENT_LOOP(cgh) \
+#define END_COMPONENT_LOOP \
} \
} \
- set_component ((cGH*)(cgh), _saved_component); \
+ if (_grouptype==CCTK_GF) set_component (_cgh, -1); \
_cl = 0; \
} while (0)
-#define BEGIN_LOCAL_COMPONENT_LOOP(cgh) \
+#define BEGIN_LOCAL_COMPONENT_LOOP(cgh, grouptype) \
do { \
int _lcl; \
- assert (reflevel>=0 && reflevel<hh->reflevels()); \
- assert (mglevel>=0 && mglevel<mglevels); \
- assert (hh->local_components(reflevel)==1 || component==-1); \
- int const _saved_component = component; \
- for (int _c=0; _c<hh->components(reflevel); ++_c) { \
- if (hh->is_local(reflevel,_c)) { \
- set_component ((cGH*)(cgh), _c); \
- {
-#define END_LOCAL_COMPONENT_LOOP(cgh) \
- } \
+ cGH * const _cgh = const_cast<cGH*>(cgh); \
+ int const _grouptype = (grouptype); \
+ int _mincl, _maxcl; \
+ if (_grouptype == CCTK_GF) { \
+ assert (reflevel>=0 && reflevel<hh->reflevels()); \
+ assert (mglevel>=0 && mglevel<mglevels); \
+ assert (hh->local_components(reflevel)==1 || component==-1); \
+ _mincl=0; \
+ _maxcl=hh->components(reflevel); \
+ } else { \
+ _mincl=component; \
+ _maxcl=component; \
+ } \
+ for (int _c=_mincl; _c<_maxcl; ++_c) { \
+ if (_grouptype==CCTK_GF || hh->is_local(reflevel,_c)) { \
+ if (_grouptype==CCTK_GF) set_component (_cgh, _c); \
+ {
+#define END_LOCAL_COMPONENT_LOOP \
+ } \
} \
} \
- set_component ((cGH*)(cgh), _saved_component); \
+ if (_grouptype==CCTK_GF) set_component (_cgh, -1); \
_lcl = 0; \
} while (0)
diff --git a/Carpet/Carpet/src/helpers.cc b/Carpet/Carpet/src/helpers.cc
index 1d9eea288..430966caf 100644
--- a/Carpet/Carpet/src/helpers.cc
+++ b/Carpet/Carpet/src/helpers.cc
@@ -5,6 +5,7 @@
#include <mpi.h>
#include "cctk.h"
+#include "cctk_FortranString.h"
#include "cctk_Parameters.h"
#include "Carpet/CarpetLib/src/defs.hh"
@@ -14,7 +15,7 @@
#include "carpet.hh"
extern "C" {
- static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/helpers.cc,v 1.35 2003/05/27 12:01:11 schnetter Exp $";
+ static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/helpers.cc,v 1.36 2003/06/18 18:24:28 schnetter Exp $";
CCTK_FILEVERSION(Carpet_Carpet_helpers_cc);
}
@@ -216,6 +217,15 @@ namespace Carpet {
assert (reflevel>=0 && reflevel<hh->reflevels());
assert (component == -1);
+ // Save
+ if (mglevel == -1) {
+ assert (cgh->cctk_time == 0xdeadbeef);
+ assert (cgh->cctk_delta_time == 0xdeadbeef);
+ } else {
+ refleveltimes[reflevel] = cgh->cctk_time;
+ delta_time = cgh->cctk_delta_time;
+ }
+
// Change
mglevel = ml;
mglevelfact = ipow(mgfact, mglevel);
@@ -227,10 +237,13 @@ namespace Carpet {
mglevelfact = 0xdeadbeef;
cgh->cctk_convlevel = 0xdeadbeef;
- cgh->cctk_delta_time = 0xdeadbeef;
+ cgh->cctk_timefac = 0;
for (int d=0; d<dim; ++d) {
- cgh->cctk_origin_space[d] = 0xdeadbeef;
+ cgh->cctk_levoff[d] = 0xdeadbeef;
+ cgh->cctk_levoffdenom[d] = 0xdeadbeef;
}
+ cgh->cctk_time = 0xdeadbeef;
+ cgh->cctk_delta_time = 0xdeadbeef;
vect<int,dim>::ref(cgh->cctk_gsh) = 0xdeadbeef;
for (int group=0; group<CCTK_NumGroups(); ++group) {
@@ -244,17 +257,18 @@ namespace Carpet {
mglevelfact = ipow(mgfact, mglevel);
cgh->cctk_convlevel = mglevel;
- // TODO: set cctk_time here as well
- cgh->cctk_delta_time = base_delta_time / reflevelfact * mglevelfact;
+ const bbox<int,dim>& baseext = dd->bases[reflevel][mglevel].exterior;
- {
- const bbox<int,dim>& baseext = dd->bases[reflevel][mglevel].exterior;
- for (int d=0; d<dim; ++d) {
- cgh->cctk_origin_space[d] = base_origin_space[d] + cgh->cctk_delta_space[d] / maxreflevelfact * baseext.lower()[d];
- }
+ assert (mglevelfact==1);
+ cgh->cctk_timefac = reflevelfact / mglevelfact;
+ cgh->cctk_time = refleveltimes[reflevel];
+ cgh->cctk_delta_time = delta_time;
+ for (int d=0; d<dim; ++d) {
+ assert (baseext.lower()[d] * reflevelfact % maxreflevelfact == 0);
+ cgh->cctk_levoff[d] = baseext.lower()[d] * reflevelfact / maxreflevelfact;
+ cgh->cctk_levoffdenom[d] = 1;
}
- const bbox<int,dim>& baseext = dd->bases[reflevel][mglevel].exterior;
vect<int,dim>::ref(cgh->cctk_gsh) = baseext.shape() / baseext.stride();
for (int group=0; group<CCTK_NumGroups(); ++group) {
if (CCTK_GroupTypeI(group) == CCTK_GF) {
@@ -355,15 +369,6 @@ namespace Carpet {
}
-#if 0
- cout << "set_component: reflevel=" << reflevel << endl;
- cout << "set_component: gsh=[" << cgh->cctk_gsh[0] << "," << cgh->cctk_gsh[1] << "," << cgh->cctk_gsh[2] << "]" << endl;
- cout << "set_component: lsh=[" << cgh->cctk_lsh[0] << "," << cgh->cctk_lsh[1] << "," << cgh->cctk_lsh[2] << "]" << endl;
- cout << "set_component: lbnd=[" << cgh->cctk_lbnd[0] << "," << cgh->cctk_lbnd[1] << "," << cgh->cctk_lbnd[2] << "]" << endl;
- cout << "set_component: origin_space=[" << cgh->cctk_origin_space[0] << "," << cgh->cctk_origin_space[1] << "," << cgh->cctk_origin_space[2] << "]" << endl;
- cout << "set_component: delta_space=[" << cgh->cctk_delta_space[0] << "," << cgh->cctk_delta_space[1] << "," << cgh->cctk_delta_space[2] << "]" << endl;
-#endif
-
for (int group=0; group<CCTK_NumGroups(); ++group) {
if (CCTK_GroupTypeI(group) == CCTK_GF) {
@@ -435,6 +440,14 @@ namespace Carpet {
return 0;
}
+ extern "C" void CCTK_FCALL CCTK_FNAME(CallScheduleGroup)
+ (int * const ierr, cGH * const cgh, ONE_FORTSTRING_ARG)
+ {
+ ONE_FORTSTRING_CREATE (group);
+ *ierr = CallScheduleGroup (cgh, group);
+ free (group);
+ }
+
// This is a temporary measure to call a local mode function from a
@@ -448,16 +461,16 @@ namespace Carpet {
// we are in global mode
BEGIN_REFLEVEL_LOOP(cgh) {
BEGIN_MGLEVEL_LOOP(cgh) {
- BEGIN_LOCAL_COMPONENT_LOOP(cgh) {
+ BEGIN_LOCAL_COMPONENT_LOOP(cgh, CCTK_GF) {
function (cgh);
- } END_LOCAL_COMPONENT_LOOP(cgh);
- } END_MGLEVEL_LOOP(cgh);
- } END_REFLEVEL_LOOP(cgh);
+ } END_LOCAL_COMPONENT_LOOP;
+ } END_MGLEVEL_LOOP;
+ } END_REFLEVEL_LOOP;
} else {
// we are in level mode
- BEGIN_LOCAL_COMPONENT_LOOP(cgh) {
+ BEGIN_LOCAL_COMPONENT_LOOP(cgh, CCTK_GF) {
function (cgh);
- } END_LOCAL_COMPONENT_LOOP(cgh);
+ } END_LOCAL_COMPONENT_LOOP;
}
return 0;
}
@@ -469,8 +482,8 @@ namespace Carpet {
BEGIN_REFLEVEL_LOOP(cgh) {
BEGIN_MGLEVEL_LOOP(cgh) {
function (cgh);
- } END_MGLEVEL_LOOP(cgh);
- } END_REFLEVEL_LOOP(cgh);
+ } END_MGLEVEL_LOOP;
+ } END_REFLEVEL_LOOP;
return 0;
}
diff --git a/Carpet/Carpet/src/variables.cc b/Carpet/Carpet/src/variables.cc
index 648d84b92..886146a10 100644
--- a/Carpet/Carpet/src/variables.cc
+++ b/Carpet/Carpet/src/variables.cc
@@ -6,7 +6,7 @@
#include "carpet.hh"
extern "C" {
- static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/variables.cc,v 1.14 2003/05/12 16:24:25 schnetter Exp $";
+ static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/variables.cc,v 1.15 2003/06/18 18:24:28 schnetter Exp $";
CCTK_FILEVERSION(Carpet_Carpet_variables_cc);
}
@@ -55,13 +55,9 @@ namespace Carpet {
// Is this the time for a global mode call?
bool do_global_mode;
-
-
- // Time step on base grid
- CCTK_REAL base_delta_time;
-
- // Spatial origin on base grid
- vect<CCTK_REAL,dim> base_origin_space;
+ // Current times on the refinement levels
+ vector<CCTK_REAL> refleveltimes;
+ CCTK_REAL delta_time;
diff --git a/Carpet/CarpetIOASCII/param.ccl b/Carpet/CarpetIOASCII/param.ccl
index 69cc43294..be609ca6a 100644
--- a/Carpet/CarpetIOASCII/param.ccl
+++ b/Carpet/CarpetIOASCII/param.ccl
@@ -1,5 +1,5 @@
# Parameter definitions for thorn CarpetIOASCII
-# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOASCII/param.ccl,v 1.7 2002/10/24 12:00:34 schnetter Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOASCII/param.ccl,v 1.8 2003/06/18 18:24:28 schnetter Exp $
@@ -51,6 +51,12 @@ BOOLEAN out3D_outer_ghosts "Output outer boundary ghost zones as well"
+BOOLEAN output_all_timelevels "Output all timelevels instead of only the current"
+{
+} "no"
+
+
+
BOOLEAN separate_grids "Separate grid levels in the output file by additional empty lines"
{
} "yes"
diff --git a/Carpet/CarpetIOASCII/src/ioascii.cc b/Carpet/CarpetIOASCII/src/ioascii.cc
index 3c2796570..247705e62 100644
--- a/Carpet/CarpetIOASCII/src/ioascii.cc
+++ b/Carpet/CarpetIOASCII/src/ioascii.cc
@@ -30,7 +30,7 @@
#include "ioascii.hh"
extern "C" {
- static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOASCII/src/ioascii.cc,v 1.48 2003/05/13 16:32:10 schnetter Exp $";
+ static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOASCII/src/ioascii.cc,v 1.49 2003/06/18 18:24:28 schnetter Exp $";
CCTK_FILEVERSION(Carpet_CarpetIOASCII_ioascii_cc);
}
@@ -166,7 +166,8 @@ namespace CarpetIOASCII {
assert (n0>=0 && n0<CCTK_NumVars());
const int var = n - n0;
assert (var>=0 && var<CCTK_NumVars());
- const int tl = 0;
+ const int num_tl = CCTK_NumTimeLevelsFromVarI(n);
+ assert (num_tl>=1);
// Check for storage
if (! CCTK_QueryGroupStorageI(cgh, group)) {
@@ -176,6 +177,9 @@ namespace CarpetIOASCII {
return 0;
}
+ const int grouptype = CCTK_GroupTypeI(group);
+ const int rl = grouptype==CCTK_GF ? reflevel : 0;
+
// Get grid hierarchy extentsion from IOUtil
const ioGH * const iogh = (const ioGH *)CCTK_GHExtension (cgh, "IO");
assert (iogh);
@@ -359,80 +363,88 @@ namespace CarpetIOASCII {
// Traverse all components on this refinement and multigrid
// level
- BEGIN_COMPONENT_LOOP(cgh) {
+ BEGIN_COMPONENT_LOOP(cgh, grouptype) {
const ggf<dim>* ff = 0;
assert (var < (int)arrdata[group].data.size());
ff = (ggf<dim>*)arrdata[group].data[var];
- const gdata<dim>* const data
- = (*ff) (tl, reflevel, component, mglevel);
- bbox<int,dim> ext = data->extent();
-
- vect<int,dim> lo = ext.lower();
- vect<int,dim> hi = ext.upper();
- vect<int,dim> str = ext.stride();
-
- // Ignore ghost zones if desired
- for (int d=0; d<dim; ++d) {
- bool output_lower_ghosts
- = cgh->cctk_bbox[2*d] ? out3D_outer_ghosts : out3D_ghosts;
- bool output_upper_ghosts
- = cgh->cctk_bbox[2*d+1] ? out3D_outer_ghosts : out3D_ghosts;
-
- if (! output_lower_ghosts) {
- lo[d] += cgh->cctk_nghostzones[d] * str[d];
- }
- if (! output_upper_ghosts) {
- hi[d] -= cgh->cctk_nghostzones[d] * str[d];
- }
- }
- ext = bbox<int,dim>(lo,hi,str);
-
- // coordinates
- const CCTK_REAL coord_time = cgh->cctk_time;
- vect<CCTK_REAL,dim> global_lower, global_upper;
- vect<CCTK_REAL,dim> coord_delta;
- if (CCTK_GroupTypeI(group) == CCTK_GF) {
+ const int mintl = output_all_timelevels ? 1-num_tl : 0;
+ const int maxtl = 0;
+ for (int tl=mintl; tl<=maxtl; ++tl) {
+
+ const gdata<dim>* const data
+ = (*ff) (tl, rl, component, mglevel);
+ bbox<int,dim> ext = data->extent();
+
+ vect<int,dim> lo = ext.lower();
+ vect<int,dim> hi = ext.upper();
+ vect<int,dim> str = ext.stride();
+
+ // Ignore ghost zones if desired
for (int d=0; d<dim; ++d) {
- const int ierr = CCTK_CoordRange
- (cgh, &global_lower[d], &global_upper[d], d+1, 0, "cart3d");
- assert (!ierr);
- coord_delta[d] = cgh->cctk_delta_space[d] / maxreflevelfact;
+ bool output_lower_ghosts
+ = cgh->cctk_bbox[2*d] ? out3D_outer_ghosts : out3D_ghosts;
+ bool output_upper_ghosts
+ = cgh->cctk_bbox[2*d+1] ? out3D_outer_ghosts : out3D_ghosts;
+
+ if (! output_lower_ghosts) {
+ lo[d] += cgh->cctk_nghostzones[d] * str[d];
+ }
+ if (! output_upper_ghosts) {
+ hi[d] -= cgh->cctk_nghostzones[d] * str[d];
+ }
+ }
+ ext = bbox<int,dim>(lo,hi,str);
+
+ // coordinates
+ const CCTK_REAL coord_time = cgh->cctk_time;
+ vect<CCTK_REAL,dim> global_lower;
+ vect<CCTK_REAL,dim> coord_delta;
+ if (grouptype == CCTK_GF) {
+ for (int d=0; d<dim; ++d) {
+ global_lower[d] = cgh->cctk_origin_space[d];
+ coord_delta[d] = cgh->cctk_delta_space[d] / maxreflevelfact;
+ }
+ } else {
+ for (int d=0; d<dim; ++d) {
+ global_lower[d] = 0.0;
+ coord_delta[d] = 1.0 / (cgh->cctk_gsh[d] - 1);
+ }
}
- } else {
+ // Note: don't permute the "coord_delta" and
+ // "data->extent().lower()"
+ // (it seems that for gcc 2.95 you then pick up the
+ // integer operator*)
+ const vect<CCTK_REAL,dim> coord_lower = global_lower + coord_delta * vect<CCTK_REAL,dim>(lo);
+ const vect<CCTK_REAL,dim> coord_upper = global_lower + coord_delta * vect<CCTK_REAL,dim>(hi);
+
+ vect<int,dim> offset1;
for (int d=0; d<dim; ++d) {
- global_lower[d] = 0;
- global_upper[d] = 1;
- coord_delta[d] = 1.0 / (cgh->cctk_gsh[d] - 1);
+ assert (cgh->cctk_levoffdenom[d]==1);
+ offset1[d] = (cgh->cctk_levoff[d] + offset[d]) * ext.stride()[d];
}
- }
- // Note: don't permute the "coord_delta" and "data->extent().lower()"
- // (it seems that for gcc 2.95 you'll then pick up the
- // integer operator*)
- const vect<CCTK_REAL,dim> coord_lower = global_lower + coord_delta * vect<CCTK_REAL,dim>(lo);
- const vect<CCTK_REAL,dim> coord_upper = global_lower + coord_delta * vect<CCTK_REAL,dim>(hi);
-
- const vect<int,dim> offset1 = ext.lower() + offset * ext.stride();
-#if 0
- cout << "IOA reflevel=" << reflevel << " global=" << global_lower << " coord=" << coord_lower << " offset=" << offset << endl;
-#endif
-
- WriteASCII (file, data, ext, n, cgh->cctk_iteration, offset1, dirs,
- tl, reflevel, component, mglevel,
- coord_time, coord_lower, coord_upper);
-
- // Append EOL after every component
- if (CCTK_MyProc(cgh)==0) {
- if (separate_components) {
- assert (file.good());
- file << endl;
- }
- }
- assert (file.good());
+ for (int d=0; d<outdim; ++d) {
+ offset1[dirs[d]] = ext.lower()[dirs[d]];
+ }
+
+ WriteASCII (file, data, ext, n, cgh->cctk_iteration, offset1, dirs,
+ tl, rl, component, mglevel,
+ coord_time, coord_lower, coord_upper);
+
+ // Append EOL after every component
+ if (CCTK_MyProc(cgh)==0) {
+ if (separate_components) {
+ assert (file.good());
+ file << endl;
+ }
+ }
+ assert (file.good());
+
+ } // for tl
- } END_COMPONENT_LOOP(cgh);
+ } END_COMPONENT_LOOP;
// Append EOL after every complete set of components
if (CCTK_MyProc(cgh)==0) {
@@ -601,29 +613,15 @@ namespace CarpetIOASCII {
{
assert (dir>=1 && dir<=dim);
-#if 0
- CCTK_REAL lower, upper;
- CCTK_CoordRange (cgh, &lower, &upper, dir, 0, "cart3d");
-
- assert (reflevel!=-1 && mglevel!=-1);
- const int npoints = (hh->baseextent.shape()[dir-1] - hh->baseextent.stride()[dir-1]) / hh->bases[reflevel][mglevel].stride()[dir-1] + 1;
-
- const CCTK_REAL rindex = (coord - lower) / (upper - lower) * (npoints-1);
- int cindex = (int)floor(rindex + 0.5 + 1e-6);
-#endif
-
assert (reflevel!=-1 && mglevel!=-1);
const int npoints = cgh->cctk_gsh[dir-1];
- const CCTK_REAL lower = cgh->cctk_origin_space[dir-1];
const CCTK_REAL delta = cgh->cctk_delta_space[dir-1] / cgh->cctk_levfac[dir-1];
+ const CCTK_REAL lower = cgh->cctk_origin_space[dir-1] + delta * cgh->cctk_levoff[dir-1] / cgh->cctk_levoffdenom[dir-1];
const CCTK_REAL upper = lower + (npoints-1) * delta;
const CCTK_REAL rindex = (coord - lower) / delta;
- int cindex = (int)floor(rindex + 0.5 + 1e-6);
-#if 0
- cout << "CTO rl=" << reflevel << " gsh=" << npoints << " l,u,d=" << lower << "," << upper << "," << delta << " coord=" << coord << " cindex=" << cindex << endl;
-#endif
+ int cindex = (int)floor(rindex + 0.75);
if (cindex<0 || cindex>=npoints) {
CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
diff --git a/Carpet/CarpetInterp/src/interp.cc b/Carpet/CarpetInterp/src/interp.cc
index 13f848729..c38bcc99f 100644
--- a/Carpet/CarpetInterp/src/interp.cc
+++ b/Carpet/CarpetInterp/src/interp.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetInterp/src/interp.cc,v 1.9 2003/05/21 16:03:32 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetInterp/src/interp.cc,v 1.10 2003/06/18 18:24:28 schnetter Exp $
#include <assert.h>
#include <math.h>
@@ -19,7 +19,7 @@
#include "interp.hh"
extern "C" {
- static char const * const rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetInterp/src/interp.cc,v 1.9 2003/05/21 16:03:32 schnetter Exp $";
+ static char const * const rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetInterp/src/interp.cc,v 1.10 2003/06/18 18:24:28 schnetter Exp $";
CCTK_FILEVERSION(Carpet_CarpetInterp_interp_cc);
}
@@ -104,12 +104,12 @@ namespace CarpetInterp {
const char * coord_system_name
= CCTK_CoordSystemName (coord_system_handle);
assert (coord_system_name);
- rvect clower, cupper, cdelta;
+ rvect lower, upper, delta;
for (int d=0; d<dim; ++d) {
ierr = CCTK_CoordRange
- (cgh, &clower[d], &cupper[d], d+1, 0, coord_system_name);
+ (cgh, &lower[d], &upper[d], d+1, 0, coord_system_name);
assert (!ierr);
- cdelta[d] = cgh->cctk_delta_space[d];
+ delta[d] = (upper[d] - lower[d]) / (hh->baseextent.shape()[d] - hh->baseextent.stride()[d]);
}
assert (interp_coords);
@@ -127,6 +127,7 @@ namespace CarpetInterp {
int const minrl = reflevel==-1 ? 0 : reflevel;
int const maxrl = reflevel==-1 ? hh->reflevels() : reflevel+1;
+ int const ml = 0;
int maxncomps = 0;
for (int rl=minrl; rl<maxrl; ++rl) {
maxncomps = max(maxncomps, hh->components(rl));
@@ -138,10 +139,9 @@ namespace CarpetInterp {
// TODO: interpolate in time
for (int rl=minrl; rl<maxrl; ++rl) {
int const tl = 0;
- int const ml = 0;
CCTK_REAL const time1 = tt->time (tl, rl, ml);
- CCTK_REAL const time2 = cgh->cctk_time / base_delta_time;
- assert (fabs((time1 - time2) / (fabs(time1) + fabs(time2) + fabs(base_delta_time))) < 1e-12);
+ CCTK_REAL const time2 = cgh->cctk_time / cgh->cctk_delta_time;
+ assert (fabs((time1 - time2) / (fabs(time1) + fabs(time2) + fabs(cgh->cctk_delta_time))) < 1e-12);
}
@@ -157,7 +157,7 @@ namespace CarpetInterp {
assert (interp_coords_type_code == CCTK_VARIABLE_REAL);
rvect pos;
for (int d=0; d<dim; ++d) {
- pos[d] = ((CCTK_REAL const *)interp_coords[d])[n];
+ pos[d] = static_cast<CCTK_REAL const *>(interp_coords[d])[n];
}
// Find the component that this grid point belongs to
@@ -165,17 +165,12 @@ namespace CarpetInterp {
home[n] = -1;
for (int rl=maxrl-1; rl>=minrl; --rl) {
- bbox<int,dim> const& baseext = dd->bases[rl][0].exterior;
- rvect const lower = clower + cdelta / maxreflevelfact * rvect(baseext.lower());
- rvect const delta = cdelta / maxreflevelfact;
-
CCTK_REAL (* const rfloor) (CCTK_REAL const) = floor;
- int const stride = maxreflevelfact / reflevelfact;
- ivect const ipos = ivect(map(rfloor, (pos - lower) / delta / stride + 0.5)) * stride;
+ ivect const ipos = ivect(map(rfloor, (pos - lower) / delta + 0.5));
// TODO: use something faster than a linear search
for (int c=0; c<hh->components(rl); ++c) {
- if (hh->extents[reflevel][c][mglevel].contains(ipos)) {
+ if (hh->extents[rl][c][ml].contains(ipos)) {
rlev[n] = rl;
home[n] = c;
goto found;
@@ -228,7 +223,7 @@ namespace CarpetInterp {
assert (dim==3);
for (int d=0; d<dim; ++d) {
allcoords[ind_prc(myproc,rl,c)][ivect(tmpcnts[ind_rc(rl,c)],d,0)]
- = ((CCTK_REAL const *)interp_coords[d])[n];
+ = static_cast<CCTK_REAL const *>(interp_coords[d])[n];
}
++ tmpcnts[c + (rl-minrl)*maxncomps];
}
@@ -287,15 +282,15 @@ namespace CarpetInterp {
if (reflevel>=minrl && reflevel<maxrl) {
BEGIN_MGLEVEL_LOOP(cgh) {
- BEGIN_LOCAL_COMPONENT_LOOP(cgh) {
+ BEGIN_LOCAL_COMPONENT_LOOP(cgh, CCTK_GF) {
// Find out about the local geometry
ivect lsh;
rvect coord_origin, coord_delta;
for (int d=0; d<dim; ++d) {
lsh[d] = cgh->cctk_lsh[d];
- coord_origin[d] = cgh->cctk_origin_space[d];
coord_delta[d] = cgh->cctk_delta_space[d] / cgh->cctk_levfac[d];
+ coord_origin[d] = cgh->cctk_origin_space[d] + (1.0 * cgh->cctk_levoff[d] / cgh->cctk_levoffdenom[d] + cgh->cctk_lbnd[d]) * coord_delta[d];
}
@@ -354,29 +349,24 @@ namespace CarpetInterp {
tmp_output_arrays[m]
= &alloutputs[ind_prc(p,reflevel,component)][ivect(0,m,0)];
}
- ierr = CCTK_InterpLocalUniform (N_dims,
- local_interp_handle,
- param_table_handle,
- &coord_origin[0],
- &coord_delta[0],
- allhomecnts[ind_prc(p,reflevel,component)],
- interp_coords_type_code,
- &tmp_interp_coords[0],
- N_input_arrays,
- &lsh[0],
- &input_array_type_codes[0],
- &input_arrays[0],
- N_output_arrays,
- output_array_type_codes,
- &tmp_output_arrays[0]);
+
+ ierr = CCTK_InterpLocalUniform
+ (N_dims, local_interp_handle, param_table_handle,
+ &coord_origin[0], &coord_delta[0],
+ allhomecnts[ind_prc(p,reflevel,component)],
+ interp_coords_type_code, &tmp_interp_coords[0],
+ N_input_arrays, &lsh[0],
+ &input_array_type_codes[0], &input_arrays[0],
+ N_output_arrays,
+ output_array_type_codes, &tmp_output_arrays[0]);
assert (!ierr);
} // for processors
- } END_LOCAL_COMPONENT_LOOP(cgh);
- } END_MGLEVEL_LOOP(cgh);
- }
- } END_REFLEVEL_LOOP(cgh);
+ } END_LOCAL_COMPONENT_LOOP;
+ } END_MGLEVEL_LOOP;
+ } // if reflevel active
+ } END_REFLEVEL_LOOP;
if (saved_reflevel!=-1) {
set_reflevel ((cGH*)cgh, saved_reflevel);
}
@@ -407,7 +397,7 @@ namespace CarpetInterp {
for (int m=0; m<N_output_arrays; ++m) {
assert (interp_coords_type_code == CCTK_VARIABLE_REAL);
assert (dim==3);
- ((CCTK_REAL *)output_arrays[m])[n] =
+ static_cast<CCTK_REAL *>(output_arrays[m])[n] =
alloutputs[ind_prc(myproc,rl,c)][ivect(tmpcnts[ind_rc(rl,c)],m,0)];
}
++ tmpcnts[ind_rc(rl,c)];
diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc
index 067542fc8..ba33f8a48 100644
--- a/Carpet/CarpetLib/src/ggf.cc
+++ b/Carpet/CarpetLib/src/ggf.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.cc,v 1.24 2003/05/02 14:23:12 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.cc,v 1.25 2003/06/18 18:24:28 schnetter Exp $
#include <assert.h>
#include <stdlib.h>
@@ -482,7 +482,7 @@ void ggf<D>::ref_restrict (int tl, int rl, int c, int ml,
// if (t.get_time(rl,ml) != t.get_time(rl+1,ml)) {
// printf("WARNING: Time on rl %d is %d, time on rl %d is %d\n",rl,t.get_time(rl,ml),rl+1,t.get_time(rl+1,ml));
// }
- //assert (t.get_time(rl,ml) == t.get_time(rl+1,ml));
+ assert (t.get_time(rl,ml) == t.get_time(rl+1,ml));
const vector<int> tl2s(1,tl);
intercat (tl ,rl ,c,ml, &dh<D>::dboxes::recv_ref_fine,
tl2s,rl+1, ml, &dh<D>::dboxes::send_ref_coarse,
diff --git a/Carpet/CarpetReduce/src/reduce.cc b/Carpet/CarpetReduce/src/reduce.cc
index 1a5bc0271..3996385f1 100644
--- a/Carpet/CarpetReduce/src/reduce.cc
+++ b/Carpet/CarpetReduce/src/reduce.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/src/reduce.cc,v 1.23 2003/05/21 14:31:03 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/src/reduce.cc,v 1.24 2003/06/18 18:24:28 schnetter Exp $
#include <assert.h>
#include <float.h>
@@ -22,7 +22,7 @@
#include "reduce.hh"
extern "C" {
- static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/src/reduce.cc,v 1.23 2003/05/21 14:31:03 schnetter Exp $";
+ static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/src/reduce.cc,v 1.24 2003/06/18 18:24:28 schnetter Exp $";
CCTK_FILEVERSION(Carpet_CarpetReduce_reduce_cc);
}
@@ -716,7 +716,11 @@ namespace CarpetReduce {
Initialise (cgh, proc, num_outvals, &myoutvals[0], outtype, &mycounts[0],
red);
- if (reduce_arrays) {
+ int const saved_component = component;
+ if (component!=-1) {
+ set_component ((cGH*)cgh, -1);
+ }
+ BEGIN_LOCAL_COMPONENT_LOOP(cgh, reduce_arrays ? CCTK_ARRAY : CCTK_GF) {
assert (grpdim<=dim);
int lsh[dim], bbox[2*dim], nghostzones[dim];
@@ -727,34 +731,34 @@ namespace CarpetReduce {
ierr = CCTK_GroupnghostzonesVI(cgh, grpdim, nghostzones, vi);
assert (!ierr);
for (int d=0; d<grpdim; ++d) {
- assert (lsh[d]>=0);
- assert (nghostzones[d]>=0 && 2*nghostzones[d]<=lsh[d]);
+ assert (lsh[d]>=0);
+ assert (nghostzones[d]>=0 && 2*nghostzones[d]<=lsh[d]);
}
vector<const void*> inarrays (num_invars);
for (int n=0; n<num_invars; ++n) {
- inarrays[n] = CCTK_VarDataPtrI(cgh, 0, invars[n]);
- assert (inarrays[n]);
+ inarrays[n] = CCTK_VarDataPtrI(cgh, 0, invars[n]);
+ assert (inarrays[n]);
}
const int intype = CCTK_VarTypeI(vi);
for (int n=0; n<num_invars; ++n) {
- assert (CCTK_VarTypeI(invars[n]) == intype);
+ assert (CCTK_VarTypeI(invars[n]) == intype);
}
vect<int,dim> mylsh, mynghostzones;
vect<vect<int,2>,dim> mybbox;
for (int d=0; d<grpdim; ++d) {
- mylsh[d] = lsh[d];
+ mylsh[d] = lsh[d];
mybbox[d][0] = bbox[2*d ];
mybbox[d][1] = bbox[2*d+1];
- mynghostzones[d] = nghostzones[d];
+ mynghostzones[d] = nghostzones[d];
}
for (int d=grpdim; d<dim; ++d) {
- mylsh[d] = 1;
+ mylsh[d] = 1;
mybbox[d][0] = 0;
mybbox[d][1] = 0;
- mynghostzones[d] = 0;
+ mynghostzones[d] = 0;
}
Reduce (cgh, proc, &mylsh[0], &mybbox[0][0], &mynghostzones[0],
@@ -762,64 +766,10 @@ namespace CarpetReduce {
num_outvals, &myoutvals[0], outtype,
&mycounts[0], red);
- } else { // ! reduce_arrays
-
- int const saved_component = component;
- if (component!=-1) {
- set_component ((cGH*)cgh, -1);
- }
- BEGIN_LOCAL_COMPONENT_LOOP(cgh) {
-
- assert (grpdim<=dim);
- int lsh[dim], bbox[2*dim], nghostzones[dim];
- ierr = CCTK_GrouplshVI(cgh, grpdim, lsh, vi);
- assert (!ierr);
- ierr = CCTK_GroupbboxVI(cgh, 2*grpdim, bbox, vi);
- assert (!ierr);
- ierr = CCTK_GroupnghostzonesVI(cgh, grpdim, nghostzones, vi);
- assert (!ierr);
- for (int d=0; d<grpdim; ++d) {
- assert (lsh[d]>=0);
- assert (nghostzones[d]>=0 && 2*nghostzones[d]<=lsh[d]);
- }
-
- vector<const void*> inarrays (num_invars);
- for (int n=0; n<num_invars; ++n) {
- inarrays[n] = CCTK_VarDataPtrI(cgh, 0, invars[n]);
- assert (inarrays[n]);
- }
-
- const int intype = CCTK_VarTypeI(vi);
- for (int n=0; n<num_invars; ++n) {
- assert (CCTK_VarTypeI(invars[n]) == intype);
- }
-
- vect<int,dim> mylsh, mynghostzones;
- vect<vect<int,2>,dim> mybbox;
- for (int d=0; d<grpdim; ++d) {
- mylsh[d] = lsh[d];
- mybbox[d][0] = bbox[2*d ];
- mybbox[d][1] = bbox[2*d+1];
- mynghostzones[d] = nghostzones[d];
- }
- for (int d=grpdim; d<dim; ++d) {
- mylsh[d] = 1;
- mybbox[d][0] = 0;
- mybbox[d][1] = 0;
- mynghostzones[d] = 0;
- }
-
- Reduce (cgh, proc, &mylsh[0], &mybbox[0][0], &mynghostzones[0],
- num_invars, &inarrays[0], intype,
- num_outvals, &myoutvals[0], outtype,
- &mycounts[0], red);
-
- } END_LOCAL_COMPONENT_LOOP(cgh);
- if (saved_component!=-1) {
- set_component ((cGH*)cgh, saved_component);
- }
-
- } // ! reduce_arrays
+ } END_LOCAL_COMPONENT_LOOP;
+ if (saved_component!=-1) {
+ set_component ((cGH*)cgh, saved_component);
+ }
Finalise (cgh, proc, num_outvals, outvals, outtype,
&myoutvals[0], &mycounts[0],
diff --git a/Carpet/CarpetSlab/src/slab.cc b/Carpet/CarpetSlab/src/slab.cc
index 3d25430e1..6f522d671 100644
--- a/Carpet/CarpetSlab/src/slab.cc
+++ b/Carpet/CarpetSlab/src/slab.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/slab.cc,v 1.8 2003/05/21 14:31:29 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/slab.cc,v 1.9 2003/06/18 18:24:28 schnetter Exp $
#include <assert.h>
#include <stdlib.h>
@@ -21,7 +21,7 @@
#include "slab.hh"
extern "C" {
- static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/slab.cc,v 1.8 2003/05/21 14:31:29 schnetter Exp $";
+ static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/slab.cc,v 1.9 2003/06/18 18:24:28 schnetter Exp $";
CCTK_FILEVERSION(Carpet_CarpetSlab_slab_cc);
}
@@ -43,20 +43,6 @@ namespace CarpetSlab {
const int stride[/*hdim*/],
const int length[/*hdim*/])
{
- if (reflevel == -1) {
- CCTK_WARN (0, "It is not possible to use hyperslabbing in global mode");
- }
-
- // Check global state
- assert (reflevel>=0);
- assert (mglevel>=0);
-
- // Save global state
- int saved_component = component;
- if (component!=-1) {
- set_component ((cGH*)cgh, -1);
- }
-
// Check Cactus grid hierarchy
assert (cgh);
@@ -82,6 +68,11 @@ namespace CarpetSlab {
const int typesize = CCTK_VarTypeSize(gp.vartype);
assert (typesize>0);
+ if (gp.grouptype==CCTK_GF && reflevel==-1) {
+ CCTK_WARN (0, "It is not possible to use hyperslabbing for a grid function in global mode");
+ }
+ const int rl = gp.grouptype==CCTK_GF ? reflevel : 0;
+
// Check dimension
assert (hdim>=0 && hdim<=gp.dim);
@@ -148,98 +139,85 @@ namespace CarpetSlab {
memset (hdata, 0, totalsize * typesize);
}
- if (hh->components(reflevel) > 0) {
-
- // Only temporarily
- component = 0;
+ // Get sample data
+ const gdata<dim>* mydata;
+ mydata = (*myff)(tl, rl, 0, 0);
+
+ // Stride of data in memory
+ const vect<int,dim> str = mydata->extent().stride();
+
+ // Stride of collected data
+ vect<int,dim> hstr = str;
+ for (int dd=0; dd<hdim; ++dd) {
+ hstr[dirs[dd]-1] *= stride[dd];
+ }
+
+ // Lower bound of collected data
+ vect<int,dim> hlb(0);
+ for (int d=0; d<gp.dim; ++d) {
+ hlb[d] = origin[d] * str[d];
+ }
+
+ // Upper bound of collected data
+ vect<int,dim> hub = hlb;
+ for (int dd=0; dd<hdim; ++dd) {
+ hub[dirs[dd]-1] += (length[dd]-1) * hstr[dirs[dd]-1];
+ }
+
+ // Calculate extent to collect
+ const bbox<int,dim> hextent (hlb, hub, hstr);
+ assert (hextent.num_points() == totalsize);
+
+ // Create collector data object
+ void* myhdata = rank==collect_proc ? hdata : 0;
+ gdata<dim>* const alldata = mydata->make_typed();
+ alldata->allocate (hextent, collect_proc, myhdata);
+
+ // Done with the temporary stuff
+ mydata = 0;
+
+ // Loop over all components, copying data from them
+ BEGIN_LOCAL_COMPONENT_LOOP (cgh, gp.grouptype) {
- // Get sample data
- const gdata<dim>* mydata;
- mydata = (*myff)(tl, reflevel, component, mglevel);
+ const int c = gp.grouptype==CCTK_GF ? component : 0;
- // Stride of data in memory
- const vect<int,dim> str = mydata->extent().stride();
+ // Get data object
+ mydata = (*myff)(tl, rl, c, mglevel);
- // Stride of collected data
- vect<int,dim> hstr = str;
- for (int dd=0; dd<hdim; ++dd) {
- hstr[dirs[dd]-1] *= stride[dd];
- }
+ // Calculate overlapping extents
+ const bboxset<int,dim> myextents
+ = ((mydd->boxes[rl][c][mglevel].sync_not
+ | mydd->boxes[rl][c][mglevel].interior)
+ & hextent);
- // Lower bound of collected data
- vect<int,dim> hlb(0);
- for (int d=0; d<gp.dim; ++d) {
- hlb[d] = origin[d] * str[d];
+ // Loop over overlapping extents
+ for (bboxset<int,dim>::const_iterator ext_iter = myextents.begin();
+ ext_iter != myextents.end();
+ ++ext_iter) {
+
+ // Copy data
+ alldata->copy_from (mydata, *ext_iter);
+
}
- // Upper bound of collected data
- vect<int,dim> hub = hlb;
- for (int dd=0; dd<hdim; ++dd) {
- hub[dirs[dd]-1] += (length[dd]-1) * hstr[dirs[dd]-1];
+ } END_LOCAL_COMPONENT_LOOP;
+
+ // Copy result to all processors
+ if (dest_proc == -1) {
+ for (int proc=0; proc<CCTK_nProcs(cgh); ++proc) {
+ if (proc != collect_proc) {
+
+ void* myhdata = rank==proc ? hdata : 0;
+ gdata<dim>* const tmpdata = mydata->make_typed();
+ tmpdata->allocate (alldata->extent(), proc, myhdata);
+ tmpdata->copy_from (alldata, alldata->extent());
+ delete tmpdata;
+
+ }
}
-
- // Calculate extent to collect
- const bbox<int,dim> hextent (hlb, hub, hstr);
- assert (hextent.num_points() == totalsize);
-
- // Create collector data object
- void* myhdata = rank==collect_proc ? hdata : 0;
- gdata<dim>* const alldata = mydata->make_typed();
- alldata->allocate (hextent, collect_proc, myhdata);
-
- // Done with the temporary stuff
- mydata = 0;
- component = -1;
-
- // Loop over all components, copying data from them
- assert (component == -1);
- for (component=0; component<hh->components(reflevel); ++component) {
-
- // Get data object
- mydata = (*myff)(tl, reflevel, component, mglevel);
-
- // Calculate overlapping extents
- const bboxset<int,dim> myextents
- = ((mydd->boxes[reflevel][component][mglevel].sync_not
- | mydd->boxes[reflevel][component][mglevel].interior)
- & hextent);
-
- // Loop over overlapping extents
- for (bboxset<int,dim>::const_iterator ext_iter = myextents.begin();
- ext_iter != myextents.end();
- ++ext_iter) {
-
- // Copy data
- alldata->copy_from (mydata, *ext_iter);
-
- }
-
- } // Loop over components
- component = -1;
-
- // Copy result to all processors
- if (dest_proc == -1) {
- for (int proc=0; proc<CCTK_nProcs(cgh); ++proc) {
- if (proc != collect_proc) {
-
- void* myhdata = rank==proc ? hdata : 0;
- gdata<dim>* const tmpdata = mydata->make_typed();
- tmpdata->allocate (alldata->extent(), proc, myhdata);
- tmpdata->copy_from (alldata, alldata->extent());
- delete tmpdata;
-
- }
- }
- } // Copy result
-
- delete alldata;
-
- } // if components>0
+ } // Copy result
- // Restore global state
- if (saved_component!=-1) {
- set_component ((cGH*)cgh, saved_component);
- }
+ delete alldata;
// Success
return hdata;
generated by cgit v1.2.3 (git 2.39.1) at 2024-07-27 22:48:24 +0000