diff options
author | Erik Schnetter <schnetter@cct.lsu.edu> | 2009-09-03 16:19:15 -0500 |
---|---|---|
committer | Barry Wardell <barry.wardell@gmail.com> | 2011-12-14 16:42:31 +0000 |
commit | 11c4d98017cbb86d08e15fd1b549180184b58a26 (patch) | |
tree | 2546a154c6f7bc0bec87de7316125ae7d1453569 /Carpet/Carpet | |
parent | f520477b1c14e02f1495cfa8d3e09f4e21ab34d0 (diff) |
Import Carpet
Ignore-this: 309b4dd613f4af2b84aa5d6743fdb6b3
Diffstat (limited to 'Carpet/Carpet')
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 |