aboutsummaryrefslogtreecommitdiff
path: root/Carpet/Carpet
diff options
context:
space:
mode:
authorErik Schnetter <schnetter@cct.lsu.edu>2009-09-03 16:19:15 -0500
committerBarry Wardell <barry.wardell@gmail.com>2011-12-14 16:42:31 +0000
commit11c4d98017cbb86d08e15fd1b549180184b58a26 (patch)
tree2546a154c6f7bc0bec87de7316125ae7d1453569 /Carpet/Carpet
parentf520477b1c14e02f1495cfa8d3e09f4e21ab34d0 (diff)
Import Carpet
Ignore-this: 309b4dd613f4af2b84aa5d6743fdb6b3
Diffstat (limited to 'Carpet/Carpet')
-rw-r--r--Carpet/Carpet/README6
-rw-r--r--Carpet/Carpet/configuration.ccl6
-rw-r--r--Carpet/Carpet/interface.ccl42
-rw-r--r--Carpet/Carpet/param.ccl52
-rw-r--r--Carpet/Carpet/schedule.ccl1
-rw-r--r--Carpet/Carpet/src/AllGatherString.cc69
-rw-r--r--Carpet/Carpet/src/CallFunction.cc28
-rw-r--r--Carpet/Carpet/src/CarpetBasegrid.cc8
-rw-r--r--Carpet/Carpet/src/CarpetParamCheck.cc8
-rw-r--r--Carpet/Carpet/src/CarpetStartup.cc12
-rw-r--r--Carpet/Carpet/src/Checksum.cc124
-rw-r--r--Carpet/Carpet/src/Comm.cc124
-rw-r--r--Carpet/Carpet/src/Cycle.cc40
-rw-r--r--Carpet/Carpet/src/Evolve.cc120
-rw-r--r--Carpet/Carpet/src/Initialise.cc112
-rw-r--r--Carpet/Carpet/src/Limits.cc93
-rw-r--r--Carpet/Carpet/src/MultiModel.cc61
-rw-r--r--Carpet/Carpet/src/OutputGH.cc18
-rw-r--r--Carpet/Carpet/src/Poison.cc152
-rw-r--r--Carpet/Carpet/src/Recompose.cc1037
-rw-r--r--Carpet/Carpet/src/Requirements.cc25
-rw-r--r--Carpet/Carpet/src/Restrict.cc24
-rw-r--r--Carpet/Carpet/src/ScheduleWrapper.cc8
-rw-r--r--Carpet/Carpet/src/SetupGH.cc437
-rw-r--r--Carpet/Carpet/src/Shutdown.cc23
-rw-r--r--Carpet/Carpet/src/Storage.cc85
-rw-r--r--Carpet/Carpet/src/Timers.cc34
-rw-r--r--Carpet/Carpet/src/Timers.hh31
-rw-r--r--Carpet/Carpet/src/Timing.cc119
-rw-r--r--Carpet/Carpet/src/carpet.hh14
-rw-r--r--Carpet/Carpet/src/carpet_public.h22
-rw-r--r--Carpet/Carpet/src/functions.hh59
-rw-r--r--Carpet/Carpet/src/helpers.cc108
-rw-r--r--Carpet/Carpet/src/make.code.defn4
-rw-r--r--Carpet/Carpet/src/modes.cc404
-rw-r--r--Carpet/Carpet/src/modes.hh87
-rw-r--r--Carpet/Carpet/src/variables.cc9
-rw-r--r--Carpet/Carpet/src/variables.hh24
-rw-r--r--Carpet/Carpet/test/small.par16
39 files changed, 2211 insertions, 1435 deletions
diff --git a/Carpet/Carpet/README b/Carpet/Carpet/README
index 03ce32952..75caae781 100644
--- a/Carpet/Carpet/README
+++ b/Carpet/Carpet/README
@@ -1,8 +1,10 @@
Cactus Code Thorn Carpet
-Authors : Erik Schnetter <schnetter@uni-tuebingen.de>
+Author(s) : Erik Schnetter <schnetter@cct.lsu.edu>
+Maintainer(s): Erik Schnetter <schnetter@cct.lsu.edu>
+Licence : GPLv2+
--------------------------------------------------------------------------
-Purpose of the thorn:
+1. Purpose
This thorn provides a parallel AMR (adaptive mesh refinement) driver
with MPI.
diff --git a/Carpet/Carpet/configuration.ccl b/Carpet/Carpet/configuration.ccl
index a80a319df..a7294fc09 100644
--- a/Carpet/Carpet/configuration.ccl
+++ b/Carpet/Carpet/configuration.ccl
@@ -2,10 +2,8 @@
PROVIDES Carpet
{
- SCRIPT
- LANG
}
-REQUIRES IOUtil CarpetLib
+REQUIRES IOUtil CarpetLib LoopControl
-REQUIRES THORNS: IOUtil CarpetLib
+REQUIRES THORNS: IOUtil CarpetLib LoopControl
diff --git a/Carpet/Carpet/interface.ccl b/Carpet/Carpet/interface.ccl
index 87c44e941..a61f22fc6 100644
--- a/Carpet/Carpet/interface.ccl
+++ b/Carpet/Carpet/interface.ccl
@@ -7,6 +7,10 @@ include header: carpet_public.h in carpet.h
include header: Timers.hh in CarpetTimers.hh
+uses include header: loopcontrol.h
+
+uses include header: mpi_string.hh
+
uses include header: defs.hh
uses include header: dist.hh
@@ -217,6 +221,21 @@ PROVIDES FUNCTION GetRefinementLevels \
+# Get pointer to grid variable for a specific map and refinement level
+CCTK_POINTER FUNCTION \
+ VarDataPtrI \
+ (CCTK_POINTER_TO_CONST IN cctkGH, \
+ CCTK_INT IN map, \
+ CCTK_INT IN reflevel, \
+ CCTK_INT IN component, \
+ CCTK_INT IN timelevel, \
+ CCTK_INT IN varindex)
+PROVIDES FUNCTION VarDataPtrI \
+ WITH Carpet_VarDataPtrI \
+ LANGUAGE C
+
+
+
# The true prototype of the routine below:
# int Carpet_Regrid (const cGH * cctkGH,
# gh::rregs * superregss,
@@ -249,11 +268,32 @@ REQUIRES FUNCTION IO_TruncateOutputFiles
+# TODO: make this somehow public, e.g. by moving it into its own thorn
+
+CCTK_INT point_classification TYPE=gf TAGS='checkpoint="no" prolongation="none"'
+{
+ point_class
+ # negative: needs to be set explicitly (e.g. boundary)
+ # zero: unused (e.g. ghost)
+ # positive: needs to be evolved
+ # -1: boundary point (needs to be set explicitly)
+ # 0: unused (e.g. ghost point, or restriction target)
+ # n=1..N: evolved, used for integrator substeps i<=n
+ # (i=N..1, counting backwards; see MoL documentation)
+ # i.e.: n=1: used between time steps (i.e., should be visualised)
+ # n>1: used only while time stepping (e.g. buffer zones)
+} "Grid point classification"
+
+
+
CCTK_REAL timing TAGS='checkpoint="no"'
{
physical_time_per_hour
+ current_physical_time_per_hour
+
+ time_total time_evolution time_computing time_communicating time_io
- time_total time_computing time_communicating time_io
+ evolution_steps_count
local_grid_points_per_second total_grid_points_per_second
local_grid_point_updates_count total_grid_point_updates_count
diff --git a/Carpet/Carpet/param.ccl b/Carpet/Carpet/param.ccl
index d4887cec7..d221d2158 100644
--- a/Carpet/Carpet/param.ccl
+++ b/Carpet/Carpet/param.ccl
@@ -171,21 +171,12 @@ STRING model "Model name for multi-model simulations -- the model name is used t
CCTK_INT prolongation_order_space "Order of prolongation operator in space" STEERABLE=recover
{
- 1 :: "first order (linear)"
- 3 :: "third order (cubic)"
- 5 :: "fifth order"
- 7 :: "seventh order"
- 9 :: "ninth order"
- 11 :: "eleventh order (one more than tenth)"
+ 0:* :: "vertex centred orders must be odd"
} 1
CCTK_INT prolongation_order_time "Order of prolongation operator in time" STEERABLE=recover
{
- 0 :: "zeroth order (constant)"
- 1 :: "first order (linear)"
- 2 :: "second order (quadratic)"
- 3 :: "third order (cubic)"
- 4 :: "fourth order (quartic)"
+ 0:* :: ""
} 1
@@ -291,21 +282,20 @@ BOOLEAN enable_all_storage "Enable storage for all grid functions" STEERABLE=rec
{
} "no"
+BOOLEAN enable_no_storage "Exit before beginning to enable storage for grid functions" STEERABLE=recover
+{
+} "no"
+
BOOLEAN poison_new_timelevels "Try to catch uninitialised grid elements by setting new timelevels to values that will catch your attention" STEERABLE=always
{
-} "no"
+} "yes"
BOOLEAN check_for_poison "Explicitely check for the poison value after every time step" STEERABLE=always
{
} "no"
-CCTK_INT poison_value "Integer value (0..255) used to poison new timelevels (with memset)" STEERABLE=always
-{
- 0:255 :: "Must fit into a byte. Use 0 for zero, 255 for nan, and e.g. 113 for a large value."
-} 255
-
CCTK_INT max_poison_locations "Maximum number of poison locations that are printed to the screen" STEERABLE=always
{
-1 :: "print all locations"
@@ -314,13 +304,6 @@ CCTK_INT max_poison_locations "Maximum number of poison locations that are print
-CCTK_INT deadbeef "A strange integer value that indicates that something has gone wrong; the integer equivalent of a nan"
-{
- *:* :: "should be large and positive"
-} 666 # 7353315
-
-
-
BOOLEAN checksum_timelevels "Try to catch unintentionally changed timelevels by taking checksums and comparing against these" STEERABLE=always
{
} "no"
@@ -357,6 +340,11 @@ BOOLEAN output_internal_data "Periodically print internal data to the screen for
{
} "no"
+REAL timing_average_window_minutes "Time interval (in wall time minutes) for calculating the current physics time per hour" STEERABLE=always
+{
+ (0.0:* :: ""
+} 10.0
+
INT print_timestats_every "Print interesting timing statistics periodically" STEERABLE=always
{
-1 :: "don't report"
@@ -377,21 +365,15 @@ STRING timer_file "File name in which detailed timing statistics are collected"
"^.+$" :: "file name"
} "carpet-timing-statistics"
-INT max_core_size_MB "Maximum size of a core file" STEERABLE=recover
+BOOLEAN timers_verbose "Output (debug) messages when a timer is started or stopped" STEERABLE=always
{
- -2 :: "unchanged"
- -1 :: "unlimited"
- 0:* :: "limited"
-} -2
+} "no"
-INT max_memory_size_MB "Maximum amount of memory per MPI process" STEERABLE=recover
-{
- -2 :: "unchanged"
- -1 :: "unlimited"
- 0:* :: "limited"
-} -2
+BOOLEAN recompose_verbose "Output debug information during recomposing" STEERABLE=ALWAYS
+{
+} "no"
KEYWORD processor_topology "How to determine the processor topology" STEERABLE=recover
{
diff --git a/Carpet/Carpet/schedule.ccl b/Carpet/Carpet/schedule.ccl
index 3c9322b5f..0f839a5c8 100644
--- a/Carpet/Carpet/schedule.ccl
+++ b/Carpet/Carpet/schedule.ccl
@@ -1,5 +1,6 @@
# Schedule definitions for thorn Carpet
+storage: point_classification
storage: timing timing2
schedule CarpetMultiModelStartup at STARTUP as MultiModel_Startup before Driver_Startup
diff --git a/Carpet/Carpet/src/AllGatherString.cc b/Carpet/Carpet/src/AllGatherString.cc
deleted file mode 100644
index 5efa96114..000000000
--- a/Carpet/Carpet/src/AllGatherString.cc
+++ /dev/null
@@ -1,69 +0,0 @@
-#include <cassert>
-#include <cstring>
-//#include <iostream>
-//#include <map>
-#include <string>
-#include <vector>
-
-#include <mpi.h>
-
-#include "cctk.h"
-
-#include "functions.hh"
-
-
-
-namespace Carpet
-{
-
- using namespace std;
-
-
-
- vector <string>
- AllGatherString (MPI_Comm const world,
- string const & data)
- {
- // Get the total number of processors
- int num_procs;
- MPI_Comm_size (world, & num_procs);
-
- // Gather the lengths of the data strings
- int const length = data.length();
- vector <int> lengths (num_procs);
-
- MPI_Allgather (const_cast <int *> (& length), 1, MPI_INT,
- & lengths.front(), 1, MPI_INT,
- world);
-
- // Allocate space for all data strings
- vector <int> offsets (num_procs + 1);
- offsets.at(0) = 0;
- for (int n = 0; n < num_procs; ++ n)
- {
- offsets.at(n + 1) = offsets.at(n) + lengths.at(n);
- }
- int const total_length = offsets.at(num_procs);
-
- // Gather all data strings
- vector <char> alldata_buffer (total_length);
-
- MPI_Allgatherv (const_cast <char *> (data.c_str()), length, MPI_CHAR,
- & alldata_buffer.front(),
- const_cast <int *> (& lengths.front()),
- const_cast <int *> (& offsets.front()),
- MPI_CHAR,
- world);
-
- // Convert data buffer with C strings to C++ strings
- vector <string> alldata (num_procs);
- for (int n = 0; n < num_procs; ++ n)
- {
- alldata.at(n) =
- string (& alldata_buffer.at (offsets.at(n)), lengths.at(n));
- }
-
- return alldata;
- }
-
-} // namespace Carpet
diff --git a/Carpet/Carpet/src/CallFunction.cc b/Carpet/Carpet/src/CallFunction.cc
index 696b5ce49..054ea6808 100644
--- a/Carpet/Carpet/src/CallFunction.cc
+++ b/Carpet/Carpet/src/CallFunction.cc
@@ -8,8 +8,8 @@
#include <gh.hh>
-#include "carpet.hh"
-#include "Timers.hh"
+#include <carpet.hh>
+#include <Timers.hh>
@@ -76,13 +76,13 @@ namespace Carpet {
BEGIN_META_MODE(cctkGH) {
BEGIN_MGLEVEL_LOOP(cctkGH) {
BEGIN_REFLEVEL_LOOP(cctkGH) {
- BEGIN_MAP_LOOP(cctkGH, CCTK_GF) {
+ BEGIN_LOCAL_MAP_LOOP(cctkGH, CCTK_GF) {
BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) {
CallScheduledFunction
("Meta time local mode",
function, attribute, data, user_timer);
} END_LOCAL_COMPONENT_LOOP;
- } END_MAP_LOOP;
+ } END_LOCAL_MAP_LOOP;
sync_timer.start();
SyncGroupsInScheduleBlock (attribute, cctkGH) ;
sync_timer.stop();
@@ -159,13 +159,13 @@ namespace Carpet {
if (attribute->loop_local) {
BEGIN_GLOBAL_MODE(cctkGH) {
BEGIN_REFLEVEL_LOOP(cctkGH) {
- BEGIN_MAP_LOOP(cctkGH, CCTK_GF) {
+ BEGIN_LOCAL_MAP_LOOP(cctkGH, CCTK_GF) {
BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) {
CallScheduledFunction
("Global time local mode",
function, attribute, data, user_timer);
} END_LOCAL_COMPONENT_LOOP;
- } END_MAP_LOOP;
+ } END_LOCAL_MAP_LOOP;
sync_timer.start();
SyncGroupsInScheduleBlock (attribute, cctkGH) ;
sync_timer.stop();
@@ -213,13 +213,13 @@ namespace Carpet {
// Level operation: call once per refinement level
if (attribute->loop_local) {
- BEGIN_MAP_LOOP(cctkGH, CCTK_GF) {
+ BEGIN_LOCAL_MAP_LOOP(cctkGH, CCTK_GF) {
BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) {
CallScheduledFunction
("Level time local mode",
function, attribute, data, user_timer);
} END_LOCAL_COMPONENT_LOOP;
- } END_MAP_LOOP;
+ } END_LOCAL_MAP_LOOP;
} else if (attribute->loop_singlemap) {
BEGIN_MAP_LOOP(cctkGH, CCTK_GF) {
CallScheduledFunction
@@ -239,13 +239,13 @@ namespace Carpet {
// Single map operation: call once per refinement level and map
if (attribute->loop_local) {
- BEGIN_MAP_LOOP(cctkGH, CCTK_GF) {
+ BEGIN_LOCAL_MAP_LOOP(cctkGH, CCTK_GF) {
BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) {
CallScheduledFunction
("Singlemap time local mode",
function, attribute, data, user_timer);
} END_LOCAL_COMPONENT_LOOP;
- } END_MAP_LOOP;
+ } END_LOCAL_MAP_LOOP;
} else {
BEGIN_MAP_LOOP(cctkGH, CCTK_GF) {
CallScheduledFunction
@@ -260,13 +260,13 @@ namespace Carpet {
} else {
// Local operation: call once per component
- BEGIN_MAP_LOOP(cctkGH, CCTK_GF) {
+ BEGIN_LOCAL_MAP_LOOP(cctkGH, CCTK_GF) {
BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) {
CallScheduledFunction
("Local mode",
function, attribute, data, user_timer);
} END_LOCAL_COMPONENT_LOOP;
- } END_MAP_LOOP;
+ } END_LOCAL_MAP_LOOP;
sync_timer.start();
SyncGroupsInScheduleBlock (attribute, cctkGH) ;
sync_timer.stop();
@@ -274,6 +274,7 @@ namespace Carpet {
}
if (schedule_barriers) {
+#if 0
static unsigned int magic = 0xe8932329UL; // a random starting value
unsigned int recv = magic;
Checkpoint ("About to Bcast %u", magic);
@@ -287,6 +288,9 @@ namespace Carpet {
Checkpoint ("About to Barrier");
MPI_Barrier (dist::comm());
Checkpoint ("Finished Barrier");
+#endif
+ static int id = 513400912; // arbitrary starting value
+ CCTK_NamedBarrier (NULL, id++);
}
total_timer.stop();
diff --git a/Carpet/Carpet/src/CarpetBasegrid.cc b/Carpet/Carpet/src/CarpetBasegrid.cc
index c2dc45f81..7e3ab4ae0 100644
--- a/Carpet/Carpet/src/CarpetBasegrid.cc
+++ b/Carpet/Carpet/src/CarpetBasegrid.cc
@@ -1,10 +1,10 @@
#include <limits>
-#include "cctk.h"
-#include "cctk_Arguments.h"
-#include "cctk_Parameters.h"
+#include <cctk.h>
+#include <cctk_Arguments.h>
+#include <cctk_Parameters.h>
-#include "carpet.hh"
+#include <carpet.hh>
diff --git a/Carpet/Carpet/src/CarpetParamCheck.cc b/Carpet/Carpet/src/CarpetParamCheck.cc
index fb9b40861..cf3d861ed 100644
--- a/Carpet/Carpet/src/CarpetParamCheck.cc
+++ b/Carpet/Carpet/src/CarpetParamCheck.cc
@@ -1,11 +1,11 @@
#include <cassert>
#include <cstdlib>
-#include "cctk.h"
-#include "cctk_Arguments.h"
-#include "cctk_Parameters.h"
+#include <cctk.h>
+#include <cctk_Arguments.h>
+#include <cctk_Parameters.h>
-#include "carpet.hh"
+#include <carpet.hh>
diff --git a/Carpet/Carpet/src/CarpetStartup.cc b/Carpet/Carpet/src/CarpetStartup.cc
index 822c03845..07677c80d 100644
--- a/Carpet/Carpet/src/CarpetStartup.cc
+++ b/Carpet/Carpet/src/CarpetStartup.cc
@@ -1,12 +1,12 @@
#include <cassert>
#include <cstdlib>
-#include "cctk.h"
-#include "cctk_Parameters.h"
+#include <cctk.h>
+#include <cctk_Parameters.h>
-#include "carpet.hh"
+#include <dist.hh>
-#include "dist.hh"
+#include <carpet.hh>
@@ -19,8 +19,11 @@ namespace Carpet {
DECLARE_CCTK_PARAMETERS;
comm_universe = MPI_COMM_WORLD;
+ // cerr << "QQQ: CarpetMultiModelStartup[1]" << endl;
SplitUniverse (comm_universe, model, comm_world, true);
+ // cerr << "QQQ: CarpetMultiModelStartup[2]" << endl;
dist::pseudoinit (comm_world);
+ // cerr << "QQQ: CarpetMultiModelStartup[3]" << endl;
return 0;
}
@@ -46,6 +49,7 @@ namespace Carpet {
CCTK_OverloadEnableGroupComm (EnableGroupComm);
CCTK_OverloadDisableGroupComm (DisableGroupComm);
CCTK_OverloadBarrier (Barrier);
+ CCTK_OverloadNamedBarrier (NamedBarrier);
CCTK_OverloadExit (Exit);
CCTK_OverloadAbort (Abort);
CCTK_OverloadMyProc (MyProc);
diff --git a/Carpet/Carpet/src/Checksum.cc b/Carpet/Carpet/src/Checksum.cc
index e65710038..139ba1555 100644
--- a/Carpet/Carpet/src/Checksum.cc
+++ b/Carpet/Carpet/src/Checksum.cc
@@ -2,12 +2,14 @@
#include <cstdlib>
#include <vector>
-#include "cctk.h"
-#include "cctk_Parameters.h"
+#include <cctk.h>
+#include <cctk_Parameters.h>
+#include <util_ErrorCodes.h>
+#include <util_Table.h>
-#include "gh.hh"
+#include <gh.hh>
-#include "carpet.hh"
+#include <carpet.hh>
@@ -20,7 +22,7 @@ namespace Carpet {
// Checksum information
struct ckdesc {
bool valid;
- unsigned int sum;
+ unsigned long sum;
};
// Helper class
@@ -33,6 +35,30 @@ namespace Carpet {
+ // Calculate the internet checksum (see RFC 1071)
+ static unsigned short
+ internet_checksum (void const * restrict const addr,
+ size_t const len)
+ {
+ unsigned long chk = 0;
+#pragma omp parallel for reduction (+: chk)
+ for (ptrdiff_t i=0; i<(ptrdiff_t)len; i+=2) {
+ unsigned long const lb = ((unsigned char const*)addr)[i];
+ unsigned long const ub = ((unsigned char const*)addr)[i+1];
+ chk += lb + (ub << 8);
+ }
+ if (len % 1) {
+ unsigned long const lb = ((unsigned char const*)addr)[len-1];
+ chk += lb;
+ }
+ while (chk >> 16) {
+ chk = (chk & 0xffffUL) + (chk >> 16);
+ }
+ return ~chk;
+ }
+
+
+
// The parameter where specifies which time levels should be
// poisoned. what specifies what kind of grid variables should be
// poisoned.
@@ -52,11 +78,11 @@ namespace Carpet {
const int grouptype = CCTK_GroupTypeI(group);
if (reflevel == 0 or grouptype == CCTK_GF) {
checksums.at(reflevel).at(mglevel).at(group).a.resize(arrdata.at(group).size());
- BEGIN_MAP_LOOP(cgh, grouptype) {
- checksums.at(reflevel).at(mglevel).at(group).a.at(map).resize(arrdata.at(group).at(map).hh->components(reflevel));
+ BEGIN_LOCAL_MAP_LOOP(cgh, grouptype) {
+ checksums.at(reflevel).at(mglevel).at(group).a.at(map).resize(arrdata.at(group).at(map).hh->local_components(reflevel));
BEGIN_LOCAL_COMPONENT_LOOP(cgh, grouptype) {
const int nvars = CCTK_NumVarsInGroupI(group);
- checksums.at(reflevel).at(mglevel).at(group).a.at(map).at(component).resize(nvars);
+ checksums.at(reflevel).at(mglevel).at(group).a.at(map).at(local_component).resize(nvars);
if (nvars > 0) {
const int n0 = CCTK_FirstVarIndexI(group);
@@ -70,30 +96,48 @@ namespace Carpet {
}
const int np = prod(size);
+ int const table = CCTK_GroupTagsTableI (group);
+ assert (table >= 0);
+ bool persistent;
+ char buf[100];
+ int const ilen =
+ Util_TableGetString (table, sizeof buf, buf, "Persistent");
+ if (ilen > 0) {
+ if (CCTK_EQUALS(buf, "yes")) {
+ persistent = true;
+ } else if (CCTK_EQUALS(buf, "no")) {
+ persistent = false;
+ } else {
+ assert (0);
+ }
+ } else if (ilen == UTIL_ERROR_TABLE_NO_SUCH_KEY) {
+ // default
+ persistent = true;
+ } else {
+ assert (0);
+ }
+
const int num_tl = CCTK_NumTimeLevelsFromVarI(n0);
assert (num_tl>0);
- const int min_tl = min_timelevel(where, num_tl);
- const int max_tl = max_timelevel(where, num_tl);
+ const int min_tl = min_timelevel(where, num_tl, persistent);
+ const int max_tl = max_timelevel(where, num_tl, persistent);
for (int var=0; var<nvars; ++var) {
- checksums.at(reflevel).at(mglevel).at(group).a.at(map).at(component).at(var).resize(num_tl);
+ checksums.at(reflevel).at(mglevel).at(group).a.at(map).at(local_component).at(var).resize(num_tl);
for (int tl=min_tl; tl<=max_tl; ++tl) {
const int n = n0 + var;
const void* data = cgh->data[n][tl];
- unsigned int chk = 0;
- for (int i=0; i<np*sz/(int)sizeof chk; ++i) {
- chk += ((const unsigned int*)data)[i];
- }
+ unsigned long const chk = internet_checksum (data, np*sz);
- checksums.at(reflevel).at(mglevel).at(group).a.at(map).at(component).at(var).at(tl).sum = chk;
- checksums.at(reflevel).at(mglevel).at(group).a.at(map).at(component).at(var).at(tl).valid = true;
+ checksums.at(reflevel).at(mglevel).at(group).a.at(map).at(local_component).at(var).at(tl).sum = chk;
+ checksums.at(reflevel).at(mglevel).at(group).a.at(map).at(local_component).at(var).at(tl).valid = true;
} // for tl
} // for var
} // if group has vars
} END_LOCAL_COMPONENT_LOOP;
- } END_MAP_LOOP;
+ } END_LOCAL_MAP_LOOP;
} // if grouptype fits
} // if storage
} // for group
@@ -118,11 +162,11 @@ namespace Carpet {
const int grouptype = CCTK_GroupTypeI(group);
if (reflevel == 0 or grouptype == CCTK_GF) {
assert (checksums.at(reflevel).at(mglevel).at(group).a.size()==arrdata.at(group).size());
- BEGIN_MAP_LOOP(cgh, grouptype) {
- assert ((int)checksums.at(reflevel).at(mglevel).at(group).a.at(map).size()==arrdata.at(group).at(map).hh->components(reflevel));
+ BEGIN_LOCAL_MAP_LOOP(cgh, grouptype) {
+ assert ((int)checksums.at(reflevel).at(mglevel).at(group).a.at(map).size()==arrdata.at(group).at(map).hh->local_components(reflevel));
BEGIN_LOCAL_COMPONENT_LOOP(cgh, grouptype) {
const int nvars = CCTK_NumVarsInGroupI(group);
- assert ((int)checksums.at(reflevel).at(mglevel).at(group).a.at(map).at(component).size()==nvars);
+ assert ((int)checksums.at(reflevel).at(mglevel).at(group).a.at(map).at(local_component).size()==nvars);
if (nvars > 0) {
const int n0 = CCTK_FirstVarIndexI(group);
@@ -136,24 +180,42 @@ namespace Carpet {
}
const int np = prod(size);
+ int const table = CCTK_GroupTagsTableI (group);
+ assert (table >= 0);
+ bool persistent;
+ char buf[100];
+ int const ilen =
+ Util_TableGetString (table, sizeof buf, buf, "Persistent");
+ if (ilen > 0) {
+ if (CCTK_EQUALS(buf, "yes")) {
+ persistent = true;
+ } else if (CCTK_EQUALS(buf, "no")) {
+ persistent = false;
+ } else {
+ assert (0);
+ }
+ } else if (ilen == UTIL_ERROR_TABLE_NO_SUCH_KEY) {
+ // default
+ persistent = true;
+ } else {
+ assert (0);
+ }
+
const int num_tl = CCTK_NumTimeLevelsFromVarI(n0);
assert (num_tl>0);
- const int min_tl = min_timelevel(where, num_tl);
- const int max_tl = max_timelevel(where, num_tl);
+ const int min_tl = min_timelevel(where, num_tl, persistent);
+ const int max_tl = max_timelevel(where, num_tl, persistent);
for (int var=0; var<nvars; ++var) {
- assert ((int)checksums.at(reflevel).at(mglevel).at(group).a.at(map).at(component).at(var).size()==num_tl);
+ assert ((int)checksums.at(reflevel).at(mglevel).at(group).a.at(map).at(local_component).at(var).size()==num_tl);
for (int tl=min_tl; tl<=max_tl; ++tl) {
- if (checksums.at(reflevel).at(mglevel).at(group).a.at(map).at(component).at(var).at(tl).valid) {
+ if (checksums.at(reflevel).at(mglevel).at(group).a.at(map).at(local_component).at(var).at(tl).valid) {
const int n = n0 + var;
const void* data = cgh->data[n][tl];
- unsigned int chk = 0;
- for (int i=0; i<np*sz/(int)sizeof chk; ++i) {
- chk += ((const unsigned int*)data)[i];
- }
+ unsigned long const chk = internet_checksum (data, np*sz);
const bool unexpected_change =
- chk != checksums.at(reflevel).at(mglevel).at(group).a.at(map).at(component).at(var).at(tl).sum;
+ chk != checksums.at(reflevel).at(mglevel).at(group).a.at(map).at(local_component).at(var).at(tl).sum;
if (unexpected_change) {
char* fullname = CCTK_FullName(n);
@@ -168,7 +230,7 @@ namespace Carpet {
} // for var
} // if group has vars
} END_LOCAL_COMPONENT_LOOP;
- } END_MAP_LOOP;
+ } END_LOCAL_MAP_LOOP;
} // if grouptype fits
} // if storage
} // for group
diff --git a/Carpet/Carpet/src/Comm.cc b/Carpet/Carpet/src/Comm.cc
index 674d7a40c..64f2f0480 100644
--- a/Carpet/Carpet/src/Comm.cc
+++ b/Carpet/Carpet/src/Comm.cc
@@ -4,13 +4,14 @@
#include <cstdlib>
#include <iostream>
-#include "cctk.h"
-#include "cctk_Parameters.h"
+#include <cctk.h>
+#include <cctk_Parameters.h>
-#include "ggf.hh"
-#include "gh.hh"
+#include <ggf.hh>
+#include <gh.hh>
-#include "carpet.hh"
+#include <carpet.hh>
+#include <Timers.hh>
@@ -75,10 +76,11 @@ namespace Carpet {
// check consistency of all groups:
// create a new set with empty and no-storage groups removed
vector<int> goodgroups;
- for (size_t g = 0; g < groups.size(); g++) {
- const int group = groups[g];
- const int grouptype = CCTK_GroupTypeI (group);
- char* groupname = CCTK_GroupName (group);
+ goodgroups.reserve (groups.size());
+ for (size_t group = 0; group < groups.size(); group++) {
+ const int g = groups.AT(group);
+ const int grouptype = CCTK_GroupTypeI (g);
+ char* const groupname = CCTK_GroupName (g);
Checkpoint ("SyncGroup \"%s\" time=%g",
groupname, (double) cctkGH->cctk_time);
@@ -102,7 +104,7 @@ namespace Carpet {
}
}
if (component != -1) {
- if (maps == 1 and vhh.at(map)->local_components(reflevel) == 1) {
+ if (maps == 1 and vhh.AT(map)->local_components(reflevel) == 1) {
CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
"Synchronising group \"%s\" in local mode",
groupname);
@@ -115,83 +117,39 @@ namespace Carpet {
}
}
- if (not CCTK_QueryGroupStorageI (cctkGH, group)) {
+ if (not CCTK_QueryGroupStorageI (cctkGH, g)) {
CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
"Cannot synchronise group \"%s\" because it has no storage",
groupname);
retval = -1;
}
- else if (CCTK_NumVarsInGroupI (group) > 0) {
- goodgroups.push_back(group);
+ else if (CCTK_NumVarsInGroupI (g) > 0) {
+ goodgroups.push_back(g);
}
free (groupname);
- }
+ } // for g
if (goodgroups.size() > 0) {
- bool local_do_prolongate;
- if (do_prolongate) {
- if (do_taper) {
- if (reflevel == 0) {
- local_do_prolongate = true;
- } else { // on a fine level
-#warning "TODO: Check iteration number instead -- this is wrong while regridding"
- CCTK_REAL mytime;
- CCTK_REAL parenttime;
- if (map == -1) {
- mytime = vtt.at(0)->time (0, reflevel, mglevel);
- parenttime = vtt.at(0)->time (0, reflevel - 1, mglevel);
- } else {
- mytime = vtt.at(map)->time (0, reflevel, mglevel);
- parenttime = vtt.at(map)->time (0, reflevel - 1, mglevel);
- }
- CCTK_REAL const eps = 1.0e-12;
- bool const in_sync =
- abs (mytime - parenttime) <= eps * abs (delta_time);
-#if 0
- int const parent_do_every =
- ipow(mgfact, mglevel) *
- (maxtimereflevelfact / timereffacts.at(reflevel-1));
- bool const parent_is_active =
- cctkGH->cctk_iteration == 0 or
- (cctkGH->cctk_iteration-1) % parent_do_every == 0;
- int const do_every =
- ipow(mgfact, mglevel) *
- (maxtimereflevelfact / timereffacts.at(reflevel));
- bool const is_active =
- cctkGH->cctk_iteration == 0 or
- (cctkGH->cctk_iteration-1) % do_every == 0;
- bool const new_in_sync = is_active and parent_is_active;
-#warning "just for testing"
-#warning "if this breaks, fix also CarpetRegrid2"
- assert (new_in_sync == in_sync);
- if (not (new_in_sync == in_sync)) {
- CCTK_WARN (CCTK_WARN_ABORT, "assert (new_in_sync == in_sync)");
- }
-#endif
- local_do_prolongate = in_sync;
- }
- } else { // no tapered grids
- local_do_prolongate = true;
- }
- } else { // not do_prolongate
- local_do_prolongate = false;
- }
-
// prolongate boundaries
+ bool const local_do_prolongate = do_prolongate and not do_taper;
if (local_do_prolongate) {
- if (reflevel > 0) {
- ProlongateGroupBoundaries (cctkGH, cctk_initial_time, goodgroups);
- }
+ static Timer timer ("Evolve::Prolongate");
+ timer.start();
+ ProlongateGroupBoundaries (cctkGH, cctk_initial_time, goodgroups);
+ timer.stop();
}
// synchronise ghostzones
if (sync_during_time_integration or local_do_prolongate) {
+ static Timer timer ("Evolve::Sync");
+ timer.start();
SyncGroups (cctkGH, goodgroups);
+ timer.stop();
}
- } // for g
+ }
return retval;
}
@@ -204,10 +162,16 @@ namespace Carpet {
DECLARE_CCTK_PARAMETERS;
const int tl = 0;
- // use the current time here (which may be modified by the user)
- const CCTK_REAL time
- = (cctkGH->cctk_time - initial_time) / delta_time;
+ if (reflevel == 0) return;
+
+ Checkpoint ("ProlongateGroups");
+ assert (groups.size() > 0);
+
+ // use the current time here (which may be modified by the user)
+ const CCTK_REAL time =
+ (cctkGH->cctk_time - initial_time) / delta_time;
+
for (comm_state state; not state.done(); state.step()) {
for (int group = 0; group < (int)groups.size(); ++group) {
const int g = groups.AT(group);
@@ -216,10 +180,10 @@ namespace Carpet {
continue;
}
assert (reflevel>=0 and reflevel<reflevels);
-
- for (int m = 0; m < (int)arrdata.at(g).size(); ++m) {
- for (int v = 0; v < (int)arrdata.at(g).at(m).data.size(); ++v) {
- ggf *const gv = arrdata.at(g).at(m).data.at(v);
+
+ for (int m = 0; m < (int)arrdata.AT(g).size(); ++m) {
+ for (int v = 0; v < (int)arrdata.AT(g).AT(m).data.size(); ++v) {
+ ggf *const gv = arrdata.AT(g).AT(m).data.AT(v);
gv->ref_bnd_prolongate_all (state, tl, reflevel, mglevel, time);
}
}
@@ -234,18 +198,20 @@ namespace Carpet {
DECLARE_CCTK_PARAMETERS;
const int tl = 0;
+ Checkpoint ("SyncGroups");
+
assert (groups.size() > 0);
for (comm_state state; not state.done(); state.step()) {
for (int group = 0; group < (int)groups.size(); ++group) {
- const int g = groups[group];
+ const int g = groups.AT(group);
const int grouptype = CCTK_GroupTypeI (g);
const int ml = grouptype == CCTK_GF ? mglevel : 0;
const int rl = grouptype == CCTK_GF ? reflevel : 0;
- for (int m = 0; m < (int)arrdata.at(g).size(); ++m) {
- for (int v = 0; v < (int)arrdata.at(g).at(m).data.size(); ++v) {
- arrdesc& array = arrdata.at(g).at(m);
- array.data.at(v)->sync_all (state, tl, rl, ml);
+ for (int m = 0; m < (int)arrdata.AT(g).size(); ++m) {
+ for (int v = 0; v < (int)arrdata.AT(g).AT(m).data.size(); ++v) {
+ arrdesc& array = arrdata.AT(g).AT(m);
+ array.data.AT(v)->sync_all (state, tl, rl, ml);
}
}
}
diff --git a/Carpet/Carpet/src/Cycle.cc b/Carpet/Carpet/src/Cycle.cc
index 1a2931dd0..889f305f1 100644
--- a/Carpet/Carpet/src/Cycle.cc
+++ b/Carpet/Carpet/src/Cycle.cc
@@ -1,12 +1,12 @@
#include <cassert>
#include <cstdlib>
-#include "cctk.h"
+#include <cctk.h>
-#include "ggf.hh"
-#include "gh.hh"
+#include <ggf.hh>
+#include <gh.hh>
-#include "carpet.hh"
+#include <carpet.hh>
@@ -27,9 +27,9 @@ namespace Carpet {
case CCTK_GF:
assert (reflevel>=0 and reflevel<reflevels);
- for (int m=0; m<(int)arrdata.at(group).size(); ++m) {
+ for (int m=0; m<(int)arrdata.AT(group).size(); ++m) {
for (int var=0; var<CCTK_NumVarsInGroupI(group); ++var) {
- arrdata.at(group).at(m).data.at(var)->
+ arrdata.AT(group).AT(m).data.AT(var)->
cycle_all (reflevel, mglevel);
}
}
@@ -41,15 +41,14 @@ namespace Carpet {
int const numtimelevels = CCTK_NumTimeLevelsI (group);
int const firstvarindex = CCTK_FirstVarIndexI (group);
for (int var=0; var<CCTK_NumVarsInGroupI(group); ++var) {
- arrdata.at(group).at(0).data.at(var)->cycle_all (0, mglevel);
+ arrdata.AT(group).AT(0).data.AT(var)->cycle_all (0, mglevel);
{
int const varindex = firstvarindex + var;
- int const c = CCTK_MyProc(cgh);
for (int tl=0; tl<numtimelevels; ++tl) {
cgh->data[varindex][tl]
- = (tl < groupdata.at(group).info.activetimelevels
- ? ((*arrdata.at(group).at(0).data.at(var))
- (tl, 0, c, 0)->storage())
+ = (tl < groupdata.AT(group).info.activetimelevels
+ ? ((*arrdata.AT(group).AT(0).data.AT(var))
+ (tl, 0, 0, 0)->storage())
: NULL);
}
}
@@ -77,9 +76,9 @@ namespace Carpet {
case CCTK_GF:
assert (reflevel>=0 and reflevel<reflevels);
- for (int m=0; m<(int)arrdata.at(group).size(); ++m) {
+ for (int m=0; m<(int)arrdata.AT(group).size(); ++m) {
for (int var=0; var<CCTK_NumVarsInGroupI(group); ++var) {
- arrdata.at(group).at(m).data.at(var)->
+ arrdata.AT(group).AT(m).data.AT(var)->
flip_all (reflevel, mglevel);
}
}
@@ -91,15 +90,14 @@ namespace Carpet {
int const numtimelevels = CCTK_NumTimeLevelsI (group);
int const firstvarindex = CCTK_FirstVarIndexI (group);
for (int var=0; var<CCTK_NumVarsInGroupI(group); ++var) {
- arrdata.at(group).at(0).data.at(var)->flip_all (0, mglevel);
+ arrdata.AT(group).AT(0).data.AT(var)->flip_all (0, mglevel);
{
int const varindex = firstvarindex + var;
- int const c = CCTK_MyProc(cgh);
for (int tl=0; tl<numtimelevels; ++tl) {
cgh->data[varindex][tl]
- = (tl < groupdata.at(group).info.activetimelevels
- ? ((*arrdata.at(group).at(0).data.at(var))
- (tl, 0, c, 0)->storage())
+ = (tl < groupdata.AT(group).info.activetimelevels
+ ? ((*arrdata.AT(group).AT(0).data.AT(var))
+ (tl, 0, 0, 0)->storage())
: NULL);
}
}
@@ -127,9 +125,9 @@ namespace Carpet {
case CCTK_GF:
assert (reflevel>=0 and reflevel<reflevels);
- for (int m=0; m<(int)arrdata.at(group).size(); ++m) {
+ for (int m=0; m<(int)arrdata.AT(group).size(); ++m) {
for (int var=0; var<CCTK_NumVarsInGroupI(group); ++var) {
- arrdata.at(group).at(m).data.at(var)->
+ arrdata.AT(group).AT(m).data.AT(var)->
fill_all (reflevel, mglevel);
}
}
@@ -139,7 +137,7 @@ namespace Carpet {
case CCTK_ARRAY:
if (do_global_mode) {
for (int var=0; var<CCTK_NumVarsInGroupI(group); ++var) {
- arrdata.at(group).at(0).data.at(var)->fill_all (0, mglevel);
+ arrdata.AT(group).AT(0).data.AT(var)->fill_all (0, mglevel);
}
}
break;
diff --git a/Carpet/Carpet/src/Evolve.cc b/Carpet/Carpet/src/Evolve.cc
index 5b9d8ff96..103221088 100644
--- a/Carpet/Carpet/src/Evolve.cc
+++ b/Carpet/Carpet/src/Evolve.cc
@@ -16,8 +16,8 @@
#include <dist.hh>
#include <th.hh>
-#include "carpet.hh"
-#include "Timers.hh"
+#include <carpet.hh>
+#include <Timers.hh>
@@ -53,27 +53,21 @@ namespace Carpet {
int const convlev = 0;
cGH* cctkGH = fc->GH[convlev];
- // Tapered grids
- do_taper = use_tapered_grids;
-
// Main loop
- BeginTiming (cctkGH);
+ BeginTimingEvolution (cctkGH);
static Timer timer ("Evolve");
timer.start();
while (not do_terminate (cctkGH)) {
AdvanceTime (cctkGH);
{
- int const do_every = maxtimereflevelfact / timereffacts.at(reflevels-1);
+ int const do_every = maxtimereflevelfact / timereffacts.AT(reflevels-1);
if ((cctkGH->cctk_iteration - 1) % do_every == 0) {
- bool const old_do_taper = do_taper;
- do_taper = false;
ENTER_GLOBAL_MODE (cctkGH, 0) {
BEGIN_REFLEVEL_LOOP (cctkGH) {
CallRegrid (cctkGH);
} END_REFLEVEL_LOOP;
} LEAVE_GLOBAL_MODE;
- do_taper = old_do_taper;
}
}
CallEvol (cctkGH);
@@ -83,7 +77,7 @@ namespace Carpet {
// Print timer values
{
- int const do_every = maxtimereflevelfact / timereffacts.at(reflevels-1);
+ int const do_every = maxtimereflevelfact / timereffacts.AT(reflevels-1);
if (output_timers_every > 0 and
cctkGH->cctk_iteration % output_timers_every == 0 and
cctkGH->cctk_iteration % do_every == 0)
@@ -99,9 +93,9 @@ namespace Carpet {
for (int ml=0; ml<mglevels; ++ml) {
for (int rl=0; rl<reflevels; ++rl) {
int const do_every =
- ipow (mgfact, ml) * (maxtimereflevelfact / timereffacts.at(rl));
+ ipow (mgfact, ml) * (maxtimereflevelfact / timereffacts.AT(rl));
if (cctkGH->cctk_iteration % do_every == 0) {
- assert (abs (leveltimes.at(ml).at(rl) - global_time) <=
+ assert (abs (leveltimes.AT(ml).AT(rl) - global_time) <=
eps * global_time);
}
}
@@ -111,8 +105,6 @@ namespace Carpet {
} // end main loop
timer.stop();
- do_taper = false;
-
Waypoint ("Done with evolution loop");
return 0;
@@ -132,7 +124,7 @@ namespace Carpet {
// Do not test on non-active reflevels to save the call to
// MPI_Allreduce below
- int const do_every = maxtimereflevelfact / timereffacts.at(reflevels-1);
+ int const do_every = maxtimereflevelfact / timereffacts.AT(reflevels-1);
if (cctkGH->cctk_iteration % do_every != 0)
{
@@ -228,7 +220,7 @@ namespace Carpet {
}
if ((cctkGH->cctk_iteration-1)
- % (maxtimereflevelfact / timereffacts.at(reflevels-1)) == 0) {
+ % (maxtimereflevelfact / timereffacts.AT(reflevels-1)) == 0) {
Waypoint ("Evolving iteration %d at t=%g",
cctkGH->cctk_iteration, (double)cctkGH->cctk_time);
}
@@ -273,25 +265,31 @@ namespace Carpet {
// Regrid
Checkpoint ("Regrid");
int const oldreflevels = reflevels;
- bool const did_regrid = Regrid (cctkGH, false);
+ bool const did_regrid = Regrid (cctkGH, false, true);
bool const did_remove_level = reflevels < oldreflevels;
assert (not did_remove_level or did_regrid);
if (did_regrid) {
+ bool did_any_recompose = false;
BEGIN_META_MODE (cctkGH) {
for (int rl=0; rl<reflevels; ++rl) {
bool const did_recompose = Recompose (cctkGH, rl, true);
+ did_any_recompose = did_any_recompose or did_recompose;
// Do not omit the global mode call when the finest level
// does not change:
// if (did_recompose or (did_remove_level and rl == reflevels - 1)) {
- if (did_recompose or rl == reflevels - 1) {
+ if (did_recompose or
+ ((did_remove_level or did_any_recompose) and
+ rl == reflevels - 1))
+ {
BEGIN_MGLEVEL_LOOP (cctkGH) {
ENTER_LEVEL_MODE (cctkGH, rl) {
do_early_global_mode = reflevel==0;
do_late_global_mode = reflevel==reflevels-1;
- do_early_meta_mode = do_early_global_mode and mglevel==mglevels-1;
+ do_early_meta_mode =
+ do_early_global_mode and mglevel==mglevels-1;
do_late_meta_mode = do_late_global_mode and mglevel==0;
do_global_mode = do_late_global_mode;
do_meta_mode = do_late_meta_mode;
@@ -309,17 +307,17 @@ namespace Carpet {
// Rewind times
for (int m=0; m<maps; ++m) {
- vtt.at(m)->set_delta
+ vtt.AT(m)->set_delta
(reflevel, mglevel,
- - vtt.at(m)->get_delta (reflevel, mglevel));
+ - vtt.AT(m)->get_delta (reflevel, mglevel));
FlipTimeLevels (cctkGH);
for (int tl=0; tl<num_tl; ++tl) {
- vtt.at(m)->advance_time (reflevel, mglevel);
+ vtt.AT(m)->advance_time (reflevel, mglevel);
CycleTimeLevels (cctkGH);
}
- vtt.at(m)->set_delta
+ vtt.AT(m)->set_delta
(reflevel, mglevel,
- - vtt.at(m)->get_delta (reflevel, mglevel));
+ - vtt.AT(m)->get_delta (reflevel, mglevel));
FlipTimeLevels (cctkGH);
} // for m
CCTK_REAL const old_cctk_time = cctkGH->cctk_time;
@@ -330,7 +328,7 @@ namespace Carpet {
// Advance times
for (int m=0; m<maps; ++m) {
- vtt.at(m)->advance_time (reflevel, mglevel);
+ vtt.AT(m)->advance_time (reflevel, mglevel);
}
CycleTimeLevels (cctkGH);
cctkGH->cctk_time +=
@@ -352,6 +350,8 @@ namespace Carpet {
} END_META_MODE;
} // if did_regrid
+ RegridFree (cctkGH, true);
+
do_global_mode = old_do_global_mode;
do_early_global_mode = old_do_early_global_mode;
do_late_global_mode = old_do_late_global_mode;
@@ -380,7 +380,7 @@ namespace Carpet {
for (int rl=0; rl<reflevels; ++rl) {
int const do_every
- = ipow(mgfact, ml) * (maxtimereflevelfact / timereffacts.at(rl));
+ = ipow(mgfact, ml) * (maxtimereflevelfact / timereffacts.AT(rl));
if ((cctkGH->cctk_iteration-1) % do_every == 0) {
ENTER_GLOBAL_MODE (cctkGH, ml) {
ENTER_LEVEL_MODE (cctkGH, rl) {
@@ -395,6 +395,15 @@ namespace Carpet {
have_done_global_mode |= do_global_mode;
have_done_anything = true;
+ if (use_tapered_grids and reflevel > 0) {
+ int const parent_do_every =
+ ipow(mgfact, mglevel) *
+ (maxtimereflevelfact / timereffacts.AT(reflevel-1));
+ bool const parent_is_active =
+ (cctkGH->cctk_iteration-1) % parent_do_every == 0;
+ do_taper = not parent_is_active;
+ }
+
// Advance times
cctkGH->cctk_time
= (global_time
@@ -402,7 +411,7 @@ namespace Carpet {
+ delta_time * mglevelfact / timereflevelfact);
CCTK_REAL const carpet_time = cctkGH->cctk_time / delta_time;
for (int m=0; m<maps; ++m) {
- vtt.at(m)->advance_time (reflevel, mglevel);
+ vtt.AT(m)->advance_time (reflevel, mglevel);
if (not adaptive_stepsize) {
#if 0
// We must not perform this check, since the
@@ -413,7 +422,7 @@ namespace Carpet {
static_assert (abs(0.1) > 0,
"Function CarpetLib::abs has wrong signature");
CCTK_REAL const level_time =
- vtt.at(m)->get_time (reflevel, mglevel);
+ vtt.AT(m)->get_time (reflevel, mglevel);
if (not (abs (level_time - carpet_time) <=
eps * max (carpet_time, 1.0))) {
int const oldprecision = cerr.precision();
@@ -429,15 +438,16 @@ namespace Carpet {
assert (abs (level_time - carpet_time) <=
eps * max (carpet_time, 1.0));
#endif
- vtt.at(m)->set_time (reflevel, mglevel, carpet_time);
+ vtt.AT(m)->set_time (reflevel, mglevel, carpet_time);
}
}
CycleTimeLevels (cctkGH);
- Waypoint ("Evolution I at iteration %d time %g%s%s",
+ Waypoint ("Evolution I at iteration %d time %g%s%s%s",
cctkGH->cctk_iteration, (double)cctkGH->cctk_time,
(do_global_mode ? " (global)" : ""),
- (do_meta_mode ? " (meta)" : ""));
+ (do_meta_mode ? " (meta)" : ""),
+ (do_taper ? " (tapering)" : ""));
// Checking
CalculateChecksums (cctkGH, allbutcurrenttime);
@@ -451,7 +461,9 @@ namespace Carpet {
PoisonCheck (cctkGH, currenttime);
// Timing statistics
- StepTiming (cctkGH);
+ StepTimingEvolution (cctkGH);
+
+ do_taper = false;
} LEAVE_LEVEL_MODE;
} LEAVE_GLOBAL_MODE;
@@ -476,7 +488,7 @@ namespace Carpet {
for (int ml=mglevels-1; ml>=0; --ml) {
for (int rl=reflevels-1; rl>=0; --rl) {
int const do_every
- = ipow(mgfact, ml) * (maxtimereflevelfact / timereffacts.at(rl));
+ = ipow(mgfact, ml) * (maxtimereflevelfact / timereffacts.AT(rl));
if (cctkGH->cctk_iteration % do_every == 0) {
ENTER_GLOBAL_MODE (cctkGH, ml) {
ENTER_LEVEL_MODE (cctkGH, rl) {
@@ -513,7 +525,7 @@ namespace Carpet {
for (int rl=0; rl<reflevels; ++rl) {
int const do_every
- = ipow(mgfact, ml) * (maxtimereflevelfact / timereffacts.at(rl));
+ = ipow(mgfact, ml) * (maxtimereflevelfact / timereffacts.AT(rl));
if (cctkGH->cctk_iteration % do_every == 0) {
ENTER_GLOBAL_MODE (cctkGH, ml) {
ENTER_LEVEL_MODE (cctkGH, rl) {
@@ -528,10 +540,20 @@ namespace Carpet {
have_done_global_mode |= do_global_mode;
have_done_anything = true;
- Waypoint ("Evolution II at iteration %d time %g%s%s",
+ if (use_tapered_grids and reflevel > 0) {
+ int const parent_do_every =
+ ipow(mgfact, mglevel) *
+ (maxtimereflevelfact / timereffacts.AT(reflevel-1));
+ bool const parent_is_active =
+ (cctkGH->cctk_iteration-1) % parent_do_every == 0;
+ do_taper = not parent_is_active;
+ }
+
+ Waypoint ("Evolution II at iteration %d time %g%s%s%s",
cctkGH->cctk_iteration, (double)cctkGH->cctk_time,
(do_global_mode ? " (global)" : ""),
- (do_meta_mode ? " (meta)" : ""));
+ (do_meta_mode ? " (meta)" : ""),
+ (do_taper ? " (tapering)" : ""));
if (reflevel < reflevels-1) {
ScheduleTraverse (where, "CCTK_POSTRESTRICT", cctkGH);
@@ -566,6 +588,8 @@ namespace Carpet {
PrintTimingStats (cctkGH);
}
+ do_taper = false;
+
} LEAVE_LEVEL_MODE;
} LEAVE_GLOBAL_MODE;
} // if do_every
@@ -607,15 +631,27 @@ namespace Carpet {
void ScheduleTraverse (char const * const where, char const * const name,
cGH * const cctkGH)
{
+ // Obtain the set of timers, creating it explicitly if it does not
+ // yet exist
+ typedef std::map <string, Timer *> timers_t;
+ // static timers_t timers;
+ static timers_t * timersp = NULL;
+ if (not timersp) timersp = new timers_t;
+ timers_t & timers = * timersp;
+
+ // Obtain timer, creating a new one if it does not yet exist
ostringstream timernamebuf;
timernamebuf << where << "::" << name;
string const timername = timernamebuf.str();
- static std::map <string, Timer *> timers;
- Timer * & mapped = timers[timername];
- if (not mapped) {
- mapped = new Timer (timername.c_str());
+ timers_t::iterator ti = timers.find (timername);
+ if (ti == timers.end()) {
+ pair <string, Timer *> const
+ newtimer (timername, new Timer (timername.c_str()));
+ ti = timers.insert(newtimer).first;
+ // It is possible to find and insert with the same function
+ // call, but this makes the code significantly more complicated
}
- Timer & timer = * mapped;
+ Timer & timer = * ti->second;
timer.start();
ostringstream infobuf;
diff --git a/Carpet/Carpet/src/Initialise.cc b/Carpet/Carpet/src/Initialise.cc
index 836be8635..94bd24521 100644
--- a/Carpet/Carpet/src/Initialise.cc
+++ b/Carpet/Carpet/src/Initialise.cc
@@ -11,8 +11,8 @@
#include <cctki_ScheduleBindings.h>
#include <cctki_WarnLevel.h>
-#include "carpet.hh"
-#include "Timers.hh"
+#include <carpet.hh>
+#include <Timers.hh>
@@ -64,8 +64,8 @@ namespace Carpet {
global_time = cctk_initial_time;
delta_time = 1.0;
for (int ml = 0; ml < mglevel; ++ ml) {
- assert (leveltimes.at(ml).size() == 1);
- leveltimes.at(ml).at(0) = global_time;
+ assert (leveltimes.AT(ml).size() == 1);
+ leveltimes.AT(ml).AT(0) = global_time;
}
cctkGH->cctk_iteration = 0;
@@ -87,8 +87,8 @@ namespace Carpet {
#if 0
// Write grid structure to file
for (int m=0; m<maps; ++m) {
- OutputGridStructure (cctkGH, m, vhh.at(m)->regions);
- OutputGridCoordinates (cctkGH, m, vhh.at(m)->regions);
+ OutputGridStructure (cctkGH, m, vhh.AT(m)->regions);
+ OutputGridCoordinates (cctkGH, m, vhh.AT(m)->regions);
} // for m
#endif
@@ -332,15 +332,15 @@ namespace Carpet {
not CCTK_EQUALS (initial_data_setup_method, "init_single_level");
for (int m=0; m<maps; ++m) {
- vtt.at(m)->set_delta
- (reflevel, mglevel, - vtt.at(m)->get_delta (reflevel, mglevel));
+ vtt.AT(m)->set_delta
+ (reflevel, mglevel, - vtt.AT(m)->get_delta (reflevel, mglevel));
FlipTimeLevels (cctkGH);
for (int tl=0; tl<num_tl; ++tl) {
- vtt.at(m)->advance_time (reflevel, mglevel);
+ vtt.AT(m)->advance_time (reflevel, mglevel);
CycleTimeLevels (cctkGH);
}
- vtt.at(m)->set_delta
- (reflevel, mglevel, - vtt.at(m)->get_delta (reflevel, mglevel));
+ vtt.AT(m)->set_delta
+ (reflevel, mglevel, - vtt.AT(m)->get_delta (reflevel, mglevel));
FlipTimeLevels (cctkGH);
} // for m
@@ -348,7 +348,7 @@ namespace Carpet {
// Advance times
for (int m=0; m<maps; ++m) {
- vtt.at(m)->advance_time (reflevel, mglevel);
+ vtt.AT(m)->advance_time (reflevel, mglevel);
}
cctkGH->cctk_time =
(+ global_time
@@ -468,7 +468,7 @@ namespace Carpet {
(do_meta_mode ? " (meta)" : ""));
int const do_every =
- ipow(mgfact, mglevel) * (maxtimereflevelfact / timereffacts.at(rl));
+ ipow(mgfact, mglevel) * (maxtimereflevelfact / timereffacts.AT(rl));
if (cctkGH->cctk_iteration % do_every == 0)
{
// Checkpoint
@@ -580,17 +580,17 @@ namespace Carpet {
// Rewind times
for (int m=0; m<maps; ++m) {
- vtt.at(m)->set_delta
+ vtt.AT(m)->set_delta
(reflevel, mglevel,
- - vtt.at(m)->get_delta (reflevel, mglevel));
+ - vtt.AT(m)->get_delta (reflevel, mglevel));
FlipTimeLevels (cctkGH);
for (int tl=0; tl<num_tl; ++tl) {
- vtt.at(m)->advance_time (reflevel, mglevel);
+ vtt.AT(m)->advance_time (reflevel, mglevel);
CycleTimeLevels (cctkGH);
}
- vtt.at(m)->set_delta
+ vtt.AT(m)->set_delta
(reflevel, mglevel,
- - vtt.at(m)->get_delta (reflevel, mglevel));
+ - vtt.AT(m)->get_delta (reflevel, mglevel));
FlipTimeLevels (cctkGH);
} // for m
CCTK_REAL const old_cctk_time = cctkGH->cctk_time;
@@ -601,7 +601,7 @@ namespace Carpet {
// Advance times
for (int m=0; m<maps; ++m) {
- vtt.at(m)->advance_time (reflevel, mglevel);
+ vtt.AT(m)->advance_time (reflevel, mglevel);
}
CycleTimeLevels (cctkGH);
cctkGH->cctk_time +=
@@ -627,6 +627,10 @@ namespace Carpet {
} END_META_MODE;
} // if did_regrid
+ if (callregrid) {
+ RegridFree (cctkGH);
+ }
+
do_global_mode = old_do_global_mode;
do_meta_mode = old_do_meta_mode;
}
@@ -670,15 +674,15 @@ namespace Carpet {
// Rewind times
for (int m=0; m<maps; ++m) {
- vtt.at(m)->set_delta
- (reflevel, mglevel, - vtt.at(m)->get_delta (reflevel, mglevel));
+ vtt.AT(m)->set_delta
+ (reflevel, mglevel, - vtt.AT(m)->get_delta (reflevel, mglevel));
FlipTimeLevels (cctkGH);
for (int tl=0; tl<num_tl; ++tl) {
- vtt.at(m)->advance_time (reflevel, mglevel);
+ vtt.AT(m)->advance_time (reflevel, mglevel);
CycleTimeLevels (cctkGH);
}
- vtt.at(m)->set_delta
- (reflevel, mglevel, - vtt.at(m)->get_delta (reflevel, mglevel));
+ vtt.AT(m)->set_delta
+ (reflevel, mglevel, - vtt.AT(m)->get_delta (reflevel, mglevel));
FlipTimeLevels (cctkGH);
} // for m
CCTK_REAL const old_cctk_time = cctkGH->cctk_time;
@@ -689,7 +693,7 @@ namespace Carpet {
// Advance times
for (int m=0; m<maps; ++m) {
- vtt.at(m)->advance_time (reflevel, mglevel);
+ vtt.AT(m)->advance_time (reflevel, mglevel);
}
CycleTimeLevels (cctkGH);
cctkGH->cctk_time +=
@@ -752,7 +756,7 @@ namespace Carpet {
// Regrid
Checkpoint ("Regrid");
- bool const did_regrid = Regrid (cctkGH, true);
+ bool const did_regrid = Regrid (cctkGH, true, prolongate_initial_data);
if (did_regrid) {
BEGIN_META_MODE (cctkGH) {
@@ -786,17 +790,17 @@ namespace Carpet {
// Rewind times
for (int m=0; m<maps; ++m) {
- vtt.at(m)->set_delta
+ vtt.AT(m)->set_delta
(reflevel, mglevel,
- - vtt.at(m)->get_delta (reflevel, mglevel));
+ - vtt.AT(m)->get_delta (reflevel, mglevel));
FlipTimeLevels (cctkGH);
for (int tl=0; tl<num_tl; ++tl) {
- vtt.at(m)->advance_time (reflevel, mglevel);
+ vtt.AT(m)->advance_time (reflevel, mglevel);
CycleTimeLevels (cctkGH);
}
- vtt.at(m)->set_delta
+ vtt.AT(m)->set_delta
(reflevel, mglevel,
- - vtt.at(m)->get_delta (reflevel, mglevel));
+ - vtt.AT(m)->get_delta (reflevel, mglevel));
FlipTimeLevels (cctkGH);
} // for m
CCTK_REAL const old_cctk_time = cctkGH->cctk_time;
@@ -807,7 +811,7 @@ namespace Carpet {
// Advance times
for (int m=0; m<maps; ++m) {
- vtt.at(m)->advance_time (reflevel, mglevel);
+ vtt.AT(m)->advance_time (reflevel, mglevel);
}
CycleTimeLevels (cctkGH);
cctkGH->cctk_time +=
@@ -829,6 +833,8 @@ namespace Carpet {
} END_META_MODE;
} // if did_regrid
+ RegridFree (cctkGH, prolongate_initial_data);
+
do_global_mode = old_do_global_mode;
do_early_global_mode = old_do_early_global_mode;
do_late_global_mode = old_do_late_global_mode;
@@ -877,14 +883,18 @@ namespace Carpet {
// Regrid
Checkpoint ("Regrid");
- bool const did_regrid = Regrid (cctkGH, true);
+ bool const did_regrid = Regrid (cctkGH, true, prolongate_initial_data);
if (did_regrid) {
for (int rl=0; rl<reflevels; ++rl) {
- Recompose (cctkGH, rl, prolongate_initial_data);
+ if (not enable_no_storage) {
+ Recompose (cctkGH, rl, prolongate_initial_data);
+ }
} // for rl
} // if did_regrid
+ RegridFree (cctkGH, prolongate_initial_data);
+
} LEAVE_LEVEL_MODE;
} LEAVE_GLOBAL_MODE;
@@ -936,7 +946,7 @@ namespace Carpet {
// Regrid
Checkpoint ("Regrid");
- bool const did_regrid = Regrid (cctkGH, true);
+ bool const did_regrid = Regrid (cctkGH, true, prolongate_initial_data);
if (did_regrid) {
BEGIN_META_MODE (cctkGH) {
@@ -971,17 +981,17 @@ namespace Carpet {
// Rewind times
for (int m=0; m<maps; ++m) {
- vtt.at(m)->set_delta
+ vtt.AT(m)->set_delta
(reflevel, mglevel,
- - vtt.at(m)->get_delta (reflevel, mglevel));
+ - vtt.AT(m)->get_delta (reflevel, mglevel));
FlipTimeLevels (cctkGH);
for (int tl=0; tl<num_tl; ++tl) {
- vtt.at(m)->advance_time (reflevel, mglevel);
+ vtt.AT(m)->advance_time (reflevel, mglevel);
CycleTimeLevels (cctkGH);
}
- vtt.at(m)->set_delta
+ vtt.AT(m)->set_delta
(reflevel, mglevel,
- - vtt.at(m)->get_delta (reflevel, mglevel));
+ - vtt.AT(m)->get_delta (reflevel, mglevel));
FlipTimeLevels (cctkGH);
} // for m
CCTK_REAL const old_cctk_time = cctkGH->cctk_time;
@@ -992,7 +1002,7 @@ namespace Carpet {
// Advance times
for (int m=0; m<maps; ++m) {
- vtt.at(m)->advance_time (reflevel, mglevel);
+ vtt.AT(m)->advance_time (reflevel, mglevel);
}
CycleTimeLevels (cctkGH);
cctkGH->cctk_time +=
@@ -1014,6 +1024,8 @@ namespace Carpet {
} END_META_MODE;
} // if did_regrid
+ RegridFree (cctkGH, prolongate_initial_data);
+
do_global_mode = old_do_global_mode;
do_early_global_mode = old_do_early_global_mode;
do_late_global_mode = old_do_late_global_mode;
@@ -1051,8 +1063,20 @@ namespace Carpet {
Waypoint ("Initialising three timelevels");
// TODO: ensure that there are 3 timelevels
+ if (not (prolongation_order_time == 2)) {
+ CCTK_WARN (CCTK_WARN_ABORT,
+ "The 3 timelevel initialisation scheme (init_3_timelevels=yes) requires 3 timelevels and a temporal prolongation order of 2 (prolongation_order_time=2)");
+ }
assert (prolongation_order_time == 2);
+ for (int rl=0; rl<int(timereffacts.size()); ++rl) {
+ if (not (timereffacts.AT(rl) == ipow (2, rl))) {
+ CCTK_WARN (CCTK_WARN_ABORT,
+ "The 3 timelevel initialisation scheme (init_3_timelevels=yes) requires temporal refinement factors of 2 for all refinement levels (time_refinement_factors[rl]=pow(2,rl))");
+ }
+ assert (timereffacts.AT(rl) == ipow (2, rl));
+ }
+
BEGIN_MGLEVEL_LOOP(cctkGH) {
BEGIN_REFLEVEL_LOOP(cctkGH) {
do_early_global_mode = reflevel==0;
@@ -1120,7 +1144,7 @@ namespace Carpet {
cctkGH->cctk_time
= global_time + delta_time * mglevelfact / timereflevelfact;
for (int m=0; m<maps; ++m) {
- vtt.at(m)->advance_time (reflevel, mglevel);
+ vtt.AT(m)->advance_time (reflevel, mglevel);
}
CycleTimeLevels (cctkGH);
@@ -1236,7 +1260,7 @@ namespace Carpet {
cctkGH->cctk_time
= global_time + 2 * delta_time * mglevelfact / timereflevelfact;
for (int m=0; m<maps; ++m) {
- vtt.at(m)->advance_time (reflevel, mglevel);
+ vtt.AT(m)->advance_time (reflevel, mglevel);
}
CycleTimeLevels (cctkGH);
@@ -1275,7 +1299,7 @@ namespace Carpet {
cctkGH->cctk_time = global_time;
for (int m=0; m<maps; ++m) {
- vtt.at(m)->set_time (reflevel, mglevel, 0);
+ vtt.AT(m)->set_time (reflevel, mglevel, 0);
}
}
diff --git a/Carpet/Carpet/src/Limits.cc b/Carpet/Carpet/src/Limits.cc
deleted file mode 100644
index 47b92f945..000000000
--- a/Carpet/Carpet/src/Limits.cc
+++ /dev/null
@@ -1,93 +0,0 @@
-#include <cctk.h>
-#include <cctk_Parameters.h>
-
-#include <algorithm>
-#include <cassert>
-#include <iostream>
-#include <sys/resource.h>
-
-#include "defs.hh"
-
-namespace Carpet {
-
- using namespace std;
-
-
-
- static
- void
- set_limit (int resource, char const * name, CCTK_INT value);
-
- static
- ostream &
- operator<< (ostream & s, struct rlimit const & limit);
-
- static
- void
- output (ostream & s, rlim_t const & value);
-
-
-
- void
- SetSystemLimits ()
- {
- DECLARE_CCTK_PARAMETERS;
- set_limit (RLIMIT_CORE, "core file size", max_core_size_MB);
- set_limit (RLIMIT_AS, "memory size", max_memory_size_MB);
- }
-
-
-
- void
- set_limit (int const resource, char const * const name, CCTK_INT const value)
- {
- struct rlimit limit;
- check (not getrlimit (resource, & limit));
-
- if (value == -2 ) {
- // Only show limit
- cout << "Current " << name << " limit: " << limit << endl;
- return;
- }
-
- cout << "Old " << name << " limit: " << limit << endl;
-
- if (value == -1) {
- limit.rlim_cur = limit.rlim_max;
- } else {
- limit.rlim_cur = min ((rlim_t) value * 1024 * 1024, limit.rlim_max);
- }
-
- check (not setrlimit (resource, & limit));
- check (not getrlimit (resource, & limit));
-
- cout << "New " << name << " limit: " << limit << endl;
- }
-
-
-
- static
- ostream &
- operator<< (ostream & s, struct rlimit const & limit)
- {
- s << "hard=";
- output (s, limit.rlim_max);
- s << ", soft=";
- output (s, limit.rlim_cur);
- return s;
- }
-
-
-
- static
- void
- output (ostream & s, rlim_t const & value)
- {
- if (value == RLIM_INFINITY) {
- s << "[unlimited]";
- } else {
- s << (value / (CCTK_REAL) (1024*1024)) << " MB";
- }
- }
-
-} // namespace Carpet
diff --git a/Carpet/Carpet/src/MultiModel.cc b/Carpet/Carpet/src/MultiModel.cc
index 975e48aec..c59eb44aa 100644
--- a/Carpet/Carpet/src/MultiModel.cc
+++ b/Carpet/Carpet/src/MultiModel.cc
@@ -7,9 +7,10 @@
#include <mpi.h>
-#include "cctk.h"
+#include <cctk.h>
-#include "functions.hh"
+#include <functions.hh>
+#include <mpi_string.hh>
@@ -22,8 +23,8 @@ namespace Carpet
vector <string> models; // Model id to model name
std::map <string, int> model_map; // Model name to model id
- vector <int> model_ids; // Processor to model id
- vector <vector <int> > model_procs; // Model id to processors
+ vector <int> model_ids; // Process to model id
+ vector <vector <int> > model_procs; // Model id to processes
@@ -65,14 +66,14 @@ namespace Carpet
SplitUniverse (MPI_Comm const world, string const model, MPI_Comm & comm,
bool const verbose)
{
- // Get the total number of processors
+ // Get the total number of processes
int num_procs;
MPI_Comm_size (world, & num_procs);
int my_proc;
MPI_Comm_rank (world, & my_proc);
// Gather all model names
- models = AllGatherString (world, model);
+ models = allgather_string (world, model);
// Map model strings to small integers
int num_models = 0;
@@ -80,81 +81,81 @@ namespace Carpet
model_map.clear ();
for (int n = 0; n < num_procs; ++ n)
{
- if (model_map.find (models.at(n)) != model_map.end())
+ if (model_map.find (models.AT(n)) != model_map.end())
{
- model_ids.at(n) = model_map[models.at(n)];
+ model_ids.AT(n) = model_map[models.AT(n)];
}
else
{
- model_map[models.at(n)] = num_models;
- model_ids.at(n) = num_models;
+ model_map[models.AT(n)] = num_models;
+ model_ids.AT(n) = num_models;
++ num_models;
}
}
- // Determine processors per model
+ // Determine processes per model
vector <int> num_model_procs (num_models, 0);
for (int n = 0; n < num_procs; ++ n)
{
- ++ num_model_procs.at (model_ids.at(n));
+ ++ num_model_procs.at (model_ids.AT(n));
}
model_procs.resize (num_models);
for (int m = 0; m < num_models; ++ m)
{
- model_procs.at(m).reserve (num_model_procs.at(m));
+ model_procs.AT(m).reserve (num_model_procs.AT(m));
}
for (int n = 0; n < num_procs; ++ n)
{
- model_procs.at (model_ids.at(n)).push_back (n);
+ model_procs.at (model_ids.AT(n)).push_back (n);
}
for (int m = 0; m < num_models; ++ m)
{
- assert (static_cast<int> (model_procs.at(m).size())
- == num_model_procs.at(m));
+ assert (static_cast<int> (model_procs.AT(m).size())
+ == num_model_procs.AT(m));
}
// Create a new communicator for each model
- MPI_Comm_split (world, model_ids.at(my_proc), my_proc, & comm);
+ MPI_Comm_split (world, model_ids.AT(my_proc), my_proc, & comm);
if (verbose)
{
CCTK_INFO ("Multi-Model listing:");
for (int m = 0; m < num_models; ++ m)
{
- cout << " model " << m << ": \"" << models.at(m) << "\"" << endl;
+ cout << " model " << m << ": \"" << models.AT(m) << "\"" << endl;
}
- CCTK_INFO ("Multi-Model processor distribution:");
+ CCTK_INFO ("Multi-Model process distribution:");
for (int n = 0; n < num_procs; ++ n)
{
- int const m = model_ids.at(n);
+ int const m = model_ids.AT(n);
bool const same_model_as_prev =
- n-1 >= 0 and model_ids.at(n-1) == m;
+ n-1 >= 0 and model_ids.AT(n-1) == m;
bool const same_model_as_next =
- n+1 < num_procs and model_ids.at(n+1) == m;
+ n+1 < num_procs and model_ids.AT(n+1) == m;
if (same_model_as_next) {
if (same_model_as_prev) {
// Output nothing
} else {
- // This processor has the same model as the next one:
+ // This process has the same model as the next one:
// output only a partial line
- cout << " processors " << n << "-";
+ cout << " processes " << n << "-";
}
} else {
if (same_model_as_prev) {
- // This processor has the same model as the previous one:
+ // This process has the same model as the previous one:
// finish a partial line
cout << n << ": "
- << "model " << m << " \"" << models.at(m) << "\"" << endl;
+ << "model " << m << " \"" << models.AT(m) << "\"" << endl;
} else {
- cout << " processor " << n << ": "
- << "model " << m << " \"" << models.at(m) << "\"" << endl;
+ cout << " process " << n << ": "
+ << "model " << m << " \"" << models.AT(m) << "\"" << endl;
}
}
}
- int const my_model = model_ids.at(my_proc);
+ int const my_model = model_ids.AT(my_proc);
CCTK_VInfo (CCTK_THORNSTRING,
- "Multi-Model: This is processor %d, model %d \"%s\"",
+ "Multi-Model: This is process %d, model %d \"%s\"",
my_proc, my_model, model.c_str());
}
}
diff --git a/Carpet/Carpet/src/OutputGH.cc b/Carpet/Carpet/src/OutputGH.cc
index 11b611754..52a877390 100644
--- a/Carpet/Carpet/src/OutputGH.cc
+++ b/Carpet/Carpet/src/OutputGH.cc
@@ -3,13 +3,13 @@
#include <cstdlib>
#include <sstream>
-#include "cctk.h"
-#include "cctk_Parameters.h"
+#include <cctk.h>
+#include <cctk_Parameters.h>
-#include "dist.hh"
+#include <dist.hh>
-#include "carpet.hh"
-#include "Timers.hh"
+#include <carpet.hh>
+#include <Timers.hh>
@@ -42,18 +42,18 @@ namespace Carpet {
IOMethod const * const method = CCTK_IOMethod (handle);
assert (method);
- if (not timers.at(handle)) {
+ if (not timers.AT(handle)) {
ostringstream buf;
buf << "OutputGH"
<< "::" << method->implementation
<< "::" << method->name
<< " [" << handle << "]";
- timers.at(handle) = new Timer (buf.str().c_str());
+ timers.AT(handle) = new Timer (buf.str().c_str());
}
- timers.at(handle)->start();
+ timers.AT(handle)->start();
num_vars += method->OutputGH (cctkGH);
- timers.at(handle)->stop();
+ timers.AT(handle)->stop();
} // for handle
diff --git a/Carpet/Carpet/src/Poison.cc b/Carpet/Carpet/src/Poison.cc
index 2811b42e2..83b4650b5 100644
--- a/Carpet/Carpet/src/Poison.cc
+++ b/Carpet/Carpet/src/Poison.cc
@@ -2,10 +2,14 @@
#include <cstdlib>
#include <cstring>
-#include "cctk.h"
-#include "cctk_Parameters.h"
+#include <cctk.h>
+#include <cctk_Parameters.h>
+#include <util_ErrorCodes.h>
+#include <util_Table.h>
-#include "carpet.hh"
+#include <defs.hh>
+
+#include <carpet.hh>
@@ -18,21 +22,25 @@ namespace Carpet {
// The parameter where specifies which time levels should be
// poisoned. what specifies what kind of grid variables should be
// poisoned.
- void Poison (const cGH* cgh, const checktimes where, const int what)
+ void
+ Poison (cGH const * const cctkGH,
+ checktimes const where, int const what)
{
DECLARE_CCTK_PARAMETERS;
- if (! poison_new_timelevels) return;
+ assert (what == 0 or what == CCTK_GF or what == CCTK_ARRAY);
+
+ if (not poison_new_timelevels) return;
for (int group=0; group<CCTK_NumGroups(); ++group) {
- if (CCTK_QueryGroupStorageI(cgh, group)) {
+ if (CCTK_QueryGroupStorageI(cctkGH, group)) {
int const grouptype = CCTK_GroupTypeI (group);
if (what == 0 or
(what == CCTK_GF and grouptype == CCTK_GF) or
(what == CCTK_ARRAY and (grouptype == CCTK_ARRAY or
grouptype == CCTK_SCALAR)))
{
- PoisonGroup (cgh, group, where);
+ PoisonGroup (cctkGH, group, where);
}
} // if has storage
} // for group
@@ -40,34 +48,56 @@ namespace Carpet {
- void PoisonGroup (const cGH* cgh, const int group, const checktimes where)
+ void
+ PoisonGroup (cGH const * const cctkGH,
+ int const group, checktimes const where)
{
DECLARE_CCTK_PARAMETERS;
assert (group>=0 and group<CCTK_NumGroups());
- if (! poison_new_timelevels) return;
+ if (not poison_new_timelevels) return;
- if (! CCTK_QueryGroupStorageI(cgh, group)) {
+ if (not CCTK_QueryGroupStorageI(cctkGH, group)) {
char * const groupname = CCTK_GroupName(group);
CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
- "Cannot poison group \"%s\" because it has no storage",
- groupname);
+ "Cannot poison group \"%s\" because it has no storage",
+ groupname);
free (groupname);
return;
}
- const int nvar = CCTK_NumVarsInGroupI(group);
+ int const nvar = CCTK_NumVarsInGroupI(group);
if (nvar == 0) return;
- const int n0 = CCTK_FirstVarIndexI(group);
+ int const n0 = CCTK_FirstVarIndexI(group);
assert (n0>=0);
- const int sz = CCTK_VarTypeSize(CCTK_VarTypeI(n0));
+ int const sz = CCTK_VarTypeSize(CCTK_VarTypeI(n0));
assert (sz>0);
- const int num_tl = CCTK_ActiveTimeLevelsVI(cgh, n0);
+ int const table = CCTK_GroupTagsTableI (group);
+ assert (table >= 0);
+ bool persistent;
+ char buf[100];
+ int const ilen = Util_TableGetString (table, sizeof buf, buf, "Persistent");
+ if (ilen > 0) {
+ if (CCTK_EQUALS(buf, "yes")) {
+ persistent = true;
+ } else if (CCTK_EQUALS(buf, "no")) {
+ persistent = false;
+ } else {
+ assert (0);
+ }
+ } else if (ilen == UTIL_ERROR_TABLE_NO_SUCH_KEY) {
+ // default
+ persistent = true;
+ } else {
+ assert (0);
+ }
+
+ int const num_tl = CCTK_ActiveTimeLevelsVI(cctkGH, n0);
assert (num_tl>0);
- const int min_tl = min_timelevel(where, num_tl);
- const int max_tl = max_timelevel(where, num_tl);
+ int const min_tl = min_timelevel(where, num_tl, persistent);
+ int const max_tl = max_timelevel(where, num_tl, persistent);
if (min_tl <= max_tl) {
@@ -77,76 +107,102 @@ namespace Carpet {
free (groupname);
}
- const int grouptype = CCTK_GroupTypeI(group);
+ CCTK_INT const poison_value = get_poison_value();
- BEGIN_MAP_LOOP(cgh, grouptype) {
- BEGIN_LOCAL_COMPONENT_LOOP(cgh, grouptype) {
+ int const grouptype = CCTK_GroupTypeI(group);
+
+ BEGIN_LOCAL_MAP_LOOP(cctkGH, grouptype) {
+ BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, grouptype) {
ivect size(1);
- const int gpdim = groupdata.at(group).info.dim;
+ int const gpdim = groupdata.AT(group).info.dim;
for (int d=0; d<gpdim; ++d) {
- size[d] = groupdata.at(group).info.lsh[d];
+ size[d] = groupdata.AT(group).info.lsh[d];
}
- const int np = prod(size);
+ int const np = prod(size);
for (int var=0; var<nvar; ++var) {
- const int n = n0 + var;
+ int const n = n0 + var;
for (int tl=min_tl; tl<=max_tl; ++tl) {
- memset (cgh->data[n][tl], poison_value, np*sz);
+ memset (cctkGH->data[n][tl], poison_value, np*sz);
} // for tl
} // for var
} END_LOCAL_COMPONENT_LOOP;
- } END_MAP_LOOP;
+ } END_LOCAL_MAP_LOOP;
} // if tl
}
- void PoisonCheck (const cGH* cgh, const checktimes where)
+ void
+ PoisonCheck (cGH const * const cctkGH, checktimes const where)
{
DECLARE_CCTK_PARAMETERS;
- if (! check_for_poison) return;
+ if (not check_for_poison) return;
Checkpoint ("PoisonCheck");
for (int group=0; group<CCTK_NumGroups(); ++group) {
- const int nvar = CCTK_NumVarsInGroupI(group);
- if (nvar > 0 && CCTK_QueryGroupStorageI(cgh, group)) {
+ int const nvar = CCTK_NumVarsInGroupI(group);
+ if (nvar > 0 and CCTK_QueryGroupStorageI(cctkGH, group)) {
- const int grouptype = CCTK_GroupTypeI(group);
- const int n0 = CCTK_FirstVarIndexI(group);
+ int const grouptype = CCTK_GroupTypeI(group);
+ int const n0 = CCTK_FirstVarIndexI(group);
assert (n0>=0);
- const int tp = CCTK_VarTypeI(n0);
- const int gpdim = groupdata.at(group).info.dim;
+ int const tp = CCTK_VarTypeI(n0);
+ int const gpdim = groupdata.AT(group).info.dim;
+
+ int const table = CCTK_GroupTagsTableI (group);
+ assert (table >= 0);
+ bool persistent;
+ char buf[100];
+ int const ilen =
+ Util_TableGetString (table, sizeof buf, buf, "Persistent");
+ if (ilen > 0) {
+ if (CCTK_EQUALS(buf, "yes")) {
+ persistent = true;
+ } else if (CCTK_EQUALS(buf, "no")) {
+ persistent = false;
+ } else {
+ assert (0);
+ }
+ } else if (ilen == UTIL_ERROR_TABLE_NO_SUCH_KEY) {
+ // default
+ persistent = true;
+ } else {
+ assert (0);
+ }
+
+ CCTK_INT const poison_value = get_poison_value();
- const int num_tl = CCTK_ActiveTimeLevelsVI(cgh, n0);
+ int const num_tl = CCTK_ActiveTimeLevelsVI(cctkGH, n0);
assert (num_tl>0);
- const int min_tl = min_timelevel(where, num_tl);
- const int max_tl = max_timelevel(where, num_tl);
+ int const min_tl = min_timelevel(where, num_tl, persistent);
+ int const max_tl = max_timelevel(where, num_tl, persistent);
- BEGIN_MAP_LOOP(cgh, grouptype) {
- BEGIN_LOCAL_COMPONENT_LOOP(cgh, grouptype) {
+ BEGIN_LOCAL_MAP_LOOP(cctkGH, grouptype) {
+ BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, grouptype) {
ivect size(1);
for (int d=0; d<gpdim; ++d) {
- size[d] = groupdata.at(group).info.lsh[d];
+ size[d] = groupdata.AT(group).info.lsh[d];
}
- const int np = prod(size);
+ int const np = prod(size);
for (int var=0; var<nvar; ++var) {
- const int n = n0 + var;
+ int const n = n0 + var;
for (int tl=min_tl; tl<=max_tl; ++tl) {
- const void* const data = cgh->data[n][tl];
+ const void* const data = cctkGH->data[n][tl];
int numpoison=0;
for (int k=0; k<size[2]; ++k) {
for (int j=0; j<size[1]; ++j) {
for (int i=0; i<size[0]; ++i) {
- const int idx = i + size[0] * (j + size[1] * k);
+ int const idx = i + size[0] * (j + size[1] * k);
bool poisoned=false;
switch (tp) {
#define TYPECASE(N,T) \
@@ -169,7 +225,7 @@ namespace Carpet {
char* fullname = CCTK_FullName(n);
CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
"At iteration %d: timelevel %d, component %d, map %d, refinement level %d of the variable \"%s\" contains poison at [%d,%d,%d]",
- cgh->cctk_iteration,
+ cctkGH->cctk_iteration,
tl, component, map, reflevel,
fullname, i,j,k);
free (fullname);
@@ -182,7 +238,7 @@ namespace Carpet {
char* fullname = CCTK_FullName(n);
CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
"At iteration %d: timelevel %d, component %d, map %d, refinement level %d of the variable \"%s\" contains poison at %d of %d locations; not all locations were printed",
- cgh->cctk_iteration,
+ cctkGH->cctk_iteration,
tl, component, map, reflevel,
fullname, numpoison, np);
free (fullname);
@@ -195,7 +251,7 @@ namespace Carpet {
} // for tl
} // for var
} END_LOCAL_COMPONENT_LOOP;
- } END_MAP_LOOP;
+ } END_LOCAL_MAP_LOOP;
} // if has storage
} // for group
diff --git a/Carpet/Carpet/src/Recompose.cc b/Carpet/Carpet/src/Recompose.cc
index 958a5d058..f1f0b8834 100644
--- a/Carpet/Carpet/src/Recompose.cc
+++ b/Carpet/Carpet/src/Recompose.cc
@@ -15,21 +15,25 @@
#include <sys/stat.h>
#include <sys/types.h>
-#include "cctk.h"
-#include "cctk_Parameters.h"
+#include <mpi.h>
-#include "bbox.hh"
-#include "bboxset.hh"
-#include "defs.hh"
-#include "dh.hh"
-#include "gh.hh"
-#include "region.hh"
-#include "vect.hh"
+#include <cctk.h>
+#include <cctk_Parameters.h>
-#include "carpet.hh"
-#include "modes.hh"
+#include <loopcontrol.h>
-#define DEBUG false // false or true
+#include <bbox.hh>
+#include <bboxset.hh>
+#include <defs.hh>
+#include <dh.hh>
+#include <gh.hh>
+#include <region.hh>
+#include <vect.hh>
+
+#include <carpet.hh>
+#include <modes.hh>
+#include <variables.hh>
+#include <Timers.hh>
@@ -61,6 +65,11 @@ namespace Carpet {
static void
+ ClassifyPoints (cGH const * cctkGH, int rl);
+
+
+
+ static void
SplitRegionsMaps_Automatic_Recursively (bvect const & dims,
int const firstproc,
int const nprocs,
@@ -76,57 +85,68 @@ namespace Carpet {
void
CheckRegions (gh::mregs const & regsss)
{
+ char const * const where = "Carpet::CheckRegions";
+ static Timer timer (where);
+ timer.start();
+
// At least one multigrid level
if (regsss.size() == 0) {
CCTK_WARN (0, "I cannot set up a grid hierarchy with zero multigrid levels.");
}
assert (regsss.size() > 0);
// At least one refinement level
- if (regsss.at(0).size() == 0) {
+ if (regsss.AT(0).size() == 0) {
CCTK_WARN (0, "I cannot set up a grid hierarchy with zero refinement levels.");
}
- assert (regsss.at(0).size() > 0);
+ assert (regsss.AT(0).size() > 0);
// At most maxreflevels levels
- if ((int)regsss.at(0).size() > maxreflevels) {
+ if ((int)regsss.AT(0).size() > maxreflevels) {
CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
"I cannot set up a grid hierarchy with more than Carpet::max_refinement_levels refinement levels. I found Carpet::max_refinement_levels=%d, while %d levels were requested.",
- (int)maxreflevels, (int)regsss.at(0).size());
+ (int)maxreflevels, (int)regsss.AT(0).size());
}
- assert ((int)regsss.at(0).size() <= maxreflevels);
+ assert ((int)regsss.AT(0).size() <= maxreflevels);
for (int ml=0; ml<(int)regsss.size(); ++ml) {
int num_regions = 0;
- for (int rl=0; rl<(int)regsss.at(0).size(); ++rl) {
+ for (int rl=0; rl<(int)regsss.AT(0).size(); ++rl) {
// No empty levels
// (but allow some empty maps)
- // assert (regsss.at(ml).at(rl).size() > 0);
- num_regions += regsss.at(ml).at(rl).size();
- for (int c=0; c<(int)regsss.at(ml).at(rl).size(); ++c) {
+ // assert (regsss.AT(ml).AT(rl).size() > 0);
+ num_regions += regsss.AT(ml).AT(rl).size();
+ for (int c=0; c<(int)regsss.AT(ml).AT(rl).size(); ++c) {
// Check sizes
// (but allow processors with zero grid points)
- // assert (all(regsss.at(rl).at(c).at(ml).extent.lower() <=
- // regsss.at(rl).at(c).at(ml).extent.upper()));
+ // assert (all(regsss.AT(rl).AT(c).AT(ml).extent.lower() <=
+ // regsss.AT(rl).AT(c).AT(ml).extent.upper()));
// Check strides
ivect const str =
- (maxspacereflevelfact / spacereffacts.at(rl) * ipow(mgfact, ml));
- assert (all(regsss.at(ml).at(rl).at(c).extent.stride() % str == 0));
+ (maxspacereflevelfact / spacereffacts.AT(rl) * ipow(mgfact, ml));
+ assert (all(regsss.AT(ml).AT(rl).AT(c).extent.stride() % str == 0));
// Check alignments
- assert (all(regsss.at(ml).at(rl).at(c).extent.lower() % str == 0));
- assert (all(regsss.at(ml).at(rl).at(c).extent.upper() % str == 0));
+ assert (all(regsss.AT(ml).AT(rl).AT(c).extent.lower() % str == 0));
+ assert (all(regsss.AT(ml).AT(rl).AT(c).extent.upper() % str == 0));
}
}
// No empty levels
assert (num_regions > 0);
}
+
+ timer.stop();
}
bool
Regrid (cGH const * const cctkGH,
- bool const force_recompose)
+ bool const force_recompose,
+ bool const do_init)
{
DECLARE_CCTK_PARAMETERS;
+ char const * const where = "Carpet::Regrid";
+ static Timer timer (where);
+ timer.start();
+
Checkpoint ("Regridding level %d...", reflevel);
assert (is_level_mode());
@@ -141,6 +161,7 @@ namespace Carpet {
CCTK_WARN (2, "The regridding routine Carpet_Regrid has not been provided. There will be no regridding. Maybe you forgot to activate a regridding thorn?");
didtell = true;
}
+ timer.stop();
return false;
}
@@ -150,7 +171,7 @@ namespace Carpet {
CCTK_WARN (2, "The regridding routine Carpet_RegridMaps has not been provided. Regridding will be performed in singlemap mode instead of level mode.");
didtell = true;
}
- }
+ }
@@ -160,7 +181,7 @@ namespace Carpet {
BEGIN_MAP_LOOP(cctkGH, CCTK_GF) {
- gh::rregs superregss = vhh.at(map)->superregions;
+ gh::rregs superregss = vhh.AT(map)->superregions;
gh::mregs regsss;
// Check whether to recompose
@@ -170,7 +191,7 @@ namespace Carpet {
did_change = did_change or do_recompose;
if (do_recompose) {
- RegridMap (cctkGH, map, superregss, regsss);
+ RegridMap (cctkGH, map, superregss, regsss, do_init);
}
} END_MAP_LOOP;
@@ -180,20 +201,28 @@ namespace Carpet {
vector<gh::rregs> superregsss (maps);
vector<gh::mregs> regssss (maps);
BEGIN_MAP_LOOP(cctkGH, CCTK_GF) {
- superregsss.at(map) = vhh.at(map)->superregions;
+ superregsss.AT(map) = vhh.AT(map)->superregions;
} END_MAP_LOOP;
// Check whether to recompose
CCTK_INT const do_recompose =
Carpet_RegridMaps (cctkGH, & superregsss, & regssss, force_recompose);
assert (do_recompose >= 0);
+#warning "TODO"
+#if 1 // #ifdef CARPET_DEBUG
+ {
+ int ival = do_recompose;
+ MPI_Bcast (& ival, 1, MPI_INT, 0, dist::comm());
+ assert (ival == do_recompose);
+ }
+#endif
did_change = did_change or do_recompose;
if (do_recompose) {
BEGIN_MAP_LOOP(cctkGH, CCTK_GF) {
- gh::rregs const & superregss = superregsss.at(map);
- gh::mregs const & regsss = regssss.at(map);
- RegridMap (cctkGH, map, superregss, regsss);
+ gh::rregs const & superregss = superregsss.AT(map);
+ gh::mregs const & regsss = regssss.AT(map);
+ RegridMap (cctkGH, map, superregss, regsss, do_init);
} END_MAP_LOOP;
}
@@ -208,6 +237,14 @@ namespace Carpet {
} // if did change
+
+ timer.stop();
+
+ if (enable_no_storage) {
+ CCTK_WARN (CCTK_WARN_ALERT,
+ "Carpet completed its internal setup, and would now normally go on to allocate memory. Since the parameter Carpet::enable_no_storage has been set, Carpet will exit instead.");
+ CCTK_Exit (cctkGH, 0);
+ }
return did_change;
}
@@ -218,12 +255,16 @@ namespace Carpet {
RegridMap (cGH const * const cctkGH,
int const m,
gh::rregs const & superregss,
- gh::mregs const & regsss)
+ gh::mregs const & regsss,
+ bool const do_init)
{
DECLARE_CCTK_PARAMETERS;
+ char const * const where = "Carpet::RegridMap";
+ static Timer timer (where);
+ timer.start();
+
Waypoint ("Regridding map %d...", m);
-
#warning "TODO: keep levels fixed here"
#if 0
@@ -239,15 +280,15 @@ namespace Carpet {
int const do_every =
use_tapered_grids ?
- maxtimereflevelfact / timereffacts.at(max(0,rl-1)):
- maxtimereflevelfact / timereffacts.at( rl );
+ maxtimereflevelfact / timereffacts.AT(max(0,rl-1)):
+ maxtimereflevelfact / timereffacts.AT( rl );
bool const regrid_this_level =
(cctkGH->cctk_iteration - 1) % do_every == 0;
if (not regrid_this_level) {
// Set regions from current grid structure
- regions.at(rl) = ...;
+ regions.AT(rl) = ...;
}
}
#endif
@@ -258,18 +299,18 @@ namespace Carpet {
// not change
// Regrid
- vhh.at(m)->regrid (superregss, regsss);
+ vhh.AT(m)->regrid (superregss, regsss, do_init);
+
+ // Output grid structure to screen
+ OutputSuperregions (cctkGH, m, * vhh.AT(m), superregss);
+ OutputGrids (cctkGH, m, * vhh.AT(m), * vdd.AT(m));
// Write grid structure to file
OutputGridStructure (cctkGH, m, regsss);
OutputGridCoordinates (cctkGH, m, regsss);
-#warning "TODO: output superregss"
-
- if (verbose or veryverbose) {
- OutputGrids (cctkGH, m, * vhh.at(m), * vdd.at(m));
- }
Waypoint ("Done regridding map %d.", m);
+ timer.stop();
}
@@ -279,13 +320,27 @@ namespace Carpet {
{
DECLARE_CCTK_PARAMETERS;
+ char const * const where = "Carpet::PostRegrid";
+ static Timer timer (where);
+ timer.start();
+
// Calculate new number of levels
int const oldreflevels = reflevels;
- reflevels = vhh.at(0)->reflevels();
+ reflevels = vhh.AT(0)->reflevels();
for (int m=0; m<maps; ++m) {
- assert (vhh.at(m)->reflevels() == reflevels);
+ assert (vhh.AT(m)->reflevels() == reflevels);
}
+#warning "TODO"
+#if 1 // #ifdef CARPET_DEBUG
+ {
+ // All processes must use the same number of levels
+ int ival = reflevels;
+ MPI_Bcast (& ival, 1, MPI_INT, 0, dist::comm());
+ assert (ival == reflevels);
+ }
+#endif
+
// One cannot switch off the current level
assert (reflevels > reflevel);
@@ -294,19 +349,23 @@ namespace Carpet {
int const grouptype = CCTK_GroupTypeI (n);
if (grouptype == CCTK_GF) {
for (int ml=0; ml<mglevels; ++ml) {
- groupdata.at(n).activetimelevels.at(ml).resize
- (reflevels, groupdata.at(n).activetimelevels.at(ml).at(0));
+ groupdata.AT(n).activetimelevels.AT(ml).resize
+ (reflevels, groupdata.AT(n).activetimelevels.AT(ml).AT(0));
}
}
}
// Set new level times
for (int ml=0; ml<mglevels; ++ml) {
- leveltimes.at(ml).resize
- (reflevels, leveltimes.at(ml).at(oldreflevels-1));
+ leveltimes.AT(ml).resize
+ (reflevels, leveltimes.AT(ml).AT(oldreflevels-1));
}
+ ++ regridding_epoch;
+
OutputGridStatistics (cctkGH);
+
+ timer.stop();
}
@@ -316,200 +375,291 @@ namespace Carpet {
int const rl,
bool const do_init)
{
+ char const * const where = "Carpet::Recompose";
+ static Timer timer (where);
+ timer.start();
+
bool did_recompose = false;
for (int m=0; m<maps; ++m) {
Waypoint ("Recomposing the grid hierarchy for map %d level %d...", m, rl);
- assert (rl>=0 and rl<vhh.at(m)->reflevels());
- did_recompose |= vhh.at(m)->recompose (rl, do_init);
+ assert (rl>=0 and rl<vhh.AT(m)->reflevels());
+ did_recompose |= vhh.AT(m)->recompose (rl, do_init);
Waypoint ("Done recomposing the grid hierarchy for map %d level %d.",
m, rl);
}
+
+ ClassifyPoints (cctkGH, rl);
+
+ timer.stop();
return did_recompose;
}
void
- OutputGrids (cGH const * const cctkGH,
- int const m,
- gh const & hh,
- dh const & dd)
+ RegridFree (cGH const * const cctkGH,
+ bool const do_init)
{
- CCTK_INFO ("Grid structure (grid points):");
- for (int ml=0; ml<hh.mglevels(); ++ml) {
- for (int rl=0; rl<hh.reflevels(); ++rl) {
- for (int c=0; c<hh.components(rl); ++c) {
- const int convfact = ipow(mgfact, ml);
- const ivect levfact = spacereffacts.at(rl);
- const ibbox ext = hh.extent(ml,rl,c);
- const ivect & lower = ext.lower();
- const ivect & upper = ext.upper();
- const ivect & stride = ext.stride();
- assert (all(lower * levfact % maxspacereflevelfact == 0));
- assert (all(upper * levfact % maxspacereflevelfact == 0));
- assert (all(((upper - lower) * levfact / maxspacereflevelfact)
- % convfact == 0));
- cout << " [" << ml << "][" << rl << "][" << m << "][" << c << "]"
- << " exterior: "
- << "proc "
- << hh.processor(rl,c)
- << " "
- << lower / stride
- << " : "
- << upper / stride
- << " ("
- << (upper - lower) / stride + 1
- << ") "
- << prod ((upper - lower) / stride + 1)
- << endl;
- }
- }
+ char const * const where = "Carpet::RegridFree";
+ static Timer timer (where);
+ timer.start();
+
+ Checkpoint ("Freeing after regridding level %d...", reflevel);
+
+ assert (is_level_mode());
+
+ // Free old grid structure
+ for (int m=0; m<maps; ++m) {
+ vhh.AT(m)->regrid_free (do_init);
}
- CCTK_INFO ("Grid structure (boundaries):");
- for (int rl=0; rl<hh.reflevels(); ++rl) {
- for (int c=0; c<hh.components(rl); ++c) {
+ timer.stop();
+ }
+
+
+
+ void
+ OutputSuperregions (cGH const * const cctkGH,
+ int const m,
+ gh const & hh,
+ gh::rregs const & superregss)
+ {
+ CCTK_INFO ("Grid structure (superregions, grid points):");
+ for (int rl=0; rl<(int)superregss.size(); ++rl) {
+ assert (rl < hh.reflevels());
+ for (int c=0; c<(int)superregss.AT(rl).size(); ++c) {
+ const ivect & levfact = spacereffacts.AT(rl);
+ const ibbox & ext = superregss.AT(rl).AT(c).extent;
+ const ivect & lower = ext.lower();
+ const ivect & upper = ext.upper();
+ const ivect & stride = ext.stride();
+ assert (all(lower * levfact % maxspacereflevelfact == 0));
+ assert (all(upper * levfact % maxspacereflevelfact == 0));
cout << " [" << rl << "][" << m << "][" << c << "]"
- << " bbox: "
- << hh.outer_boundaries(rl,c)
- << endl;
+ << " exterior: "
+ << lower / stride
+ << " : "
+ << upper / stride
+ << " ("
+ << (upper - lower) / stride + 1
+ << ") "
+ << prod ((upper - lower) / stride + 1)
+ << eol;
}
}
- CCTK_INFO ("Grid structure (coordinates):");
- for (int ml=0; ml<hh.mglevels(); ++ml) {
- for (int rl=0; rl<hh.reflevels(); ++rl) {
- for (int c=0; c<hh.components(rl); ++c) {
- const rvect origin = domainspecs.at(m).exterior_min;
- const rvect delta = (domainspecs.at(m).exterior_max - domainspecs.at(m).exterior_min) / rvect (domainspecs.at(m).npoints - 1);
- const ibbox ext = hh.extent(ml,rl,c);
- const ivect & lower = ext.lower();
- const ivect & upper = ext.upper();
- const int convfact = ipow(mgfact, ml);
- const ivect levfact = spacereffacts.at(rl);
- cout << " [" << ml << "][" << rl << "][" << m << "][" << c << "]"
- << " exterior: "
- << origin + delta * rvect(lower) / rvect(maxspacereflevelfact)
- << " : "
- << origin + delta * rvect(upper) / rvect(maxspacereflevelfact)
- << " : "
- << delta * rvect(convfact) / rvect(levfact) << endl;
- }
+ CCTK_INFO ("Grid structure (superregions, coordinates):");
+ for (int rl=0; rl<(int)superregss.size(); ++rl) {
+ assert (rl < hh.reflevels());
+ for (int c=0; c<(int)superregss.AT(rl).size(); ++c) {
+ const rvect origin = domainspecs.AT(m).exterior_min;
+ const rvect delta = (domainspecs.AT(m).exterior_max - domainspecs.AT(m).exterior_min) / rvect (domainspecs.AT(m).npoints - 1);
+ const ibbox & ext = superregss.AT(rl).AT(c).extent;
+ const ivect & lower = ext.lower();
+ const ivect & upper = ext.upper();
+ const ivect & levfact = spacereffacts.AT(rl);
+ cout << " [" << rl << "][" << m << "][" << c << "]"
+ << " exterior: "
+ << origin + delta * rvect(lower) / rvect(maxspacereflevelfact)
+ << " : "
+ << origin + delta * rvect(upper) / rvect(maxspacereflevelfact)
+ << " : "
+ << delta / rvect(levfact) << eol;
}
}
- CCTK_INFO ("Grid structure (coordinates, including ghosts):");
- for (int ml=0; ml<hh.mglevels(); ++ml) {
- for (int rl=0; rl<hh.reflevels(); ++rl) {
- for (int c=0; c<hh.components(rl); ++c) {
- const rvect origin = domainspecs.at(m).exterior_min;
- const rvect delta = (domainspecs.at(m).exterior_max - domainspecs.at(m).exterior_min) / rvect (domainspecs.at(m).npoints - 1);
- const ivect lower = dd.boxes.at(ml).at(rl).at(c).exterior.lower();
- const ivect upper = dd.boxes.at(ml).at(rl).at(c).exterior.upper();
- const int convfact = ipow(mgfact, ml);
- const ivect levfact = spacereffacts.at(rl);
- cout << " [" << ml << "][" << rl << "][" << m << "][" << c << "]"
- << " exterior: "
- << origin + delta * rvect(lower) / rvect(maxspacereflevelfact)
- << " : "
- << origin + delta * rvect(upper) / rvect(maxspacereflevelfact)
- << " : "
- << delta * rvect(convfact) / rvect(levfact) << endl;
+ fflush (stdout);
+ }
+
+
+
+ void
+ OutputGrids (cGH const * const cctkGH,
+ int const m,
+ gh const & hh,
+ dh const & dd)
+ {
+ DECLARE_CCTK_PARAMETERS;
+
+ if (verbose or veryverbose) {
+
+ CCTK_INFO ("Grid structure (grid points):");
+ for (int ml=0; ml<hh.mglevels(); ++ml) {
+ for (int rl=0; rl<hh.reflevels(); ++rl) {
+ for (int c=0; c<hh.components(rl); ++c) {
+ const int convfact = ipow(mgfact, ml);
+ const ivect levfact = spacereffacts.AT(rl);
+ const ibbox ext = hh.extent(ml,rl,c);
+ const ivect & lower = ext.lower();
+ const ivect & upper = ext.upper();
+ const ivect & stride = ext.stride();
+ assert (all(lower * levfact % maxspacereflevelfact == 0));
+ assert (all(upper * levfact % maxspacereflevelfact == 0));
+ assert (all(((upper - lower) * levfact / maxspacereflevelfact)
+ % convfact == 0));
+ cout << " [" << ml << "][" << rl << "][" << m << "][" << c << "]"
+ << " exterior: "
+ << "proc "
+ << hh.processor(rl,c)
+ << " "
+ << lower / stride
+ << " : "
+ << upper / stride
+ << " ("
+ << (upper - lower) / stride + 1
+ << ") "
+ << prod ((upper - lower) / stride + 1)
+ << eol;
+ }
}
}
- }
-
- CCTK_INFO ("Grid statistics:");
- const int oldprecision = cout.precision();
- const ios_base::fmtflags oldflags = cout.flags();
- cout.setf (ios::fixed);
- for (int ml=0; ml<hh.mglevels(); ++ml) {
- CCTK_REAL coarsevolume = 0;
+
+ CCTK_INFO ("Grid structure (boundaries):");
for (int rl=0; rl<hh.reflevels(); ++rl) {
-
- const CCTK_REAL basevolume = hh.baseextents.AT(0).AT(0).size();
- CCTK_REAL countvolume = 0;
- CCTK_REAL totalvolume = 0;
- CCTK_REAL totalvolume2 = 0;
-
for (int c=0; c<hh.components(rl); ++c) {
- const CCTK_REAL volume = hh.extent(ml,rl,c).size();
- ++ countvolume;
- totalvolume += volume;
- totalvolume2 += ipow(volume, 2);
+ cout << " [" << rl << "][" << m << "][" << c << "]"
+ << " bbox: "
+ << hh.outer_boundaries(rl,c)
+ << eol;
}
-
- const CCTK_REAL avgvolume = totalvolume / countvolume;
- const CCTK_REAL stddevvolume
- = sqrt (max (CCTK_REAL(0),
- totalvolume2 / countvolume - ipow (avgvolume, 2)));
-
- for (int c=0; c<hh.components(rl); ++c) {
- const CCTK_REAL volume = hh.extent(ml,rl,c).size();
- cout << " [" << ml << "][" << rl << "][" << m << "][" << c << "]"
- << " volume: " << setprecision(0) << volume
- << " of parent: " << setprecision(1) << 100 * volume / totalvolume << "%"
- << " of domain: " << setprecision(1) << 100 * volume / basevolume << "%"
- << endl;
+ }
+
+ CCTK_INFO ("Grid structure (coordinates):");
+ for (int ml=0; ml<hh.mglevels(); ++ml) {
+ for (int rl=0; rl<hh.reflevels(); ++rl) {
+ for (int c=0; c<hh.components(rl); ++c) {
+ const rvect origin = domainspecs.AT(m).exterior_min;
+ const rvect delta = (domainspecs.AT(m).exterior_max - domainspecs.AT(m).exterior_min) / rvect (domainspecs.AT(m).npoints - 1);
+ const ibbox ext = hh.extent(ml,rl,c);
+ const ivect & lower = ext.lower();
+ const ivect & upper = ext.upper();
+ const int convfact = ipow(mgfact, ml);
+ const ivect levfact = spacereffacts.AT(rl);
+ cout << " [" << ml << "][" << rl << "][" << m << "][" << c << "]"
+ << " exterior: "
+ << origin + delta * rvect(lower) / rvect(maxspacereflevelfact)
+ << " : "
+ << origin + delta * rvect(upper) / rvect(maxspacereflevelfact)
+ << " : "
+ << delta * rvect(convfact) / rvect(levfact) << eol;
+ }
}
-
- cout << " [" << ml << "][" << rl << "][" << m << "]"
- << " average volume: " << setprecision(0) << avgvolume
- << " of parent: " << setprecision(1) << 100 * avgvolume / totalvolume << "%"
- << " of domain: " << setprecision(1) << 100 * avgvolume / basevolume << "%"
- << endl;
- cout << " [" << ml << "][" << rl << "][" << m << "]"
- << " standard deviation: " << setprecision(0) << stddevvolume
- << " of parent: " << setprecision(1) << 100 * stddevvolume / totalvolume << "%"
- << " of domain: " << setprecision(1) << 100 * stddevvolume / basevolume << "%"
- << endl;
-
- // TODO: ghost points vs. interior points (and boundary
- // points)
-
- CCTK_REAL countquadrupole = 0;
- CCTK_REAL minquadrupole = 1;
- CCTK_REAL totalquadrupole = 0;
- CCTK_REAL totalquadrupole2 = 0;
-
- for (int c=0; c<hh.components(rl); ++c) {
- const ibbox ext = hh.extent(ml,rl,c);
- const ivect shape = ext.shape();
- const ivect stride = ext.stride();
- const CCTK_REAL minlength = minval (rvect (shape) / rvect (stride));
- const CCTK_REAL maxlength = maxval (rvect (shape) / rvect (stride));
- const CCTK_REAL quadrupole = minlength / maxlength;
- ++ countquadrupole;
- minquadrupole = min (minquadrupole, quadrupole);
- totalquadrupole += quadrupole;
- totalquadrupole2 += ipow (quadrupole, 2);
- cout << " [" << ml << "][" << rl << "][" << m << "][" << c << "]"
- << " length ratio: " << setprecision(3) << quadrupole
- << endl;
+ }
+
+ CCTK_INFO ("Grid structure (coordinates, including ghosts):");
+ for (int ml=0; ml<hh.mglevels(); ++ml) {
+ for (int rl=0; rl<hh.reflevels(); ++rl) {
+ for (int c=0; c<hh.components(rl); ++c) {
+ const rvect origin = domainspecs.AT(m).exterior_min;
+ const rvect delta = (domainspecs.AT(m).exterior_max - domainspecs.AT(m).exterior_min) / rvect (domainspecs.AT(m).npoints - 1);
+ const ivect lower = dd.boxes.AT(ml).AT(rl).AT(c).exterior.lower();
+ const ivect upper = dd.boxes.AT(ml).AT(rl).AT(c).exterior.upper();
+ const int convfact = ipow(mgfact, ml);
+ const ivect levfact = spacereffacts.AT(rl);
+ cout << " [" << ml << "][" << rl << "][" << m << "][" << c << "]"
+ << " exterior: "
+ << origin + delta * rvect(lower) / rvect(maxspacereflevelfact)
+ << " : "
+ << origin + delta * rvect(upper) / rvect(maxspacereflevelfact)
+ << " : "
+ << delta * rvect(convfact) / rvect(levfact) << eol;
+ }
+ }
+ }
+
+ CCTK_INFO ("Grid statistics:");
+ const int oldprecision = cout.precision();
+ const ios_base::fmtflags oldflags = cout.flags();
+ cout.setf (ios::fixed);
+ for (int ml=0; ml<hh.mglevels(); ++ml) {
+ CCTK_REAL coarsevolume = 0;
+ for (int rl=0; rl<hh.reflevels(); ++rl) {
+
+ const CCTK_REAL basevolume = hh.baseextents.AT(0).AT(0).size();
+ CCTK_REAL countvolume = 0;
+ CCTK_REAL totalvolume = 0;
+ CCTK_REAL totalvolume2 = 0;
+
+ for (int c=0; c<hh.components(rl); ++c) {
+ const CCTK_REAL volume = hh.extent(ml,rl,c).size();
+ ++ countvolume;
+ totalvolume += volume;
+ totalvolume2 += ipow(volume, 2);
+ }
+
+ const CCTK_REAL avgvolume = totalvolume / countvolume;
+ const CCTK_REAL stddevvolume
+ = sqrt (max (CCTK_REAL(0),
+ totalvolume2 / countvolume - ipow (avgvolume, 2)));
+
+ for (int c=0; c<hh.components(rl); ++c) {
+ const CCTK_REAL volume = hh.extent(ml,rl,c).size();
+ cout << " [" << ml << "][" << rl << "][" << m << "][" << c << "]"
+ << " volume: " << setprecision(0) << volume
+ << " of parent: " << setprecision(1) << 100 * volume / totalvolume << "%"
+ << " of domain: " << setprecision(1) << 100 * volume / basevolume << "%"
+ << eol;
+ }
+
+ cout << " [" << ml << "][" << rl << "][" << m << "]"
+ << " average volume: " << setprecision(0) << avgvolume
+ << " of parent: " << setprecision(1) << 100 * avgvolume / totalvolume << "%"
+ << " of domain: " << setprecision(1) << 100 * avgvolume / basevolume << "%"
+ << eol;
+ cout << " [" << ml << "][" << rl << "][" << m << "]"
+ << " standard deviation: " << setprecision(0) << stddevvolume
+ << " of parent: " << setprecision(1) << 100 * stddevvolume / totalvolume << "%"
+ << " of domain: " << setprecision(1) << 100 * stddevvolume / basevolume << "%"
+ << eol;
+
+ // TODO: ghost points vs. interior points (and boundary
+ // points)
+
+ CCTK_REAL countquadrupole = 0;
+ CCTK_REAL minquadrupole = 1;
+ CCTK_REAL totalquadrupole = 0;
+ CCTK_REAL totalquadrupole2 = 0;
+
+ for (int c=0; c<hh.components(rl); ++c) {
+ const ibbox ext = hh.extent(ml,rl,c);
+ const ivect shape = ext.shape();
+ const ivect stride = ext.stride();
+ const CCTK_REAL minlength = minval (rvect (shape) / rvect (stride));
+ const CCTK_REAL maxlength = maxval (rvect (shape) / rvect (stride));
+ const CCTK_REAL quadrupole = minlength / maxlength;
+ ++ countquadrupole;
+ minquadrupole = min (minquadrupole, quadrupole);
+ totalquadrupole += quadrupole;
+ totalquadrupole2 += ipow (quadrupole, 2);
+ cout << " [" << ml << "][" << rl << "][" << m << "][" << c << "]"
+ << " length ratio: " << setprecision(3) << quadrupole
+ << eol;
+ }
+
+ const CCTK_REAL avgquadrupole = totalquadrupole / countquadrupole;
+ const CCTK_REAL stddevquadrupole
+ = sqrt (max (CCTK_REAL(0),
+ (totalquadrupole2 / countquadrupole
+ - ipow (avgquadrupole, 2))));
+
+ cout << " [" << ml << "][" << rl << "][" << m << "]"
+ << " average length ratio: " << setprecision(3) << avgquadrupole
+ << " standard deviation: " << setprecision(3) << stddevquadrupole
+ << eol;
+
+ // TODO: processor distribution, average load, std deviation
+
+ coarsevolume = totalvolume * prod (rvect (spacereflevelfact));
}
-
- const CCTK_REAL avgquadrupole = totalquadrupole / countquadrupole;
- const CCTK_REAL stddevquadrupole
- = sqrt (max (CCTK_REAL(0),
- (totalquadrupole2 / countquadrupole
- - ipow (avgquadrupole, 2))));
-
- cout << " [" << ml << "][" << rl << "][" << m << "]"
- << " average length ratio: " << setprecision(3) << avgquadrupole
- << " standard deviation: " << setprecision(3) << stddevquadrupole
- << endl;
-
- // TODO: processor distribution, average load, std deviation
-
- coarsevolume = totalvolume * prod (rvect (spacereflevelfact));
}
+ cout.precision (oldprecision);
+ cout.setf (oldflags);
+
+ fflush (stdout);
}
- cout.precision (oldprecision);
- cout.setf (oldflags);
-
}
@@ -561,14 +711,14 @@ namespace Carpet {
file << "maps " << maps << eol;
file << m << " mglevels " << regsss.size() << eol;
for (int ml=0; ml<(int)regsss.size(); ++ml) {
- file << m << " " << ml << " reflevels " << regsss.at(ml).size() << eol;
- for (int rl=0; rl<(int)regsss.at(ml).size(); ++rl) {
- file << m << " " << ml << " " << rl << " components " << regsss.at(ml).at(rl).size() << eol;
- for (int c=0; c<(int)regsss.at(ml).at(rl).size(); ++c) {
+ file << m << " " << ml << " reflevels " << regsss.AT(ml).size() << eol;
+ for (int rl=0; rl<(int)regsss.AT(ml).size(); ++rl) {
+ file << m << " " << ml << " " << rl << " components " << regsss.AT(ml).AT(rl).size() << eol;
+ for (int c=0; c<(int)regsss.AT(ml).AT(rl).size(); ++c) {
file << m << " " << ml << " " << rl << " " << c << " "
- << regsss.at(ml).at(rl).at(c).processor << " "
- << regsss.at(ml).at(rl).at(c).extent << " "
- << regsss.at(ml).at(rl).at(c).outer_boundaries << eol;
+ << regsss.AT(ml).AT(rl).AT(c).processor << " "
+ << regsss.AT(ml).AT(rl).AT(c).extent << " "
+ << regsss.AT(ml).AT(rl).AT(c).outer_boundaries << eol;
}
}
}
@@ -668,7 +818,7 @@ namespace Carpet {
// Affine transformation between index space and coordinate space
rvect const origin = exterior_lower;
rvect const scale =
- rvect (vhh.at(m)->baseextents.at(0).at(0).stride()) / spacing;
+ rvect (vhh.AT(m)->baseextents.AT(0).AT(0).stride()) / spacing;
@@ -676,11 +826,11 @@ namespace Carpet {
file << "maps " << maps << eol;
file << m << " mglevels " << regsss.size() << eol;
for (int ml=0; ml<(int)regsss.size(); ++ml) {
- file << m << " " << ml << " reflevels " << regsss.at(ml).size() << eol;
- for (int rl=0; rl<(int)regsss.at(ml).size(); ++rl) {
+ file << m << " " << ml << " reflevels " << regsss.AT(ml).size() << eol;
+ for (int rl=0; rl<(int)regsss.AT(ml).size(); ++rl) {
ibset extents;
- for (int c=0; c<(int)regsss.at(ml).at(rl).size(); ++c) {
- extents += regsss.at(ml).at(rl).at(c).extent;
+ for (int c=0; c<(int)regsss.AT(ml).AT(rl).size(); ++c) {
+ extents += regsss.AT(ml).AT(rl).AT(c).extent;
}
extents.normalize();
file << m << " " << ml << " " << rl << " regions " << extents.setsize() << eol;
@@ -690,12 +840,12 @@ namespace Carpet {
{
#if 0
ibbox const & ext = * bi;
- ibbox const & baseext = vhh.at(m)->baseextents.at(ml).at(rl);
- ibbox const & coarseext = vhh.at(m)->baseextents.at(ml).at(0 );
+ ibbox const & baseext = vhh.AT(m)->baseextents.AT(ml).AT(rl);
+ ibbox const & coarseext = vhh.AT(m)->baseextents.AT(ml).AT(0 );
// This is nice, but is wrong since CartGrid3D has not yet
// initialised the coordinates
- ivect const cctk_levfac = spacereffacts.at(rl);
+ ivect const cctk_levfac = spacereffacts.AT(rl);
ivect const cctk_levoff = baseext.lower() - coarseext.lower();
ivect const cctk_levoffdenom = baseext.stride();
@@ -705,9 +855,9 @@ namespace Carpet {
(ext.upper() - baseext.lower()) / ext.stride();
rvect const cctk_origin_space =
- origin_space.at(m).at(ml);
+ origin_space.AT(m).AT(ml);
rvect const cctk_delta_space =
- delta_space.at(m) * rvect (ipow (mgfact, ml));
+ delta_space.AT(m) * rvect (ipow (mgfact, ml));
rvect const CCTK_ORIGIN_SPACE =
cctk_origin_space +
@@ -773,28 +923,35 @@ namespace Carpet {
OutputGridStatistics (cGH const * const cctkGH)
{
// Grid array statistics
- int num_gfs = 0;
int num_arrays = 0;
+ int num_gfs = 0;
+ int size_gfs = 0;
CCTK_REAL num_active_array_points = 0;
CCTK_REAL num_total_array_points = 0;
+ CCTK_REAL size_total_array_points = 0;
for (int g=0; g<CCTK_NumGroups(); ++g) {
+ cGroup gdata;
+ check (not CCTK_GroupData (g, &gdata));
int const num_tl = CCTK_ActiveTimeLevelsGI (cctkGH, g);
- switch (CCTK_GroupTypeI(g)) {
+ int const num_vars = gdata.numvars;
+ int const size_vars = gdata.numvars * CCTK_VarTypeSize (gdata.vartype);
+ switch (gdata.grouptype) {
case CCTK_SCALAR:
case CCTK_ARRAY: {
- int const num_vars = CCTK_NumVarsInGroupI (g);
num_arrays += num_tl * num_vars;
gh const * const hh = arrdata.AT(g).AT(0).hh;
dh const * const dd = arrdata.AT(g).AT(0).dd;
for (int c=0; c<hh->components(0); ++c) {
dh::dboxes const & b = dd->boxes.AT(0).AT(0).AT(c);
- num_active_array_points += num_tl * num_vars * b.active.size();
- num_total_array_points += num_tl * num_vars * b.exterior.size();
+ num_active_array_points += num_tl * num_vars * b.active_size;
+ num_total_array_points += num_tl * num_vars * b.exterior_size;
+ size_total_array_points += num_tl * size_vars * b.exterior_size;
}
break;
}
case CCTK_GF:
- num_gfs += num_tl;
+ num_gfs += num_tl * num_vars;
+ size_gfs += num_tl * size_vars;
break;
default:
assert (0);
@@ -807,6 +964,7 @@ namespace Carpet {
CCTK_REAL num_active_mem_points = 0;
CCTK_REAL num_owned_mem_points = 0;
CCTK_REAL num_total_mem_points = 0;
+ CCTK_REAL size_total_mem_points = 0;
CCTK_REAL num_active_cpu_points = 0;
CCTK_REAL num_owned_cpu_points = 0;
CCTK_REAL num_total_cpu_points = 0;
@@ -822,12 +980,13 @@ namespace Carpet {
for (int c=0; c<hh->components(rl); ++c) {
++ num_comps;
dh::dboxes const & b = dd->boxes.AT(ml).AT(rl).AT(c);
- num_active_mem_points += num_gfs * b.active.size();
- num_owned_mem_points += num_gfs * b.owned.size();
- num_total_mem_points += num_gfs * b.exterior.size();
- num_active_cpu_points += trf * b.active.size();
- num_owned_cpu_points += trf * b.owned.size();
- num_total_cpu_points += trf * b.exterior.size();
+ num_active_mem_points += num_gfs * b.active_size;
+ num_owned_mem_points += num_gfs * b.owned_size;
+ num_total_mem_points += num_gfs * b.exterior_size;
+ size_total_mem_points += size_gfs * b.exterior_size;
+ num_active_cpu_points += trf * b.active_size;
+ num_owned_cpu_points += trf * b.owned_size;
+ num_total_cpu_points += trf * b.exterior_size;
}
}
}
@@ -841,27 +1000,30 @@ namespace Carpet {
"Grid structure statistics:");
CCTK_VInfo (CCTK_THORNSTRING,
"GF: rhs: %.0fk active, %.0fk owned (+%.0f%%), %.0fk total (+%.0f%%), %.3g steps/time",
- double (num_active_cpu_points / 1000),
- double (num_owned_cpu_points / 1000),
+ double (num_active_cpu_points / 1e+3),
+ double (num_owned_cpu_points / 1e+3),
double (num_owned_cpu_points / num_active_cpu_points * 100 - 100),
- double (num_total_cpu_points / 1000),
+ double (num_total_cpu_points / 1e+3),
double (num_total_cpu_points / num_owned_cpu_points * 100 - 100),
double (num_steps / delta_time));
CCTK_VInfo (CCTK_THORNSTRING,
"GF: vars: %d, pts: %.0fM active, %.0fM owned (+%.0f%%), %.0fM total (+%.0f%%), %.1f comp/proc",
num_gfs,
- double (num_active_mem_points / 1000000),
- double (num_owned_mem_points / 1000000),
+ double (num_active_mem_points / 1e+6),
+ double (num_owned_mem_points / 1e+6),
double (num_owned_mem_points / num_active_mem_points * 100 - 100),
- double (num_total_mem_points / 1000000),
+ double (num_total_mem_points / 1e+6),
double (num_total_mem_points / num_owned_mem_points * 100 - 100),
double (num_comps / (reflevels * CCTK_nProcs (cctkGH))));
CCTK_VInfo (CCTK_THORNSTRING,
"GA: vars: %d, pts: %.0fM active, %.0fM total (+%.0f%%)",
num_arrays,
- double (num_active_array_points / 1000000),
- double (num_total_array_points / 1000000),
+ double (num_active_array_points / 1e+6),
+ double (num_total_array_points / 1e+6),
double (num_total_array_points / num_active_array_points * 100 - 100));
+ CCTK_VInfo (CCTK_THORNSTRING,
+ "Total required memory: %.3f GByte (for GAs and currently active GFs)",
+ double ((size_total_array_points + size_total_mem_points) / 1e+9));
// After this, we will begin to allocate memory for the grid
// structure. If we run out of memory, ensure that this output
@@ -902,16 +1064,16 @@ namespace Carpet {
regs = superregs;
if (nprocs == 1) {
- regs.at(0).processor = 0;
- pseudoregion_t const preg (regs.at(0).extent, regs.at(0).processor);
- assert (superregs.at(0).processors == NULL);
- superregs.at(0).processors = new ipfulltree (preg);
+ regs.AT(0).processor = 0;
+ pseudoregion_t const preg (regs.AT(0).extent, regs.AT(0).processor);
+ assert (superregs.AT(0).processors == NULL);
+ superregs.AT(0).processors = new ipfulltree (preg);
return;
}
assert (dir>=0 and dir<dim);
- region_t const & reg0 = regs.at(0);
+ region_t const & reg0 = regs.AT(0);
const ivect rstr0 = reg0.extent.stride();
const ivect rlb0 = reg0.extent.lower();
const ivect rub0 = reg0.extent.upper() + rstr0;
@@ -933,7 +1095,7 @@ namespace Carpet {
if (cub[dir] > rub0[dir]) cub[dir] = rub0[dir];
assert (clb[dir] <= cub[dir]);
assert (cub[dir] <= rub0[dir]);
- region_t & reg = regs.at(c);
+ region_t & reg = regs.AT(c);
ibbox & ext = reg.extent;
b2vect & obnd = reg.outer_boundaries;
int & proc = reg.processor;
@@ -943,16 +1105,16 @@ namespace Carpet {
obnd[1][dir] &= cub[dir] == rub0[dir];
proc = c;
pseudoregion_t const preg (reg.extent, c);
- bounds.at(c) = reg.extent.lower()[dir];
- subtrees.at(c) = new ipfulltree (preg);
+ bounds.AT(c) = reg.extent.lower()[dir];
+ subtrees.AT(c) = new ipfulltree (preg);
}
- bounds.at(nprocs) = rub0[dir] + rstr0[dir];
+ bounds.AT(nprocs) = rub0[dir] + rstr0[dir];
- assert (superregs.at(0).processors == NULL);
- superregs.at(0).processors = new ipfulltree (dir, bounds, subtrees);
+ assert (superregs.AT(0).processors == NULL);
+ superregs.AT(0).processors = new ipfulltree (dir, bounds, subtrees);
for (int c=0; c<(int)regs.size(); ++c) {
- assert (regs.at(c).processor == c);
+ assert (regs.AT(c).processor == c);
}
}
@@ -968,9 +1130,9 @@ namespace Carpet {
vector<vector<region_t> > regss (1);
SplitRegionsMaps_Automatic (cctkGH, superregss, regss);
assert (superregss.size() == 1);
- superregs = regss.at(0);
+ superregs = superregss.AT(0);
assert (regss.size() == 1);
- regs = regss.at(0);
+ regs = regss.AT(0);
}
@@ -993,7 +1155,7 @@ namespace Carpet {
assert (superregs.size() == 1);
- region_t const & reg0 = superregs.at(0);
+ region_t const & reg0 = superregs.AT(0);
const ivect rstr0 = reg0.extent.stride();
const ivect rlb0 = reg0.extent.lower();
const ivect rub0 = reg0.extent.upper() + rstr0;
@@ -1052,7 +1214,7 @@ namespace Carpet {
assert (all (cub <= rub0));
assert (all (not (ipos==0) or clb==rlb0));
assert (all (not (ipos==nprocs_dir-1) or cub==rub0));
- region_t & reg = regs.at(c);
+ region_t & reg = regs.AT(c);
ibbox & ext = reg.extent;
b2vect & obnd = reg.outer_boundaries;
int & proc = reg.processor;
@@ -1063,21 +1225,21 @@ namespace Carpet {
proc = c;
pseudoregion_t preg (reg.extent, c);
- subtreesx.at(i) = new ipfulltree (preg);
+ subtreesx.AT(i) = new ipfulltree (preg);
}
- boundsx.at(nprocs_dir[0]) = rub0[0] + rstr0[0];
- subtreesy.at(j) = new ipfulltree (0, boundsx, subtreesx);
+ boundsx.AT(nprocs_dir[0]) = rub0[0] + rstr0[0];
+ subtreesy.AT(j) = new ipfulltree (0, boundsx, subtreesx);
}
- boundsy.at(nprocs_dir[1]) = rub0[1] + rstr0[1];
- subtreesz.at(k) = new ipfulltree (1, boundsy, subtreesy);
+ boundsy.AT(nprocs_dir[1]) = rub0[1] + rstr0[1];
+ subtreesz.AT(k) = new ipfulltree (1, boundsy, subtreesy);
}
- boundsz.at(nprocs_dir[2]) = rub0[2] + rstr0[2];
+ boundsz.AT(nprocs_dir[2]) = rub0[2] + rstr0[2];
- assert (superregs.at(0).processors == NULL);
- superregs.at(0).processors = new ipfulltree (2, boundsz, subtreesz);
+ assert (superregs.AT(0).processors == NULL);
+ superregs.AT(0).processors = new ipfulltree (2, boundsz, subtreesz);
for (int c=0; c<(int)regs.size(); ++c) {
- assert (regs.at(c).processor == c);
+ assert (regs.AT(c).processor == c);
}
}
@@ -1121,7 +1283,9 @@ namespace Carpet {
region_t & superreg,
vector<region_t> & newregs)
{
- if (DEBUG) cout << "SRMAR enter" << endl;
+ DECLARE_CCTK_PARAMETERS;
+
+ if (recompose_verbose) cout << "SRMAR enter" << endl;
// Check preconditions
assert (nprocs >= 1);
@@ -1130,7 +1294,7 @@ namespace Carpet {
// Are we done?
if (all(dims)) {
- if (DEBUG) cout << "SRMAR bottom" << endl;
+ if (recompose_verbose) cout << "SRMAR bottom" << endl;
// Check precondition
assert (nprocs == 1);
@@ -1145,30 +1309,57 @@ namespace Carpet {
// Check postcondition
assert (newregs.size() == oldsize + nprocs);
- if (DEBUG) cout << "SRMAR exit" << endl;
+ if (recompose_verbose) cout << "SRMAR exit" << endl;
return;
}
// Special case
if (superreg.extent.empty()) {
- if (DEBUG) cout << "SRMAR empty" << endl;
+ if (recompose_verbose) cout << "SRMAR empty" << endl;
+
+ // Choose a direction
+ int mydim = -1;
+ for (int d=0; d<dim; ++d) {
+ if (not dims[d]) {
+ mydim = d;
+ break;
+ }
+ }
+ assert (mydim>=0 and mydim<dim);
// Create a new region
region_t newreg (superreg);
newreg.outer_boundaries = b2vect(false);
- if (DEBUG) cout << "SRMAR newreg " << newreg << endl;
+ if (recompose_verbose) cout << "SRMAR newreg " << newreg << endl;
+#if 0
// Store
for (int p=0; p<nprocs; ++p) {
newreg.processor = firstproc + p;
newregs.push_back (newreg);
}
+ assert (superreg.processors == NULL);
superreg.processors = new ipfulltree ();
+#endif
+
+ // Store
+ vector<int> bounds (nprocs+1);
+ vector<ipfulltree *> subtrees (nprocs);
+ for (int p=0; p<nprocs; ++p) {
+ newreg.processor = firstproc + p;
+ newregs.push_back (newreg);
+ bounds.AT(p) = newreg.extent.lower()[mydim];
+ pseudoregion_t const preg (newreg.extent, newreg.processor);
+ subtrees.AT(p) = new ipfulltree (preg);
+ }
+ bounds.AT(nprocs) = newreg.extent.lower()[mydim];
+ assert (superreg.processors == NULL);
+ superreg.processors = new ipfulltree (mydim, bounds, subtrees);
// Check postcondition
assert (newregs.size() == oldsize + nprocs);
- if (DEBUG) cout << "SRMAR exit" << endl;
+ if (recompose_verbose) cout << "SRMAR exit" << endl;
return;
}
@@ -1177,7 +1368,7 @@ namespace Carpet {
CCTK_REAL const rfact = pow (nprocs / prod(rcost), CCTK_REAL(1)/dim);
rcost *= rfact;
assert (abs (prod (rcost) - nprocs) <= 1.0e-6);
- if (DEBUG) cout << "SRMA shapes " << rcost << endl;
+ if (recompose_verbose) cout << "SRMA shapes " << rcost << endl;
// Choose a direction
int mydim = -1;
@@ -1198,22 +1389,20 @@ namespace Carpet {
}
assert (mydim>=0 and mydim<dim);
assert (mycost>=0);
- if (DEBUG) cout << "SRMAR mydim " << mydim << endl;
- if (DEBUG) cout << "SRMAR mycost " << mycost << endl;
+ if (recompose_verbose) cout << "SRMAR mydim " << mydim << endl;
+ if (recompose_verbose) cout << "SRMAR mycost " << mycost << endl;
// Mark this direction as done
assert (not dims[mydim]);
bvect const newdims = dims.replace(mydim, true);
// Choose a number of slices for this direction
- int const nslices
- = min (nprocs,
- int (floor (mycost * pow(nprocs / totalcost,
- CCTK_REAL(1) / alldims) +
- CCTK_REAL(0.5))));
+ CCTK_REAL const mycost1 =
+ mycost * pow(nprocs / totalcost, CCTK_REAL(1) / alldims);
+ int const nslices = min (nprocs, int (floor (mycost1 + CCTK_REAL(0.5))));
assert (nslices <= nprocs);
- if (DEBUG) cout << "SRMAR " << mydim << " nprocs " << nprocs << endl;
- if (DEBUG) cout << "SRMAR " << mydim << " nslices " << nslices << endl;
+ if (recompose_verbose) cout << "SRMAR " << mydim << " nprocs " << nprocs << endl;
+ if (recompose_verbose) cout << "SRMAR " << mydim << " nslices " << nslices << endl;
// Split the remaining processors
vector<int> mynprocs(nslices);
@@ -1221,11 +1410,11 @@ namespace Carpet {
int const mynprocs_left = nprocs - nslices * mynprocs_base;
int sum_mynprocs = 0;
for (int n=0; n<nslices; ++n) {
- mynprocs.at(n) = mynprocs_base + int (n < mynprocs_left);
- sum_mynprocs += mynprocs.at(n);
+ mynprocs.AT(n) = mynprocs_base + int (n < mynprocs_left);
+ sum_mynprocs += mynprocs.AT(n);
}
assert (sum_mynprocs == nprocs);
- if (DEBUG) cout << "SRMAR " << mydim << " mynprocs " << mynprocs << endl;
+ if (recompose_verbose) cout << "SRMAR " << mydim << " mynprocs " << mynprocs << endl;
// Split the region
vector<int> mynpoints(nslices);
@@ -1237,20 +1426,22 @@ namespace Carpet {
int npoints_left = npoints;
int nprocs_left = nprocs;
for (int n=0; n<nslices; ++n) {
- mynpoints.at(n) = int (floor (CCTK_REAL(1) * npoints_left * mynprocs.at(n)
- / nprocs_left + CCTK_REAL(0.5)));
- assert (mynprocs .at(n) > 0);
- npoints_left -= mynpoints.at(n);
- nprocs_left -= mynprocs.at(n);
+ assert (nprocs_left > 0);
+ CCTK_REAL const npoints1 =
+ CCTK_REAL(1) * npoints_left * mynprocs.AT(n) / nprocs_left;
+ mynpoints.AT(n) = int (floor (npoints1 + CCTK_REAL(0.5)));
+ assert (mynpoints.AT(n) >= 0 and mynpoints.AT(n) <= npoints_left);
+ npoints_left -= mynpoints.AT(n);
+ nprocs_left -= mynprocs.AT(n);
assert (npoints_left >= 0);
assert (nprocs_left >= 0);
}
assert (npoints_left == 0);
assert (nprocs_left == 0);
- if (DEBUG) cout << "SRMAR " << mydim << " mynpoints " << mynpoints << endl;
+ if (recompose_verbose) cout << "SRMAR " << mydim << " mynpoints " << mynpoints << endl;
// Create the regions and recurse
- if (DEBUG) cout << "SRMAR " << mydim << ": create bboxes" << endl;
+ if (recompose_verbose) cout << "SRMAR " << mydim << ": create bboxes" << endl;
ivect lo = superreg.extent.lower();
ivect up = superreg.extent.upper();
ivect const str = superreg.extent.stride();
@@ -1258,41 +1449,41 @@ namespace Carpet {
vector<int> bounds (nslices+1);
vector<ipfulltree *> subtrees (nslices);
for (int n=0; n<nslices; ++n) {
- if (DEBUG) cout << "SRMAR " << mydim << " n " << n << endl;
+ if (recompose_verbose) cout << "SRMAR " << mydim << " n " << n << endl;
// Create a new region
- up[mydim] = lo[mydim] + (mynpoints.at(n) - 1) * str[mydim];
+ up[mydim] = lo[mydim] + (mynpoints.AT(n) - 1) * str[mydim];
b2vect newob (superreg.outer_boundaries);
if (n > 0) {
newob[0][mydim] = false;
}
if (n < nslices-1) {
- up[mydim] = lo[mydim] + (mynpoints.at(n) - 1) * str[mydim];
+ up[mydim] = lo[mydim] + (mynpoints.AT(n) - 1) * str[mydim];
newob[1][mydim] = false;
}
region_t newreg (superreg);
newreg.extent = ibbox(lo, up, str);
newreg.outer_boundaries = newob;
- if (DEBUG) cout << "SRMAR " << mydim << " newreg " << newreg << endl;
+ if (recompose_verbose) cout << "SRMAR " << mydim << " newreg " << newreg << endl;
// Recurse
- bounds.at(n) = lo[mydim];
+ bounds.AT(n) = lo[mydim];
SplitRegionsMaps_Automatic_Recursively
- (newdims, p, mynprocs.at(n), newreg, newregs);
- if (DEBUG) cout << "SRMAR " << mydim << " newregs " << newregs << endl;
+ (newdims, p, mynprocs.AT(n), newreg, newregs);
+ if (recompose_verbose) cout << "SRMAR " << mydim << " newregs " << newregs << endl;
assert (newreg.processors != NULL);
- subtrees.at(n) = newreg.processors;
+ subtrees.AT(n) = newreg.processors;
newreg.processors = NULL;
newreg.processor = p;
// Next slice
lo[mydim] = up[mydim] + str[mydim];
- p += mynprocs.at(n);
+ p += mynprocs.AT(n);
}
assert (up[mydim] == superreg.extent.upper()[mydim]);
assert (p == firstproc + nprocs);
- bounds.at(nslices) = up[mydim] + str[mydim];
+ bounds.AT(nslices) = up[mydim] + str[mydim];
// Create tree
assert (superreg.processors == NULL);
@@ -1301,7 +1492,7 @@ namespace Carpet {
// Check postcondition
assert (newregs.size() == oldsize + nprocs);
- if (DEBUG) cout << "SRMAR exit" << endl;
+ if (recompose_verbose) cout << "SRMAR exit" << endl;
}
@@ -1313,20 +1504,20 @@ namespace Carpet {
{
DECLARE_CCTK_PARAMETERS;
- if (DEBUG) cout << "SRMA enter" << endl;
+ if (recompose_verbose) cout << "SRMA enter" << endl;
int const nmaps = superregss.size();
int minmap = 1000000000;
for (int m=0; m<nmaps; ++m) {
- for (int r=0; r<int(superregss.at(m).size()); ++r) {
- minmap = min (minmap, superregss.at(m).at(r).map);
+ for (int r=0; r<int(superregss.AT(m).size()); ++r) {
+ minmap = min (minmap, superregss.AT(m).AT(r).map);
}
}
int nregs = 0;
for (int m=0; m<nmaps; ++m) {
- nregs += superregss.at(m).size();
+ nregs += superregss.AT(m).size();
}
- if (DEBUG) cout << "SRMA nregs " << nregs << endl;
+ if (recompose_verbose) cout << "SRMA nregs " << nregs << endl;
// Something to do?
if (nregs == 0) {
@@ -1337,7 +1528,7 @@ namespace Carpet {
vector<region_t> superregs;
{
for (int m=0; m<nmaps; ++m) {
- combine_regions (superregss.at(m), superregs);
+ combine_regions (superregss.AT(m), superregs);
}
nregs = superregs.size();
@@ -1351,14 +1542,13 @@ namespace Carpet {
reg.extent = ibbox (ivect (0), ivect (-1), ivect (1));
reg.outer_boundaries = b2vect (bvect (true), bvect (true));
reg.map = 0;
- reg.processor = -1;
superregs.push_back (reg);
nregs = superregs.size();
}
}
int const real_nprocs = CCTK_nProcs (cctkGH);
- if (DEBUG) cout << "SRMA real_nprocs " << real_nprocs << endl;
+ if (recompose_verbose) cout << "SRMA real_nprocs " << real_nprocs << endl;
// Deactivate some processors if there are too many
int nprocs;
@@ -1367,109 +1557,155 @@ namespace Carpet {
} else {
CCTK_REAL mycost = 0;
for (int r=0; r<nregs; ++r) {
- mycost += prod (cost (superregs.at(r)));
+ mycost += prod (cost (superregs.AT(r)));
}
int const goodnprocs = int (floor (mycost / min_points_per_proc));
nprocs = max (1, min (real_nprocs, goodnprocs));
}
- if (DEBUG) cout << "SRMA nprocs " << nprocs << endl;
+ if (recompose_verbose) cout << "SRMA nprocs " << nprocs << endl;
// ncomps: number of components per processor
int const ncomps = (nregs + nprocs - 1) / nprocs;
- if (DEBUG) cout << "SRMA ncomps " << ncomps << endl;
+ if (recompose_verbose) cout << "SRMA ncomps " << ncomps << endl;
assert (ncomps > 0);
int const newnregs = nprocs * ncomps;
- if (DEBUG) cout << "SRMA newnregs " << newnregs << endl;
+ if (recompose_verbose) cout << "SRMA newnregs " << newnregs << endl;
- if (DEBUG) cout << "SRMA: distributing processors to regions" << endl;
+ if (recompose_verbose) cout << "SRMA: distributing processors to regions" << endl;
vector<CCTK_REAL> mycosts(nregs);
for (int r=0; r<nregs; ++r) {
- mycosts.at(r) = prod (cost (superregs.at(r)));
+ mycosts.AT(r) = prod (cost (superregs.AT(r)));
}
int nregs_left = newnregs;
vector<int> mynprocs(nregs);
for (int r=0; r<nregs; ++r) {
- mynprocs.at(r) = 1;
+ mynprocs.AT(r) = 1;
-- nregs_left;
}
#warning "TODO: split regions if necessary"
while (nregs_left > 0) {
- if (DEBUG) cout << "SRMA nregs_left " << nregs_left << endl;
+ if (recompose_verbose) cout << "SRMA nregs_left " << nregs_left << endl;
int maxr = -1;
CCTK_REAL maxratio = -1;
for (int r=0; r<nregs; ++r) {
- CCTK_REAL const ratio = mycosts.at(r) / mynprocs.at(r);
+ CCTK_REAL const ratio = mycosts.AT(r) / mynprocs.AT(r);
if (ratio > maxratio) { maxr=r; maxratio=ratio; }
}
assert (maxr>=0 and maxr<nregs);
- ++ mynprocs.at(maxr);
- if (DEBUG) cout << "SRMA maxr " << maxr << endl;
- if (DEBUG) cout << "SRMA mynprocs[maxr] " << mynprocs.at(maxr) << endl;
+ ++ mynprocs.AT(maxr);
+ if (recompose_verbose) cout << "SRMA maxr " << maxr << endl;
+ if (recompose_verbose) cout << "SRMA mynprocs[maxr] " << mynprocs.AT(maxr) << endl;
-- nregs_left;
}
assert (nregs_left == 0);
int sum_nprocs = 0;
for (int r=0; r<nregs; ++r) {
- sum_nprocs += mynprocs.at(r);
+ sum_nprocs += mynprocs.AT(r);
}
assert (sum_nprocs == newnregs);
- if (DEBUG) cout << "SRMA mynprocs " << mynprocs << endl;
+ if (recompose_verbose) cout << "SRMA mynprocs " << mynprocs << endl;
- if (DEBUG) cout << "SRMA: splitting work units" << endl;
+ if (recompose_verbose) cout << "SRMA: splitting work units" << endl;
#warning "TODO: rename newregs to regs"
#warning "TODO: rename nregs to nsuperregs"
#warning "TODO: rename newnregs to nregs"
vector<region_t> newregs;
newregs.reserve (newnregs);
- for (int r=0, p=0; r<nregs; p+=mynprocs.at(r), ++r) {
- if (DEBUG) cout << "SRMA superreg[" << r << "] " << superregs.at(r) << endl;
+ for (int r=0, p=0; r<nregs; p+=mynprocs.AT(r), ++r) {
+ if (recompose_verbose) cout << "SRMA superreg[" << r << "] " << superregs.AT(r) << endl;
bvect const dims = false;
SplitRegionsMaps_Automatic_Recursively
- (dims, p, mynprocs.at(r), superregs.at(r), newregs);
+ (dims, p, mynprocs.AT(r), superregs.AT(r), newregs);
} // for r
- if (DEBUG) cout << "SRMA newregs " << newregs << endl;
+ if (recompose_verbose) cout << "SRMA newregs " << newregs << endl;
assert (int(newregs.size()) == newnregs);
// Count components per map
vector<int> myncomps(nmaps, 0);
+#if 0
vector<int> empty_comps(nmaps, 0);
+#endif
for (int r=0; r<newnregs; ++r) {
- int const m = newregs.at(r).map - minmap;
+ int const m = newregs.AT(r).map - minmap;
assert (m>=0 and m<nmaps);
- if (not newregs.at(r).extent.empty()) {
+#if 0
+ if (not newregs.AT(r).extent.empty()) {
// Ignore empty regions, which may come from empty grid arrays
- ++ myncomps.at(m);
+ ++ myncomps.AT(m);
} else {
- ++ empty_comps.at(m);
+ ++ empty_comps.AT(m);
}
+#endif
+ ++ myncomps.AT(m);
}
vector<int> mynregs(nmaps, 0);
+#if 0
+ vector<int> empty_regs(nmaps, 0);
+#endif
for (int r=0; r<nregs; ++r) {
- int const m = superregs.at(r).map - minmap;
+ int const m = superregs.AT(r).map - minmap;
assert (m>=0 and m<nmaps);
- ++ mynregs.at(m);
+ ++ mynregs.AT(m);
+#if 0
+ if (not superregs.AT(r).extent.empty()) {
+ // Ignore empty regions, which may come from empty grid arrays
+ ++ mynregs.AT(m);
+ } else {
+ ++ empty_regs.AT(m);
+ }
+#endif
}
// Convert virtual to real processors
for (int r=0; r<newnregs; ++r) {
- newregs.at(r).processor /= ncomps;
- assert (newregs.at(r).processor >= 0 and
- newregs.at(r).processor < nprocs);
+ newregs.AT(r).processor /= ncomps;
+ assert (newregs.AT(r).processor >= 0 and
+ newregs.AT(r).processor < nprocs);
}
{
vector<int> tmpncomps(nmaps, 0);
for (int r=0; r<nregs; ++r) {
- int const m = superregs.at(r).map - minmap;
- ipfulltree * const regf = superregs.at(r).processors;
+ int const m = superregs.AT(r).map - minmap;
+ ipfulltree * const regf = superregs.AT(r).processors;
assert (regf != NULL);
for (ipfulltree::iterator fti (* regf); not fti.done(); ++ fti) {
pseudoregion_t & preg = (* fti).payload();
- preg.component = tmpncomps.at(m)++;
+ preg.component = tmpncomps.AT(m)++;
}
}
for (int m=0; m<nmaps; ++m) {
- assert (tmpncomps.at(m) == myncomps.at(m));
+#warning "TODO"
+ if (not (tmpncomps.AT(m) == myncomps.AT(m))) {
+ cout << "Recompose.cc" << endl
+ << "superregss=" << superregss << endl
+ << "regss=" << regss << endl
+ << "nregs=" << nregs << endl
+ << "newregs=" << newregs << endl
+ << "newnregs=" << newnregs << endl
+ << "m=" << m << endl
+ << "tmpncomps.AT(m)=" << tmpncomps.AT(m) << endl
+ << "myncomps.AT(m)=" << myncomps.AT(m) << endl;
+ cout << "newregs:" << endl;
+ for (int r=0; r<newnregs; ++r) {
+ int const m = newregs.AT(r).map - minmap;
+ assert (m>=0 and m<nmaps);
+ cout << "r=" << r << endl
+ << "newregs=" << newregs.at(r) << endl
+ << "newregs.extent.size()=" << newregs.at(r).extent.size() << endl
+ << "newregs.extent.empty()=" << newregs.at(r).extent.empty() << endl;
+ }
+ cout << "*regf:" << endl;
+ for (int r=0; r<nregs; ++r) {
+ int const m = superregs.AT(r).map - minmap;
+ ipfulltree * const regf = superregs.AT(r).processors;
+ cout << "r=" << r << endl
+ << "m=" << m << endl;
+ cout << "*regf=" << *regf << endl;
+ }
+ cout.flush();
+ }
+ assert (tmpncomps.AT(m) == myncomps.AT(m));
}
}
@@ -1477,34 +1713,37 @@ namespace Carpet {
// Allocate regions
assert ((int)regss.size() == nmaps);
for (int m=0; m<nmaps; ++m) {
- assert (regss.at(m).empty());
- regss.at(m).reserve (myncomps.at(m));
- superregss.at(m).clear();
- superregss.at(m).reserve (mynregs.at(m));
+ assert (regss.AT(m).empty());
+ regss.AT(m).reserve (myncomps.AT(m));
+ superregss.AT(m).clear();
+ superregss.AT(m).reserve (mynregs.AT(m));
}
// Assign regions
for (int r=0; r<newnregs; ++r) {
- int const m = newregs.at(r).map - minmap;
+ int const m = newregs.AT(r).map - minmap;
assert (m>=0 and m<nmaps);
- regss.at(m).push_back (newregs.at(r));
+ regss.AT(m).push_back (newregs.AT(r));
}
for (int r=0; r<nregs; ++r) {
- int const m = superregs.at(r).map - minmap;
+ int const m = superregs.AT(r).map - minmap;
assert (m>=0 and m<nmaps);
- superregss.at(m).push_back (superregs.at(r));
+ superregss.AT(m).push_back (superregs.AT(r));
}
// Output regions
- if (DEBUG) {
+ if (recompose_verbose) {
cout << "SRMA superregss " << superregss << endl;
cout << "SRMA regss " << regss << endl;
}
// Check sizes
for (int m=0; m<nmaps; ++m) {
- assert (int(regss.at(m).size()) == myncomps.at(m) + empty_comps.at(m));
- assert (int(superregss.at(m).size()) == mynregs.at(m));
+#if 0
+ assert (int(regss.AT(m).size()) == myncomps.AT(m) + empty_comps.AT(m));
+#endif
+ assert (int(regss.AT(m).size()) == myncomps.AT(m));
+ assert (int(superregss.AT(m).size()) == mynregs.AT(m));
}
- if (DEBUG) cout << "SRMA exit" << endl;
+ if (recompose_verbose) cout << "SRMA exit" << endl;
}
@@ -1547,30 +1786,30 @@ namespace Carpet {
}
}
vector<ibbox> bases(mglevels);
- bases.at(0) = base;
+ bases.AT(0) = base;
for (int ml=1; ml<mglevels; ++ml) {
// next finer base
- ivect const fbaselo = bases.at(ml-1).lower();
- ivect const fbasehi = bases.at(ml-1).upper();
- ivect const fbasestr = bases.at(ml-1).stride();
+ ivect const fbaselo = bases.AT(ml-1).lower();
+ ivect const fbasehi = bases.AT(ml-1).upper();
+ ivect const fbasestr = bases.AT(ml-1).stride();
// this base
ivect const basestr = fbasestr * mgfact;
ivect const baselo = fbaselo + (xpose(offset)[0] - ivect(mgfact) * xpose(offset)[0]) * fbasestr;
ivect const basehi = fbasehi + (xpose(offset)[1] - ivect(mgfact) * xpose(offset)[1]) * fbasestr;
ivect const baselo1 = baselo;
ivect const basehi1 = baselo1 + (basehi - baselo1) / basestr * basestr;
- bases.at(ml) = ibbox(baselo1, basehi1, basestr);
+ bases.AT(ml) = ibbox(baselo1, basehi1, basestr);
// next finer grid
- ivect const flo = regs.at(ml-1).extent.lower();
- ivect const fhi = regs.at(ml-1).extent.upper();
- ivect const fstr = regs.at(ml-1).extent.stride();
+ ivect const flo = regs.AT(ml-1).extent.lower();
+ ivect const fhi = regs.AT(ml-1).extent.upper();
+ ivect const fstr = regs.AT(ml-1).extent.stride();
// this grid
ivect const str = fstr * mgfact;
ivect const lo = flo + either (reg.outer_boundaries[0], (xpose(offset)[0] - ivect(mgfact) * xpose(offset)[0]) * fstr, ivect(0));
ivect const hi = fhi + either (reg.outer_boundaries[1], - (xpose(offset)[1] - ivect(mgfact) * xpose(offset)[1]) * fstr, ivect(0));
ivect const lo1 = baselo1 + (lo - baselo1 + str - 1) / str * str;
ivect const hi1 = lo1 + (hi - lo1) / str * str;
- regs.at(ml).extent = ibbox(lo1, hi1, str);
+ regs.AT(ml).extent = ibbox(lo1, hi1, str);
}
}
}
@@ -1583,25 +1822,25 @@ namespace Carpet {
{
regsss.resize (mglevels);
for (int ml=0; ml<mglevels; ++ml) {
- regsss.at(ml).resize (regss.size());
+ regsss.AT(ml).resize (regss.size());
for (int rl=0; rl<(int)regss.size(); ++rl) {
- regsss.at(ml).at(rl).resize (regss.at(rl).size());
+ regsss.AT(ml).AT(rl).resize (regss.AT(rl).size());
}
}
for (int rl=0; rl<(int)regss.size(); ++rl) {
ibbox base;
- for (int c=0; c<(int)regss.at(rl).size(); ++c) {
- base = base.expanded_containing(regss.at(rl).at(c).extent);
+ for (int c=0; c<(int)regss.AT(rl).size(); ++c) {
+ base = base.expanded_containing(regss.AT(rl).AT(c).extent);
}
- for (int c=0; c<(int)regss.at(rl).size(); ++c) {
+ for (int c=0; c<(int)regss.AT(rl).size(); ++c) {
vector<region_t> mg_regs;
- MakeMultigridBoxes (cctkGH, m, base, regss.at(rl).at(c), mg_regs);
+ MakeMultigridBoxes (cctkGH, m, base, regss.AT(rl).AT(c), mg_regs);
assert ((int)mg_regs.size() == mglevels);
for (int ml=0; ml<mglevels; ++ml) {
- regsss.at(ml).at(rl).at(c) = mg_regs.at(ml);
+ regsss.AT(ml).AT(rl).AT(c) = mg_regs.AT(ml);
}
} // for c
@@ -1615,8 +1854,46 @@ namespace Carpet {
vector<vector<vector<vector<region_t> > > > & regssss)
{
for (int m = 0; m < maps; ++m) {
- MakeMultigridBoxes (cctkGH, m, regsss.at(m), regssss.at(m));
+ MakeMultigridBoxes (cctkGH, m, regsss.AT(m), regssss.AT(m));
} // for m
}
+
+
+ static
+ void
+ ClassifyPoints (cGH const * const cctkGH, int const rl)
+ {
+ // negative: needs to be set explicitly (e.g. boundary)
+ // zero: unused (e.g. ghost)
+ // positive: needs to be evolved
+ // -1: boundary point (needs to be set explicitly)
+ // 0: unused (e.g. ghost point, or restriction target)
+ // n=1..N: evolved, used for integrator substeps i<=n
+ // (i=N..1, counting backwards; see MoL documentation)
+ // i.e.: n=1: used between time steps (i.e., should be visualised)
+ // n>1: used only while time stepping (e.g. buffer zones)
+
+ BEGIN_META_MODE(cctkGH) {
+ BEGIN_MGLEVEL_LOOP (cctkGH) {
+ ENTER_LEVEL_MODE (cctkGH, rl) {
+ BEGIN_LOCAL_MAP_LOOP(cctkGH, CCTK_GF) {
+ BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) {
+ DECLARE_CCTK_ARGUMENTS;
+#pragma omp parallel
+ LC_LOOP3 (CarpetClassifyPoints,
+ i,j,k,
+ 0,0,0, cctk_lsh[0],cctk_lsh[1],cctk_lsh[2],
+ cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
+ {
+ int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
+ point_class[ind] = 1;
+ } LC_ENDLOOP3(CarpetClassifyPoints);
+ } END_LOCAL_COMPONENT_LOOP;
+ } END_LOCAL_MAP_LOOP;
+ } LEAVE_LEVEL_MODE;
+ } END_MGLEVEL_LOOP;
+ } END_META_MODE;
+ }
+
} // namespace Carpet
diff --git a/Carpet/Carpet/src/Requirements.cc b/Carpet/Carpet/src/Requirements.cc
index 6a4a83ebc..da0beaf6b 100644
--- a/Carpet/Carpet/src/Requirements.cc
+++ b/Carpet/Carpet/src/Requirements.cc
@@ -8,7 +8,7 @@
#include <cctki_GHExtensions.h>
#include <cctki_Schedule.h>
-#include "carpet.hh"
+#include <carpet.hh>
using namespace std;
@@ -21,6 +21,14 @@ using namespace std;
typedef enum {sched_none, sched_group, sched_function} iSchedType;
typedef enum {schedpoint_misc, schedpoint_analysis} iSchedPoint;
+typedef struct t_timer
+{
+ struct t_timer *next;
+ int timer_handle;
+ char *schedule_bin;
+ int has_been_output;
+} t_timer;
+
typedef struct
{
/* Static data */
@@ -41,8 +49,7 @@ typedef struct
int *comm_groups;
/* Timer data */
-
- int timer_handle;
+ t_timer *timers;
/* Dynamic data */
int *CommOnEntry;
@@ -68,6 +75,10 @@ namespace Carpet {
// 2. Things can be provided only once, not multiple times.
// Except when they are also provided.
+ // Keep track of which time levels contain good data; modify this
+ // while time level cycling; routine specify how many time levels
+ // they require/provide
+
int CheckEntry (void * attribute, void * data);
@@ -306,8 +317,8 @@ namespace Carpet {
CheckOneGroup (cctkGH, "CCTK_TERMINATE");
}
- }; // namespace Carpet
-}; // namespace Requirements
+ } // namespace Carpet
+} // namespace Requirements
@@ -324,8 +335,8 @@ namespace Carpet {
Checkpoint ("Skipping check of schedule requirements (no flesh support)");
// do nothing
}
- }; // namespace Carpet
-}; // namespace Requirements
+ } // namespace Carpet
+} // namespace Requirements
diff --git a/Carpet/Carpet/src/Restrict.cc b/Carpet/Carpet/src/Restrict.cc
index a0a94ddd4..bad055641 100644
--- a/Carpet/Carpet/src/Restrict.cc
+++ b/Carpet/Carpet/src/Restrict.cc
@@ -2,13 +2,13 @@
#include <cmath>
#include <cstdlib>
-#include "cctk.h"
-#include "cctk_Parameters.h"
+#include <cctk.h>
+#include <cctk_Parameters.h>
-#include "ggf.hh"
-#include "gh.hh"
+#include <ggf.hh>
+#include <gh.hh>
-#include "carpet.hh"
+#include <carpet.hh>
@@ -58,7 +58,7 @@ namespace Carpet {
}
- // restricts a set of groups which all have the same vartype
+ // restrict a set of groups
static void RestrictGroups (const cGH* cctkGH, const vector<int>& groups) {
DECLARE_CCTK_PARAMETERS;
@@ -66,14 +66,14 @@ namespace Carpet {
for (comm_state state; not state.done(); state.step()) {
for (int group = 0; group < (int)groups.size(); ++group) {
- const int g = groups[group];
- for (int m=0; m<(int)arrdata.at(g).size(); ++m) {
+ const int g = groups.AT(group);
+ for (int m=0; m<(int)arrdata.AT(g).size(); ++m) {
// use background time here (which may not be modified
// by the user)
- const CCTK_REAL time = vtt.at(m)->time (tl, reflevel, mglevel);
+ const CCTK_REAL time = vtt.AT(m)->time (tl, reflevel, mglevel);
- const CCTK_REAL time1 = vtt.at(m)->time (0, reflevel, mglevel);
+ const CCTK_REAL time1 = vtt.AT(m)->time (0, reflevel, mglevel);
const CCTK_REAL time2 =
(cctkGH->cctk_time - cctk_initial_time) / delta_time;
const CCTK_REAL time0 =
@@ -81,8 +81,8 @@ namespace Carpet {
const CCTK_REAL eps = 1.0e-12;
assert (abs(time1 - time2) <= eps * time0);
- for (int v = 0; v < (int)arrdata.at(g).at(m).data.size(); ++v) {
- ggf *const gv = arrdata.at(g).at(m).data.at(v);
+ for (int v = 0; v < (int)arrdata.AT(g).AT(m).data.size(); ++v) {
+ ggf *const gv = arrdata.AT(g).AT(m).data.AT(v);
gv->ref_restrict_all (state, tl, reflevel, mglevel, time);
}
}
diff --git a/Carpet/Carpet/src/ScheduleWrapper.cc b/Carpet/Carpet/src/ScheduleWrapper.cc
index c782c0f91..624dacbf9 100644
--- a/Carpet/Carpet/src/ScheduleWrapper.cc
+++ b/Carpet/Carpet/src/ScheduleWrapper.cc
@@ -1,10 +1,12 @@
#include <cassert>
#include <list>
-#include "cctk.h"
-#include "cctk_Arguments.h"
+#include <cctk.h>
+#include <cctk_Arguments.h>
+
+#include <carpet.hh>
+
-#include "carpet.hh"
namespace Carpet {
diff --git a/Carpet/Carpet/src/SetupGH.cc b/Carpet/Carpet/src/SetupGH.cc
index e9f0df32f..6be1c7008 100644
--- a/Carpet/Carpet/src/SetupGH.cc
+++ b/Carpet/Carpet/src/SetupGH.cc
@@ -4,6 +4,7 @@
#include <cstdlib>
#include <cstring>
#include <iostream>
+#include <limits>
#include <sstream>
#include <stack>
#include <string>
@@ -13,7 +14,6 @@
#include <cctk.h>
#include <cctk_Parameters.h>
-
#include <util_ErrorCodes.h>
#include <util_Network.h>
#include <util_Table.h>
@@ -23,10 +23,11 @@
#include <dist.hh>
#include <ggf.hh>
#include <gh.hh>
+#include <mpi_string.hh>
#include <region.hh>
#include <vect.hh>
-#include "carpet.hh"
+#include <carpet.hh>
@@ -181,9 +182,6 @@ namespace Carpet {
{
DECLARE_CCTK_PARAMETERS;
- // Check and/or modify system limits
- SetSystemLimits ();
-
// Some statistics
{
int const nprocs = dist::size();
@@ -191,14 +189,42 @@ namespace Carpet {
dist::set_num_threads (num_threads);
int const mynthreads = dist::num_threads();
int const nthreads_total = dist::total_num_threads();
+#ifdef CCTK_MPI
+ CCTK_VInfo (CCTK_THORNSTRING,
+ "MPI is enabled");
CCTK_VInfo (CCTK_THORNSTRING,
"Carpet is running on %d processes", nprocs);
CCTK_VInfo (CCTK_THORNSTRING,
"This is process %d", myproc);
+#else
+ CCTK_VInfo (CCTK_THORNSTRING,
+ "MPI is disabled");
+#endif
+ char const * const OMP_NUM_THREADS = getenv ("OMP_NUM_THREADS");
+#ifdef _OPENMP
+ CCTK_VInfo (CCTK_THORNSTRING,
+ "OpenMP is enabled");
+ if (not OMP_NUM_THREADS and num_threads == -1) {
+ CCTK_VWarn (CCTK_WARN_COMPLAIN, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "Although OpenMP is enabled, neither the environment variable OMP_NUM_THREADS nor the parameter Carpet::num_threads are set. A system-specific default value is used instead.");
+ }
CCTK_VInfo (CCTK_THORNSTRING,
"This process contains %d threads", mynthreads);
CCTK_VInfo (CCTK_THORNSTRING,
"There are %d threads in total", nthreads_total);
+#else
+ CCTK_VInfo (CCTK_THORNSTRING,
+ "OpenMP is disabled");
+ int const omp_num_threads = OMP_NUM_THREADS ? atoi (OMP_NUM_THREADS) : 0;
+ if (omp_num_threads > 0) {
+ CCTK_VWarn (CCTK_WARN_ALERT, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "Although OpenMP is disabled, the environment variable OMP_NUM_THREADS is set to %d. It will be ignored.", omp_num_threads);
+ }
+ if (num_threads > 0) {
+ CCTK_VWarn (CCTK_WARN_ALERT, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "Although OpenMP is disabled, the parameter Carpet::num_threads is set to %d. It will be ignored.", num_threads);
+ }
+#endif
#if 0
// Do not call Util_GetHostName. Certain InfiniBand libraries
@@ -209,7 +235,7 @@ namespace Carpet {
char hostnamebuf[1000];
Util_GetHostName (hostnamebuf, sizeof hostnamebuf);
string const hostname (hostnamebuf);
- vector <string> hostnames = AllGatherString (dist::comm(), hostname);
+ vector <string> hostnames = allgather_string (dist::comm(), hostname);
// Collect process ids
int const mypid = getpid ();
vector <int> pids (nprocs);
@@ -227,7 +253,7 @@ namespace Carpet {
for (int n = 0; n < nprocs; ++ n) {
CCTK_VInfo (CCTK_THORNSTRING,
" %6d: %s, pid=%d, num_threads=%d",
- n, hostnames.at(n).c_str(), pids.at(n), nthreads.at(n));
+ n, hostnames.AT(n).c_str(), pids.AT(n), nthreads.AT(n));
}
}
#endif
@@ -240,6 +266,7 @@ namespace Carpet {
mc_grouptype = -1;
map = -1;
component = -1;
+ cctkGH->cctk_mode = CCTK_MODE_META;
// Say hello
Waypoint ("Setting up the grid hierarchy");
@@ -275,7 +302,9 @@ namespace Carpet {
do_warn_about_storage = false; // This is enabled later
if (enable_all_storage) {
- enable_storage_for_all_groups (cctkGH);
+ if (not enable_no_storage) {
+ enable_storage_for_all_groups (cctkGH);
+ }
}
Waypoint ("Done with setting up the grid hierarchy");
@@ -327,8 +356,11 @@ namespace Carpet {
if (CCTK_EQUALS (space_refinement_factors, "")) {
// Calculate them from the default refinement factor
spacereffacts.resize (maxreflevels);
- for (int n=0; n<maxreflevels; ++n) {
- spacereffacts.at(n) = ivect (ipow (int(refinement_factor), n));
+ assert (ipow (2, 0) == 1);
+ assert (ipow (3, 1) == 3);
+ assert (ipow (4, 2) == 16);
+ for (int rl=0; rl<maxreflevels; ++rl) {
+ spacereffacts.AT(rl) = ivect (ipow (int(refinement_factor), rl));
}
} else {
// Read them from the parameter
@@ -342,10 +374,10 @@ namespace Carpet {
}
// TODO: turn these into real error messages
assert (int(spacereffacts.size()) >= maxreflevels);
- assert (all (spacereffacts.front() == 1));
- for (int n=1; n<maxreflevels; ++n) {
- assert (all (spacereffacts.at(n) >= spacereffacts.at(n-1)));
- assert (all (spacereffacts.at(n) % spacereffacts.at(n-1) == 0));
+ assert (all (spacereffacts.AT(0) == 1));
+ for (int rl=1; rl<maxreflevels; ++rl) {
+ assert (all (spacereffacts.AT(rl) >= spacereffacts.AT(rl-1)));
+ assert (all (spacereffacts.AT(rl) % spacereffacts.AT(rl-1) == 0));
}
spacereffacts.resize (maxreflevels);
@@ -353,8 +385,8 @@ namespace Carpet {
if (CCTK_EQUALS (time_refinement_factors, "")) {
// Calculate them from the default refinement factor
timereffacts.resize (maxreflevels);
- for (int n=0; n<maxreflevels; ++n) {
- timereffacts.at(n) = ipow (int(refinement_factor), n);
+ for (int rl=0; rl<maxreflevels; ++rl) {
+ timereffacts.AT(rl) = ipow (int(refinement_factor), rl);
}
} else {
// Read them from the parameter
@@ -367,16 +399,16 @@ namespace Carpet {
}
// TODO: turn these into real error messages
assert (int(timereffacts.size()) >= maxreflevels);
- assert (timereffacts.front() == 1);
- for (int n=1; n<maxreflevels; ++n) {
- assert (timereffacts.at(n) >= timereffacts.at(n-1));
- assert (timereffacts.at(n) % timereffacts.at(n-1) == 0);
+ assert (timereffacts.AT(0) == 1);
+ for (int rl=1; rl<maxreflevels; ++rl) {
+ assert (timereffacts.AT(rl) >= timereffacts.AT(rl-1));
+ assert (timereffacts.AT(rl) % timereffacts.AT(rl-1) == 0);
}
timereffacts.resize (maxreflevels);
// Calculate the maximum refinement factors
- maxtimereflevelfact = timereffacts.at (maxreflevels-1);
- maxspacereflevelfact = spacereffacts.at (maxreflevels-1);
+ maxtimereflevelfact = timereffacts.AT(maxreflevels-1);
+ maxspacereflevelfact = spacereffacts.AT(maxreflevels-1);
}
@@ -502,7 +534,7 @@ namespace Carpet {
// Allocate grid hierarchy
vhh.resize(maps);
- vhh.at(m) = new gh (spacereffacts, refcentering,
+ vhh.AT(m) = new gh (spacereffacts, refcentering,
convergence_factor, mgcentering,
baseexts, nboundaryzones);
}
@@ -541,7 +573,7 @@ namespace Carpet {
assert (all (all (buffers >= 0)));
vdd.resize(maps);
- vdd.at(m) = new dh (* vhh.at(m),
+ vdd.AT(m) = new dh (* vhh.AT(m),
ghosts, buffers,
prolongation_order_space);
@@ -559,7 +591,7 @@ namespace Carpet {
DECLARE_CCTK_PARAMETERS;
vtt.resize (maps);
- vtt.at(m) = new th (* vhh.at(m),
+ vtt.AT(m) = new th (* vhh.AT(m),
timereffacts,
1.0);
}
@@ -573,7 +605,7 @@ namespace Carpet {
vector<vector<vector<region_t> > > superregsss (maps);
for (int m=0; m<maps; ++m) {
- set_base_extent (m, superregsss.at(m));
+ set_base_extent (m, superregsss.AT(m));
}
vector<vector<vector<vector<region_t> > > > regssss;
@@ -582,40 +614,45 @@ namespace Carpet {
for (int m=0; m<maps; ++m) {
// Check the regions
- CheckRegions (regssss.at(m));
+ CheckRegions (regssss.AT(m));
// Recompose grid hierarchy
- vhh.at(m)->regrid (superregsss.at(m), regssss.at(m));
+ vhh.AT(m)->regrid (superregsss.AT(m), regssss.AT(m), false);
int const rl = 0;
- vhh.at(m)->recompose (rl, false);
+ vhh.AT(m)->recompose (rl, false);
+ vhh.AT(m)->regrid_free (false);
} // for m
- CCTK_INFO ("Grid structure (grid points):");
- for (int ml=0; ml<mglevels; ++ml) {
- int const rl = 0;
- for (int m=0; m<maps; ++m) {
- for (int c=0; c<vhh.at(m)->components(rl); ++c) {
- ibbox const ext = vhh.at(m)->extent(ml,rl,c);
- ivect const lower = ext.lower();
- ivect const upper = ext.upper();
- int const convfact = ipow(mgfact, ml);
- assert (all(lower % maxspacereflevelfact == 0));
- assert (all(upper % maxspacereflevelfact == 0));
- assert (all(((upper - lower) / maxspacereflevelfact) % convfact == 0));
- cout << " [" << ml << "][" << rl << "][" << m << "][" << c << "]"
- << " exterior: "
- << "proc "
- << vhh.at(m)->processor(rl,c)
- << " "
- << lower / maxspacereflevelfact
- << " : "
- << upper / maxspacereflevelfact
- << " ("
- << (upper - lower) / maxspacereflevelfact / convfact + 1
- << ") "
- << prod ((upper - lower) / maxspacereflevelfact / convfact + 1)
- << endl;
+ regridding_epoch = 0;
+
+ if (verbose or veryverbose) {
+ CCTK_INFO ("Grid structure (grid points):");
+ for (int ml=0; ml<mglevels; ++ml) {
+ int const rl = 0;
+ for (int m=0; m<maps; ++m) {
+ for (int c=0; c<vhh.AT(m)->components(rl); ++c) {
+ ibbox const ext = vhh.AT(m)->extent(ml,rl,c);
+ ivect const lower = ext.lower();
+ ivect const upper = ext.upper();
+ int const convfact = ipow(mgfact, ml);
+ assert (all(lower % maxspacereflevelfact == 0));
+ assert (all(upper % maxspacereflevelfact == 0));
+ assert (all(((upper - lower) / maxspacereflevelfact) % convfact == 0));
+ cout << " [" << ml << "][" << rl << "][" << m << "][" << c << "]"
+ << " exterior: "
+ << "proc "
+ << vhh.AT(m)->processor(rl,c)
+ << " "
+ << lower / maxspacereflevelfact
+ << " : "
+ << upper / maxspacereflevelfact
+ << " ("
+ << (upper - lower) / maxspacereflevelfact / convfact + 1
+ << ") "
+ << prod ((upper - lower) / maxspacereflevelfact / convfact + 1)
+ << endl;
+ }
}
}
}
@@ -623,7 +660,7 @@ namespace Carpet {
// Assert that all maps have one refinement level
reflevels = 1;
for (int m=0; m<maps; ++m) {
- assert (vhh.at(m)->reflevels() == reflevels);
+ assert (vhh.AT(m)->reflevels() == reflevels);
}
}
@@ -638,13 +675,13 @@ namespace Carpet {
// Create one refinement level
int const rl = 0;
regss.resize(1);
- vector<region_t> & regs = regss.at(rl);
+ vector<region_t> & regs = regss.AT(rl);
if (CCTK_EQUALS (base_extents, "")) {
// Default: one grid component covering everything
region_t reg;
- reg.extent = vhh.at(m)->baseextents.at(0).at(0);
+ reg.extent = vhh.AT(m)->baseextents.AT(0).AT(0);
reg.outer_boundaries = b2vect (bvect (true));
reg.map = m;
regs.push_back (reg);
@@ -677,8 +714,8 @@ namespace Carpet {
for (size_t n=0; n<exts.size(); ++n) {
region_t reg;
- reg.extent = exts.at(n);
- reg.outer_boundaries = xpose(obs.at(n));
+ reg.extent = exts.AT(n);
+ reg.outer_boundaries = xpose(obs.AT(n));
reg.map = m;
regs.push_back (reg);
}
@@ -711,17 +748,17 @@ namespace Carpet {
assert (gdata.dim == dim);
// Set up one refinement level
- groupdata.at(group).activetimelevels.resize(mglevels);
+ groupdata.AT(group).activetimelevels.resize(mglevels);
for (int ml=0; ml<mglevels; ++ml) {
- groupdata.at(group).activetimelevels.at(ml).resize(1);
+ groupdata.AT(group).activetimelevels.AT(ml).resize(1);
}
// Grid function groups use the global grid descriptors
- arrdata.at(group).resize(maps);
+ arrdata.AT(group).resize(maps);
for (int m=0; m<maps; ++m) {
- arrdata.at(group).at(m).hh = vhh.at(m);
- arrdata.at(group).at(m).dd = vdd.at(m);
- arrdata.at(group).at(m).tt = vtt.at(m);
+ arrdata.AT(group).AT(m).hh = vhh.AT(m);
+ arrdata.AT(group).AT(m).dd = vdd.AT(m);
+ arrdata.AT(group).AT(m).tt = vtt.AT(m);
}
break;
@@ -735,13 +772,13 @@ namespace Carpet {
assert (gdata.dim >= 0 and gdata.dim <= dim);
// Use only one refinement level for grid arrays
- groupdata.at(group).activetimelevels.resize(mglevels);
+ groupdata.AT(group).activetimelevels.resize(mglevels);
for (int ml=0; ml<mglevels; ++ml) {
- groupdata.at(group).activetimelevels.at(ml).resize(1);
+ groupdata.AT(group).activetimelevels.AT(ml).resize(1);
}
// Use only one map for grid arrays
- arrdata.at(group).resize(1);
+ arrdata.AT(group).resize(1);
ivect sizes;
i2vect ghosts;
@@ -788,33 +825,33 @@ namespace Carpet {
ivect const str(1);
ibbox const baseext(lb, ub, str);
vector <vector <ibbox> > baseexts (1);
- baseexts.at(0).resize (1);
- baseexts.at(0).at(0) = baseext;
+ baseexts.AT(0).resize (1);
+ baseexts.AT(0).AT(0) = baseext;
i2vect const nboundaryzones(0);
// One refinement level
vector<int> grouptimereffacts(1);
- grouptimereffacts.at(0) = 1;
+ grouptimereffacts.AT(0) = 1;
vector<ivect> groupspacereffacts(1);
- groupspacereffacts.at(0) = ivect(1);
+ groupspacereffacts.AT(0) = ivect(1);
// There is only one map
- int const m=0;
+ int const m = 0;
- arrdata.at(group).at(m).hh =
+ arrdata.AT(group).AT(m).hh =
new gh (groupspacereffacts, vertex_centered,
convergence_factor, vertex_centered,
baseexts, nboundaryzones);
i2vect const buffers = i2vect (0);
int const my_prolongation_order_space = 0;
- arrdata.at (group).at(m).dd =
- new dh (*arrdata.at (group).at(m).hh,
+ arrdata.AT(group).AT(m).dd =
+ new dh (*arrdata.AT(group).AT(m).hh,
ghosts, buffers, my_prolongation_order_space);
CCTK_REAL const basedelta = 1.0;
- arrdata.at (group).at(m).tt =
- new th (*arrdata.at (group).at(m).hh, grouptimereffacts, basedelta);
+ arrdata.AT(group).AT(m).tt =
+ new th (*arrdata.AT(group).AT(m).hh, grouptimereffacts, basedelta);
}
@@ -826,15 +863,26 @@ namespace Carpet {
ivect const & convpowers,
ivect const & convoffsets)
{
+#if 0
+ // Do not set up anything for groups that have zero variables
+ // TODO: To do this, need to change modes.cc to provide illegal
+ // entries for GroupDynamicData
+ int const nvars = CCTK_NumVarsInGroupI (group);
+ assert (nvars >= 0);
+ if (nvars == 0) return;
+#endif
+
// Set refinement structure for scalars and arrays
vector<region_t> superregs(1);
int const m = 0;
int const rl = 0;
- int const c = 0;
- superregs.at(c).extent =
- arrdata.at(group).at(m).hh->baseextents.at(rl).at(c);
- superregs.at(c).outer_boundaries = b2vect(true);
- superregs.at(c).map = m;
+ {
+ int const c = 0;
+ superregs.AT(c).extent =
+ arrdata.AT(group).AT(m).hh->baseextents.AT(rl).AT(c);
+ superregs.AT(c).outer_boundaries = b2vect(true);
+ superregs.AT(c).map = m;
+ }
vector<region_t> regs;
// Split it into components, one for each processor
@@ -871,7 +919,7 @@ namespace Carpet {
region_t empty;
empty.extent = ebox;
empty.processor = c;
- regs.at(c) = empty;
+ regs.AT(c) = empty;
}
}
}
@@ -881,21 +929,21 @@ namespace Carpet {
int const nprocs = CCTK_nProcs (cctkGH);
vector<bool> used (nprocs, false);
for (int c=0; c<nprocs; ++c) {
- int const p = regs.at(c).processor;
+ int const p = regs.AT(c).processor;
assert (p >= 0 and p < nprocs);
- assert (not used.at(p));
- used.at(p) = true;
+ assert (not used.AT(p));
+ used.AT(p) = true;
}
for (int c=0; c<nprocs; ++c) {
- assert (used.at(c));
+ assert (used.AT(c));
}
}
// Only one refinement level
vector<vector<region_t> > superregss(1);
- superregss.at(rl) = superregs;
+ superregss.AT(rl) = superregs;
vector<vector<region_t> > regss(1);
- regss.at(rl) = regs;
+ regss.AT(rl) = regs;
// Create all multigrid levels
vector<vector<vector<region_t> > > regsss (mglevels);
@@ -906,32 +954,31 @@ namespace Carpet {
offset[0][d] = 0;
offset[1][d] = convoffsets[d];
}
- regsss.at(0) = regss;
+ regsss.AT(0) = regss;
for (int ml=1; ml<mglevels; ++ml) {
- int const rl = 0;
- for (int c=0; c<int(regss.at(rl).size()); ++c) {
+ for (int c=0; c<int(regss.AT(rl).size()); ++c) {
// this base
ivect const baselo = ivect(0);
ivect const baselo1 = baselo;
// next finer grid
- ivect const flo = regsss.at(ml-1).at(rl).at(c).extent.lower();
- ivect const fhi = regsss.at(ml-1).at(rl).at(c).extent.upper();
- ivect const fstr = regsss.at(ml-1).at(rl).at(c).extent.stride();
+ ivect const flo = regsss.AT(ml-1).AT(rl).AT(c).extent.lower();
+ ivect const fhi = regsss.AT(ml-1).AT(rl).AT(c).extent.upper();
+ ivect const fstr = regsss.AT(ml-1).AT(rl).AT(c).extent.stride();
// this grid
ivect const str = fstr * mgfact1;
ivect const lo = flo +
- either (regsss.at(ml-1).at(rl).at(c).outer_boundaries[0],
+ either (regsss.AT(ml-1).AT(rl).AT(c).outer_boundaries[0],
+ (offset[0] - mgfact1 * offset[0]) * fstr,
ivect(0));
ivect const hi = fhi +
- either (regsss.at(ml-1).at(rl).at(c).outer_boundaries[1],
+ either (regsss.AT(ml-1).AT(rl).AT(c).outer_boundaries[1],
- (offset[1] - mgfact1 * offset[1]) * fstr,
ivect(0));
ivect const lo1 = baselo1 + (lo - baselo1 + str - 1) / str * str;
ivect const hi1 = lo1 + (hi - lo1) / str * str;
- regsss.at(ml).at(rl).at(c) = regsss.at(ml-1).at(rl).at(c);
- regsss.at(ml).at(rl).at(c).extent = ibbox(lo1, hi1, str);
+ regsss.AT(ml).AT(rl).AT(c) = regsss.AT(ml-1).AT(rl).AT(c);
+ regsss.AT(ml).AT(rl).AT(c).extent = ibbox(lo1, hi1, str);
}
}
@@ -940,8 +987,9 @@ namespace Carpet {
char * const groupname = CCTK_GroupName (group);
assert (groupname);
Checkpoint ("Recomposing grid array group \"%s\"...", groupname);
- arrdata.at(group).at(0).hh->regrid (superregss, regsss);
- arrdata.at(group).at(0).hh->recompose (0, false);
+ arrdata.AT(group).AT(0).hh->regrid (superregss, regsss, false);
+ arrdata.AT(group).AT(0).hh->recompose (0, false);
+ arrdata.AT(group).AT(0).hh->regrid_free (false);
Checkpoint ("Done recomposing grid array group \"%s\".", groupname);
free (groupname);
}
@@ -957,25 +1005,25 @@ namespace Carpet {
DECLARE_CCTK_PARAMETERS;
// Initialise group information
- groupdata.at(group).info.dim = gdata.dim;
- groupdata.at(group).info.gsh = new int [dim];
- groupdata.at(group).info.lsh = new int [dim];
- groupdata.at(group).info.lbnd = new int [dim];
- groupdata.at(group).info.ubnd = new int [dim];
- groupdata.at(group).info.bbox = new int [2*dim];
- groupdata.at(group).info.nghostzones = new int [dim];
-
- groupdata.at(group).transport_operator =
+ groupdata.AT(group).info.dim = gdata.dim;
+ groupdata.AT(group).info.gsh = new int [dim];
+ groupdata.AT(group).info.lsh = new int [dim];
+ groupdata.AT(group).info.lbnd = new int [dim];
+ groupdata.AT(group).info.ubnd = new int [dim];
+ groupdata.AT(group).info.bbox = new int [2*dim];
+ groupdata.AT(group).info.nghostzones = new int [dim];
+
+ groupdata.AT(group).transport_operator =
get_transport_operator (cctkGH, group, gdata);
- groupdata.at(group).info.activetimelevels = 0;
+ groupdata.AT(group).info.activetimelevels = 0;
// Initialise group variables
- for (size_t m=0; m<arrdata.at(group).size(); ++m) {
+ for (size_t m=0; m<arrdata.AT(group).size(); ++m) {
- arrdata.at(group).at(m).data.resize(CCTK_NumVarsInGroupI(group));
- for (size_t var=0; var<arrdata.at(group).at(m).data.size(); ++var) {
- arrdata.at(group).at(m).data.at(var) = 0;
+ arrdata.AT(group).AT(m).data.resize(CCTK_NumVarsInGroupI(group));
+ for (size_t var=0; var<arrdata.AT(group).AT(m).data.size(); ++var) {
+ arrdata.AT(group).AT(m).data.AT(var) = 0;
}
}
@@ -989,14 +1037,14 @@ namespace Carpet {
// Allocate level times
leveltimes.resize (mglevels);
for (int ml=0; ml<mglevels; ++ml) {
- leveltimes.at(ml).resize (1);
+ leveltimes.AT(ml).resize (1);
}
// Allocate orgin and spacings
origin_space.resize (maps);
delta_space.resize (maps);
for (int m=0; m<maps; ++m) {
- origin_space.at(m).resize (mglevels);
+ origin_space.AT(m).resize (mglevels);
}
// Current state
@@ -1009,11 +1057,12 @@ namespace Carpet {
}
// Set up things as if in local mode
- mglevel = 0;
- reflevel = 0;
- mc_grouptype = CCTK_GF;
- map = 0;
- component = 0;
+ mglevel = 0;
+ reflevel = 0;
+ mc_grouptype = CCTK_GF;
+ map = 0;
+ component = 0;
+ local_component = -1;
// Leave everything, so that everything is set up correctly
leave_local_mode (cctkGH);
@@ -1402,17 +1451,36 @@ namespace Carpet {
// Sanity check
assert (all (npoints <= INT_MAX));
+#if 0
int max = INT_MAX;
for (int d=0; d<dim; ++d) {
assert (npoints[d] <= max);
max /= npoints[d];
}
+#endif
+ {
+ CCTK_REAL const total_npoints = prod(rvect(npoints));
+ CCTK_REAL const size_max = numeric_limits<ibbox::size_type>::max();
+ if (total_npoints > size_max) {
+ CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "The domain for map %d contains %g grid points. This number is larger than the maximum number supported by Carpet (%g).",
+ m, double(total_npoints), double(size_max));
+ }
+ CCTK_REAL const int_max = numeric_limits<int>::max();
+ if (total_npoints > int_max) {
+ if (dist::rank() == 0) {
+ CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "The domain for map %d contains %g grid points. This number is larger than the maximum number that can be represented as an integer (%g). This can lead to strange problems in thorns which try to calculate the total number of grid points.",
+ m, double(total_npoints), double(int_max));
+ }
+ }
+ }
// Save domain specification
domainspecs.resize(maps);
- domainspecs.at(m).exterior_min = exterior_min;
- domainspecs.at(m).exterior_max = exterior_max;
- domainspecs.at(m).npoints = npoints;
+ domainspecs.AT(m).exterior_min = exterior_min;
+ domainspecs.AT(m).exterior_max = exterior_max;
+ domainspecs.AT(m).npoints = npoints;
}
@@ -1432,26 +1500,53 @@ namespace Carpet {
assert (baseextents.empty());
baseextents.resize (num_convergence_levels);
for (int ml=0; ml<num_convergence_levels; ++ml) {
- baseextents.at(ml).resize (maxreflevels);
+ baseextents.AT(ml).resize (maxreflevels);
for (int rl=0; rl<maxreflevels; ++rl) {
if (ml == 0) {
if (rl == 0) {
// Use base extent
- baseextents.at(ml).at(rl) = baseextent;
+ baseextents.AT(ml).AT(rl) = baseextent;
} else {
// Refine next coarser refinement level
- assert (refcentering == vertex_centered);
- assert (not any (any (is_staggered)));
- i2vect const bnd_shift =
- (nboundaryzones - 1) * i2vect (not is_internal) + shiftout;
- ibbox const & cbox = baseextents.at(ml).at(rl-1);
- ibbox const cbox_phys = cbox.expand (- bnd_shift);
- assert (all (baseextent.stride() % spacereffacts.at(rl-1) == 0));
- ivect const fstride = baseextent.stride() / spacereffacts.at(rl);
- ibbox const fbox_phys =
- ibbox (cbox_phys.lower(), cbox_phys.upper(), fstride);
- ibbox const fbox = fbox_phys.expand (bnd_shift);
- baseextents.at(ml).at(rl) = fbox;
+ if (refcentering == vertex_centered) {
+ ibbox const & cbox = baseextents.AT(ml).AT(rl-1);
+ assert (not any (any (is_staggered)));
+ i2vect const bnd_shift =
+ (nboundaryzones - 1) * i2vect (not is_internal) + shiftout;
+ ibbox const cbox_phys = cbox.expand (- bnd_shift);
+ assert (all (baseextent.stride() % spacereffacts.AT(rl-1) == 0));
+ ivect const fstride = baseextent.stride() / spacereffacts.AT(rl);
+ ibbox const fbox_phys =
+ ibbox (cbox_phys.lower(), cbox_phys.upper(), fstride);
+ ibbox const fbox = fbox_phys.expand (bnd_shift);
+ baseextents.AT(ml).AT(rl) = fbox;
+ } else {
+ ibbox const & cbox = baseextents.AT(ml).AT(rl-1);
+ assert (all (all (is_staggered)));
+ ivect const cstride = cbox.stride();
+ assert (all (cstride % 2 == 0));
+ i2vect const bnd_shift_cstride =
+ + (cstride / 2 - i2vect (is_internal) * cstride)
+ + (nboundaryzones - 1) * i2vect (not is_internal) * cstride
+ + shiftout * cstride;
+ ibbox const cbox_phys
+ (cbox.lower() - (- bnd_shift_cstride[0]),
+ cbox.upper() + (- bnd_shift_cstride[1]),
+ cbox.stride());
+ assert (all (baseextent.stride() % spacereffacts.AT(rl-1) == 0));
+ ivect const fstride = baseextent.stride() / spacereffacts.AT(rl);
+ ibbox const fbox_phys =
+ ibbox (cbox_phys.lower(), cbox_phys.upper(), fstride);
+ assert (all (cstride % fstride == 0));
+ assert (all (all (bnd_shift_cstride % (cstride / fstride) == 0)));
+ i2vect const bnd_shift_fstride =
+ bnd_shift_cstride / (cstride / fstride);
+ ibbox const fbox
+ (fbox_phys.lower() - bnd_shift_fstride[0],
+ fbox_phys.upper() + bnd_shift_fstride[1],
+ fbox_phys.stride());
+ baseextents.AT(ml).AT(rl) = fbox;
+ }
}
} else {
// Coarsen next finer convergence level
@@ -1459,13 +1554,13 @@ namespace Carpet {
assert (not any (any (is_staggered)));
i2vect const bnd_shift =
(nboundaryzones - 1) * i2vect (not is_internal) + shiftout;
- ibbox const & fbox = baseextents.at(ml-1).at(rl);
+ ibbox const & fbox = baseextents.AT(ml-1).AT(rl);
ibbox const fbox_phys = fbox.expand (- bnd_shift);
ivect const cstride = fbox.stride() * int (convergence_factor);
ibbox const cbox_phys =
ibbox (fbox_phys.lower(), fbox_phys.upper(), cstride);
ibbox const cbox = cbox_phys.expand (bnd_shift);
- baseextents.at(ml).at(rl) = cbox;
+ baseextents.AT(ml).AT(rl) = cbox;
}
}
}
@@ -1493,10 +1588,10 @@ namespace Carpet {
vector<vector<region_t> > regss(1);
// Distribute onto the processors
- SplitRegions (cctkGH, superregsss.at(m).at(rl), regss.at(rl));
+ SplitRegions (cctkGH, superregsss.AT(m).AT(rl), regss.AT(rl));
// Create all multigrid levels
- MakeMultigridBoxes (cctkGH, m, regss, regssss.at(m));
+ MakeMultigridBoxes (cctkGH, m, regss, regssss.AT(m));
} // for m
} else {
@@ -1506,7 +1601,7 @@ namespace Carpet {
vector<vector<region_t> > superregss(maps);
for (int m=0; m<maps; ++m) {
- superregss.at(m) = superregsss.at(m).at(rl);
+ superregss.AT(m) = superregsss.AT(m).AT(rl);
}
vector<vector<region_t> > regss(maps);
@@ -1514,9 +1609,9 @@ namespace Carpet {
vector<vector<vector<region_t> > > regsss(maps);
for (int m=0; m<maps; ++m) {
- superregsss.at(m).at(rl) = superregss.at(m);
- regsss.at(m).resize(1);
- regsss.at(m).at(rl) = regss.at(m);
+ superregsss.AT(m).AT(rl) = superregss.AT(m);
+ regsss.AT(m).resize(1);
+ regsss.AT(m).AT(rl) = regss.AT(m);
}
// Create all multigrid levels
@@ -1583,8 +1678,7 @@ namespace Carpet {
ivect baseconvpowers = convpowers * int(basemglevel);
rvect real_sizes =
- (((rvect (sizes)
- - rvect (convoffsets))
+ (((rvect (sizes) - rvect (convoffsets))
/ ipow (rvect(convergence_factor), baseconvpowers))
+ rvect (convoffsets));
// Do not modify extra dimensions
@@ -1720,8 +1814,8 @@ namespace Carpet {
int num_gf_vars = 0;
vector<int> num_array_groups(dim+1), num_array_vars(dim+1);
for (int d=0; d<=dim; ++d) {
- num_array_groups.at(d) = 0;
- num_array_vars.at(d) = 0;
+ num_array_groups.AT(d) = 0;
+ num_array_vars.AT(d) = 0;
}
for (int group=0; group<CCTK_NumGroups(); ++group) {
@@ -1737,8 +1831,8 @@ namespace Carpet {
case CCTK_SCALAR:
case CCTK_ARRAY:
assert (gdata.dim<=dim);
- num_array_groups.at(gdata.dim) += 1;
- num_array_vars.at(gdata.dim) += gdata.numvars * gdata.numtimelevels;
+ num_array_groups.AT(gdata.dim) += 1;
+ num_array_vars.AT(gdata.dim) += gdata.numvars * gdata.numtimelevels;
break;
default:
assert (0);
@@ -1751,11 +1845,11 @@ namespace Carpet {
num_gf_vars, num_gf_groups);
CCTK_VInfo (CCTK_THORNSTRING,
" There are %d grid scalars in %d groups",
- num_array_vars.at(0), num_array_groups.at(0));
+ num_array_vars.AT(0), num_array_groups.AT(0));
for (int d=1; d<=dim; ++d) {
CCTK_VInfo (CCTK_THORNSTRING,
" There are %d %d-dimensional grid arrays in %d groups",
- num_array_vars.at(d), d, num_array_groups.at(d));
+ num_array_vars.AT(d), d, num_array_groups.AT(d));
}
CCTK_VInfo (CCTK_THORNSTRING,
" (The number of variables counts all time levels)");
@@ -1770,13 +1864,31 @@ namespace Carpet {
{
assert (group>=0 and group<CCTK_NumGroups());
- if (CCTK_GroupTypeI(group) != CCTK_GF) {
+ if (gdata.grouptype != CCTK_GF) {
// Ignore everything but true grid functions
- return op_error;
+ return op_copy;
}
bool const can_transfer = can_transfer_variable_type (cctkGH, group, gdata);
+ // Catch a common error (using the tag "prolongate" instead of
+ // "prolongation")
+ {
+ int const prolong_length = Util_TableGetString
+ (gdata.tagstable, 0, NULL, "Prolongate");
+ if (prolong_length >= 0) {
+ char * const groupname = CCTK_GroupName (group);
+ CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "Group \"%s\" contains the illegal tag \"Prolongate\". (Did you mean \"Prolongation\" instead?)",
+ groupname);
+ free (groupname);
+ } else if (prolong_length == UTIL_ERROR_TABLE_NO_SUCH_KEY) {
+ // good -- do nothing
+ } else {
+ assert (0);
+ }
+ }
+
// Get prolongation method
char prolong_string[1000];
bool have_prolong_string = false;
@@ -1861,6 +1973,7 @@ namespace Carpet {
if (not have_prolong_string) {
if (can_transfer) {
// Use the default
+#if 0
if (gdata.numtimelevels == 1) {
// Only one time level:
char * const groupname = CCTK_GroupName (group);
@@ -1873,6 +1986,8 @@ namespace Carpet {
// Several time levels: use the default
return op_Lagrange;
}
+#endif
+ return op_Lagrange;
} else {
if (gdata.grouptype == CCTK_GF) {
char * const groupname = CCTK_GroupName (group);
@@ -1895,6 +2010,8 @@ namespace Carpet {
return op_sync;
} else if (CCTK_Equals(prolong_string, "sync")) {
return op_sync;
+ } else if (CCTK_Equals(prolong_string, "restrict")) {
+ return op_restrict;
} else if (CCTK_Equals(prolong_string, "copy")) {
return op_copy;
} else if (CCTK_Equals(prolong_string, "Lagrange")) {
@@ -1903,6 +2020,8 @@ namespace Carpet {
return op_ENO;
} else if (CCTK_Equals(prolong_string, "WENO")) {
return op_WENO;
+ } else if (CCTK_Equals(prolong_string, "Lagrange_monotone")) {
+ return op_Lagrange_monotone;
} else {
char * const groupname = CCTK_GroupName (group);
CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
@@ -2170,7 +2289,7 @@ namespace Carpet {
DECLARE_CCTK_PARAMETERS;
int const prolongation_stencil_size
- = vdd.at(m)->prolongation_stencil_size();
+ = vdd.AT(m)->prolongation_stencil_size();
int const min_nghosts
= ((prolongation_stencil_size + refinement_factor - 1)
/ (refinement_factor - 1));
diff --git a/Carpet/Carpet/src/Shutdown.cc b/Carpet/Carpet/src/Shutdown.cc
index 4edf04374..6da0daf3e 100644
--- a/Carpet/Carpet/src/Shutdown.cc
+++ b/Carpet/Carpet/src/Shutdown.cc
@@ -1,13 +1,13 @@
#include <cstdio>
#include <cstdlib>
-#include "cctk.h"
-#include "cctk_Parameters.h"
+#include <cctk.h>
+#include <cctk_Parameters.h>
-#include "dist.hh"
+#include <dist.hh>
-#include "carpet.hh"
-#include "Timers.hh"
+#include <carpet.hh>
+#include <Timers.hh>
@@ -52,6 +52,15 @@ namespace Carpet {
} END_REVERSE_MGLEVEL_LOOP;
} // for rl
+ // Stop all timers before shutdown, since timers may rely on data
+ // structures which are destroyed during shutdown
+ int const ierr = CCTK_TimerStop ("CCTK total time");
+ assert (not ierr);
+ timer.stop();
+ if (output_timers_every > 0) {
+ TimerSet::writeData (cctkGH, timer_file);
+ }
+
BEGIN_REVERSE_MGLEVEL_LOOP(cctkGH) {
do_early_global_mode = true;
do_late_global_mode = true;
@@ -65,10 +74,6 @@ namespace Carpet {
CCTK_ScheduleTraverse ("CCTK_SHUTDOWN", cctkGH, CallFunction);
} END_REVERSE_MGLEVEL_LOOP;
- timer.stop();
- if (output_timers_every > 0) {
- TimerSet::writeData (cctkGH, timer_file);
- }
// earlier checkpoint before finalising MPI
Waypoint ("Done with shutdown");
diff --git a/Carpet/Carpet/src/Storage.cc b/Carpet/Carpet/src/Storage.cc
index e038a261f..636247d20 100644
--- a/Carpet/Carpet/src/Storage.cc
+++ b/Carpet/Carpet/src/Storage.cc
@@ -1,14 +1,15 @@
#include <cassert>
#include <cstdlib>
-#include "cctk.h"
-#include "cctk_Parameters.h"
+#include <cctk.h>
+#include <cctk_Parameters.h>
-#include "dh.hh"
-#include "gf.hh"
-#include "operators.hh"
+#include <defs.hh>
+#include <dh.hh>
+#include <gf.hh>
+#include <operators.hh>
-#include "carpet.hh"
+#include <carpet.hh>
@@ -94,16 +95,16 @@ namespace Carpet {
// Record previous number of allocated time levels
if (status) {
// Note: This remembers only the last level
- status[n] = groupdata.at(group).activetimelevels.at(ml).at(rl);
+ status[n] = groupdata.AT(group).activetimelevels.AT(ml).AT(rl);
}
// Only do something if the number of time levels actually
// needs to be changed -- do nothing otherwise
const bool do_increase
- = inc and timelevels[n] > groupdata.at(group).activetimelevels.at(ml).at(rl);
+ = inc and timelevels[n] > groupdata.AT(group).activetimelevels.AT(ml).AT(rl);
const bool do_decrease
- = not inc and timelevels[n] < groupdata.at(group).activetimelevels.at(ml).at(rl);
+ = not inc and timelevels[n] < groupdata.AT(group).activetimelevels.AT(ml).AT(rl);
if (do_increase or do_decrease) {
if (not can_do) {
@@ -127,13 +128,13 @@ namespace Carpet {
if (gp.grouptype == CCTK_GF) {
assert ((map == -1 or maps == 1)
and (component == -1
- or vhh.at(0)->local_components(reflevel) == 1));
+ or vhh.AT(0)->local_components(reflevel) == 1));
}
// Set the new number of active time levels
- groupdata.at(group).activetimelevels.at(ml).at(rl) = timelevels[n];
+ groupdata.AT(group).activetimelevels.AT(ml).AT(rl) = timelevels[n];
- for (int m=0; m<(int)arrdata.at(group).size(); ++m) {
+ for (int m=0; m<(int)arrdata.AT(group).size(); ++m) {
for (int var=0; var<gp.numvars; ++var) {
#ifdef CCTK_HAVE_CONTIGUOUS_GROUPS
bool const contiguous = gp.contiguous;
@@ -152,19 +153,19 @@ namespace Carpet {
assert (vectorlength>0 and vectorlength<=gp.numvars);
ggf* const vectorleader
= (vectorindex>0
- ? arrdata.at(group).at(m).data.at(var - vectorindex)
+ ? arrdata.AT(group).AT(m).data.AT(var - vectorindex)
: NULL);
const int varindex = firstvarindex + var;
#warning "TODO: allocate these in SetupGH, and after recomposing"
- if (not arrdata.at(group).at(m).data.at(var)) {
+ if (not arrdata.AT(group).AT(m).data.AT(var)) {
switch (gp.vartype) {
#define TYPECASE(N,T) \
case N: \
- arrdata.at(group).at(m).data.at(var) = new gf<T> \
+ arrdata.AT(group).AT(m).data.AT(var) = new gf<T> \
(varindex, \
- groupdata.at(group).transport_operator, \
- *arrdata.at(group).at(m).tt, \
- *arrdata.at(group).at(m).dd, \
+ groupdata.AT(group).transport_operator, \
+ *arrdata.AT(group).AT(m).tt, \
+ *arrdata.AT(group).AT(m).dd, \
prolongation_order_time, \
vectorlength, vectorindex, (gf<T>*)vectorleader); \
break;
@@ -175,18 +176,17 @@ namespace Carpet {
} // switch gp.vartype
} // if not allocated
- arrdata.at(group).at(m).data.at(var)->set_timelevels
+ arrdata.AT(group).AT(m).data.AT(var)->set_timelevels
(ml, rl, timelevels[n]);
// Set the data pointers for grid arrays
if (gp.grouptype != CCTK_GF) {
assert (rl==0 and m==0);
- int const c = CCTK_MyProc(cgh);
for (int tl=0; tl<gp.numtimelevels; ++tl) {
cgh->data[varindex][tl]
- = (tl < groupdata.at(group).activetimelevels.at(ml).at(rl)
- ? ((*arrdata.at(group).at(m).data.at(var))
- (tl, 0, c, 0)->storage())
+ = (tl < groupdata.AT(group).activetimelevels.AT(ml).AT(rl)
+ ? ((*arrdata.AT(group).AT(m).data.AT(var))
+ (tl, 0, 0, 0)->storage())
: NULL);
}
} // if grouptype != GF
@@ -202,7 +202,7 @@ namespace Carpet {
// Record current number of time levels
// Note: This adds the time levels of all refinement levels
total_num_timelevels
- += groupdata.at(group).activetimelevels.at(ml).at(rl);
+ += groupdata.AT(group).activetimelevels.AT(ml).AT(rl);
} // for rl
} // for ml
@@ -292,8 +292,9 @@ namespace Carpet {
if (groupdata.size() == 0) return -3;
int const rl = grouptype == CCTK_GF ? reflevel : 0;
// Return whether storage is allocated
- assert (groupdata.at(group).activetimelevels.at(mglevel).at(rl) != deadbeef);
- return groupdata.at(group).activetimelevels.at(mglevel).at(rl) > 0;
+ CCTK_INT const deadbeef = get_deadbeef();
+ assert (groupdata.AT(group).activetimelevels.AT(mglevel).AT(rl) != deadbeef);
+ return groupdata.AT(group).activetimelevels.AT(mglevel).AT(rl) > 0;
}
@@ -318,12 +319,12 @@ namespace Carpet {
return &error; // global or level mode for a GF
}
- const int gpdim = groupdata.at(group).info.dim;
+ const int gpdim = groupdata.AT(group).info.dim;
assert (dir>=0 and dir<gpdim);
if (CCTK_QueryGroupStorageI(cgh, group)) {
- return &groupdata.at(group).info.lsh[dir];
+ return &groupdata.AT(group).info.lsh[dir];
} else {
@@ -337,8 +338,19 @@ namespace Carpet {
int GroupDynamicData (const cGH* cgh, int group, cGroupDynamicData* data)
{
+ // Return values:
+ // 0 for success
+ // -1 if given pointer to data structure is NULL
+ // -3 if given GH pointer is invalid
+ // (-77 if group has zero variables)
+ // -78 if group does not exist
+ if (not cgh) return -3;
+ if (not (group>=0 and group<CCTK_NumGroups())) return -78;
+ if (not data) return -1;
+ assert (cgh);
assert (group>=0 and group<CCTK_NumGroups());
- *data = groupdata.at(group).info;
+ assert (data);
+ *data = groupdata.AT(group).info;
return 0;
}
@@ -358,12 +370,13 @@ namespace Carpet {
check (not CCTK_GroupData (group, & gp));
if (gp.grouptype == CCTK_GF) {
- if (groupdata.at(group).transport_operator != op_none and
- groupdata.at(group).transport_operator != op_sync and
- groupdata.at(group).transport_operator != op_copy)
+ if (groupdata.AT(group).transport_operator != op_none and
+ groupdata.AT(group).transport_operator != op_sync and
+ groupdata.AT(group).transport_operator != op_restrict and
+ groupdata.AT(group).transport_operator != op_copy)
{
- if (groupdata.at(group).activetimelevels.at(ml).at(rl) != 0 and
- (groupdata.at(group).activetimelevels.at(ml).at(rl) <
+ if (groupdata.AT(group).activetimelevels.AT(ml).AT(rl) != 0 and
+ (groupdata.AT(group).activetimelevels.AT(ml).AT(rl) <
prolongation_order_time+1))
{
static vector<bool> didwarn;
@@ -371,9 +384,9 @@ namespace Carpet {
if ((int)didwarn.size() < numgroups) {
didwarn.resize (numgroups, false);
}
- if (not didwarn.at(group)) {
+ if (not didwarn.AT(group)) {
// Warn only once per group
- didwarn.at(group) = true;
+ didwarn.AT(group) = true;
char * const groupname = CCTK_GroupName (group);
CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
"There are not enough time levels for the desired temporal prolongation order in the grid function group \"%s\". With Carpet::prolongation_order_time=%d, you need at least %d time levels.",
diff --git a/Carpet/Carpet/src/Timers.cc b/Carpet/Carpet/src/Timers.cc
index 265c2cd4c..c279bf0cc 100644
--- a/Carpet/Carpet/src/Timers.cc
+++ b/Carpet/Carpet/src/Timers.cc
@@ -3,10 +3,9 @@
#include <cstring>
#include <list>
-#include "cctk.h"
-#include "cctk_Parameters.h"
-
-#include "util_String.h"
+#include <cctk.h>
+#include <cctk_Parameters.h>
+#include <util_String.h>
#if HAVE_UNISTD_H
# include <fcntl.h>
@@ -15,7 +14,7 @@
#include <defs.hh>
-#include "Timers.hh"
+#include <Timers.hh>
@@ -236,6 +235,8 @@ namespace Carpet {
case val_double:
printf (" %g", timer->vals[i].val.d);
break;
+ case val_none:
+ break;
default:
assert (0);
}
@@ -245,4 +246,27 @@ namespace Carpet {
if (was_running) start();
}
+
+
+ // Output (debug) messages that a timer is starting or stopping
+ void
+ Timer::msgStart ()
+ const
+ {
+ DECLARE_CCTK_PARAMETERS;
+ if (timers_verbose) {
+ CCTK_VInfo (CCTK_THORNSTRING, "Timer \"%s\" starting", name());
+ }
+ }
+
+ void
+ Timer::msgStop ()
+ const
+ {
+ DECLARE_CCTK_PARAMETERS;
+ if (timers_verbose) {
+ CCTK_VInfo (CCTK_THORNSTRING, "Timer \"%s\" stopping", name());
+ }
+ }
+
} // namespace Carpet
diff --git a/Carpet/Carpet/src/Timers.hh b/Carpet/Carpet/src/Timers.hh
index 2c2543168..af887cb69 100644
--- a/Carpet/Carpet/src/Timers.hh
+++ b/Carpet/Carpet/src/Timers.hh
@@ -1,6 +1,7 @@
+#include <iostream>
#include <list>
-#include "cctk.h"
+#include <cctk.h>
@@ -79,6 +80,7 @@ namespace Carpet {
void
start ()
{
+ msgStart ();
running = true;
CCTK_TimerStartI (handle);
}
@@ -89,6 +91,7 @@ namespace Carpet {
{
CCTK_TimerStopI (handle);
running = false;
+ msgStop ();
}
// Reset the timer
@@ -107,6 +110,32 @@ namespace Carpet {
void
printData ();
+ private:
+
+ // Output (debug) messages that a timer is starting or stopping
+ void
+ msgStart ()
+ const;
+
+ void
+ msgStop ()
+ const;
+
};
+
+
+ // Macros for using timers in a convenient manner
+
+#define TIMING_BEGIN(name) \
+ { \
+ static Carpet::Timer timer (name); \
+ timer.start(); \
+ {
+
+#define TIMING_END \
+ } \
+ timer.stop(); \
+} while (0)
+
} // namespace Carpet
diff --git a/Carpet/Carpet/src/Timing.cc b/Carpet/Carpet/src/Timing.cc
index 6c4102733..2fd24f3bd 100644
--- a/Carpet/Carpet/src/Timing.cc
+++ b/Carpet/Carpet/src/Timing.cc
@@ -3,9 +3,9 @@
#include <cmath>
#include <cstdlib>
-#include "cctk.h"
-#include "cctk_Arguments.h"
-#include "cctk_Parameters.h"
+#include <cctk.h>
+#include <cctk_Arguments.h>
+#include <cctk_Parameters.h>
// IRIX wants this before <time.h>
#if HAVE_SYS_TYPES_H
@@ -27,10 +27,10 @@
# include <unistd.h>
#endif
-#include "defs.hh"
-#include "dist.hh"
+#include <defs.hh>
+#include <dist.hh>
-#include "carpet.hh"
+#include <carpet.hh>
@@ -77,17 +77,17 @@ namespace Carpet {
for (int m = 0; m < maps; ++ m) {
assert (reflevel >= 0);
int const rl = reflevel;
- for (int c = 0; c < vhh.at(m)->components(rl); ++ c) {
+ for (int c = 0; c < vhh.AT(m)->components(rl); ++ c) {
assert (mglevel >= 0);
int const ml = mglevel;
// Base region
- ibbox const ext = vhh.at(m)->extent(ml,rl,c);
+ ibbox const ext = vhh.AT(m)->extent(ml,rl,c);
// Count the grid points
CCTK_REAL const domainsize = ext.size();
- if (vhh.at(m)->is_local (rl, c)) {
+ if (vhh.AT(m)->is_local (rl, c)) {
local_num_grid_points += domainsize;
}
global_num_grid_points += domainsize;
@@ -117,7 +117,11 @@ namespace Carpet {
+ // Time at which the simulation started
+ CCTK_REAL startup_walltime; // in seconds
+
// Time at which the evolution started
+ bool in_evolution = false;
CCTK_REAL initial_walltime; // in seconds
CCTK_REAL initial_phystime;
@@ -136,13 +140,19 @@ namespace Carpet {
{
DECLARE_CCTK_ARGUMENTS;
- * physical_time_per_hour = 0.0;
+ startup_walltime = get_walltime();
+
+ * physical_time_per_hour = 0.0;
+ * current_physical_time_per_hour = 0.0;
* time_total = 0.0;
+ * time_evolution = 0.0;
* time_computing = 0.0;
* time_communicating = 0.0;
* time_io = 0.0;
+ * evolution_steps_count = 0.0;
+
* local_grid_points_per_second = 0.0;
* total_grid_points_per_second = 0.0;
* local_grid_point_updates_count = 0.0;
@@ -171,10 +181,11 @@ namespace Carpet {
// Begin timing (to be called after initialisation, just before the
// main evolution begins)
void
- BeginTiming (cGH const * const cctkGH)
+ BeginTimingEvolution (cGH const * const cctkGH)
{
DECLARE_CCTK_ARGUMENTS;
+ in_evolution = true;
initial_walltime = get_walltime();
initial_phystime = cctkGH->cctk_time;
}
@@ -184,15 +195,18 @@ namespace Carpet {
// Take a step on the current refinement and multigrid level (to be
// called when EVOL is scheduled)
void
- StepTiming (cGH const * const cctkGH)
+ StepTimingEvolution (cGH const * const cctkGH)
{
DECLARE_CCTK_ARGUMENTS;
+ assert (in_evolution);
assert (timing_state == state_computing);
CCTK_REAL local_updates, global_updates;
current_level_updates (cctkGH, local_updates, global_updates);
+ ++ * evolution_steps_count;
+
* local_grid_point_updates_count += local_updates;
* total_grid_point_updates_count += global_updates;
@@ -266,7 +280,9 @@ namespace Carpet {
assert (timing_state == state_computing);
// Measure the elapsed time
- * time_total = get_walltime() - initial_walltime;
+ double const walltime = get_walltime();
+ * time_total = walltime - startup_walltime;
+ * time_evolution = in_evolution ? walltime - initial_walltime : 0.0;
* time_computing = * time_total - (* time_communicating + * time_io);
}
@@ -327,12 +343,52 @@ namespace Carpet {
UpdatePhysicalTimePerHour (cGH const * const cctkGH)
{
DECLARE_CCTK_ARGUMENTS;
+ DECLARE_CCTK_PARAMETERS;
+
+ if (not in_evolution) {
+ * physical_time_per_hour = 0.0;
+ * current_physical_time_per_hour = 0.0;
+ return;
+ }
+
+ static int last_iteration = -1;
+ static size_t num_samples = 0;
+ static CCTK_REAL last_physical_time;
+ static CCTK_REAL last_time_evolution;
+ assert (cctk_iteration > last_iteration); // expect progress
// Calculate elapsed physical time
CCTK_REAL const physical_time = cctkGH->cctk_time - initial_phystime;
- // Calculate physical time per hour
- * physical_time_per_hour = 3600.0 * physical_time / max (* time_total, eps);
+ // Calculate average physical time per hour
+ * physical_time_per_hour =
+ 3600.0 * physical_time / max (* time_evolution, eps);
+
+ // Calculate current physical time per hour as moving average
+ if (last_iteration < 0) {
+ // No past data are available
+ * current_physical_time_per_hour = * physical_time_per_hour;
+ } else if (num_samples < 3 or
+ * time_evolution < 0.01 * timing_average_window_minutes * 60.0)
+ {
+ // Less than three previous samples are available, or less thatn
+ // one percent of a window of past data are available
+ * current_physical_time_per_hour = * physical_time_per_hour;
+ } else {
+ CCTK_REAL const window =
+ min (* time_evolution, timing_average_window_minutes * 60.0);
+ CCTK_REAL const alpha =
+ exp (- (* time_evolution - last_time_evolution) / window);
+ * current_physical_time_per_hour =
+ (1.0 - alpha) * * physical_time_per_hour +
+ ( alpha) * * current_physical_time_per_hour;
+ }
+
+ // Remember last iteration
+ last_iteration = cctk_iteration;
+ ++ num_samples;
+ last_physical_time = physical_time;
+ last_time_evolution = * time_evolution;
}
@@ -357,7 +413,13 @@ namespace Carpet {
DECLARE_CCTK_ARGUMENTS;
CCTK_VInfo (CCTK_THORNSTRING,
- "Total run time: %g h", double (* time_total / 3600.0));
+ "Total run time: %g h (iteration %d)",
+ double (* time_total / 3600.0),
+ cctk_iteration);
+ CCTK_VInfo (CCTK_THORNSTRING,
+ "Total evolution time: %g h (%d steps)",
+ double (* time_evolution / 3600.0),
+ int (* evolution_steps_count));
CCTK_VInfo (CCTK_THORNSTRING,
"(Comp, Comm, I/O) fractions = (%3.1f%%, %3.1f%%, %3.1f%%)",
double (100.0 * * time_computing /
@@ -461,7 +523,10 @@ namespace Carpet {
DECLARE_CCTK_ARGUMENTS;
CCTK_VInfo (CCTK_THORNSTRING,
- "Physical time per hour: %g",
+ "Current physical time per hour: %g",
+ double (* current_physical_time_per_hour));
+ CCTK_VInfo (CCTK_THORNSTRING,
+ "Average physical time per hour: %g",
double (* physical_time_per_hour));
}
@@ -479,25 +544,25 @@ namespace Carpet {
<< "Memory statistics:" << eol
<< " Grid hierarchy:" << eol;
for (int m = 0; m < Carpet::maps; ++ m) {
- cout << " gh[" << m << "]: " << PRINTMEM(*vhh.at(m)) << eol
- << " dh[" << m << "]: " << PRINTMEM(*vdd.at(m)) << eol
- << " th[" << m << "]: " << PRINTMEM(*vtt.at(m)) << eol;
+ cout << " gh[" << m << "]: " << PRINTMEM(*vhh.AT(m)) << eol
+ << " dh[" << m << "]: " << PRINTMEM(*vdd.AT(m)) << eol
+ << " th[" << m << "]: " << PRINTMEM(*vtt.AT(m)) << eol;
}
#if 0
for (int g = 0; g < (int)arrdata.size(); ++ g) {
if (CCTK_GroupTypeI(g) != CCTK_GF) {
char * const groupname = CCTK_GroupName(g);
- for (int m = 0; m < (int)arrdata.at(g).size(); ++ m) {
+ for (int m = 0; m < (int)arrdata.AT(g).size(); ++ m) {
cout << " Group " << groupname << ":" << eol
- << " gh[" << m << "]: " << PRINTMEM(*arrdata.at(g).at(m).hh) << eol
- << " dh[" << m << "]: " << PRINTMEM(*arrdata.at(g).at(m).dd) << eol
- << " th[" << m << "]: " << PRINTMEM(*arrdata.at(g).at(m).tt) << eol;
- for (int v = 0; v < (int)arrdata.at(g).at(m).data.size(); ++ v) {
+ << " gh[" << m << "]: " << PRINTMEM(*arrdata.AT(g).AT(m).hh) << eol
+ << " dh[" << m << "]: " << PRINTMEM(*arrdata.AT(g).AT(m).dd) << eol
+ << " th[" << m << "]: " << PRINTMEM(*arrdata.AT(g).AT(m).tt) << eol;
+ for (int v = 0; v < (int)arrdata.AT(g).AT(m).data.size(); ++ v) {
char * const fullname = CCTK_FullName(CCTK_FirstVarIndexI(g)+v);
cout << " Variable " << fullname << ":" << eol
<< " ggf[" << m << "]: ";
- if (arrdata.at(g).at(m).data.at(v)) {
- cout << PRINTMEM(*arrdata.at(g).at(m).data.at(v));
+ if (arrdata.AT(g).AT(m).data.AT(v)) {
+ cout << PRINTMEM(*arrdata.AT(g).AT(m).data.AT(v));
} else {
cout << "<null>";
}
diff --git a/Carpet/Carpet/src/carpet.hh b/Carpet/Carpet/src/carpet.hh
index 1a6ee0f3b..bbcb4296a 100644
--- a/Carpet/Carpet/src/carpet.hh
+++ b/Carpet/Carpet/src/carpet.hh
@@ -3,11 +3,11 @@
#include <vector>
-#include "cctk.h"
-#include "cctk_Arguments.h"
-#include "cctk_Functions.h"
+#include <cctk.h>
+#include <cctk_Arguments.h>
+#include <cctk_Functions.h>
-#include "gh.hh"
+#include <gh.hh>
#include "carpet_public.hh"
@@ -36,7 +36,7 @@ namespace Carpet {
int CallFunction (void* function, cFunctionData* attribute, void* data);
// Other functions
- bool Regrid (cGH const * cctkGH, bool force_recompose);
+ bool Regrid (cGH const * cctkGH, bool force_recompose, bool do_init);
void CycleTimeLevels (const cGH* cgh);
void FlipTimeLevels (const cGH* cgh);
@@ -52,8 +52,8 @@ namespace Carpet {
allbutcurrenttime,
alltimes };
- int min_timelevel (checktimes where, int num_tl);
- int max_timelevel (checktimes where, int num_tl);
+ int min_timelevel (checktimes where, int num_tl, bool persistent) CCTK_ATTRIBUTE_CONST;
+ int max_timelevel (checktimes where, int num_tl, bool persistent) CCTK_ATTRIBUTE_CONST;
void Poison (const cGH* cgh, checktimes where, int what = 0);
void PoisonGroup (const cGH* cgh, int group, checktimes where);
diff --git a/Carpet/Carpet/src/carpet_public.h b/Carpet/Carpet/src/carpet_public.h
index ce2faaad5..20a265848 100644
--- a/Carpet/Carpet/src/carpet_public.h
+++ b/Carpet/Carpet/src/carpet_public.h
@@ -3,7 +3,7 @@
#include <mpi.h>
-#include "cctk.h"
+#include <cctk.h>
@@ -76,7 +76,7 @@ namespace Carpet {
};
- struct CarpetGH const * GetCarpetGH (const cGH * const cgh);
+ struct CarpetGH const * GetCarpetGH (const cGH * const cgh) CCTK_ATTRIBUTE_CONST;
@@ -86,6 +86,16 @@ namespace Carpet {
+ /* Grid function access */
+ CCTK_POINTER Carpet_VarDataPtrI (CCTK_POINTER_TO_CONST const cctkGH,
+ CCTK_INT const m,
+ CCTK_INT const rl,
+ CCTK_INT const c,
+ CCTK_INT const tl,
+ CCTK_INT const varindex);
+
+
+
/* Call a schedule group */
int CallScheduleGroup (cGH * const cgh, const char * const group);
@@ -104,10 +114,10 @@ namespace Carpet {
/* Helper functions */
- MPI_Comm CarpetMPIComm (void);
- MPI_Datatype CarpetMPIDatatype (int vartype);
- MPI_Datatype CarpetSimpleMPIDatatype (int vartype);
- int CarpetSimpleMPIDatatypeLength (int vartype);
+ MPI_Comm CarpetMPIComm (void) CCTK_ATTRIBUTE_CONST;
+ MPI_Datatype CarpetMPIDatatype (int vartype) CCTK_ATTRIBUTE_CONST;
+ MPI_Datatype CarpetSimpleMPIDatatype (int vartype) CCTK_ATTRIBUTE_CONST;
+ int CarpetSimpleMPIDatatypeLength (int vartype) CCTK_ATTRIBUTE_CONST;
diff --git a/Carpet/Carpet/src/functions.hh b/Carpet/Carpet/src/functions.hh
index d8b508441..e1fd104b8 100644
--- a/Carpet/Carpet/src/functions.hh
+++ b/Carpet/Carpet/src/functions.hh
@@ -7,13 +7,13 @@
#include <mpi.h>
-#include "cctk.h"
-#include "cctk_Schedule.h"
+#include <cctk.h>
+#include <cctk_Schedule.h>
-#include "bbox.hh"
-#include "dh.hh"
-#include "gh.hh"
-#include "vect.hh"
+#include <bbox.hh>
+#include <dh.hh>
+#include <gh.hh>
+#include <vect.hh>
@@ -33,24 +33,20 @@ namespace Carpet {
int GroupStorageDecrease (const cGH* cgh, int n_groups, const int* groups,
const int* timelevels, int* status);
int Barrier (const cGH* cgh);
- int Exit (cGH* cgh, int retval);
- int Abort (cGH* cgh, int retval);
- int MyProc (const cGH* cgh);
- int nProcs (const cGH* cgh);
+ int NamedBarrier (const cGH* cgh, int id);
+ int Exit (const cGH* cgh, int retval);
+ int Abort (const cGH* cgh, int retval);
+ int MyProc (const cGH* cgh) CCTK_ATTRIBUTE_CONST;
+ int nProcs (const cGH* cgh) CCTK_ATTRIBUTE_CONST;
const int* ArrayGroupSizeB (const cGH* cgh, int dir, int group,
- const char* groupname);
- int QueryGroupStorageB (const cGH* cgh, int group, const char* groupname);
+ const char* groupname) CCTK_ATTRIBUTE_PURE;
+ int QueryGroupStorageB (const cGH* cgh, int group, const char* groupname) CCTK_ATTRIBUTE_PURE;
int GroupDynamicData (const cGH* cgh, int group, cGroupDynamicData* data);
void Restrict (const cGH* cgh);
- // Strings
- vector <string>
- AllGatherString (MPI_Comm const world,
- string const & data);
-
// Multi-Model
void
SplitUniverse (MPI_Comm const world, string const model, MPI_Comm & comm,
@@ -58,19 +54,19 @@ namespace Carpet {
// Model id to model name
vector <string> Models ();
- string Model (int id);
+ string Model (int id) CCTK_ATTRIBUTE_PURE;
// Model name to model id
std::map <string, int> ModelMap ();
- int ModelMap (string name);
+ int ModelMap (string name) CCTK_ATTRIBUTE_PURE;
// Processor to model id
vector <int> ModelIds ();
- int ModelId (int proc);
+ int ModelId (int proc) CCTK_ATTRIBUTE_PURE;
// Model id to processors
- vector <vector <int> > ModelProcs ();
- vector <int> ModelProcs (int proc);
+ vector <vector <int> > ModelProcs () CCTK_ATTRIBUTE_PURE;
+ vector <int> ModelProcs (int proc) CCTK_ATTRIBUTE_PURE;
extern "C" {
CCTK_POINTER_TO_CONST
@@ -90,10 +86,6 @@ namespace Carpet {
- void SetSystemLimits ();
-
-
-
// Helpers for storage
void GroupsStorageCheck (cGH const * const cctkGH);
@@ -102,18 +94,27 @@ namespace Carpet {
RegridMap (cGH const * cctkGH,
int m,
gh::rregs const & supeerregss,
- gh::mregs const & regsss);
+ gh::mregs const & regsss,
+ bool do_init);
void
PostRegrid (cGH const * cctkGH);
bool
Recompose (cGH const * cctkGH,
int rl,
bool do_init);
+ void
+ RegridFree (cGH const * cctkGH,
+ bool do_init);
void
CheckRegions (gh::mregs const & regsss);
void
+ OutputSuperregions (cGH const * cctkGH,
+ int m,
+ gh const & hh,
+ gh::rregs const & superregss);
+ void
OutputGrids (cGH const * cctkGH,
int const m,
gh const & hh,
@@ -177,8 +178,8 @@ namespace Carpet {
// Timing statistics functions
void InitTimingStats (cGH const * cctkGH);
- void BeginTiming (cGH const * cctkGH);
- void StepTiming (cGH const * cctkGH);
+ void BeginTimingEvolution (cGH const * cctkGH);
+ void StepTimingEvolution (cGH const * cctkGH);
void BeginTimingIO (cGH const * cctkGH);
void EndTimingIO (cGH const * cctkGH,
CCTK_REAL files, CCTK_REAL bytes, bool is_binary);
diff --git a/Carpet/Carpet/src/helpers.cc b/Carpet/Carpet/src/helpers.cc
index c799e07f4..d5b49fc61 100644
--- a/Carpet/Carpet/src/helpers.cc
+++ b/Carpet/Carpet/src/helpers.cc
@@ -6,15 +6,15 @@
#include <mpi.h>
-#include "cctk.h"
-#include "cctk_FortranString.h"
-#include "cctk_Parameters.h"
+#include <cctk.h>
+#include <cctk_FortranString.h>
+#include <cctk_Parameters.h>
-#include "defs.hh"
-#include "dist.hh"
-#include "ggf.hh"
+#include <defs.hh>
+#include <dist.hh>
+#include <ggf.hh>
-#include "carpet.hh"
+#include <carpet.hh>
@@ -64,9 +64,51 @@ namespace Carpet {
{
return do_prolongate;
}
-
-
-
+
+
+
+ // Get pointer to grid variable for a specific map and refinement level
+ CCTK_POINTER
+ Carpet_VarDataPtrI (CCTK_POINTER_TO_CONST const cctkGH,
+ CCTK_INT const m,
+ CCTK_INT const rl,
+ CCTK_INT const c,
+ CCTK_INT const tl,
+ CCTK_INT const varindex)
+ {
+ assert (cctkGH);
+ assert (varindex >= 0 and varindex < CCTK_NumVars());
+ int const groupindex = CCTK_GroupIndexFromVarI (varindex);
+ assert (groupindex >= 0);
+ int const grouptype = CCTK_GroupTypeI (groupindex);
+ assert (mglevel >= 0);
+ if (grouptype == CCTK_GF) {
+ assert (m >= 0 and m < maps);
+ assert (rl >= 0 and rl < reflevels);
+ assert (c >= 0 and c < arrdata.AT(groupindex).AT(m).hh->components(reflevel));
+ assert (arrdata.AT(groupindex).AT(m).hh->is_local(rl, c));
+ } else {
+ assert (m == 0);
+ assert (rl == 0);
+ assert (c == CCTK_MyProc (NULL));
+ }
+ int const maxtimelevels = CCTK_MaxTimeLevelsGI (groupindex);
+ assert (tl >= 0 and tl < maxtimelevels);
+
+ int const activetimelevels =
+ groupdata.AT(groupindex).activetimelevels.AT(mglevel).AT(rl);
+ if (tl < activetimelevels) {
+ int const var = varindex - CCTK_FirstVarIndexI (groupindex);
+ ggf * const ff = arrdata.AT(groupindex).AT(m).data.AT(var);
+ gdata * const data = (*ff) (tl, rl, c, mglevel);
+ return data->storage();
+ } else {
+ return NULL;
+ }
+ }
+
+
+
// Multi-Model
CCTK_POINTER_TO_CONST
Carpet_GetMPICommUniverse (CCTK_POINTER_TO_CONST const cctkGH)
@@ -106,13 +148,13 @@ namespace Carpet {
assert (size == dim);
- ibbox const & baseext = vhh.at(m)->baseextents.at(ml).at(0);
+ ibbox const & baseext = vhh.AT(m)->baseextents.AT(ml).AT(0);
ivect const igsh = baseext.shape() / baseext.stride();
for (int d = 0; d < dim; ++ d) {
gsh[d] = igsh[d];
- lower[d] = origin_space.at(m).at(ml)[d];
- delta[d] = delta_space.at(m)[d] * mglevelfact;
+ lower[d] = origin_space.AT(m).AT(ml)[d];
+ delta[d] = delta_space.AT(m)[d] * mglevelfact;
upper[d] = lower[d] + delta[d] * (gsh[d] - 1);
}
@@ -125,16 +167,37 @@ namespace Carpet {
int Barrier (const cGH* cgh)
{
- void *dummy = &dummy;
+ const void *dummy = &dummy;
dummy = &cgh;
MPI_Barrier (dist::comm());
return 0;
}
+ int NamedBarrier (const cGH* const cgh, const int id)
+ {
+ const void *dummy = &dummy;
+ dummy = &cgh;
+
+ assert (id >= 0);
+ const int root = 0;
+ int my_id = dist::rank()==root ? id : -1;
+ Checkpoint ("About to Bcast %d", id);
+ MPI_Bcast (&my_id, 1, MPI_INT, root, dist::comm());
+ Checkpoint ("Finished Bcast");
+ if (my_id != id) {
+ CCTK_VWarn (CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "Wrong Barrier name: expected %d, found %d", id, my_id);
+ }
+ Checkpoint ("About to Barrier %d", id);
+ MPI_Barrier (dist::comm());
+ Checkpoint ("Finished Barrier");
+ return 0;
+ }
+
- int Exit (cGH* cgh, int retval)
+ int Exit (const cGH* cgh, int retval)
{
CCTK_Barrier (cgh);
dist::finalize();
@@ -142,7 +205,7 @@ namespace Carpet {
return -1;
}
- int Abort (cGH* cgh, int retval)
+ int Abort (const cGH* cgh, int retval)
{
void *dummy = &dummy;
dummy = &cgh;
@@ -189,7 +252,7 @@ namespace Carpet {
#define TYPECASE(N,T) \
case N: { \
T dummy; \
- return dist::datatype(dummy); \
+ return dist::mpi_datatype(dummy); \
}
#include "typecase"
#undef TYPECASE
@@ -266,15 +329,17 @@ namespace Carpet {
// Timelevels
- int min_timelevel (const checktimes where, const int num_tl)
+ int min_timelevel (const checktimes where, const int num_tl,
+ const bool persistent)
{
assert (num_tl>0);
switch (where) {
case currenttime:
return 0;
case currenttimebutnotifonly:
- // don't include current time if there is only one time level
- return num_tl>1 ? 0 : 1;
+ // don't include current time if there is only one (persistent)
+ // time level
+ return (not persistent or num_tl>1) ? 0 : 1;
case previoustime:
return 1;
case allbutlasttime:
@@ -290,7 +355,8 @@ namespace Carpet {
return -999;
}
- int max_timelevel (const checktimes where, const int num_tl)
+ int max_timelevel (const checktimes where, const int num_tl,
+ const bool persistent)
{
assert (num_tl>0);
switch (where) {
diff --git a/Carpet/Carpet/src/make.code.defn b/Carpet/Carpet/src/make.code.defn
index db03ff6d7..8edcaac69 100644
--- a/Carpet/Carpet/src/make.code.defn
+++ b/Carpet/Carpet/src/make.code.defn
@@ -1,8 +1,7 @@
# Main make.code.defn file for thorn Carpet -*-Makefile-*-
# Source files in this directory
-SRCS = AllGatherString.cc \
- CallFunction.cc \
+SRCS = CallFunction.cc \
CarpetBasegrid.cc \
CarpetParamCheck.cc \
CarpetStartup.cc \
@@ -11,7 +10,6 @@ SRCS = AllGatherString.cc \
Cycle.cc \
Evolve.cc \
Initialise.cc \
- Limits.cc \
MultiModel.cc \
OutputGH.cc \
Poison.cc \
diff --git a/Carpet/Carpet/src/modes.cc b/Carpet/Carpet/src/modes.cc
index 998f593ef..138842da4 100644
--- a/Carpet/Carpet/src/modes.cc
+++ b/Carpet/Carpet/src/modes.cc
@@ -11,7 +11,7 @@
#include <defs.hh>
#include <ggf.hh>
-#include "carpet.hh"
+#include <carpet.hh>
@@ -54,7 +54,7 @@ namespace Carpet {
// assert (mglevel>=0 and mglevel<mglevels);
// assert (reflevel>=0 and reflevel<reflevels);
// assert (map>=0 and map<maps);
-// assert (vhh.at(map)->local_components(reflevel)==1 or component==-1);
+// assert (vhh.AT(map)->local_components(reflevel)==1 or component==-1);
}
@@ -76,13 +76,16 @@ namespace Carpet {
// TODO: this could also just be "mglevel" instead
cctkGH->cctk_convlevel = basemglevel + mglevel;
+ // Set mode
+ cctkGH->cctk_mode = CCTK_MODE_GLOBAL;
+
// Set time delta
cctkGH->cctk_delta_time = delta_time * mglevelfact;
if (maps == 1) {
// Set space delta
for (int d=0; d<dim; ++d) {
- cctkGH->cctk_origin_space[d] = origin_space.at(0).at(mglevel)[d];
- cctkGH->cctk_delta_space[d] = delta_space.at(0)[d] * mglevelfact;
+ cctkGH->cctk_origin_space[d] = origin_space.AT(0).AT(mglevel)[d];
+ cctkGH->cctk_delta_space[d] = delta_space.AT(0)[d] * mglevelfact;
}
}
@@ -94,48 +97,48 @@ namespace Carpet {
const int rl = 0;
const int m = 0;
- const int c = CCTK_MyProc(cctkGH);
+ const int c = dist::rank();
- const ibbox& baseext = arrdata.at(group).at(m).hh->baseextents.at(ml).at(rl);
- const ibbox& ext = arrdata.at(group).at(m).dd->boxes.at(ml).at(rl).at(c).exterior;
- const b2vect& obnds = arrdata.at(group).at(m).hh->outer_boundaries(rl,c);
+ const ibbox& baseext = arrdata.AT(group).AT(m).hh->baseextents.AT(ml).AT(rl);
+ const ibbox& ext = arrdata.AT(group).AT(m).dd->boxes.AT(ml).AT(rl).AT(c).exterior;
+ const b2vect& obnds = arrdata.AT(group).AT(m).hh->outer_boundaries(rl,c);
- ivect::ref(const_cast<int*>(groupdata.at(group).info.nghostzones))
- = arrdata.at(group).at(m).dd->ghost_width[0];
- ivect::ref(const_cast<int*>(groupdata.at(group).info.gsh))
+ ivect::ref(const_cast<int*>(groupdata.AT(group).info.nghostzones))
+ = arrdata.AT(group).AT(m).dd->ghost_width[0];
+ ivect::ref(const_cast<int*>(groupdata.AT(group).info.gsh))
= baseext.shape() / baseext.stride();
- ivect::ref(const_cast<int*>(groupdata.at(group).info.lsh))
+ ivect::ref(const_cast<int*>(groupdata.AT(group).info.lsh))
= ext.shape() / ext.stride();
- ivect::ref(const_cast<int*>(groupdata.at(group).info.lbnd))
+ ivect::ref(const_cast<int*>(groupdata.AT(group).info.lbnd))
= (ext.lower() - baseext.lower()) / ext.stride();
- ivect::ref(const_cast<int*>(groupdata.at(group).info.ubnd))
+ ivect::ref(const_cast<int*>(groupdata.AT(group).info.ubnd))
= (ext.upper() - baseext.lower()) / ext.stride();
if (gp.disttype == CCTK_DISTRIB_CONSTANT) {
int const d = gp.dim==0 ? 0 : gp.dim-1;
- ivect & gsh = ivect::ref(const_cast<int*>(groupdata.at(group).info.gsh));
- ivect & lsh = ivect::ref(const_cast<int*>(groupdata.at(group).info.lsh));
- ivect & lbnd = ivect::ref(const_cast<int*>(groupdata.at(group).info.lbnd));
- ivect & ubnd = ivect::ref(const_cast<int*>(groupdata.at(group).info.ubnd));
+ ivect & gsh = ivect::ref(const_cast<int*>(groupdata.AT(group).info.gsh));
+ ivect & lsh = ivect::ref(const_cast<int*>(groupdata.AT(group).info.lsh));
+ ivect & lbnd = ivect::ref(const_cast<int*>(groupdata.AT(group).info.lbnd));
+ ivect & ubnd = ivect::ref(const_cast<int*>(groupdata.AT(group).info.ubnd));
gsh[d] = lsh[d];
lbnd[d] = 0;
ubnd[d] = lsh[d] - 1;
}
for (int d=0; d<dim; ++d) {
- const_cast<int*>(groupdata.at(group).info.bbox)[2*d ] = obnds[0][d];
- const_cast<int*>(groupdata.at(group).info.bbox)[2*d+1] = obnds[1][d];
+ const_cast<int*>(groupdata.AT(group).info.bbox)[2*d ] = obnds[0][d];
+ const_cast<int*>(groupdata.AT(group).info.bbox)[2*d+1] = obnds[1][d];
}
- groupdata.at(group).info.activetimelevels
- = groupdata.at(group).activetimelevels.at(mglevel).at(0);
+ groupdata.AT(group).info.activetimelevels
+ = groupdata.AT(group).activetimelevels.AT(mglevel).AT(0);
for (int d=0; d<dim; ++d) {
- assert (groupdata.at(group).info.lsh[d]>=0);
- assert (groupdata.at(group).info.lsh[d]<=groupdata.at(group).info.gsh[d]);
- assert (groupdata.at(group).info.lbnd[d]>=0);
- assert (groupdata.at(group).info.lbnd[d]<=groupdata.at(group).info.ubnd[d]+1);
- assert (groupdata.at(group).info.ubnd[d]<groupdata.at(group).info.gsh[d]);
- assert (groupdata.at(group).info.lbnd[d] + groupdata.at(group).info.lsh[d] - 1
- == groupdata.at(group).info.ubnd[d]);
- assert (groupdata.at(group).info.lbnd[d]<=groupdata.at(group).info.ubnd[d]+1);
+ assert (groupdata.AT(group).info.lsh[d]>=0);
+ assert (groupdata.AT(group).info.lsh[d]<=groupdata.AT(group).info.gsh[d]);
+ assert (groupdata.AT(group).info.lbnd[d]>=0);
+ assert (groupdata.AT(group).info.lbnd[d]<=groupdata.AT(group).info.ubnd[d]+1);
+ assert (groupdata.AT(group).info.ubnd[d]<groupdata.AT(group).info.gsh[d]);
+ assert (groupdata.AT(group).info.lbnd[d] + groupdata.AT(group).info.lsh[d] - 1
+ == groupdata.AT(group).info.ubnd[d]);
+ assert (groupdata.AT(group).info.lbnd[d]<=groupdata.AT(group).info.ubnd[d]+1);
}
const int numvars = CCTK_NumVarsInGroupI (group);
@@ -144,17 +147,18 @@ namespace Carpet {
assert (firstvar>=0);
const int max_tl = CCTK_MaxTimeLevelsGI (group);
assert (max_tl>=0);
- const int active_tl = groupdata.at(group).info.activetimelevels;
+ const int active_tl = groupdata.AT(group).info.activetimelevels;
assert (active_tl>=0 and active_tl<=max_tl);
- assert (arrdata.at(group).at(m).hh->is_local(rl,c));
+ assert (arrdata.AT(group).AT(m).hh->is_local(rl,c));
for (int var=0; var<numvars; ++var) {
assert (firstvar+var<CCTK_NumVars());
- ggf * const ff = arrdata.at(group).at(m).data.at(var);
+ ggf * const ff = arrdata.AT(group).AT(m).data.AT(var);
for (int tl=0; tl<max_tl; ++tl) {
if (ff and tl<active_tl) {
- gdata * const data = (*ff) (tl, rl, c, ml);
+ int const lc = 0;
+ gdata * const data = (*ff) (tl, rl, lc, ml);
assert (data);
cctkGH->data[firstvar+var][tl] = data->storage();
} else {
@@ -186,36 +190,38 @@ namespace Carpet {
if (maps == 1) {
// Save and unset space delta
for (int d=0; d<dim; ++d) {
- origin_space.at(0).at(mglevel)[d] = cctkGH->cctk_origin_space[d];
- delta_space.at(mglevel)[d] = cctkGH->cctk_delta_space[d] / mglevelfact;
+ origin_space.AT(0).AT(mglevel)[d] = cctkGH->cctk_origin_space[d];
+ delta_space.AT(mglevel)[d] = cctkGH->cctk_delta_space[d] / mglevelfact;
cctkGH->cctk_origin_space[d] = -424242.0;
cctkGH->cctk_delta_space[d] = -424242.0;
}
}
+ CCTK_INT const deadbeef = get_deadbeef();
+
// Set array information
for (int group=0; group<CCTK_NumGroups(); ++group) {
if (CCTK_GroupTypeI(group) != CCTK_GF) {
const int m = 0;
-// ivect::ref(const_cast<int*>(groupdata.at(group).info.nghostzones))
+// ivect::ref(const_cast<int*>(groupdata.AT(group).info.nghostzones))
// = deadbeef;
- ivect::ref(const_cast<int*>(groupdata.at(group).info.nghostzones))
- = arrdata.at(group).at(m).dd->ghost_width[0];
- ivect::ref(const_cast<int*>(groupdata.at(group).info.gsh))
+ ivect::ref(const_cast<int*>(groupdata.AT(group).info.nghostzones))
+ = arrdata.AT(group).AT(m).dd->ghost_width[0];
+ ivect::ref(const_cast<int*>(groupdata.AT(group).info.gsh))
= deadbeef;
- ivect::ref(const_cast<int*>(groupdata.at(group).info.lsh))
+ ivect::ref(const_cast<int*>(groupdata.AT(group).info.lsh))
= deadbeef;
- ivect::ref(const_cast<int*>(groupdata.at(group).info.lbnd))
+ ivect::ref(const_cast<int*>(groupdata.AT(group).info.lbnd))
= -deadbeef;
- ivect::ref(const_cast<int*>(groupdata.at(group).info.ubnd))
+ ivect::ref(const_cast<int*>(groupdata.AT(group).info.ubnd))
= deadbeef;
for (int d=0; d<dim; ++d) {
- const_cast<int*>(groupdata.at(group).info.bbox)[2*d ] = deadbeef;
- const_cast<int*>(groupdata.at(group).info.bbox)[2*d+1] = deadbeef;
+ const_cast<int*>(groupdata.AT(group).info.bbox)[2*d ] = deadbeef;
+ const_cast<int*>(groupdata.AT(group).info.bbox)[2*d+1] = deadbeef;
}
- groupdata.at(group).info.activetimelevels = deadbeef;
+ groupdata.AT(group).info.activetimelevels = deadbeef;
const int numvars = CCTK_NumVarsInGroupI (group);
if (numvars>0) {
@@ -236,6 +242,9 @@ namespace Carpet {
} // if grouptype
} // for group
+ // Set mode
+ cctkGH->cctk_mode = CCTK_MODE_META;
+
mglevel = -1;
mglevelfact = -deadbeef;
cctkGH->cctk_convlevel = -deadbeef;
@@ -255,27 +264,30 @@ namespace Carpet {
assert (rl>=0 and rl<reflevels);
Checkpoint ("Entering level mode");
+ // Set mode
+ cctkGH->cctk_mode = CCTK_MODE_LEVEL;
+
reflevel = rl;
- timereflevelfact = timereffacts.at (reflevel);
- spacereflevelfact = spacereffacts.at (reflevel);
+ timereflevelfact = timereffacts.AT (reflevel);
+ spacereflevelfact = spacereffacts.AT (reflevel);
ivect::ref(cctkGH->cctk_levfac) = spacereflevelfact;
cctkGH->cctk_timefac = timereflevelfact;
// Set number of time levels
for (int group=0; group<CCTK_NumGroups(); ++group) {
if (CCTK_GroupTypeI(group) == CCTK_GF) {
- groupdata.at(group).info.activetimelevels
- = groupdata.at(group).activetimelevels.at(mglevel).at(reflevel);
+ groupdata.AT(group).info.activetimelevels
+ = groupdata.AT(group).activetimelevels.AT(mglevel).AT(reflevel);
}
}
// Set current time
assert (mglevel>=0 and mglevel<(int)leveltimes.size());
- assert (reflevel>=0 and reflevel<(int)leveltimes.at(mglevel).size());
+ assert (reflevel>=0 and reflevel<(int)leveltimes.AT(mglevel).size());
if (not adaptive_stepsize) {
- cctkGH->cctk_time = leveltimes.at(mglevel).at(reflevel);
+ cctkGH->cctk_time = leveltimes.AT(mglevel).AT(reflevel);
} else {
- leveltimes.at(mglevel).at(reflevel) = cctkGH->cctk_time;
+ leveltimes.AT(mglevel).AT(reflevel) = cctkGH->cctk_time;
}
assert (is_level_mode());
@@ -290,11 +302,13 @@ namespace Carpet {
if (reflevel == -1) return; // early return
Checkpoint ("Leaving level mode");
-
+
+ CCTK_INT const deadbeef = get_deadbeef();
+
// Save and unset current time
assert (mglevel>=0 and mglevel<(int)leveltimes.size());
- assert (reflevel>=0 and reflevel<(int)leveltimes.at(mglevel).size());
- leveltimes.at(mglevel).at(reflevel) = cctkGH->cctk_time;
+ assert (reflevel>=0 and reflevel<(int)leveltimes.AT(mglevel).size());
+ leveltimes.AT(mglevel).AT(reflevel) = cctkGH->cctk_time;
if (not adaptive_stepsize) {
cctkGH->cctk_time = global_time;
} else {
@@ -304,13 +318,16 @@ namespace Carpet {
// Unset number of time levels
for (int group=0; group<CCTK_NumGroups(); ++group) {
if (CCTK_GroupTypeI(group) == CCTK_GF) {
- groupdata.at(group).info.activetimelevels = deadbeef;
+ groupdata.AT(group).info.activetimelevels = deadbeef;
}
}
+ // Set mode
+ cctkGH->cctk_mode = CCTK_MODE_GLOBAL;
+
reflevel = -1;
- timereflevelfact = timereffacts.at (reflevels - 1);
- // TODO: use spacereffacts.at (reflevel - 1) instead?
+ timereflevelfact = timereffacts.AT (reflevels - 1);
+ // TODO: use spacereffacts.AT (reflevel - 1) instead?
spacereflevelfact = ivect(-deadbeef);
ivect::ref(cctkGH->cctk_levfac) = spacereflevelfact;
cctkGH->cctk_timefac = timereflevelfact;
@@ -332,6 +349,9 @@ namespace Carpet {
or grouptype == CCTK_GF);
Checkpoint ("Entering singlemap mode");
+ // Set mode
+ cctkGH->cctk_mode = CCTK_MODE_SINGLEMAP;
+
assert (mc_grouptype == -1);
mc_grouptype = grouptype;
carpetGH.map = map = m;
@@ -341,26 +361,26 @@ namespace Carpet {
if (maps > 1) {
// Set space delta
for (int d=0; d<dim; ++d) {
- cctkGH->cctk_origin_space[d] = origin_space.at(map).at(mglevel)[d];
- cctkGH->cctk_delta_space[d] = delta_space.at(map)[d] * mglevelfact;
+ cctkGH->cctk_origin_space[d] = origin_space.AT(map).AT(mglevel)[d];
+ cctkGH->cctk_delta_space[d] = delta_space.AT(map)[d] * mglevelfact;
}
}
// Set grid shape
- const ibbox& coarseext = vhh.at(map)->baseextents.at(mglevel).at(0 );
- const ibbox& baseext = vhh.at(map)->baseextents.at(mglevel).at(reflevel);
+ const ibbox& coarseext = vhh.AT(map)->baseextents.AT(mglevel).AT(0 );
+ const ibbox& baseext = vhh.AT(map)->baseextents.AT(mglevel).AT(reflevel);
// assert (all (baseext.lower() % baseext.stride() == 0));
ivect::ref(cctkGH->cctk_levoff) = baseext.lower() - coarseext.lower();
ivect::ref(cctkGH->cctk_levoffdenom) = baseext.stride();
ivect::ref(cctkGH->cctk_gsh) = baseext.shape() / baseext.stride();
- assert (all (vdd.at(map)->ghost_width[0] == vdd.at(map)->ghost_width[1]));
- ivect::ref(cctkGH->cctk_nghostzones) = vdd.at(map)->ghost_width[0];
+ assert (all (vdd.AT(map)->ghost_width[0] == vdd.AT(map)->ghost_width[1]));
+ ivect::ref(cctkGH->cctk_nghostzones) = vdd.AT(map)->ghost_width[0];
for (int group=0; group<CCTK_NumGroups(); ++group) {
if (CCTK_GroupTypeI(group) == CCTK_GF) {
- ivect::ref(const_cast<int*>(groupdata.at(group).info.gsh))
+ ivect::ref(const_cast<int*>(groupdata.AT(group).info.gsh))
= ivect::ref(cctkGH->cctk_gsh);
- ivect::ref(const_cast<int*>(groupdata.at(group).info.nghostzones))
+ ivect::ref(const_cast<int*>(groupdata.AT(group).info.nghostzones))
= ivect::ref(cctkGH->cctk_nghostzones);
}
}
@@ -385,12 +405,14 @@ namespace Carpet {
if (mc_grouptype == CCTK_GF) {
+ CCTK_INT const deadbeef = get_deadbeef();
+
// Save space delta
// (Do this early and often, so that interpolation has access to
// the correct values right away.)
for (int d=0; d<dim; ++d) {
- origin_space.at(map).at(mglevel)[d] = cctkGH->cctk_origin_space[d];
- delta_space.at(map)[d] = cctkGH->cctk_delta_space[d] / mglevelfact;
+ origin_space.AT(map).AT(mglevel)[d] = cctkGH->cctk_origin_space[d];
+ delta_space.AT(map)[d] = cctkGH->cctk_delta_space[d] / mglevelfact;
}
if (maps > 1) {
// Unset space delta
@@ -405,19 +427,22 @@ namespace Carpet {
ivect::ref(cctkGH->cctk_levoffdenom) = 0;
ivect::ref(cctkGH->cctk_gsh) = deadbeef;
// ivect::ref(cctkGH->cctk_nghostzones) = deadbeef;
- ivect::ref(cctkGH->cctk_nghostzones) = vdd.at(map)->ghost_width[0];
+ ivect::ref(cctkGH->cctk_nghostzones) = vdd.AT(map)->ghost_width[0];
for (int group=0; group<CCTK_NumGroups(); ++group) {
if (CCTK_GroupTypeI(group) == CCTK_GF) {
- ivect::ref(const_cast<int*>(groupdata.at(group).info.gsh))
+ ivect::ref(const_cast<int*>(groupdata.AT(group).info.gsh))
= ivect::ref(cctkGH->cctk_gsh);
- ivect::ref(const_cast<int*>(groupdata.at(group).info.nghostzones))
+ ivect::ref(const_cast<int*>(groupdata.AT(group).info.nghostzones))
= ivect::ref(cctkGH->cctk_nghostzones);
}
}
} // if mc_grouptype
+ // Set mode
+ cctkGH->cctk_mode = CCTK_MODE_LEVEL;
+
mc_grouptype = -1;
carpetGH.map = map = -1;
@@ -428,25 +453,32 @@ namespace Carpet {
// Set local mode
- void enter_local_mode (cGH * const cctkGH, int const c, int const grouptype)
+ void enter_local_mode (cGH * const cctkGH,
+ int const c, int const lc, int const grouptype)
{
assert (is_singlemap_mode());
if (mc_grouptype == CCTK_GF) {
- assert (c>=0 and c<vhh.at(map)->components(reflevel));
+ assert (c>=0 and c<vhh.AT(map)->components(reflevel));
+ assert (lc==-1 or (lc>=0 and lc<vhh.AT(map)->local_components(reflevel)));
} else {
- assert (c>=0 and c<CCTK_nProcs(cctkGH));
+ assert (c>=0 and c<dist::size());
+ assert (lc==-1 or lc==0);
}
Checkpoint ("Entering local mode");
+ // Set mode
+ cctkGH->cctk_mode = CCTK_MODE_LOCAL;
+
assert (grouptype == mc_grouptype);
component = c;
+ local_component = lc;
if (mc_grouptype == CCTK_GF) {
// Set cGH fields
- const ibbox& baseext = vhh.at(map)->baseextents.at(mglevel).at(reflevel);
- const ibbox& ext = vdd.at(map)->boxes.at(mglevel).at(reflevel).at(component).exterior;
- const b2vect& obnds = vhh.at(map)->outer_boundaries(reflevel,component);
+ const ibbox& baseext = vhh.AT(map)->baseextents.AT(mglevel).AT(reflevel);
+ const ibbox& ext = vdd.AT(map)->boxes.AT(mglevel).AT(reflevel).AT(component).exterior;
+ const b2vect& obnds = vhh.AT(map)->outer_boundaries(reflevel,component);
ivect::ref(cctkGH->cctk_lsh) = ext.shape() / ext.stride();
ivect::ref(cctkGH->cctk_lbnd)
@@ -484,44 +516,47 @@ namespace Carpet {
for (int group=0; group<CCTK_NumGroups(); ++group) {
if (CCTK_GroupTypeI(group) == CCTK_GF) {
- ivect::ref(const_cast<int*>(groupdata.at(group).info.lsh))
+ ivect::ref(const_cast<int*>(groupdata.AT(group).info.lsh))
= ivect::ref(cctkGH->cctk_lsh);
- ivect::ref(const_cast<int*>(groupdata.at(group).info.lbnd))
+ ivect::ref(const_cast<int*>(groupdata.AT(group).info.lbnd))
= ivect::ref(cctkGH->cctk_lbnd);
- ivect::ref(const_cast<int*>(groupdata.at(group).info.ubnd))
+ ivect::ref(const_cast<int*>(groupdata.AT(group).info.ubnd))
= ivect::ref(cctkGH->cctk_ubnd);
for (int d=0; d<dim; ++d) {
- const_cast<int*>(groupdata.at(group).info.bbox)[2*d ]
+ const_cast<int*>(groupdata.AT(group).info.bbox)[2*d ]
= cctkGH->cctk_bbox[2*d ];
- const_cast<int*>(groupdata.at(group).info.bbox)[2*d+1]
+ const_cast<int*>(groupdata.AT(group).info.bbox)[2*d+1]
= cctkGH->cctk_bbox[2*d+1];
}
- const int numvars = CCTK_NumVarsInGroupI (group);
- if (numvars>0) {
- const int firstvar = CCTK_FirstVarIndexI (group);
- assert (firstvar>=0);
- const int max_tl = CCTK_MaxTimeLevelsGI (group);
- assert (max_tl>=0);
- const int active_tl = CCTK_ActiveTimeLevelsGI (cctkGH, group);
- assert (active_tl>=0 and active_tl<=max_tl);
- const int available_tl =
- do_allow_past_timelevels ? active_tl : min (1, active_tl);
-
- // assert (vhh.at(map)->is_local(reflevel,component));
-
- assert (group<(int)arrdata.size());
- for (int var=0; var<numvars; ++var) {
- assert (firstvar+var<CCTK_NumVars());
- ggf * const ff = arrdata.at(group).at(map).data.at(var);
- for (int tl=0; tl<max_tl; ++tl) {
- if (ff and tl<available_tl) {
- gdata * const data = (*ff) (tl, reflevel, component, mglevel);
- assert (data);
- cctkGH->data[firstvar+var][tl] = data->storage();
- } else {
- cctkGH->data[firstvar+var][tl] = NULL;
+ if (local_component != -1) {
+ const int numvars = CCTK_NumVarsInGroupI (group);
+ if (numvars>0) {
+ const int firstvar = CCTK_FirstVarIndexI (group);
+ assert (firstvar>=0);
+ const int max_tl = CCTK_MaxTimeLevelsGI (group);
+ assert (max_tl>=0);
+ const int active_tl = CCTK_ActiveTimeLevelsGI (cctkGH, group);
+ assert (active_tl>=0 and active_tl<=max_tl);
+ const int available_tl =
+ do_allow_past_timelevels ? active_tl : min (1, active_tl);
+
+ // assert (vhh.AT(map)->is_local(reflevel,component));
+
+ assert (group<(int)arrdata.size());
+ for (int var=0; var<numvars; ++var) {
+ assert (firstvar+var<CCTK_NumVars());
+ ggf * const ff = arrdata.AT(group).AT(map).data.AT(var);
+ for (int tl=0; tl<max_tl; ++tl) {
+ if (ff and tl<available_tl) {
+ gdata * const data =
+ (*ff) (tl, reflevel, local_component, mglevel);
+ assert (data);
+ cctkGH->data[firstvar+var][tl] = data->storage();
+ } else {
+ cctkGH->data[firstvar+var][tl] = NULL;
+ }
}
}
}
@@ -547,6 +582,8 @@ namespace Carpet {
if (mc_grouptype == CCTK_GF) {
+ CCTK_INT const deadbeef = get_deadbeef();
+
// Unset cGH fields
ivect::ref(cctkGH->cctk_lsh) = deadbeef;
ivect::ref(cctkGH->cctk_lbnd) = -deadbeef;
@@ -569,32 +606,34 @@ namespace Carpet {
for (int group=0; group<CCTK_NumGroups(); ++group) {
if (CCTK_GroupTypeI(group) == CCTK_GF) {
- ivect::ref(const_cast<int*>(groupdata.at(group).info.lsh))
+ ivect::ref(const_cast<int*>(groupdata.AT(group).info.lsh))
= ivect::ref(cctkGH->cctk_lsh);
- ivect::ref(const_cast<int*>(groupdata.at(group).info.lbnd))
+ ivect::ref(const_cast<int*>(groupdata.AT(group).info.lbnd))
= ivect::ref(cctkGH->cctk_lbnd);
- ivect::ref(const_cast<int*>(groupdata.at(group).info.ubnd))
+ ivect::ref(const_cast<int*>(groupdata.AT(group).info.ubnd))
= ivect::ref(cctkGH->cctk_ubnd);
for (int d=0; d<dim; ++d) {
- const_cast<int*>(groupdata.at(group).info.bbox)[2*d ]
+ const_cast<int*>(groupdata.AT(group).info.bbox)[2*d ]
= cctkGH->cctk_bbox[2*d ];
- const_cast<int*>(groupdata.at(group).info.bbox)[2*d+1]
+ const_cast<int*>(groupdata.AT(group).info.bbox)[2*d+1]
= cctkGH->cctk_bbox[2*d+1];
}
- const int numvars = CCTK_NumVarsInGroupI (group);
- if (numvars>0) {
- const int firstvar = CCTK_FirstVarIndexI (group);
- assert (firstvar>=0);
- const int max_tl = CCTK_MaxTimeLevelsGI (group);
- assert (max_tl>=0);
-
- assert (group<(int)arrdata.size());
- for (int var=0; var<numvars; ++var) {
- assert (firstvar+var<CCTK_NumVars());
- for (int tl=0; tl<max_tl; ++tl) {
- cctkGH->data[firstvar+var][tl] = NULL;
+ if (local_component != -1) {
+ const int numvars = CCTK_NumVarsInGroupI (group);
+ if (numvars>0) {
+ const int firstvar = CCTK_FirstVarIndexI (group);
+ assert (firstvar>=0);
+ const int max_tl = CCTK_MaxTimeLevelsGI (group);
+ assert (max_tl>=0);
+
+ assert (group<(int)arrdata.size());
+ for (int var=0; var<numvars; ++var) {
+ assert (firstvar+var<CCTK_NumVars());
+ for (int tl=0; tl<max_tl; ++tl) {
+ cctkGH->data[firstvar+var][tl] = NULL;
+ }
}
}
}
@@ -604,7 +643,11 @@ namespace Carpet {
} // if mc_grouptype
+ // Set mode
+ cctkGH->cctk_mode = CCTK_MODE_SINGLEMAP;
+
component = -1;
+ local_component = -1;
assert (is_singlemap_mode());
}
@@ -749,7 +792,8 @@ namespace Carpet {
bool map_iterator::done () const
{
- return grouptype == CCTK_GF ? m >= maps : m >= 1;
+ int const maxm = grouptype == CCTK_GF ? maps : 1;
+ return m >= maxm;
}
void map_iterator::step ()
@@ -763,7 +807,46 @@ namespace Carpet {
- // Local mode iterator
+ // processor local map iterator
+
+ local_map_iterator::local_map_iterator (cGH const * const cctkGH_,
+ int const grouptype_)
+ : cctkGH(const_cast<cGH*>(cctkGH_)), grouptype(grouptype_), m(0)
+ {
+ assert (grouptype == CCTK_GF
+ or grouptype == CCTK_ARRAY or grouptype == CCTK_SCALAR);
+ enter_singlemap_mode (cctkGH, m, grouptype);
+ }
+
+ local_map_iterator::~local_map_iterator ()
+ {
+ leave_singlemap_mode (cctkGH);
+ }
+
+ bool local_map_iterator::done () const
+ {
+ int const maxm = grouptype == CCTK_GF ? maps : 1;
+ return m >= maxm;
+ }
+
+ void local_map_iterator::step ()
+ {
+ while (true) {
+ ++ m;
+ if (done()) break;
+ int const maxlc =
+ grouptype == CCTK_GF ? vhh.AT(m)->local_components(reflevel) : 1;
+ if (maxlc > 0) break;
+ }
+ if (not done()) {
+ leave_singlemap_mode (cctkGH);
+ enter_singlemap_mode (cctkGH, m, grouptype);
+ }
+ }
+
+
+
+ // local mode iterator
component_iterator::component_iterator (cGH const * const cctkGH_,
int const grouptype_)
@@ -771,7 +854,6 @@ namespace Carpet {
{
assert (grouptype == CCTK_GF
or grouptype == CCTK_ARRAY or grouptype == CCTK_SCALAR);
- assert (is_singlemap_mode());
step ();
}
@@ -793,17 +875,21 @@ namespace Carpet {
++ c;
if (not done()) {
leave_local_mode (cctkGH);
- enter_local_mode (cctkGH, c, grouptype);
+ int const lc =
+ (grouptype == CCTK_GF
+ ? vhh.AT(map)->get_local_component(reflevel, c)
+ : (c == dist::rank() ? 0 : -1));
+ enter_local_mode (cctkGH, c, lc, grouptype);
}
}
- // Processor local mode iterator
+ // processor local mode iterator
local_component_iterator::local_component_iterator (cGH const * const cctkGH_,
int const grouptype_)
- : cctkGH(const_cast<cGH*>(cctkGH_)), grouptype(grouptype_), c(-1)
+ : cctkGH(const_cast<cGH*>(cctkGH_)), grouptype(grouptype_), lc(-1)
{
assert (grouptype == CCTK_GF
or grouptype == CCTK_ARRAY or grouptype == CCTK_SCALAR);
@@ -818,23 +904,21 @@ namespace Carpet {
bool local_component_iterator::done () const
{
- int const maxc = (grouptype == CCTK_GF
- ? vhh.AT(map)->components(reflevel)
- : dist::size());
- return c >= maxc;
+ int const maxlc = (grouptype == CCTK_GF
+ ? vhh.AT(map)->local_components(reflevel)
+ : 1);
+ return lc >= maxlc;
}
void local_component_iterator::step ()
{
- do {
- ++ c;
- } while (not done() and not (grouptype == CCTK_GF
- ? vhh.at(map)->is_local(reflevel, c)
- : c == CCTK_MyProc(cctkGH)));
-
+ ++ lc;
if (not done()) {
+ int const c = (grouptype == CCTK_GF
+ ? vhh.AT(map)->get_component(reflevel,lc)
+ : dist::rank());
leave_local_mode (cctkGH);
- enter_local_mode (cctkGH, c, grouptype);
+ enter_local_mode (cctkGH, c, lc, grouptype);
}
}
@@ -847,7 +931,7 @@ namespace Carpet {
// Singlemap escape
singlemap_escape::singlemap_escape (cGH const * const cctkGH_)
- : cctkGH(const_cast<cGH*>(cctkGH_)), c(component)
+ : cctkGH(const_cast<cGH*>(cctkGH_)), c(component), lc(local_component)
{
assert (not is_meta_mode());
assert (not is_global_mode());
@@ -861,7 +945,7 @@ namespace Carpet {
{
assert (is_singlemap_mode());
if (c != -1) {
- enter_local_mode (cctkGH, c, mc_grouptype);
+ enter_local_mode (cctkGH, c, lc, mc_grouptype);
}
}
@@ -871,7 +955,7 @@ namespace Carpet {
level_escape::level_escape (cGH const * const cctkGH_)
: cctkGH(const_cast<cGH*>(cctkGH_)),
- grouptype(mc_grouptype), m(map), c(component)
+ grouptype(mc_grouptype), m(map), c(component), lc(local_component)
{
assert (not is_meta_mode());
assert (not is_global_mode());
@@ -889,7 +973,7 @@ namespace Carpet {
if (m != -1) {
enter_singlemap_mode (cctkGH, m, grouptype);
if (c != -1) {
- enter_local_mode (cctkGH, c, grouptype);
+ enter_local_mode (cctkGH, c, lc, grouptype);
}
}
}
@@ -900,7 +984,8 @@ namespace Carpet {
global_escape::global_escape (cGH const * const cctkGH_)
: cctkGH(const_cast<cGH*>(cctkGH_)),
- rl(reflevel), grouptype(mc_grouptype), m(map), c(component)
+ rl(reflevel),
+ grouptype(mc_grouptype), m(map), c(component), lc(local_component)
{
assert (not is_meta_mode());
if (not is_global_mode()) {
@@ -922,7 +1007,7 @@ namespace Carpet {
if (m != -1) {
enter_singlemap_mode (cctkGH, m, grouptype);
if (c != -1) {
- enter_local_mode (cctkGH, c, grouptype);
+ enter_local_mode (cctkGH, c, lc, grouptype);
}
}
}
@@ -934,7 +1019,8 @@ namespace Carpet {
meta_escape::meta_escape (cGH const * const cctkGH_)
: cctkGH(const_cast<cGH*>(cctkGH_)),
- ml(mglevel), rl(reflevel), grouptype(mc_grouptype), m(map), c(component)
+ ml(mglevel), rl(reflevel),
+ grouptype(mc_grouptype), m(map), c(component), lc(local_component)
{
if (not is_meta_mode()) {
if (not is_global_mode()) {
@@ -960,7 +1046,7 @@ namespace Carpet {
if (m != -1) {
enter_singlemap_mode (cctkGH, m, grouptype);
if (c != -1) {
- enter_local_mode (cctkGH, c, grouptype);
+ enter_local_mode (cctkGH, c, lc, grouptype);
}
}
}
@@ -979,27 +1065,27 @@ namespace Carpet {
if (is_meta_mode()) {
BEGIN_MGLEVEL_LOOP(cctkGH) {
BEGIN_REFLEVEL_LOOP(cctkGH) {
- BEGIN_MAP_LOOP(cctkGH, CCTK_GF) {
+ BEGIN_LOCAL_MAP_LOOP(cctkGH, CCTK_GF) {
BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) {
function (cctkGH);
} END_LOCAL_COMPONENT_LOOP;
- } END_MAP_LOOP;
+ } END_LOCAL_MAP_LOOP;
} END_REFLEVEL_LOOP;
} END_MGLEVEL_LOOP;
} else if (is_global_mode()) {
BEGIN_REFLEVEL_LOOP(cctkGH) {
- BEGIN_MAP_LOOP(cctkGH, CCTK_GF) {
+ BEGIN_LOCAL_MAP_LOOP(cctkGH, CCTK_GF) {
BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) {
function (cctkGH);
} END_LOCAL_COMPONENT_LOOP;
- } END_MAP_LOOP;
+ } END_LOCAL_MAP_LOOP;
} END_REFLEVEL_LOOP;
} else if (is_level_mode()) {
- BEGIN_MAP_LOOP(cctkGH, CCTK_GF) {
+ BEGIN_LOCAL_MAP_LOOP(cctkGH, CCTK_GF) {
BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) {
function (cctkGH);
} END_LOCAL_COMPONENT_LOOP;
- } END_MAP_LOOP;
+ } END_LOCAL_MAP_LOOP;
} else if (is_singlemap_mode()) {
BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) {
function (cctkGH);
@@ -1140,7 +1226,13 @@ namespace Carpet {
: cctkGH(const_cast<cGH*>(cctkGH_))
{
assert (is_singlemap_mode());
- enter_local_mode (cctkGH, c, grouptype);
+ int const lc =
+ (c != -1
+ ? (grouptype == CCTK_GF
+ ? vhh.AT(map)->get_local_component(reflevel, c)
+ : (c == dist::rank() ? 0 : -1))
+ : -1);
+ enter_local_mode (cctkGH, c, lc, grouptype);
}
component_setter::~component_setter ()
diff --git a/Carpet/Carpet/src/modes.hh b/Carpet/Carpet/src/modes.hh
index 55a5cfd35..c2631d952 100644
--- a/Carpet/Carpet/src/modes.hh
+++ b/Carpet/Carpet/src/modes.hh
@@ -27,27 +27,27 @@ namespace Carpet {
// Mode indicators
- bool is_meta_mode ();
- bool is_global_mode ();
- bool is_level_mode ();
- bool is_singlemap_mode ();
- bool is_local_mode ();
+ bool is_meta_mode () CCTK_ATTRIBUTE_PURE;
+ bool is_global_mode () CCTK_ATTRIBUTE_PURE;
+ bool is_level_mode () CCTK_ATTRIBUTE_PURE;
+ bool is_singlemap_mode () CCTK_ATTRIBUTE_PURE;
+ bool is_local_mode () CCTK_ATTRIBUTE_PURE;
// Mode setting
- void enter_global_mode (cGH * const cctkGH, int const ml);
- void leave_global_mode (cGH * const cctkGH);
+ void enter_global_mode (cGH * cctkGH, int ml);
+ void leave_global_mode (cGH * cctkGH);
- void enter_level_mode (cGH * const cctkGH, int const rl);
- void leave_level_mode (cGH * const cctkGH);
+ void enter_level_mode (cGH * cctkGH, int rl);
+ void leave_level_mode (cGH * cctkGH);
- void enter_singlemap_mode (cGH * const cctkGH, int const m, int const grouptype);
- void leave_singlemap_mode (cGH * const cctkGH);
+ void enter_singlemap_mode (cGH * cctkGH, int m, int grouptype);
+ void leave_singlemap_mode (cGH * cctkGH);
- void enter_local_mode (cGH * const cctkGH, int const c, int const grouptype);
- void leave_local_mode (cGH * const cctkGH);
+ void enter_local_mode (cGH * cctkGH, int c, int lc, int grouptype);
+ void leave_local_mode (cGH * cctkGH);
@@ -57,7 +57,7 @@ namespace Carpet {
cGH * cctkGH;
int ml;
public:
- mglevel_iterator (cGH const * const cctkGH);
+ mglevel_iterator (cGH const * cctkGH);
~mglevel_iterator ();
bool done () const;
void step ();
@@ -67,7 +67,7 @@ namespace Carpet {
cGH * cctkGH;
int ml;
public:
- reverse_mglevel_iterator (cGH const * const cctkGH);
+ reverse_mglevel_iterator (cGH const * cctkGH);
~reverse_mglevel_iterator ();
bool done () const;
void step ();
@@ -77,7 +77,7 @@ namespace Carpet {
cGH * cctkGH;
int rl;
public:
- reflevel_iterator (cGH const * const cctkGH);
+ reflevel_iterator (cGH const * cctkGH);
~reflevel_iterator ();
bool done () const;
void step ();
@@ -87,7 +87,7 @@ namespace Carpet {
cGH * cctkGH;
int rl;
public:
- reverse_reflevel_iterator (cGH const * const cctkGH);
+ reverse_reflevel_iterator (cGH const * cctkGH);
~reverse_reflevel_iterator ();
bool done () const;
void step ();
@@ -103,12 +103,23 @@ namespace Carpet {
int grouptype;
int m;
public:
- map_iterator (cGH const * const cctkGH, int const grouptype);
+ map_iterator (cGH const * cctkGH, int grouptype);
~map_iterator ();
bool done () const;
void step ();
};
+ class local_map_iterator {
+ cGH * cctkGH;
+ int grouptype;
+ int m;
+ public:
+ local_map_iterator (cGH const * cctkGH, int grouptype);
+ ~local_map_iterator ();
+ bool done () const;
+ void step ();
+ };
+
// Loop over all components. If grouptype is CCTK_GF, then loop
// over grid function components. If grouptype is CCTK_ARRAY (or
// CCTK_SCALAR), then loop over grid array (or grid scalar)
@@ -120,7 +131,7 @@ namespace Carpet {
int grouptype;
int c;
public:
- component_iterator (cGH const * const cctkGH, int const grouptype);
+ component_iterator (cGH const * cctkGH, int grouptype);
~component_iterator ();
bool done () const;
void step ();
@@ -129,9 +140,9 @@ namespace Carpet {
class local_component_iterator {
cGH * cctkGH;
int grouptype;
- int c;
+ int lc;
public:
- local_component_iterator (cGH const * const cctkGH, int const grouptype);
+ local_component_iterator (cGH const * cctkGH, int grouptype);
~local_component_iterator ();
bool done () const;
void step ();
@@ -201,6 +212,18 @@ namespace Carpet {
map_loop_ = false; \
} while (false)
+#define BEGIN_LOCAL_MAP_LOOP(cctkGH, grouptype) \
+ do { \
+ bool local_map_loop_ = true; \
+ for (Carpet::local_map_iterator local_map_iter_(cctkGH, grouptype); \
+ not local_map_iter_.done(); \
+ local_map_iter_.step()) {
+#define END_LOCAL_MAP_LOOP \
+ } \
+ assert (local_map_loop_); \
+ local_map_loop_ = false; \
+ } while (false)
+
#define BEGIN_COMPONENT_LOOP(cctkGH, grouptype) \
do { \
bool component_loop_ = true; \
@@ -232,8 +255,9 @@ namespace Carpet {
class singlemap_escape {
cGH * cctkGH;
int c;
+ int lc;
public:
- singlemap_escape (cGH const * const cctkGH);
+ singlemap_escape (cGH const * cctkGH);
~singlemap_escape ();
};
@@ -242,8 +266,9 @@ namespace Carpet {
int grouptype;
int m;
int c;
+ int lc;
public:
- level_escape (cGH const * const cctkGH);
+ level_escape (cGH const * cctkGH);
~level_escape ();
};
@@ -253,8 +278,9 @@ namespace Carpet {
int grouptype;
int m;
int c;
+ int lc;
public:
- global_escape (cGH const * const cctkGH);
+ global_escape (cGH const * cctkGH);
~global_escape ();
};
@@ -265,8 +291,9 @@ namespace Carpet {
int grouptype;
int m;
int c;
+ int lc;
public:
- meta_escape (cGH const * const cctkGH);
+ meta_escape (cGH const * cctkGH);
~meta_escape ();
};
@@ -325,30 +352,28 @@ namespace Carpet {
class mglevel_setter {
cGH * cctkGH;
public:
- mglevel_setter (cGH const * const cctkGH, int const ml);
+ mglevel_setter (cGH const * cctkGH, int ml);
~mglevel_setter ();
};
class reflevel_setter {
cGH * cctkGH;
public:
- reflevel_setter (cGH const * const cctkGH, int const rl);
+ reflevel_setter (cGH const * cctkGH, int rl);
~reflevel_setter ();
};
class map_setter {
cGH * cctkGH;
public:
- map_setter (cGH const * const cctkGH,
- int const m, int const grouptype);
+ map_setter (cGH const * cctkGH, int m, int grouptype);
~map_setter ();
};
class component_setter {
cGH * cctkGH;
public:
- component_setter (cGH const * const cctkGH,
- int const c, int const grouptype);
+ component_setter (cGH const * cctkGH, int c, int grouptype);
~component_setter ();
};
diff --git a/Carpet/Carpet/src/variables.cc b/Carpet/Carpet/src/variables.cc
index a27a5bacb..0bf834ff4 100644
--- a/Carpet/Carpet/src/variables.cc
+++ b/Carpet/Carpet/src/variables.cc
@@ -1,5 +1,4 @@
-
-#include "variables.hh"
+#include <variables.hh>
@@ -51,6 +50,7 @@ namespace Carpet {
int mc_grouptype;
int map;
int component;
+ int local_component;
// Current refinement factors
int timereflevelfact;
@@ -90,9 +90,13 @@ namespace Carpet {
bool do_allow_past_timelevels;
// Is prolongation enabled?
+ // (This flag disables prolongation during MoL integration
+ // substeps.)
bool do_prolongate;
// Is tapering enabled?
+ // (This flag disables prolongation while the current refinement
+ // level is not aligned with the parent.)
bool do_taper;
// Should we warn about groups with insufficiently many time levels?
@@ -106,6 +110,7 @@ namespace Carpet {
vector<gh*> vhh; // [map]
vector<dh*> vdd; // [map]
vector<th*> vtt; // [map]
+ int regridding_epoch;
// Data for the groups
vector<groupdesc> groupdata; // [group]
diff --git a/Carpet/Carpet/src/variables.hh b/Carpet/Carpet/src/variables.hh
index 892664264..2746ec295 100644
--- a/Carpet/Carpet/src/variables.hh
+++ b/Carpet/Carpet/src/variables.hh
@@ -16,16 +16,16 @@
#include <mpi.h>
-#include "cctk.h"
+#include <cctk.h>
-#include "bbox.hh"
-#include "data.hh"
-#include "dh.hh"
-#include "ggf.hh"
-#include "gh.hh"
-#include "operators.hh"
-#include "th.hh"
-#include "vect.hh"
+#include <bbox.hh>
+#include <data.hh>
+#include <dh.hh>
+#include <ggf.hh>
+#include <gh.hh>
+#include <operators.hh>
+#include <th.hh>
+#include <vect.hh>
#include "carpet_public.h"
@@ -80,6 +80,7 @@ namespace Carpet {
extern int mc_grouptype; // -1, CCTK_SCALAR/CCTK_ARRAY, CCTK_GF
extern int map;
extern int component;
+ extern int local_component; // -1 for non-local
// Current refinement factors
extern int timereflevelfact;
@@ -126,9 +127,13 @@ namespace Carpet {
extern bool do_allow_past_timelevels;
// Is prolongation enabled?
+ // (This flag disables prolongation during MoL integration
+ // substeps.)
extern bool do_prolongate;
// Is tapering enabled?
+ // (This flag disables prolongation while the current refinement
+ // level is not aligned with the parent.)
extern bool do_taper;
// Should we warn about groups with insufficiently many time levels?
@@ -142,6 +147,7 @@ namespace Carpet {
extern vector<gh*> vhh; // [map]
extern vector<dh*> vdd; // [map]
extern vector<th*> vtt; // [map]
+ extern int regridding_epoch;
// Data for the groups
struct groupdesc {
diff --git a/Carpet/Carpet/test/small.par b/Carpet/Carpet/test/small.par
new file mode 100644
index 000000000..67b4aef5c
--- /dev/null
+++ b/Carpet/Carpet/test/small.par
@@ -0,0 +1,16 @@
+ActiveThorns = "
+ Boundary
+ CarpetLib
+ Carpet
+ CoordBase
+ IOUtil
+ InitBase
+ InterpToArray
+ LoopControl
+ SymBase
+"
+
+IO::out_dir = $parfile
+
+InterpToArray::nparrays1d = 1
+InterpToArray::parray1d_npoints_i = 1