From f1aa84c14ac110f5802603b4bf9ceb7923c5b471 Mon Sep 17 00:00:00 2001 From: schnetter <> Date: Wed, 4 Jul 2001 10:29:00 +0000 Subject: Broke the file "carpet.cc" into several files. Broke the file "carpet.cc" into several files. Simplified support for arrays with a dimension different from 3. Added a testing thorn. darcs-hash:20010704102946-07bb3-82132f0c643e91f0de33bbaef93f9c693ce7227f.gz --- Carpet/Carpet/src/CallFunction.cc | 292 +-- Carpet/Carpet/src/CarpetStartup.cc | 10 +- Carpet/Carpet/src/Checksum.cc | 313 +-- Carpet/Carpet/src/Comm.cc | 162 +- Carpet/Carpet/src/Cycle.cc | 135 +- Carpet/Carpet/src/Evolve.cc | 390 +--- Carpet/Carpet/src/Initialise.cc | 547 +---- Carpet/Carpet/src/Poison.cc | 269 +-- Carpet/Carpet/src/Recompose.cc | 908 +-------- Carpet/Carpet/src/Restrict.cc | 132 +- Carpet/Carpet/src/SetupGH.cc | 1041 ++-------- Carpet/Carpet/src/Shutdown.cc | 48 +- Carpet/Carpet/src/Storage.cc | 350 ++-- Carpet/Carpet/src/carpet.cc | 2138 -------------------- Carpet/Carpet/src/carpet.hh | 62 +- Carpet/Carpet/src/helpers.cc | 564 +++--- Carpet/Carpet/src/make.code.defn | 18 +- Carpet/Carpet/src/variables.cc | 90 +- Carpet/CarpetIOASCII/src/ioascii.cc | 288 +-- Carpet/CarpetLib/src/bboxset.cc | 4 +- Carpet/CarpetLib/src/data.cc | 4 +- Carpet/CarpetLib/src/defs.cc | 16 +- Carpet/CarpetLib/src/dh.cc | 4 +- Carpet/CarpetLib/src/gdata.cc | 36 +- Carpet/CarpetLib/src/gf.cc | 4 +- Carpet/CarpetLib/src/ggf.cc | 4 +- Carpet/CarpetLib/src/gh.cc | 4 +- Carpet/CarpetSlab/src/carpetslab.cc | 212 +- Carpet/CarpetTest/interface.ccl | 7 +- Carpet/CarpetTest/schedule.ccl | 6 +- .../CarpetTest/src/carpettest_check_arguments.F77 | 38 +- Carpet/CarpetTest/src/carpettest_check_sizes.c | 181 +- Carpet/CarpetTest/test/arraysizes.par | 27 +- CarpetAttic/CarpetIOFlexIO/src/ioflexio.cc | 228 +-- 34 files changed, 1865 insertions(+), 6667 deletions(-) delete mode 100644 Carpet/Carpet/src/carpet.cc diff --git a/Carpet/Carpet/src/CallFunction.cc b/Carpet/Carpet/src/CallFunction.cc index 4511c10ef..5f3b4aac2 100644 --- a/Carpet/Carpet/src/CallFunction.cc +++ b/Carpet/Carpet/src/CallFunction.cc @@ -1,19 +1,14 @@ #include #include -#include - #include "cctk.h" #include "cctki_GHExtensions.h" -#include "gh.hh" +#include "Carpet/CarpetLib/src/gh.hh" #include "carpet.hh" - -extern "C" { - static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/CallFunction.cc,v 1.19 2004/08/02 11:43:15 schnetter Exp $"; - CCTK_FILEVERSION(Carpet_Carpet_CallFunction_cc); -} + +static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/CallFunction.cc,v 1.1 2001/07/04 12:29:46 schnetter Exp $"; @@ -21,278 +16,47 @@ namespace Carpet { using namespace std; - /// Traverse one function on all components of one refinement level - /// of one multigrid level. - int CallFunction (void* function, ///< the function to call - cFunctionData* attribute, ///< attributes of the function - void* data) ///< ??? + int CallFunction (void* function, cFunctionData* attribute, void* data) { -// Checkpoint ("Starting CallFunction..."); + // Traverse one function on all components of one refinement level + // of one multigrid level - cGH* cgh = (cGH*)data; + assert (mglevel>=0); + assert (reflevel>=0); - assert (!! attribute->meta - + !! attribute->global - + !! attribute->level - + !! attribute->singlemap - + !! attribute->local <= 1); +// Checkpoint ("%*sStarting CallFunction...", 2*reflevel, ""); - assert (!! attribute->loop_global - + !! attribute->loop_level - + !! attribute->loop_singlemap - + !! attribute->loop_local <= 1); + cGH* cgh = (cGH*)data; - if (attribute->meta || is_meta_mode()) { - // Convtest operation - - if (do_meta_mode) { - if (attribute->loop_local) { - BEGIN_META_MODE(cgh) { - BEGIN_MGLEVEL_LOOP(cgh) { - BEGIN_REFLEVEL_LOOP(cgh) { - BEGIN_MAP_LOOP(cgh, CCTK_GF) { - BEGIN_LOCAL_COMPONENT_LOOP(cgh, CCTK_GF) { - Checkpoint ("Meta time local mode call at %s to %s::%s", - attribute->where, attribute->thorn, attribute->routine); - const int res = CCTK_CallFunction (function, attribute, data); - assert (res==0); - } END_LOCAL_COMPONENT_LOOP; - } END_MAP_LOOP; - for (int n=0; nn_SyncGroups; ++n) { - char * const groupname = CCTK_GroupName (attribute->SyncGroups[n]); - SyncGroup (cgh, groupname); - free (groupname); - } - } END_REFLEVEL_LOOP; - } END_MGLEVEL_LOOP; - } END_META_MODE; - } else if (attribute->loop_singlemap) { - BEGIN_META_MODE(cgh) { - BEGIN_MGLEVEL_LOOP(cgh) { - BEGIN_REFLEVEL_LOOP(cgh) { - BEGIN_MAP_LOOP(cgh, CCTK_GF) { - Checkpoint ("Meta time singlemap mode call at %s to %s::%s", - attribute->where, attribute->thorn, attribute->routine); - const int res = CCTK_CallFunction (function, attribute, data); - assert (res==0); - } END_MAP_LOOP; - for (int n=0; nn_SyncGroups; ++n) { - char * const groupname = CCTK_GroupName (attribute->SyncGroups[n]); - SyncGroup (cgh, groupname); - free (groupname); - } - } END_REFLEVEL_LOOP; - } END_MGLEVEL_LOOP; - } END_META_MODE; - } else if (attribute->loop_level) { - BEGIN_META_MODE(cgh) { - BEGIN_MGLEVEL_LOOP(cgh) { - BEGIN_REFLEVEL_LOOP(cgh) { - Checkpoint ("Meta time level mode call at %s to %s::%s", - attribute->where, attribute->thorn, attribute->routine); - const int res = CCTK_CallFunction (function, attribute, data); - assert (res==0); - for (int n=0; nn_SyncGroups; ++n) { - char * const groupname = CCTK_GroupName (attribute->SyncGroups[n]); - SyncGroup (cgh, groupname); - free (groupname); - } - } END_REFLEVEL_LOOP; - } END_MGLEVEL_LOOP; - } END_META_MODE; - } else if (attribute->loop_global) { - BEGIN_META_MODE(cgh) { - BEGIN_MGLEVEL_LOOP(cgh) { - Checkpoint ("Meta time global mode call at %s to %s::%s", - attribute->where, attribute->thorn, attribute->routine); - const int res = CCTK_CallFunction (function, attribute, data); - assert (res==0); - for (int n=0; nn_SyncGroups; ++n) { - char * const groupname = CCTK_GroupName (attribute->SyncGroups[n]); - SyncGroup (cgh, groupname); - free (groupname); - } - } END_MGLEVEL_LOOP; - } END_META_MODE; - } else { - BEGIN_META_MODE(cgh) { - Checkpoint ("Meta mode call at %s to %s::%s", - attribute->where, attribute->thorn, attribute->routine); - const int res = CCTK_CallFunction (function, attribute, data); - assert (res==0); - for (int n=0; nn_SyncGroups; ++n) { - char * const groupname = CCTK_GroupName (attribute->SyncGroups[n]); - SyncGroup (cgh, groupname); - free (groupname); - } - } END_META_MODE; - } - } - - } else if (attribute->global || is_global_mode()) { - // Global operation: call once - - assert (! attribute->loop_meta); + if (attribute->global) { + // Global operation: call once per refinement level - if (do_global_mode) { - if (attribute->loop_local) { - BEGIN_GLOBAL_MODE(cgh) { - BEGIN_REFLEVEL_LOOP(cgh) { - BEGIN_MAP_LOOP(cgh, CCTK_GF) { - BEGIN_LOCAL_COMPONENT_LOOP(cgh, CCTK_GF) { - Checkpoint ("Global time local mode call at %s to %s::%s", - attribute->where, attribute->thorn, attribute->routine); - const int res = CCTK_CallFunction (function, attribute, data); - assert (res==0); - } END_LOCAL_COMPONENT_LOOP; - } END_MAP_LOOP; - for (int n=0; nn_SyncGroups; ++n) { - char * const groupname = CCTK_GroupName (attribute->SyncGroups[n]); - SyncGroup (cgh, groupname); - free (groupname); - } - } END_REFLEVEL_LOOP; - } END_GLOBAL_MODE; - } else if (attribute->loop_singlemap) { - BEGIN_GLOBAL_MODE(cgh) { - BEGIN_REFLEVEL_LOOP(cgh) { - BEGIN_MAP_LOOP(cgh, CCTK_GF) { - Checkpoint ("Global time singlemap mode call at %s to %s::%s", - attribute->where, attribute->thorn, attribute->routine); - const int res = CCTK_CallFunction (function, attribute, data); - assert (res==0); - } END_MAP_LOOP; - for (int n=0; nn_SyncGroups; ++n) { - char * const groupname = CCTK_GroupName (attribute->SyncGroups[n]); - SyncGroup (cgh, groupname); - free (groupname); - } - } END_REFLEVEL_LOOP; - } END_GLOBAL_MODE; - } else if (attribute->loop_level) { - BEGIN_GLOBAL_MODE(cgh) { - BEGIN_REFLEVEL_LOOP(cgh) { - Checkpoint ("Global time level mode call at %s to %s::%s", - attribute->where, attribute->thorn, attribute->routine); - const int res = CCTK_CallFunction (function, attribute, data); - assert (res==0); - for (int n=0; nn_SyncGroups; ++n) { - char * const groupname = CCTK_GroupName (attribute->SyncGroups[n]); - SyncGroup (cgh, groupname); - free (groupname); - } - } END_REFLEVEL_LOOP; - } END_GLOBAL_MODE; - } else { - BEGIN_GLOBAL_MODE(cgh) { - Checkpoint ("Global mode call at %s to %s::%s", - attribute->where, attribute->thorn, attribute->routine); - const int res = CCTK_CallFunction (function, attribute, data); - assert (res==0); - for (int n=0; nn_SyncGroups; ++n) { - char * const groupname = CCTK_GroupName (attribute->SyncGroups[n]); - SyncGroup (cgh, groupname); - free (groupname); - } - } END_GLOBAL_MODE; - } - } - - } else if (attribute->level) { - // Level operation: call once per refinement level - - assert (! attribute->loop_meta); - assert (! attribute->loop_global); - - if (attribute->loop_local) { - BEGIN_MAP_LOOP(cgh, CCTK_GF) { - BEGIN_LOCAL_COMPONENT_LOOP(cgh, CCTK_GF) { - Checkpoint ("Level time local mode call at %s to %s::%s", - attribute->where, attribute->thorn, attribute->routine); - const int res = CCTK_CallFunction (function, attribute, data); - assert (res==0); - } END_LOCAL_COMPONENT_LOOP; - } END_MAP_LOOP; - } else if (attribute->loop_singlemap) { - BEGIN_MAP_LOOP(cgh, CCTK_GF) { - Checkpoint ("Level time singlemap mode call at %s to %s::%s", - attribute->where, attribute->thorn, attribute->routine); - const int res = CCTK_CallFunction (function, attribute, data); - assert (res==0); - } END_MAP_LOOP; - } else { - Checkpoint ("Level mode call at %s to %s::%s", - attribute->where, attribute->thorn, attribute->routine); - const int res = CCTK_CallFunction (function, attribute, data); - assert (res==0); - } - for (int n=0; nn_SyncGroups; ++n) { - char * const groupname = CCTK_GroupName (attribute->SyncGroups[n]); - SyncGroup (cgh, groupname); - free (groupname); - } - - } else if (attribute->singlemap) { - // Single map operation: call once per refinement level and map - - assert (! attribute->loop_meta); - assert (! attribute->loop_global); - assert (! attribute->loop_level); - - if (attribute->loop_local) { - BEGIN_MAP_LOOP(cgh, CCTK_GF) { - BEGIN_LOCAL_COMPONENT_LOOP(cgh, CCTK_GF) { - Checkpoint ("Singlemap time local mode call at %s to %s::%s", - attribute->where, attribute->thorn, attribute->routine); - const int res = CCTK_CallFunction (function, attribute, data); - assert (res==0); - } END_LOCAL_COMPONENT_LOOP; - } END_MAP_LOOP; - } else { - BEGIN_MAP_LOOP(cgh, CCTK_GF) { - Checkpoint ("Singlemap mode call at %s to %s::%s", - attribute->where, attribute->thorn, attribute->routine); - const int res = CCTK_CallFunction (function, attribute, data); - assert (res==0); - } END_MAP_LOOP; - } - for (int n=0; nn_SyncGroups; ++n) { - char * const groupname = CCTK_GroupName (attribute->SyncGroups[n]); - SyncGroup (cgh, groupname); - free (groupname); - } + const int res = CCTK_CallFunction (function, attribute, data); + assert (res==0); } else { // Local operation: call once per component - assert (! attribute->loop_meta); - assert (! attribute->loop_global); - assert (! attribute->loop_level); - assert (! attribute->loop_singlemap); - - BEGIN_MAP_LOOP(cgh, CCTK_GF) { - BEGIN_LOCAL_COMPONENT_LOOP(cgh, CCTK_GF) { - Checkpoint ("Local mode call at %s to %s::%s", - attribute->where, attribute->thorn, attribute->routine); - const int res = CCTK_CallFunction (function, attribute, data); - assert (res==0); - } END_LOCAL_COMPONENT_LOOP; - } END_MAP_LOOP; - for (int n=0; nn_SyncGroups; ++n) { - char * const groupname = CCTK_GroupName (attribute->SyncGroups[n]); - SyncGroup (cgh, groupname); - free (groupname); - } + BEGIN_COMPONENT_LOOP(cgh) { + // This requires that all processors have the same number of + // local components + if (hh->is_local(reflevel, component)) { + + const int res = CCTK_CallFunction (function, attribute, data); + assert (res==0); + + } // if is_local + + } END_COMPONENT_LOOP(cgh); } -// Checkpoint ("done with CallFunction."); +// Checkpoint ("%*sdone with CallFunction.", 2*reflevel, ""); // The return value indicates whether the grid functions have been // synchronised. - // 0: let the flesh do the synchronisation - // 1: we did the synchronisation - return 1; + // return 0: let the flesh do the synchronisation, if necessary + return 0; } } // namespace Carpet diff --git a/Carpet/Carpet/src/CarpetStartup.cc b/Carpet/Carpet/src/CarpetStartup.cc index c55c99d9a..b82589fe5 100644 --- a/Carpet/Carpet/src/CarpetStartup.cc +++ b/Carpet/Carpet/src/CarpetStartup.cc @@ -5,10 +5,7 @@ #include "carpet.hh" -extern "C" { - static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/CarpetStartup.cc,v 1.5 2003/09/19 16:04:31 schnetter Exp $"; - CCTK_FILEVERSION(Carpet_Carpet_CarpetStartup_cc); -} +static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/CarpetStartup.cc,v 1.1 2001/07/04 12:29:46 schnetter Exp $"; @@ -16,7 +13,7 @@ namespace Carpet { using namespace std; - void CarpetStartup() + int CarpetStartup() { CCTK_RegisterBanner ("AMR driver provided by Carpet"); @@ -39,7 +36,8 @@ namespace Carpet { CCTK_OverloadnProcs (nProcs); CCTK_OverloadArrayGroupSizeB (ArrayGroupSizeB); CCTK_OverloadQueryGroupStorageB (QueryGroupStorageB); - CCTK_OverloadGroupDynamicData (GroupDynamicData); + + return 0; } } // namespace Carpet diff --git a/Carpet/Carpet/src/Checksum.cc b/Carpet/Carpet/src/Checksum.cc index 6943e681a..a24461318 100644 --- a/Carpet/Carpet/src/Checksum.cc +++ b/Carpet/Carpet/src/Checksum.cc @@ -1,19 +1,14 @@ #include #include -#include - #include "cctk.h" #include "cctk_Parameters.h" -#include "gh.hh" +#include "Carpet/CarpetLib/src/gh.hh" #include "carpet.hh" -extern "C" { - static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Checksum.cc,v 1.15 2004/03/24 17:44:51 schnetter Exp $"; - CCTK_FILEVERSION(Carpet_Carpet_Checksum_cc); -} +static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Checksum.cc,v 1.1 2001/07/04 12:29:46 schnetter Exp $"; @@ -23,159 +18,201 @@ namespace Carpet { - // Checksum information - struct ckdesc { - bool valid; - unsigned int sum; - }; - - // Helper class - struct ckdesc4 { - vector > > > a; // [m][c][var][tl] - }; - - // Checksum information - vector > > checksums; // [rl][ml][group] - - - - void CalculateChecksums (const cGH* cgh, const checktimes where) + void CalculateChecksums (cGH* cgh, const checktimes where) { DECLARE_CCTK_PARAMETERS; if (! checksum_timelevels) return; - Checkpoint ("CalculateChecksums"); + Checkpoint ("%*sCalculateChecksums", 2*reflevel, ""); + + checksums.resize(CCTK_NumVars()); + for (int group=0; group0); + const int min_tl = mintl(where, num_tl); + const int max_tl = maxtl(where, num_tl); + + checksums[n].resize(maxreflevels); + checksums[n][reflevel].resize(num_tl); + for (int tl=min_tl; tl<=max_tl; ++tl) { + switch (CCTK_GroupTypeFromVarI(n)) { + + case CCTK_SCALAR: { + checksums[n][reflevel][tl].resize(1); + const int c=0; + checksums[n][reflevel][tl][c].valid = false; + break; + } + + case CCTK_ARRAY: + case CCTK_GF: { + checksums[n][reflevel][tl].resize(hh->components(reflevel)); + BEGIN_COMPONENT_LOOP(cgh) { + checksums[n][reflevel][tl][component].valid = false; + } END_COMPONENT_LOOP(cgh); + break; + } + + default: + abort(); + } + } + } + } - checksums.resize(maxreflevels); - checksums.at(reflevel).resize(mglevels); - checksums.at(reflevel).at(mglevel).resize(CCTK_NumGroups()); for (int group=0; groupcomponents(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); - if (nvars > 0) { - - const int n0 = CCTK_FirstVarIndexI(group); - const int sz = CCTK_VarTypeSize(CCTK_VarTypeI(n0)); - assert (sz>0); - - ivect size(1); - const int gpdim = groupdata.at(group).info.dim; - for (int d=0; d0); - const int min_tl = mintl(where, num_tl); - const int max_tl = maxtl(where, num_tl); - - for (int var=0; vardata[n][tl]; - unsigned int chk = 0; - for (int i=0; i0); + const int num_tl = CCTK_NumTimeLevelsFromVarI(n); + assert (num_tl>0); + const int min_tl = mintl(where, num_tl); + const int max_tl = maxtl(where, num_tl); + + for (int tl=min_tl; tl<=max_tl; ++tl) { + switch (CCTK_GroupTypeFromVarI(n)) { + + case CCTK_SCALAR: { + const int c=0; + int chk = 0; + const void* data = cgh->data[n][tl]; + for (int i=0; iis_local(reflevel,component)) { + const int gpdim = CCTK_GroupDimI(group); + int np = 1; + for (int d=0; ddata[n][tl]; + int chk = 0; + for (int i=0; icomponents(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); - if (nvars > 0) { - - const int n0 = CCTK_FirstVarIndexI(group); - const int sz = CCTK_VarTypeSize(CCTK_VarTypeI(n0)); - assert (sz>0); - - ivect size(1); - const int gpdim = groupdata.at(group).info.dim; - for (int d=0; d0); - const int min_tl = mintl(where, num_tl); - const int max_tl = maxtl(where, num_tl); - - for (int var=0; vardata[n][tl]; - unsigned int chk = 0; - for (int i=0; i0); + const int num_tl = CCTK_NumTimeLevelsFromVarI(n); + assert (num_tl>0); + const int min_tl = mintl(where, num_tl); + const int max_tl = maxtl(where, num_tl); + + assert ((int)checksums[n].size()==maxreflevels); + assert ((int)checksums[n][reflevel].size()==num_tl); + + for (int tl=min_tl; tl<=max_tl; ++tl) { + + bool unexpected_change = false; + + switch (CCTK_GroupTypeFromVarI(n)) { + + case CCTK_SCALAR: { + assert ((int)checksums[n][reflevel][tl].size()==1); + const int c=0; + if (checksums[n][reflevel][tl][c].valid) { + int chk = 0; + const void* data = cgh->data[n][tl]; + for (int i=0; icomponents(reflevel)); + BEGIN_COMPONENT_LOOP(cgh) { + if (checksums[n][reflevel][tl][component].valid) { + if (hh->is_local(reflevel,component)) { + const int gpdim = CCTK_GroupDimI(group); + int np = 1; + for (int d=0; ddata[n][tl]; + int chk = 0; + for (int i=0; i=0 && groupcomponents(reflevel) > 1) assert (component == -1); - Checkpoint ("SyncGroup \"%s\" time=%g", groupname, (double)cgh->cctk_time); + Checkpoint ("%*sSyncGroup %s", 2*reflevel, "", groupname); - const int grouptype = CCTK_GroupTypeI(group); - - if (grouptype == CCTK_GF) { - if (reflevel == -1) { - CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, - "Cannot synchronise in global mode " - "(Tried to synchronise group \"%s\")", - groupname); - } - if (map != -1 && component == -1) { - if (maps == 1) { - CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING, - "Synchronising group \"%s\" in singlemap mode", - groupname); - } else { - CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, - "Cannot synchronise in singlemap mode " - "(Tried to synchronise group \"%s\")", - groupname); - } - } - if (component != -1) { - if (maps == 1 && vhh.at(map)->local_components(reflevel) == 1) { - CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING, - "Synchronising group \"%s\" in local mode", - groupname); - } else { - CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, - "Cannot synchronise in local mode " - "(Tried to synchronise group \"%s\")", - groupname); - } - } - } + const int group = CCTK_GroupIndex(groupname); + assert (group>=0 && group=0); const int num_tl = CCTK_NumTimeLevelsFromVarI(n0); assert (num_tl>0); const int tl = 0; - // Prolongate the boundaries - if (do_prolongate) { - switch (grouptype) { - - case CCTK_GF: - assert (reflevel>=0 && reflevel 0) { - - // use the current time here (which may be modified by the - // user) - const CCTK_REAL time - = (cgh->cctk_time - cctk_initial_time) / delta_time; - - for (comm_state state; !state.done(); state.step()) { - for (int m=0; m<(int)arrdata.at(group).size(); ++m) { - for (int var=0; varcomponents(reflevel); ++c) { - arrdata.at(group).at(m).data.at(var)->ref_bnd_prolongate - (state, tl, reflevel, c, mglevel, time); - } - } - } - } // for state - } // if reflevel>0 - break; - - case CCTK_SCALAR: - case CCTK_ARRAY: - // do nothing - break; - - default: - assert (0); - } // switch grouptype - } // if do_prolongate - - // Sync - for (comm_state state; !state.done(); state.step()) { - switch (CCTK_GroupTypeI(group)) { - - case CCTK_GF: - for (int m=0; m<(int)arrdata.at(group).size(); ++m) { - for (int var=0; varcomponents(reflevel); ++c) { - arrdata.at(group).at(m).data.at(var)->sync - (state, tl, reflevel, c, mglevel); - } - } - } - break; - - case CCTK_SCALAR: - case CCTK_ARRAY: - for (int var=0; var<(int)arrdata.at(group).at(0).data.size(); ++var) { - arrdata.at(group).at(0).data.at(var)->sync (state, 0, 0, 0, 0); - } - break; - - default: - assert (0); - } // switch grouptype - } // for state + switch (CCTK_GroupTypeI(group)) { + + case CCTK_SCALAR: + // TODO: Check whether the local values are consistent + break; + + case CCTK_ARRAY: + assert (group<(int)arrdata.size()); + for (int var=0; var<(int)arrdata[group].data.size(); ++var) { + if (num_tl>1 && reflevel>0) { + for (int c=0; ccomponents(reflevel); ++c) { + arrdata[group].data[var]->ref_bnd_prolongate + (tl, reflevel, c, mglevel, + prolongation_order_space, prolongation_order_time); + } + } + for (int c=0; ccomponents(reflevel); ++c) { + arrdata[group].data[var]->sync (tl, reflevel, c, mglevel); + } + } + break; + + case CCTK_GF: + assert (group<(int)gfdata.size()); + for (int var=0; var<(int)gfdata[group].data.size(); ++var) { + if (num_tl>1 && reflevel>0) { + for (int c=0; ccomponents(reflevel); ++c) { + gfdata[group].data[var]->ref_bnd_prolongate + (tl, reflevel, c, mglevel, + prolongation_order_space, prolongation_order_time); + } + } + for (int c=0; ccomponents(reflevel); ++c) { + gfdata[group].data[var]->sync (tl, reflevel, c, mglevel); + } + } + break; + + default: + abort(); + } return 0; } - int EnableGroupComm (const cGH* cgh, const char* groupname) + int EnableGroupComm (cGH* cgh, const char* groupname) { // Communication is always enabled return 0; } - int DisableGroupComm (const cGH* cgh, const char* groupname) + int DisableGroupComm (cGH* cgh, const char* groupname) { // Communication is always enabled return -1; diff --git a/Carpet/Carpet/src/Cycle.cc b/Carpet/Carpet/src/Cycle.cc index 69a195563..d0c841cae 100644 --- a/Carpet/Carpet/src/Cycle.cc +++ b/Carpet/Carpet/src/Cycle.cc @@ -3,15 +3,12 @@ #include "cctk.h" -#include "ggf.hh" -#include "gh.hh" +#include "Carpet/CarpetLib/src/ggf.hh" +#include "Carpet/CarpetLib/src/gh.hh" #include "carpet.hh" -extern "C" { - static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Cycle.cc,v 1.19 2004/05/21 18:16:23 schnetter Exp $"; - CCTK_FILEVERSION(Carpet_Carpet_Cycle_cc); -} +static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Cycle.cc,v 1.1 2001/07/04 12:29:46 schnetter Exp $"; @@ -21,95 +18,51 @@ namespace Carpet { - void CycleTimeLevels (const cGH* cgh) + void CycleTimeLevels (cGH* cgh) { - Checkpoint ("CycleTimeLevels"); - assert (is_level_mode()); + Checkpoint ("%*sCycleTimeLevels", 2*reflevel, ""); for (int group=0; group=0 && reflevelcomponents(reflevel); ++c) { - arrdata.at(group).at(m).data.at(var)->cycle (reflevel, c, mglevel); - } - } - } - break; - - case CCTK_SCALAR: - case CCTK_ARRAY: - if (do_global_mode) { - for (int var=0; varcomponents(0); ++c) { - arrdata.at(group).at(0).data.at(var)->cycle (0, c, mglevel); - } - } - } - break; - - default: - assert (0); - } // switch grouptype - } // if storage - } // for group - } - - - - void FlipTimeLevels (const cGH* cgh) - { - Checkpoint ("FlipTimeLevels"); - assert (is_level_mode()); - - for (int group=0; group0) { - const int var0 = CCTK_FirstVarIndexI(group); - assert (var0>=0); - - switch (CCTK_GroupTypeI(group)) { - - case CCTK_GF: - for (int m=0; m<(int)arrdata.at(group).size(); ++m) { - for (int var=0; varcomponents(reflevel); ++c) { - arrdata.at(group).at(m).data.at(var)->flip (reflevel, c, mglevel); - } - } - } - break; - - case CCTK_SCALAR: - case CCTK_ARRAY: - if (do_global_mode) { - for (int var=0; varcomponents(0); ++c) { - arrdata.at(group).at(0).data.at(var)->flip (0, c, mglevel); - } - } - } - break; - - default: - assert (0); - } // switch grouptype - - } // if num_vars>0 - } // if storage - } // for group + for (int var=0; var0); + void* tmpdata = scdata[group].data[var][reflevel][0]; + for (int tl=1; tlcomponents(reflevel); ++c) { + arrdata[group].data[var]->cycle (reflevel, c, mglevel); + } + break; + } + case CCTK_GF: { + assert (group<(int)gfdata.size()); + assert (var<(int)gfdata[group].data.size()); + for (int c=0; ccomponents(reflevel); ++c) { + gfdata[group].data[var]->cycle (reflevel, c, mglevel); + } + break; + } + default: + abort(); + } + } + } + } } } // namespace Carpet diff --git a/Carpet/Carpet/src/Evolve.cc b/Carpet/Carpet/src/Evolve.cc index e226a8330..a9e464c44 100644 --- a/Carpet/Carpet/src/Evolve.cc +++ b/Carpet/Carpet/src/Evolve.cc @@ -3,37 +3,12 @@ #include "cctk.h" #include "cctk_Parameters.h" -#include "cctk_Termination.h" -// IRIX wants this before -#if HAVE_SYS_TYPES_H -# include -#endif - -#if TIME_WITH_SYS_TIME -# include -# include -#else -# if HAVE_SYS_TIME_H -# include -# elif HAVE_TIME_H -# include -# endif -#endif - -#if HAVE_UNISTD_H -# include -#endif - -#include "dist.hh" -#include "th.hh" +#include "Carpet/CarpetLib/src/th.hh" #include "carpet.hh" -extern "C" { - static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Evolve.cc,v 1.54 2004/08/19 15:38:20 schnetter Exp $"; - CCTK_FILEVERSION(Carpet_Carpet_Evolve_cc); -} +static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Evolve.cc,v 1.1 2001/07/04 12:29:46 schnetter Exp $"; @@ -43,307 +18,96 @@ namespace Carpet { - static bool do_terminate (const cGH *cgh, - const CCTK_REAL time, const int iteration) - { - DECLARE_CCTK_PARAMETERS; - - bool term; - - // Early shortcut - if (terminate_next || CCTK_TerminationReached(cgh)) { - - term = true; - - } else { - - const bool term_iter = iteration >= cctk_itlast; - const bool term_time - = (delta_time > 0 - ? time >= cctk_final_time - 1.0e-8 * cgh->cctk_delta_time - : time <= cctk_final_time - 1.0e-8 * cgh->cctk_delta_time); -#ifdef HAVE_TIME_GETTIMEOFDAY - // get the current time - struct timeval tv; - gettimeofday (&tv, 0); - const double thetime = tv.tv_sec + tv.tv_usec / 1e6; - - static bool firsttime = true; - static double initial_runtime; - if (firsttime) { - firsttime = false; - initial_runtime = thetime; - } - - const double runtime = thetime - initial_runtime; - const bool term_runtime = (max_runtime > 0 - && runtime >= 60.0 * max_runtime); -#else - const bool term_runtime = false; -#endif - - if (CCTK_Equals(terminate, "never")) { - term = false; - } else if (CCTK_Equals(terminate, "iteration")) { - term = term_iter; - } else if (CCTK_Equals(terminate, "time")) { - term = term_time; - } else if (CCTK_Equals(terminate, "runtime")) { - term = term_runtime; - } else if (CCTK_Equals(terminate, "any")) { - term = term_iter || term_time || term_runtime; - } else if (CCTK_Equals(terminate, "all")) { - term = term_iter && term_time && term_runtime; - } else if (CCTK_Equals(terminate, "either")) { - term = term_iter || term_time; - } else if (CCTK_Equals(terminate, "both")) { - term = term_iter && term_time; - } else if (CCTK_Equals(terminate, "immediately")) { - term = true; - } else { - CCTK_WARN (0, "Unsupported termination condition"); - } - - } - - { - int local, global; - local = term; - MPI_Allreduce (&local, &global, 1, MPI_INT, MPI_LOR, dist::comm); - term = global; - } - - return term; - } - - - int Evolve (tFleshConfig* fc) { DECLARE_CCTK_PARAMETERS; - Waypoint ("Starting evolution loop"); + Checkpoint ("starting Evolve..."); const int convlev = 0; cGH* cgh = fc->GH[convlev]; // Main loop - while (! do_terminate(cgh, cgh->cctk_time, cgh->cctk_iteration)) { - - + while (cgh->cctk_iteration < cctk_itlast + || (cctk_final_time >= cctk_initial_time + && cgh->cctk_time < cctk_final_time)) { // Advance time ++cgh->cctk_iteration; - global_time = cctk_initial_time - + cgh->cctk_iteration * delta_time / maxreflevelfact; - cgh->cctk_time = global_time; - if ((cgh->cctk_iteration-1) - % (maxreflevelfact / ipow(reffact, reflevels-1)) == 0) { - Waypoint ("Evolving iteration %d at t=%g", - cgh->cctk_iteration, (double)cgh->cctk_time); - } - - - - // Regrid - { - bool did_regrid = false; - for (int rl=0; rlcctk_iteration-1) % do_every == 0) { - enter_global_mode (cgh, ml); - enter_level_mode (cgh, rl); - - Checkpoint ("Regrid"); - did_regrid |= Regrid (cgh, false, true); - - leave_level_mode (cgh); - leave_global_mode (cgh); - } // if do_every - } // ml - } // for rl - - if (did_regrid) { - for (int rl=0; rl=0; --ml) { - enter_global_mode (cgh, ml); - enter_level_mode (cgh, rl); - - do_global_mode = reflevel==0; - do_meta_mode = do_global_mode && mglevel==mglevels-1; - - Waypoint ("Postregrid at iteration %d time %g%s%s", - cgh->cctk_iteration, (double)cgh->cctk_time, - (do_global_mode ? " (global)" : ""), - (do_meta_mode ? " (meta)" : "")); - - // Postregrid - Checkpoint ("Scheduling POSTREGRID"); - CCTK_ScheduleTraverse ("CCTK_POSTREGRID", cgh, CallFunction); - - leave_level_mode (cgh); - leave_global_mode (cgh); - } // for ml - } // for rl - } // if did_regrid - } - - - - for (int ml=mglevels-1; ml>=0; --ml) { - - bool have_done_global_mode = false; - bool have_done_anything = false; - - for (int rl=0; rlcctk_iteration-1) % do_every == 0) { - enter_global_mode (cgh, ml); - enter_level_mode (cgh, rl); - - do_global_mode = ! have_done_global_mode; - do_meta_mode = do_global_mode && mglevel==mglevels-1; - assert (! (have_done_global_mode && do_global_mode)); - have_done_global_mode |= do_global_mode; - have_done_anything = true; - - // Advance times - for (int m=0; madvance_time (reflevel, mglevel); - } - cgh->cctk_time = (global_time - - delta_time / maxreflevelfact - + delta_time * mglevelfact / reflevelfact); - CycleTimeLevels (cgh); - - Waypoint ("Evolution I at iteration %d time %g%s%s", - cgh->cctk_iteration, (double)cgh->cctk_time, - (do_global_mode ? " (global)" : ""), - (do_meta_mode ? " (meta)" : "")); - - // Checking - CalculateChecksums (cgh, allbutcurrenttime); - Poison (cgh, currenttimebutnotifonly); - - // Evolve - Checkpoint ("Scheduling PRESTEP"); - CCTK_ScheduleTraverse ("CCTK_PRESTEP", cgh, CallFunction); - Checkpoint ("Scheduling EVOL"); - CCTK_ScheduleTraverse ("CCTK_EVOL", cgh, CallFunction); - - // Checking - PoisonCheck (cgh, currenttime); - - leave_level_mode (cgh); - leave_global_mode (cgh); - } // if do_every - } // for rl - - if (have_done_anything) assert (have_done_global_mode); - - } // for ml - - - - for (int ml=mglevels-1; ml>=0; --ml) { - for (int rl=reflevels-1; rl>=0; --rl) { - const int do_every - = ipow(mgfact, ml) * (maxreflevelfact / ipow(reffact, rl)); - if (cgh->cctk_iteration % do_every == 0) { - enter_global_mode (cgh, ml); - enter_level_mode (cgh, rl); - - Waypoint ("Evolution/Restrict at iteration %d time %g", - cgh->cctk_iteration, (double)cgh->cctk_time); - - // Restrict - Restrict (cgh); - - leave_level_mode (cgh); - leave_global_mode (cgh); - } // if do_every - } // for rl - } // for ml - - - - for (int ml=mglevels-1; ml>=0; --ml) { - - bool have_done_global_mode = false; - bool have_done_anything = false; - - for (int rl=0; rlcctk_iteration % do_every == 0) { - enter_global_mode (cgh, ml); - enter_level_mode (cgh, rl); - - int finest_active_reflevel = -1; - { - for (int rl_=0; rl_cctk_iteration % do_every_ == 0) { - finest_active_reflevel = rl_; - } - } - assert (finest_active_reflevel >= 0); - } - do_global_mode = rl == finest_active_reflevel; - do_meta_mode = do_global_mode && mglevel==mglevels-1; - assert (! (have_done_global_mode && do_global_mode)); - have_done_global_mode |= do_global_mode; - have_done_anything = true; - - Waypoint ("Evolution II at iteration %d time %g%s%s", - cgh->cctk_iteration, (double)cgh->cctk_time, - (do_global_mode ? " (global)" : ""), - (do_meta_mode ? " (meta)" : "")); - - Checkpoint ("Scheduling POSTRESTRICT"); - CCTK_ScheduleTraverse ("CCTK_POSTRESTRICT", cgh, CallFunction); - - // Poststep - Checkpoint ("Scheduling POSTSTEP"); - CCTK_ScheduleTraverse ("CCTK_POSTSTEP", cgh, CallFunction); - - // Checking - PoisonCheck (cgh, currenttime); - CalculateChecksums (cgh, currenttime); - - // Checkpoint - Checkpoint ("Scheduling CHECKPOINT"); - CCTK_ScheduleTraverse ("CCTK_CHECKPOINT", cgh, CallFunction); - - // Analysis - Checkpoint ("Scheduling ANALYSIS"); - CCTK_ScheduleTraverse ("CCTK_ANALYSIS", cgh, CallFunction); - - // Output - Checkpoint ("OutputGH"); - CCTK_OutputGH (cgh); - - // Checking - CheckChecksums (cgh, alltimes); - - leave_level_mode (cgh); - leave_global_mode (cgh); - } // if do_every - } // for rl - - if (have_done_anything) assert (have_done_global_mode); - - } // for ml - - + cgh->cctk_time += base_delta_time / maxreflevelfact; + + Checkpoint ("Evolving iteration %d...", cgh->cctk_iteration); + + BEGIN_REFLEVEL_LOOP(cgh) { + if ((cgh->cctk_iteration-1) % (maxreflevelfact/reflevelfact) == 0) { + + // Cycle time levels + CycleTimeLevels (cgh); + + // Advance level times + tt->advance_time (reflevel, mglevel); + for (int group=0; groupadvance_time (reflevel, mglevel); + break; + case CCTK_GF: + break; + default: + abort(); + } + } + + // Checking + CalculateChecksums (cgh, allbutcurrenttime); + Poison (cgh, currenttimebutnotifonly); + + // Evolve + Checkpoint ("%*sScheduling PRESTEP", 2*reflevel, ""); + CCTK_ScheduleTraverse ("CCTK_PRESTEP", cgh, CallFunction); + Checkpoint ("%*sScheduling EVOL", 2*reflevel, ""); + CCTK_ScheduleTraverse ("CCTK_EVOL", cgh, CallFunction); + Checkpoint ("%*sScheduling POSTSTEP", 2*reflevel, ""); + CCTK_ScheduleTraverse ("CCTK_POSTSTEP", cgh, CallFunction); + + // Checking + PoisonCheck (cgh, currenttimebutnotifonly); + + } + } END_REFLEVEL_LOOP(cgh); + + BEGIN_REVERSE_REFLEVEL_LOOP(cgh) { + if (cgh->cctk_iteration % (maxreflevelfact/reflevelfact) == 0) { + + // Restrict + Restrict (cgh); + + // Checking + CalculateChecksums (cgh, currenttime); + + // Checkpoint + Checkpoint ("%*sScheduling CHECKPOINT", 2*reflevel, ""); + CCTK_ScheduleTraverse ("CCTK_CHECKPOINT", cgh, CallFunction); + + // Analysis + Checkpoint ("%*sScheduling ANALYSIS", 2*reflevel, ""); + CCTK_ScheduleTraverse ("CCTK_ANALYSIS", cgh, CallFunction); + + // Output + Checkpoint ("%*sOutputGH", 2*reflevel, ""); + CCTK_OutputGH (cgh); + + // Checking + CheckChecksums (cgh, alltimes); + + } + } END_REVERSE_REFLEVEL_LOOP(cgh); } // main loop - Waypoint ("Done with evolution loop"); + Checkpoint ("done with Evolve."); return 0; } diff --git a/Carpet/Carpet/src/Initialise.cc b/Carpet/Carpet/src/Initialise.cc index fde570065..45c6af58d 100644 --- a/Carpet/Carpet/src/Initialise.cc +++ b/Carpet/Carpet/src/Initialise.cc @@ -1,8 +1,6 @@ #include #include -#include - #include "cctk.h" #include "cctk_Parameters.h" #include "cctki_GHExtensions.h" @@ -11,10 +9,7 @@ #include "carpet.hh" -extern "C" { - static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Initialise.cc,v 1.52 2004/08/19 15:38:20 schnetter Exp $"; - CCTK_FILEVERSION(Carpet_Carpet_Initialise_cc); -} +static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Initialise.cc,v 1.1 2001/07/04 12:29:47 schnetter Exp $"; @@ -34,17 +29,11 @@ namespace Carpet { CCTKi_AddGH (fc, convlev, cgh); // Delay checkpoint until MPI has been initialised - Waypoint ("Starting initialisation"); + Checkpoint ("starting Initialise..."); // Initialise stuff cgh->cctk_iteration = 0; - global_time = cctk_initial_time; - delta_time = 1.0; - cgh->cctk_time = global_time; - cgh->cctk_delta_time = delta_time; - - do_global_mode = true; - do_meta_mode = true; + cgh->cctk_time = cctk_initial_time; // Enable storage and communtication CCTKi_ScheduleGHInit (cgh); @@ -52,503 +41,71 @@ namespace Carpet { // Initialise stuff CCTKi_InitGHExtensions (cgh); + // Check parameters + Checkpoint ("PARAMCHECK"); + CCTK_ScheduleTraverse ("CCTK_PARAMCHECK", cgh, CallFunction); + CCTKi_FinaliseParamWarn(); - - // Output the grid structure - { - // Loop over maps - for (int m=0; mextents, vhh.at(m)->outer_boundaries, - vhh.at(m)->processors); - } // loop over maps - } - - - - BEGIN_MGLEVEL_LOOP(cgh) { - do_global_mode = true; - do_meta_mode = mglevel==mglevels-1; - - // Register coordinates - Checkpoint ("Scheduling CCTK_WRAGH"); - CCTK_ScheduleTraverse ("CCTK_WRAGH", cgh, CallFunction); - - // Check parameters - Checkpoint ("Scheduling PARAMCHECK"); - CCTK_ScheduleTraverse ("CCTK_PARAMCHECK", cgh, CallFunction); - CCTKi_FinaliseParamWarn(); - } END_MGLEVEL_LOOP; - - - - if (fc->recovered) { - // if recovering - - - - for (int rl=0; rlcctk_time = global_time; - - Waypoint ("Recovering I at iteration %d time %g%s%s", - cgh->cctk_iteration, (double)cgh->cctk_time, - (do_global_mode ? " (global)" : ""), - (do_meta_mode ? " (meta)" : "")); - - // Set up the grids - Checkpoint ("Scheduling BASEGRID"); - CCTK_ScheduleTraverse ("CCTK_BASEGRID", cgh, CallFunction); - - // Recover - Checkpoint ("Scheduling RECOVER_VARIABLES"); - CCTK_ScheduleTraverse ("CCTK_RECOVER_VARIABLES", cgh, CallFunction); - - leave_level_mode (cgh); - } END_MGLEVEL_LOOP; - - // Regrid - { - bool did_regrid = false; - { - const int ml=0; - enter_global_mode (cgh, ml); - enter_level_mode (cgh, rl); - - // Regrid - Checkpoint ("Regrid"); - did_regrid |= Regrid (cgh, true, false); - - leave_level_mode (cgh); - leave_global_mode (cgh); - } // ml - - if (did_regrid) { - BEGIN_MGLEVEL_LOOP(cgh) { - enter_level_mode (cgh, rl); - do_global_mode = true; - do_meta_mode = do_global_mode && mglevel==mglevels-1; - - Waypoint ("Postregrid at iteration %d time %g%s%s", - cgh->cctk_iteration, (double)cgh->cctk_time, - (do_global_mode ? " (global)" : ""), - (do_meta_mode ? " (meta)" : "")); - - // Postregrid - Checkpoint ("Scheduling POSTREGRID"); - CCTK_ScheduleTraverse ("CCTK_POSTREGRID", cgh, CallFunction); - - leave_level_mode (cgh); - } END_MGLEVEL_LOOP; - } // if did_regrid - } - - } // for rl - - + BEGIN_REFLEVEL_LOOP(cgh) { - for (int rl=0; rlcctk_iteration, (double)cgh->cctk_time, - (do_global_mode ? " (global)" : ""), - (do_meta_mode ? " (meta)" : "")); - - // Post recover - Checkpoint ("Scheduling POST_RECOVER_VARIABLES"); - CCTK_ScheduleTraverse - ("CCTK_POST_RECOVER_VARIABLES", cgh, CallFunction); - - // Checking - PoisonCheck (cgh, alltimes); - CheckChecksums (cgh, allbutcurrenttime); - - leave_level_mode (cgh); - } END_MGLEVEL_LOOP; - } // for rl + // Checking + Poison (cgh, allbutlasttime); + // Set up the grid + Checkpoint ("%*sScheduling BASEGRID", 2*reflevel, ""); + CCTK_ScheduleTraverse ("CCTK_BASEGRID", cgh, CallFunction); + if (reflevel==0) { + base_delta_time = cgh->cctk_delta_time; + } else { +// assert (abs(cgh->cctk_delta_time - base_delta_time / reflevelfactor) +// < 1e-6 * base_delta_time); + // This circumvents a bug in CactusBase/Time + cgh->cctk_delta_time = base_delta_time / reflevelfact; + } + // Set up the initial data + Checkpoint ("%*sScheduling INITIAL", 2*reflevel, ""); + CCTK_ScheduleTraverse ("CCTK_INITIAL", cgh, CallFunction); + Checkpoint ("%*sScheduling POSTINITIAL", 2*reflevel, ""); + CCTK_ScheduleTraverse ("CCTK_POSTINITIAL", cgh, CallFunction); - } else { - // if not recovering + // Recover + Checkpoint ("%*sScheduling RECOVER_VARIABLES", 2*reflevel, ""); + CCTK_ScheduleTraverse ("CCTK_RECOVER_VARIABLES", cgh, CallFunction); + Checkpoint ("%*sScheduling CPINITIAL", 2*reflevel, ""); + CCTK_ScheduleTraverse ("CCTK_CPINITIAL", cgh, CallFunction); + // Poststep + Checkpoint ("%*sScheduling POSTSTEP", 2*reflevel, ""); + CCTK_ScheduleTraverse ("CCTK_POSTSTEP", cgh, CallFunction); + // Checking + PoisonCheck (cgh, allbutlasttime); - for (int rl=0; rlcctk_time = global_time; - - Waypoint ("Initialisation I at iteration %d time %g%s%s", - cgh->cctk_iteration, (double)cgh->cctk_time, - (do_global_mode ? " (global)" : ""), - (do_meta_mode ? " (meta)" : "")); - - // Checking - Poison (cgh, alltimes); - - // Set up the grids - Checkpoint ("Scheduling BASEGRID"); - CCTK_ScheduleTraverse ("CCTK_BASEGRID", cgh, CallFunction); - - const int num_tl = init_each_timelevel ? 3 : 1; - - // Rewind - for (int m=0; mset_delta - (reflevel, mglevel, - vtt.at(m)->get_delta (reflevel, mglevel)); - FlipTimeLevels (cgh); - for (int tl=0; tladvance_time (reflevel, mglevel); - CycleTimeLevels (cgh); - } - vtt.at(m)->set_delta - (reflevel, mglevel, - vtt.at(m)->get_delta (reflevel, mglevel)); - FlipTimeLevels (cgh); - } - - const bool outer_do_global_mode = do_global_mode; - for (int tl=num_tl-1; tl>=0; --tl) { - do_global_mode = outer_do_global_mode && tl==0; - - // Advance times - for (int m=0; madvance_time (reflevel, mglevel); - } - cgh->cctk_time - = global_time - tl * delta_time * mglevelfact / reflevelfact; - CycleTimeLevels (cgh); - - // Set up the initial data - Checkpoint ("Scheduling INITIAL"); - CCTK_ScheduleTraverse ("CCTK_INITIAL", cgh, CallFunction); - - } // for tl - do_global_mode = outer_do_global_mode; - - // Checking - PoisonCheck (cgh, currenttime); - - leave_level_mode (cgh); - } END_MGLEVEL_LOOP; - - // Regrid - { - bool did_regrid = false; - { - const int ml=0; - enter_global_mode (cgh, ml); - enter_level_mode (cgh, rl); - - // Regrid - Checkpoint ("Regrid"); - did_regrid |= Regrid (cgh, false, prolongate_initial_data); - - leave_level_mode (cgh); - leave_global_mode (cgh); - } // ml - - if (did_regrid) { - BEGIN_MGLEVEL_LOOP(cgh) { - enter_level_mode (cgh, rl); - do_global_mode = true; - do_meta_mode = do_global_mode && mglevel==mglevels-1; - - Waypoint ("Postregrid at iteration %d time %g%s%s", - cgh->cctk_iteration, (double)cgh->cctk_time, - (do_global_mode ? " (global)" : ""), - (do_meta_mode ? " (meta)" : "")); - - // Postregrid - Checkpoint ("Scheduling POSTREGRID"); - CCTK_ScheduleTraverse ("CCTK_POSTREGRID", cgh, CallFunction); - - leave_level_mode (cgh); - } END_MGLEVEL_LOOP; - } // if did_regrid - } - - } // for rl - - - - for (int rl=reflevels-1; rl>=0; --rl) { - BEGIN_MGLEVEL_LOOP(cgh) { - enter_level_mode (cgh, rl); - - Waypoint ("Initialisation/Restrict at iteration %d time %g", - cgh->cctk_iteration, (double)cgh->cctk_time); - - // Restrict - Restrict (cgh); - - leave_level_mode (cgh); - } END_MGLEVEL_LOOP; - } // for rl - - - - for (int rl=0; rlcctk_iteration, (double)cgh->cctk_time, - (do_global_mode ? " (global)" : ""), - (do_meta_mode ? " (meta)" : "")); - - Checkpoint ("Scheduling POSTRESTRICTINITIAL"); - CCTK_ScheduleTraverse - ("CCTK_POSTRESTRICTINITIAL", cgh, CallFunction); - - // Postinitial - Checkpoint ("Scheduling POSTINITIAL"); - CCTK_ScheduleTraverse ("CCTK_POSTINITIAL", cgh, CallFunction); - - // Poststep - Checkpoint ("Scheduling POSTSTEP"); - CCTK_ScheduleTraverse ("CCTK_POSTSTEP", cgh, CallFunction); - - // Checking - PoisonCheck (cgh, alltimes); - CheckChecksums (cgh, allbutcurrenttime); - - leave_level_mode (cgh); - } END_MGLEVEL_LOOP; - } // for rl - - + } END_REFLEVEL_LOOP(cgh); - if (init_3_timelevels) { - // Use Scott Hawley's algorithm for getting two extra - // timelevels of data - Waypoint ("Initialising three timelevels"); + BEGIN_REVERSE_REFLEVEL_LOOP(cgh) { - for (int rl=0; rladvance_time (reflevel, mglevel); - } - cgh->cctk_time - = global_time + delta_time * mglevelfact / reflevelfact; - CycleTimeLevels (cgh); - - Waypoint ("Initialisation 3TL evolution I (a) (forwards) at iteration %d time %g%s%s", - cgh->cctk_iteration, (double)cgh->cctk_time, - (do_global_mode ? " (global)" : ""), - (do_meta_mode ? " (meta)" : "")); - - // Checking - CalculateChecksums (cgh, allbutcurrenttime); - Poison (cgh, currenttimebutnotifonly); - - // Evolve forward - Checkpoint ("Scheduling PRESTEP"); - CCTK_ScheduleTraverse ("CCTK_PRESTEP", cgh, CallFunction); - Checkpoint ("Scheduling EVOL"); - CCTK_ScheduleTraverse ("CCTK_EVOL", cgh, CallFunction); - - // Checking - PoisonCheck (cgh, currenttime); - - leave_level_mode (cgh); - } END_MGLEVEL_LOOP; - } // for rl + // Restrict + Restrict (cgh); - delta_time *= -1; - for (int rl=0; rlcctk_time - = global_time + delta_time * mglevelfact / reflevelfact; - - leave_level_mode (cgh); - } END_MGLEVEL_LOOP; - } // for rl + // Checking + CalculateChecksums (cgh, allbutlasttime); - for (int rl=0; rlcctk_iteration, (double)cgh->cctk_time, - (do_global_mode ? " (global)" : ""), - (do_meta_mode ? " (meta)" : "")); - - // Checking - CalculateChecksums (cgh, allbutcurrenttime); - Poison (cgh, currenttimebutnotifonly); - - // Evolve backward - Checkpoint ("Scheduling PRESTEP"); - CCTK_ScheduleTraverse ("CCTK_PRESTEP", cgh, CallFunction); - Checkpoint ("Scheduling EVOL"); - CCTK_ScheduleTraverse ("CCTK_EVOL", cgh, CallFunction); - - // Checking - PoisonCheck (cgh, alltimes); - - leave_level_mode (cgh); - } END_MGLEVEL_LOOP; - } // for rl + // Analysis + Checkpoint ("%*sScheduling ANALYSIS", 2*reflevel, ""); + CCTK_ScheduleTraverse ("CCTK_ANALYSIS", cgh, CallFunction); - Waypoint ("Hourglass structure in place"); + // Output + Checkpoint ("%*sOutputGH", 2*reflevel, ""); + CCTK_OutputGH (cgh); - // Evolve each level "backwards" one more timestep - // Starting with the finest level and proceeding to the coarsest - for (int rl=reflevels-1; rl>=0; --rl) { - BEGIN_MGLEVEL_LOOP(cgh) { - enter_level_mode (cgh, rl); - do_global_mode = reflevel==0; - do_meta_mode = do_global_mode && mglevel==mglevels-1; - - Waypoint ("Initialisation 3TL evolution II (b) (backwards) at iteration %d time %g%s%s", - cgh->cctk_iteration, (double)cgh->cctk_time, - (do_global_mode ? " (global)" : ""), - (do_meta_mode ? " (meta)" : "")); - - // Restrict - Restrict (cgh); - - Checkpoint ("Scheduling POSTRESTRICT"); - CCTK_ScheduleTraverse ("CCTK_POSTRESTRICT", cgh, CallFunction); - - // Poststep - Checkpoint ("Scheduling POSTSTEP"); - CCTK_ScheduleTraverse ("CCTK_POSTSTEP", cgh, CallFunction); - - // Checking - PoisonCheck (cgh, alltimes); - - // Advance times - for (int m=0; madvance_time (reflevel, mglevel); - } - cgh->cctk_time - = global_time + 2 * delta_time * mglevelfact / reflevelfact; - CycleTimeLevels (cgh); - - Waypoint ("Initialisation 3TL evolution I (c) (backwards) at iteration %d time %g%s%s", - cgh->cctk_iteration, (double)cgh->cctk_time, - (do_global_mode ? " (global)" : ""), - (do_meta_mode ? " (meta)" : "")); - - // Checking - CalculateChecksums (cgh, allbutcurrenttime); - Poison (cgh, currenttimebutnotifonly); - - // Evolve backward - Checkpoint ("Scheduling PRESTEP"); - CCTK_ScheduleTraverse ("CCTK_PRESTEP", cgh, CallFunction); - Checkpoint ("Scheduling EVOL"); - CCTK_ScheduleTraverse ("CCTK_EVOL", cgh, CallFunction); - Checkpoint ("Scheduling POSTSTEP"); - CCTK_ScheduleTraverse ("CCTK_POSTSTEP", cgh, CallFunction); - - // Checking - PoisonCheck (cgh, alltimes); - - leave_level_mode (cgh); - } END_MGLEVEL_LOOP; - } // for rl + // Checking + CheckChecksums (cgh, allbutlasttime); - delta_time *= -1; - for (int rl=0; rlset_delta - (reflevel, mglevel, - vtt.at(m)->get_delta (reflevel, mglevel)); - vtt.at(m)->advance_time (reflevel, mglevel); - vtt.at(m)->advance_time (reflevel, mglevel); - vtt.at(m)->set_delta - (reflevel, mglevel, - vtt.at(m)->get_delta (reflevel, mglevel)); - } - cgh->cctk_time = global_time; - - leave_level_mode (cgh); - } END_MGLEVEL_LOOP; - } // for rl - - Waypoint ("Finished initialising three timelevels"); - - } // if init_3_timelevels - - } // if not recovering - + } END_REVERSE_REFLEVEL_LOOP(cgh); - - for (int rl=0; rlcctk_iteration, (double)cgh->cctk_time, - (do_global_mode ? " (global)" : ""), - (do_meta_mode ? " (meta)" : "")); - - // Checkpoint - Checkpoint ("Scheduling CPINITIAL"); - CCTK_ScheduleTraverse ("CCTK_CPINITIAL", cgh, CallFunction); - - // Analysis - Checkpoint ("Scheduling ANALYSIS"); - CCTK_ScheduleTraverse ("CCTK_ANALYSIS", cgh, CallFunction); - - // Output - Checkpoint ("OutputGH"); - CCTK_OutputGH (cgh); - - // Checking - PoisonCheck (cgh, alltimes); - CheckChecksums (cgh, allbutcurrenttime); - - leave_level_mode (cgh); - } END_MGLEVEL_LOOP; - } // for rl - - - - Waypoint ("Done with initialisation"); + Checkpoint ("done with Initialise."); return 0; } diff --git a/Carpet/Carpet/src/Poison.cc b/Carpet/Carpet/src/Poison.cc index 7f1a00f2e..4b0194bf3 100644 --- a/Carpet/Carpet/src/Poison.cc +++ b/Carpet/Carpet/src/Poison.cc @@ -7,10 +7,7 @@ #include "carpet.hh" -extern "C" { - static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Poison.cc,v 1.17 2004/03/23 19:30:14 schnetter Exp $"; - CCTK_FILEVERSION(Carpet_Carpet_Poison_cc); -} +static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Poison.cc,v 1.1 2001/07/04 12:29:47 schnetter Exp $"; @@ -20,7 +17,7 @@ namespace Carpet { - void Poison (const cGH* cgh, const checktimes where) + void Poison (cGH* cgh, const checktimes where) { DECLARE_CCTK_PARAMETERS; @@ -35,7 +32,7 @@ namespace Carpet { - void PoisonGroup (const cGH* cgh, const int group, const checktimes where) + void PoisonGroup (cGH* cgh, const int group, const checktimes where) { DECLARE_CCTK_PARAMETERS; @@ -43,150 +40,170 @@ namespace Carpet { if (! poison_new_timelevels) return; + Checkpoint ("%*sPoisonGroup %s", 2*reflevel, "", CCTK_GroupName(group)); + if (! CCTK_QueryGroupStorageI(cgh, group)) { - char * const groupname = CCTK_GroupName(group); CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING, "Cannot poison group \"%s\" because it has no storage", - groupname); - free (groupname); + CCTK_GroupName(group)); return; } - const int n0 = CCTK_FirstVarIndexI(group); - assert (n0>=0); - const int nvar = CCTK_NumVarsInGroupI(group); - const int sz = CCTK_VarTypeSize(CCTK_VarTypeI(n0)); - assert (sz>0); - - const int num_tl = CCTK_NumTimeLevelsFromVarI(n0); - assert (num_tl>0); - const int min_tl = mintl(where, num_tl); - const int max_tl = maxtl(where, num_tl); - - if (min_tl <= max_tl) { + for (int var=0; var0); - const int grouptype = CCTK_GroupTypeI(group); + const int num_tl = CCTK_NumTimeLevelsFromVarI(n); + assert (num_tl>0); + const int min_tl = mintl(where, num_tl); + const int max_tl = maxtl(where, num_tl); - BEGIN_MAP_LOOP(cgh, grouptype) { - BEGIN_LOCAL_COMPONENT_LOOP(cgh, grouptype) { - - ivect size(1); - const int gpdim = groupdata.at(group).info.dim; - for (int d=0; ddata[n][tl], poison_value, np*sz); - } // for tl - } // for var - - } END_LOCAL_COMPONENT_LOOP; - } END_MAP_LOOP; + for (int tl=min_tl; tl<=max_tl; ++tl) { + + switch (CCTK_GroupTypeFromVarI(n)) { + case CCTK_SCALAR: { + assert (group<(int)scdata.size()); + assert (var<(int)scdata[group].data.size()); + memset (cgh->data[n][tl], poison_value, sz); + break; + } + case CCTK_ARRAY: + case CCTK_GF: { + BEGIN_COMPONENT_LOOP(cgh) { + if (hh->is_local(reflevel,component)) { + int np = 1; + const int gpdim = CCTK_GroupDimI(group); + for (int d=0; ddata[n][tl], poison_value, np*sz); + } + } END_COMPONENT_LOOP(cgh); + break; + } + default: + abort(); + } + + } // for tl - } // if tl + } // for var } - void PoisonCheck (const cGH* cgh, const checktimes where) + void PoisonCheck (cGH* cgh, const checktimes where) { DECLARE_CCTK_PARAMETERS; if (! check_for_poison) return; - Checkpoint ("PoisonCheck"); + Checkpoint ("%*sPoisonCheck", 2*reflevel, ""); for (int group=0; group=0); - const int nvar = CCTK_NumVarsInGroupI(group); - const int tp = CCTK_VarTypeI(n0); - const int gpdim = groupdata.at(group).info.dim; - - const int num_tl = CCTK_NumTimeLevelsFromVarI(n0); - assert (num_tl>0); - const int min_tl = mintl(where, num_tl); - const int max_tl = maxtl(where, num_tl); - - BEGIN_MAP_LOOP(cgh, grouptype) { - BEGIN_LOCAL_COMPONENT_LOOP(cgh, grouptype) { - - ivect size(1); - for (int d=0; ddata[n][tl]; - int numpoison=0; - for (int k=0; k0); + const int min_tl = mintl(where, num_tl); + const int max_tl = maxtl(where, num_tl); + + for (int tl=min_tl; tl<=max_tl; ++tl) { + + switch (CCTK_GroupTypeFromVarI(n)) { + + case CCTK_SCALAR: { + bool poisoned=false; + switch (CCTK_VarTypeI(n)) { +#define TYPECASE(N,T) \ + case N: { \ + T worm; \ + memset (&worm, poison_value, sizeof(worm)); \ + poisoned = *(const T*)cgh->data[n][tl] == worm; \ + break; \ + } +#include "typecase" +#undef TYPECASE + default: + UnsupportedVarType(n); + } + if (poisoned) { + char* fullname = CCTK_FullName(n); + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "The variable \"%s\" contains poison in timelevel %d", + fullname, tl); + free (fullname); + } + break; + } + + case CCTK_ARRAY: + case CCTK_GF: { + BEGIN_COMPONENT_LOOP(cgh) { + if (hh->is_local(reflevel,component)) { + vect size(1); + const int gpdim = CCTK_GroupDimI(group); + for (int d=0; ddata[n][tl]; + int numpoison=0; + for (int k=0; kmax_poison_locations) { - char* fullname = CCTK_FullName(n); - CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, - "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", - tl, component, map, reflevel, fullname, numpoison, np); - free (fullname); - } else if (numpoison>0) { - CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, - "Found poison at %d of %d locations", - numpoison, np); - } - - } // for tl - } // for var - } END_LOCAL_COMPONENT_LOOP; - } END_MAP_LOOP; - + default: + UnsupportedVarType(n); + } + if (poisoned) { + ++numpoison; + if (numpoison<=10) { + char* fullname = CCTK_FullName(n); + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "The variable \"%s\" contains poison at [%d,%d,%d] in timelevel %d", + fullname, i,j,k, tl); + free (fullname); + } + } // if poisoned + } // for i + } // for j + } // for k + if (numpoison>10) { + char* fullname = CCTK_FullName(n); + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "The variable \"%s\" contains poison at %d locations in timelevel %d; not all locations were printed.", + fullname, numpoison, tl); + free (fullname); + } + } // if is local + } END_COMPONENT_LOOP(cgh); + break; + } + + default: + abort(); + } + + } // for tl + + } // for var } // if has storage } // for group } diff --git a/Carpet/Carpet/src/Recompose.cc b/Carpet/Carpet/src/Recompose.cc index e86761021..e1c9afba8 100644 --- a/Carpet/Carpet/src/Recompose.cc +++ b/Carpet/Carpet/src/Recompose.cc @@ -1,37 +1,18 @@ #include -#include #include -#include -#include -#include -#include -#include -#include -#include -#include #include #include "cctk.h" #include "cctk_Parameters.h" -#include "CactusBase/IOUtil/src/ioGH.h" - -#include "bbox.hh" -#include "bboxset.hh" -#include "defs.hh" -#include "gh.hh" -#include "vect.hh" +#include "Carpet/CarpetLib/src/bbox.hh" +#include "Carpet/CarpetLib/src/gh.hh" +#include "Carpet/CarpetLib/src/vect.hh" #include "carpet.hh" -#include "modes.hh" - -extern "C" { - static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Recompose.cc,v 1.72 2004/08/19 15:38:20 schnetter Exp $"; - CCTK_FILEVERSION(Carpet_Carpet_Recompose_cc); -} -#define DEBUG false // false or true +static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Recompose.cc,v 1.1 2001/07/04 12:29:47 schnetter Exp $"; @@ -41,821 +22,94 @@ namespace Carpet { - // Reduction operator - template - static typename func::result_type - reduce (iter const first, iter const last, - typename func::result_type const & init) - { - typename func::result_type res (init); - for (iter it (first); it != last; ++it) { - res = func::operator() (res, *it); - } - return res; - } - + static void Recompose_gh (cGH* cgh, gh* hh); - static void SplitRegions_Automatic_Recursively (bvect const & dims, - int const nprocs, - rvect const rshape, - ibbox const & bb, - bbvect const & ob, - int const & p, - vector & bbs, - vector & obs, - vector & ps); - static void SplitRegions_AsSpecified (const cGH* cgh, - vector& bbs, - vector& obs, - vector& ps); - - - void CheckRegions (const gh::rexts & bbsss, - const gh::rbnds & obss, - const gh::rprocs& pss) - { - // At least one level - if (bbsss.size() == 0) { - CCTK_WARN (0, "I cannot set up a grid hierarchy with zero refinement levels."); - } - assert (bbsss.size() > 0); - // At most maxreflevels levels - if ((int)bbsss.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)bbsss.size()); - } - assert ((int)bbsss.size() <= maxreflevels); - for (int rl=0; rl<(int)bbsss.size(); ++rl) { - // No empty levels - assert (bbsss.at(rl).size() > 0); - for (int c=0; c<(int)bbsss.at(rl).size(); ++c) { - // At least one multigrid level - assert (bbsss.at(rl).at(c).size() > 0); - for (int ml=0; ml<(int)bbsss.at(rl).at(c).size(); ++ml) { - // Check sizes - // Do allow processors with zero grid points -// assert (all(bbsss.at(rl).at(c).at(ml).lower() <= bbsssi.at(rl).at(c).at(ml).upper())); - // Check strides - const int str = ipow(reffact, maxreflevels-rl-1) * ipow(mgfact, ml); - assert (all(bbsss.at(rl).at(c).at(ml).stride() == str)); - // Check alignments - assert (all(bbsss.at(rl).at(c).at(ml).lower() % str == 0)); - assert (all(bbsss.at(rl).at(c).at(ml).upper() % str == 0)); - } - } - } - - assert (pss.size() == bbsss.size()); - assert (obss.size() == bbsss.size()); - for (int rl=0; rl<(int)bbsss.size(); ++rl) { - assert (obss.at(rl).size() == bbsss.at(rl).size()); - assert (pss.at(rl).size() == bbsss.at(rl).size()); - } - - } - - - - bool Regrid (const cGH* cgh, const bool force_recompose, const bool do_init) - { - assert (is_level_mode()); - - if (! CCTK_IsFunctionAliased ("Carpet_Regrid")) { - static bool didtell = false; - if (!didtell) { - CCTK_WARN (1, "No regridding routine has been provided. There will be no regridding. Maybe you forgot to activate a regridding thorn?"); - didtell = true; - } - return false; - } - - bool did_change = false; - BEGIN_MAP_LOOP(cgh, CCTK_GF) { - - gh::rexts bbsss = vhh.at(map)->extents; - gh::rbnds obss = vhh.at(map)->outer_boundaries; - gh::rprocs pss = vhh.at(map)->processors; - - // Check whether to recompose - CCTK_INT const do_recompose - = Carpet_Regrid (cgh, &bbsss, &obss, &pss, force_recompose); - assert (do_recompose >= 0); - did_change = did_change || do_recompose; - - if (do_recompose) { - - CCTK_INFO ("Recomposing the grid hierarchy"); - - // Check the regions - CheckRegions (bbsss, obss, pss); - // TODO: check also that the current and all coarser levels - // did not change - - // Write grid structure to file - OutputGridStructure (cgh, map, bbsss, obss, pss); - - // Recompose - vhh.at(map)->recompose (bbsss, obss, pss, do_init); - - OutputGrids (cgh, map, *vhh.at(map)); - - } - - } END_MAP_LOOP; - - // Calculate new number of levels - reflevels = vhh.at(0)->reflevels(); - for (int m=0; mreflevels() == reflevels); - } - - // One cannot switch off the current level - assert (reflevels > reflevel); - - return did_change; - } - - - - void OutputGrids (const cGH* cgh, const int m, const gh& hh) - { - CCTK_INFO ("New grid structure (grid points):"); - cout << " Refinement level " << reflevel << ", map " << map << endl; - for (int rl=0; rl::rexts & bbsss, - const gh::rbnds & obss, - const gh::rprocs& pss) + void Recompose (cGH* cgh) { DECLARE_CCTK_PARAMETERS; - // Output only on the root processor - if (CCTK_MyProc(cgh) != 0) return; - - // Output only if output is desired - if (strcmp(grid_structure_filename, "") == 0) return; - - // Get grid hierarchy extention from IOUtil - const ioGH * const iogh = (const ioGH *)CCTK_GHExtension (cgh, "IO"); - - // Output only if IO exists and has been initialised - if (! iogh) return; - - assert (iogh); - - // Create the output directory - CCTK_CreateDirectory (0755, out_dir); - - ostringstream filenamebuf; - filenamebuf << out_dir << "/" << grid_structure_filename; - // we need a persistent temporary here - string filenamestr = filenamebuf.str(); - const char * filename = filenamestr.c_str(); - - ofstream file; - - static bool do_truncate = true; - - if (do_truncate) { - do_truncate = false; - struct stat fileinfo; - if (! iogh->recovered - || stat(filename, &fileinfo)!=0) { - file.open (filename, ios::out | ios::trunc); - assert (file.good()); - file << "# grid structure" << endl - << "# format: map reflevel component mglevel processor bounding-box is-outer-boundary" << endl; - assert (file.good()); - } - } - if (! file.is_open()) { - file.open (filename, ios::app); - assert (file.good()); - } - - file << "iteration " << cgh->cctk_iteration << endl; - file << "maps " << maps << endl; - file << m << " reflevels " << bbsss.size() << endl; - for (int rl=0; rl<(int)bbsss.size(); ++rl) { - file << m << " " << rl << " components " << bbsss.at(rl).size() << endl; - for (int c=0; c<(int)bbsss.at(rl).size(); ++c) { - file << m << " " << rl << " " << c << " mglevels " << bbsss.at(rl).at(c).size() << endl; - for (int ml=0; ml<(int)bbsss.at(rl).at(c).size(); ++ml) { - file << m << " " << rl << " " << c << " " << ml << " " << pss.at(rl).at(c) << " " << bbsss.at(rl).at(c).at(ml) << obss.at(rl).at(c) << endl; - } - } - } - file << endl; - - file.close(); - assert (file.good()); + assert (component == -1); + Checkpoint ("%*sRecompose", 2*reflevel, ""); + + Recompose_gh (cgh, hh); + + for (int group=0; group& bbs, vector& obs, - vector& ps) + static void Recompose_gh (cGH* cgh, gh* hh) { DECLARE_CCTK_PARAMETERS; - if (CCTK_EQUALS (processor_topology, "along-z")) { - SplitRegions_AlongZ (cgh, bbs, obs, ps); - } else if (CCTK_EQUALS (processor_topology, "along-dir")) { - SplitRegions_AlongDir (cgh, bbs, obs, ps, split_direction); - } else if (CCTK_EQUALS (processor_topology, "automatic")) { - SplitRegions_Automatic (cgh, bbs, obs, ps); - } else if (CCTK_EQUALS (processor_topology, "manual")) { - SplitRegions_AsSpecified (cgh, bbs, obs, ps); - } else { - assert (0); - } - } - - - - void SplitRegions_AlongZ (const cGH* cgh, vector& bbs, - vector& obs, vector& ps) - { - SplitRegions_AlongDir (cgh, bbs, obs, ps, 2); - } - - - - void SplitRegions_AlongDir (const cGH* cgh, vector& bbs, - vector& obs, vector& ps, - const int dir) - { - // Something to do? - if (bbs.size() == 0) { - ps.resize(0); - return; - } - - const int nprocs = CCTK_nProcs(cgh); - - if (nprocs==1) { - ps.resize(1); - ps.at(0) = 0; - return; - } - - assert (bbs.size() == 1); - - assert (dir>=0 && dir rub[dir]) clb[dir] = rub[dir]; - if (cub[dir] > rub[dir]) cub[dir] = rub[dir]; - assert (clb[dir] <= cub[dir]); - assert (cub[dir] <= rub[dir]); - bbs.at(c) = ibbox(clb, cub-cstr, cstr); - obs.at(c) = obnd; - ps.at(c) = c; - if (c>0) obs.at(c)[dir][0] = false; - if (c & bbs, - vector & obs, - vector & ps) - { - if (DEBUG) cout << "SRAR enter" << endl; - // check preconditions - assert (nprocs >= 1); - - // are we done? - if (all(dims)) { - if (DEBUG) cout << "SRAR bottom" << endl; - - // check precondition - assert (nprocs == 1); - - // return arguments - bbs.assign (1, bb); - obs.assign (1, ob); - ps.assign (1, p); - - // return - if (DEBUG) cout << "SRAR exit" << endl; - return; - } - - // choose a direction - int mydim = -1; - CCTK_REAL mysize = 0; - int alldims = 0; - CCTK_REAL allsizes = 1; - for (int d=0; d= mysize) { - mydim = d; - mysize = rshape[d]; - } - } - } - assert (mydim>=0 && mydim=0); - if (DEBUG) cout << "SRAR mydim " << mydim << endl; - if (DEBUG) cout << "SRAR mysize " << mysize << endl; - - if (mysize == 0) { - // the bbox is empty - if (DEBUG) cout << "SRAR empty" << endl; - - // create the bboxes - bbs.clear(); - obs.clear(); - ps.clear(); - bbs.reserve(nprocs); - obs.reserve(nprocs); - ps.reserve(nprocs); - - // create a new bbox - assert (bb.empty()); - bbvect const newob (false); - ibbox const newbb (bb); - int const newp (p); - if (DEBUG) cout << "SRAR " << mydim << " newbb " << newbb << endl; - if (DEBUG) cout << "SRAR " << mydim << " newob " << newob << endl; - if (DEBUG) cout << "SRAR " << mydim << " newp " << newp << endl; - - // store - bbs.insert (bbs.end(), nprocs, newbb); - obs.insert (obs.end(), nprocs, newob); - for (int pp=0; pp mynprocs(nslices); - int const mynprocs_base = nprocs / nslices; - int const mynprocs_left = nprocs - nslices * mynprocs_base; - for (int n=0; n myslice(nslices); - int slice_left = ((bb.upper() - bb.lower()) / bb.stride())[mydim] + 1; - int nprocs_left = nprocs; - for (int n=0; n= 0); - slice_left -= myslice.at(n); - nprocs_left -= mynprocs.at(n); - } - assert (slice_left == 0); - assert (nprocs_left == 0); - if (DEBUG) cout << "SRAR " << mydim << " myslice " << myslice << endl; - - // create the bboxes and recurse - if (DEBUG) cout << "SRAR " << mydim << ": create bboxes" << endl; - bbs.clear(); - obs.clear(); - ps.clear(); - bbs.reserve(nprocs); - obs.reserve(nprocs); - ps.reserve(nprocs); - ivect last_up; - for (int n=0; n 0) { - lo[mydim] = last_up[mydim] + str[mydim]; - newob[mydim][0] = false; - } - if (n < nslices-1) { - up[mydim] = lo[mydim] + (myslice.at(n)-1) * str[mydim]; - newob[mydim][1] = false; - last_up = up; - } - ibbox newbb(lo, up, str); - int newp(p + n * mynprocs_base + (n < mynprocs_left ? n : mynprocs_left)); - if (DEBUG) cout << "SRAR " << mydim << " newbb " << newbb << endl; - if (DEBUG) cout << "SRAR " << mydim << " newob " << newob << endl; - if (DEBUG) cout << "SRAR " << mydim << " newp " << newp << endl; - - // recurse - vector newbbs; - vector newobs; - vector newps; - SplitRegions_Automatic_Recursively - (newdims, mynprocs.at(n), rshape, - newbb, newob, newp, newbbs, newobs, newps); - if (DEBUG) cout << "SRAR " << mydim << " newbbs " << newbbs << endl; - if (DEBUG) cout << "SRAR " << mydim << " newobs " << newobs << endl; - if (DEBUG) cout << "SRAR " << mydim << " newps " << newps << endl; - - // store - assert ((int)newbbs.size() == mynprocs.at(n)); - assert ((int)newobs.size() == mynprocs.at(n)); - assert ((int)newps.size() == mynprocs.at(n)); - bbs.insert (bbs.end(), newbbs.begin(), newbbs.end()); - obs.insert (obs.end(), newobs.begin(), newobs.end()); - ps.insert (ps.end(), newps.begin(), newps.end()); - } - - // check postconditions - assert ((int)bbs.size() == nprocs); - assert ((int)obs.size() == nprocs); - assert ((int)ps.size() == nprocs); - for (int n=0; n<(int)ps.size(); ++n) { - assert ((int)ps.at(n) == p+n); - } - if (DEBUG) cout << "SRAR exit" << endl; - } - - - - void SplitRegions_Automatic (const cGH* cgh, vector& bbs, - vector& obs, vector& ps) - { - if (DEBUG) cout << "SRA enter" << endl; - // Something to do? - if (bbs.size() == 0) { - ps.resize(0); - return; - } - - const int nprocs = CCTK_nProcs(cgh); - if (DEBUG) cout << "SRA nprocs " << nprocs << endl; - - // nslices: number of disjoint bboxes - int const nslices = bbs.size(); - if (DEBUG) cout << "SRA nslices " << nslices << endl; - // ncomps: number of components per processor - int const ncomps = (nslices + nprocs - 1) / nprocs; - if (DEBUG) cout << "SRA ncomps " << ncomps << endl; - assert (ncomps > 0); - vector mysize(nslices); - for (int c=0; c mynprocs(nslices); - { - if (DEBUG) cout << "SRA: distributing processors to slices" << endl; - int ncomps_left = nprocs * ncomps; - for (int c=0; c 0) { - if (DEBUG) cout << "SRA ncomps_left " << ncomps_left << endl; - int maxc = -1; - CCTK_REAL maxratio = -1; - for (int c=0; c maxratio) { maxc=c; maxratio=ratio; } - } - assert (maxc>=0 && maxc allbbs; - vector allobs; - vector allps; - - if (DEBUG) cout << "SRA: splitting regions" << endl; - for (int c=0, p=0; c thebbs; - vector theobs; - vector theps; - - SplitRegions_Automatic_Recursively - (dims, mynprocs.at(c), rshape, bb, ob, p, thebbs, theobs, theps); - if (DEBUG) cout << "SRA thebbs " << thebbs << endl; - if (DEBUG) cout << "SRA theobs " << theobs << endl; - if (DEBUG) cout << "SRA theps " << theps << endl; - - allbbs.insert(allbbs.end(), thebbs.begin(), thebbs.end()); - allobs.insert(allobs.end(), theobs.begin(), theobs.end()); - allps.insert(allps.end(), theps.begin(), theps.end()); - - } // for c - - bbs = allbbs; - obs = allobs; - ps = allps; - for (int n=0; n<(int)ps.size(); ++n) { - ps.at(n) /= ncomps; - assert (ps.at(n) >= 0 && ps.at(n) < nprocs); - } - - if (DEBUG) cout << "SRA exit" << endl; - } - - - - static void SplitRegions_AsSpecified (const cGH* cgh, - vector& bbs, - vector& obs, - vector& ps) - { - DECLARE_CCTK_PARAMETERS; - - // Something to do? - if (bbs.size() == 0) { - ps.resize(0); - return; - } - - const int nprocs = CCTK_nProcs(cgh); - - assert (bbs.size() == 1); - - const ivect rstr = bbs.at(0).stride(); - const ivect rlb = bbs.at(0).lower(); - const ivect rub = bbs.at(0).upper() + rstr; - const bbvect obnd = obs.at(0); - - const ivect nprocs_dir - (processor_topology_3d_x, processor_topology_3d_y, - processor_topology_3d_z); - assert (all (nprocs_dir > 0)); - if (prod(nprocs_dir) != nprocs) { - CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, - "The specified processor topology [%d,%d,%d] requires %d processors, but there are %d processors", nprocs_dir[0], nprocs_dir[1], nprocs_dir[2], prod(nprocs_dir), nprocs); - } - assert (prod(nprocs_dir) == nprocs); - - bbs.resize(nprocs); - obs.resize(nprocs); - ps.resize(nprocs); - const ivect cstr = rstr; - const ivect glonp = (rub - rlb) / cstr; -// const ivect locnp = (glonp + nprocs_dir - 1) / nprocs_dir; - const ivect locnp = glonp / nprocs_dir; - const ivect rem = glonp % nprocs_dir; - const ivect step = locnp * cstr; - assert (dim==3); - for (int k=0; k= 0)); - assert (all (clb <= cub)); - assert (all (cub <= rub)); - assert (all (! (ipos==0) || clb==rlb)); - assert (all (! (ipos==nprocs_dir-1) || cub==rub)); - bbs.at(c) = ibbox(clb, cub-cstr, cstr); - obs.at(c) = obnd; - ps.at(c) = c; - if (i>0) obs.at(c)[0][0] = false; - if (j>0) obs.at(c)[1][0] = false; - if (k>0) obs.at(c)[2][0] = false; - if (i& bbs) - { - bbs.resize (mglevels); - bbs.at(0) = bb; - // boundary offsets - jjvect nboundaryzones, is_internal, is_staggered, shiftout; - const int ierr = GetBoundarySpecification - (2*dim, &nboundaryzones[0][0], &is_internal[0][0], - &is_staggered[0][0], &shiftout[0][0]); - assert (!ierr); - // (distance in grid points between the exterior and the physical boundary) - iivect offset; - for (int d=0; d bases(mglevels); - bases.at(0) = base; - for (int ml=1; ml const & bbs, - vector const & obs, - vector >& bbss) - { - assert (bbs.size() == obs.size()); - ibbox base; - for (int c=0; c<(int)bbs.size(); ++c) { - base = base.expanded_containing(bbs.at(c)); - } - bbss.resize(bbs.size()); - for (int c=0; c<(int)bbs.size(); ++c) { - MakeMultigridBoxes (cgh, base, bbs.at(c), obs.at(c), bbss.at(c)); - } + const int nprocs = CCTK_nProcs(cgh); + const int reflevels = max_refinement_levels; // arbitrary value + const int mglevels = 1; // arbitrary value + vector > > bbss(reflevels); + // note: what this routine calls "ub" is "ub+str" elsewhere + vect rstr = hh->baseextent.stride(); + vect rlb = hh->baseextent.lower(); + vect rub = hh->baseextent.upper() + rstr; + for (int rl=0; rl0) { + // save old values + const vect oldrlb = rlb; + const vect oldrub = rub; + // calculate extent and centre + const vect rextent = rub - rlb + rstr; + const vect rcentre = rlb + (rextent / 2 / rstr) * rstr; + // calculate new extent + assert (all(rextent % hh->reffact == 0)); + const vect newrextent = rextent / hh->reffact; + // refined boxes have smaller stride + assert (all(rstr%hh->reffact == 0)); + rstr /= hh->reffact; + // refine (arbitrarily) around the center only + rlb = rcentre - (newrextent/2 / rstr) * rstr; + rub = rlb + newrextent - rstr; +// // refine (arbitrarily) around the lower boundary only +// rlb = rlb; +// rub = rlb + newrextent - rstr; + assert (all(rlb >= oldrlb && rub <= oldrub)); + } + vector > bbs(nprocs); + for (int c=0; c cstr = rstr; + vect clb = rlb; + vect cub = rub; + const int glonpz = (rub[dim-1] - rlb[dim-1]) / cstr[dim-1]; + const int locnpz = (glonpz + nprocs - 1) / nprocs; + const int zstep = locnpz * cstr[dim-1]; + clb[dim-1] = rlb[dim-1] + zstep * c; + cub[dim-1] = rlb[dim-1] + zstep * (c+1); + if (c == nprocs-1) cub[dim-1] = rub[dim-1]; + bbs[c] = bbox(clb, cub-cstr, cstr); + } + bbss[rl] = bbs; + } + vector > > > bbsss + = hh->make_multigrid_boxes(bbss, mglevels); + vector > pss(bbss.size()); + for (int rl=0; rl(bbss[rl].size()); + // make sure all processors have the same number of components + assert (bbss[rl].size() % nprocs == 0); + for (int c=0; c<(int)bbss[rl].size(); ++c) { + pss[rl][c] = c % nprocs; // distribute among processors + } + } + hh->recompose(bbsss, pss); } } // namespace Carpet diff --git a/Carpet/Carpet/src/Restrict.cc b/Carpet/Carpet/src/Restrict.cc index c44540dba..c935b9336 100644 --- a/Carpet/Carpet/src/Restrict.cc +++ b/Carpet/Carpet/src/Restrict.cc @@ -1,19 +1,14 @@ #include -#include #include #include "cctk.h" -#include "cctk_Parameters.h" -#include "ggf.hh" -#include "gh.hh" +#include "Carpet/CarpetLib/src/ggf.hh" +#include "Carpet/CarpetLib/src/gh.hh" #include "carpet.hh" -extern "C" { - static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Restrict.cc,v 1.27 2004/05/21 18:16:23 schnetter Exp $"; - CCTK_FILEVERSION(Carpet_Carpet_Restrict_cc); -} +static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Restrict.cc,v 1.1 2001/07/04 12:29:47 schnetter Exp $"; @@ -23,84 +18,57 @@ namespace Carpet { - void Restrict (const cGH* cgh) + void Restrict (cGH* cgh) { - DECLARE_CCTK_PARAMETERS; + assert (component == -1); - assert (is_level_mode()); + Checkpoint ("%*sRestrict", 2*reflevel, ""); - if (suppress_restriction) { - Checkpoint ("Restriction suppressed"); - return; - } - - Checkpoint ("Restrict"); - - // Restrict - if (reflevel < reflevels-1) { - for (comm_state state; !state.done(); state.step()) { - for (int group=0; grouptime (tl, reflevel, mglevel); - - const CCTK_REAL time1 = vtt.at(m)->time (0, reflevel, mglevel); - const CCTK_REAL time2 - = (cgh->cctk_time - cctk_initial_time) / delta_time; - assert (fabs(time1 - time2) / (fabs(time1) + fabs(time2) + fabs(cgh->cctk_delta_time)) < 1e-12); - - for (int var=0; varcomponents(reflevel); ++c) { - arrdata.at(group).at(m).data.at(var)->ref_restrict - (state, tl, reflevel, c, mglevel, time); - } - } - } - - } // if group has storage - } // if grouptype == CCTK_GF - } // loop over groups - } // for state - } // if not finest refinement level - - - - // Sync - if (reflevel < reflevels-1) { - for (comm_state state; !state.done(); state.step()) { - for (int group=0; groupcomponents(reflevel); ++c) { - arrdata.at(group).at(m).data.at(var)->sync - (state, tl, reflevel, c, mglevel); - } - } - } - - } // if group has storage - } // if grouptype == CCTK_GF - } // loop over groups - } // for state - } // if not finest refinement level + if (reflevel == hh->reflevels()-1) return; + for (int group=0; groupcomponents(reflevel); ++c) { + arrdata[group].data[var]->ref_restrict + (tl, reflevel, c, mglevel); + } + for (int c=0; ccomponents(reflevel); ++c) { + arrdata[group].data[var]->sync (tl, reflevel, c, mglevel); + } + } + break; + + case CCTK_GF: + for (int var=0; var<(int)gfdata[group].data.size(); ++var) { + for (int c=0; ccomponents(reflevel); ++c) { + gfdata[group].data[var]->ref_restrict + (tl, reflevel, c, mglevel); + } + for (int c=0; ccomponents(reflevel); ++c) { + gfdata[group].data[var]->sync (tl, reflevel, c, mglevel); + } + } + break; + + default: + abort(); + } + + } // if group has storage + + } // loop over groups } } // namespace Carpet - diff --git a/Carpet/Carpet/src/SetupGH.cc b/Carpet/Carpet/src/SetupGH.cc index d19ab3dc3..c7837edc9 100644 --- a/Carpet/Carpet/src/SetupGH.cc +++ b/Carpet/Carpet/src/SetupGH.cc @@ -1,32 +1,16 @@ #include -#include -#include #include -#include - -#include -#include -#include #include "cctk.h" #include "cctk_Parameters.h" -#include "util_ErrorCodes.h" -#include "util_Table.h" - -#include "bbox.hh" -#include "defs.hh" -#include "dist.hh" -#include "ggf.hh" -#include "gh.hh" -#include "vect.hh" +#include "Carpet/CarpetLib/src/dist.hh" +#include "Carpet/CarpetLib/src/ggf.hh" +#include "Carpet/CarpetLib/src/vect.hh" #include "carpet.hh" -extern "C" { - static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/SetupGH.cc,v 1.85 2004/08/07 20:07:27 schnetter Exp $"; - CCTK_FILEVERSION(Carpet_Carpet_SetupGH_cc); -} +static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/SetupGH.cc,v 1.1 2001/07/04 12:29:47 schnetter Exp $"; @@ -36,925 +20,220 @@ namespace Carpet { - static bool CanTransferVariableType (const cGH * const cgh, const int group) - { - // Find out which types correspond to the default types -#if CCTK_INTEGER_PRECISION_1 -# define CCTK_DEFAULT_INTEGER_TYPE CCTK_VARIABLE_INT1 -#elif CCTK_INTEGER_PRECISION_2 -# define CCTK_DEFAULT_INTEGER_TYPE CCTK_VARIABLE_INT2 -#elif CCTK_INTEGER_PRECISION_4 -# define CCTK_DEFAULT_INTEGER_TYPE CCTK_VARIABLE_INT4 -#elif CCTK_INTEGER_PRECISION_8 -# define CCTK_DEFAULT_INTEGER_TYPE CCTK_VARIABLE_INT8 -#else -# error "Unsupported default integer type" -#endif - -#if CCTK_REAL_PRECISION_4 -# define CCTK_DEFAULT_REAL_TYPE CCTK_VARIABLE_REAL4 -# define CCTK_DEFAULT_COMPLEX_TYPE CCTK_VARIABLE_COMPLEX8 -#elif CCTK_REAL_PRECISION_8 -# define CCTK_DEFAULT_REAL_TYPE CCTK_VARIABLE_REAL8 -# define CCTK_DEFAULT_COMPLEX_TYPE CCTK_VARIABLE_COMPLEX16 -#elif CCTK_REAL_PRECISION_16 -# define CCTK_DEFAULT_REAL_TYPE CCTK_VARIABLE_REAL16 -# define CCTK_DEFAULT_COMPLEX_TYPE CCTK_VARIABLE_COMPLEX32 -#else -# error "Unsupported default real type" -#endif - - if (CCTK_NumVarsInGroupI(group) == 0) return true; - - const int var0 = CCTK_FirstVarIndexI(group); - const int type0 = CCTK_VarTypeI(var0); - int type1; - switch (type0) { - case CCTK_VARIABLE_INT: - type1 = CCTK_DEFAULT_INTEGER_TYPE; - break; - case CCTK_VARIABLE_REAL: - type1 = CCTK_DEFAULT_REAL_TYPE; - break; - case CCTK_VARIABLE_COMPLEX: - type1 = CCTK_DEFAULT_COMPLEX_TYPE; - break; - default: - type1 = type0; - } - switch (type1) { - -#ifdef CCTK_REAL8 - case CCTK_VARIABLE_REAL8: - // This type is supported. - return true; -#endif - -#ifdef CCTK_REAL4 - case CCTK_VARIABLE_REAL4: -#endif -#ifdef CCTK_REAL16 - case CCTK_VARIABLE_REAL16: -#endif -#ifdef CCTK_REAL4 /* CCTK_COMPLEX8 */ - case CCTK_VARIABLE_COMPLEX8: -#endif -#ifdef CCTK_REAL8 /* CCTK_COMPLEX16 */ - case CCTK_VARIABLE_COMPLEX16: -#endif -#ifdef CCTK_REAL16 /* CCTK_COMPLEX32 */ - case CCTK_VARIABLE_COMPLEX32: -#endif - // This type is not supported, but could be. - return false; - - case CCTK_VARIABLE_BYTE: -#ifdef CCTK_INT1 - case CCTK_VARIABLE_INT1: -#endif -#ifdef CCTK_INT2 - case CCTK_VARIABLE_INT2: -#endif -#ifdef CCTK_INT4 - case CCTK_VARIABLE_INT4: -#endif -#ifdef CCTK_INT8 - case CCTK_VARIABLE_INT8: -#endif - // This type is not supported, and cannot be. - return false; - - default: - { - CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, - "Internal error: encountered variable type %d (%s) for group %d (%s)", - type1, CCTK_VarTypeName(type1), - group, CCTK_GroupName(group)); - } - } - - // not reached - return false; - } - - - - static operator_type GetTransportOperator (const cGH * const cgh, - const int group) - { - assert (group>=0 && group= 0) { - have_prolong_string = true; - } else if (prolong_length == UTIL_ERROR_TABLE_NO_SUCH_KEY) { - // do nothing - } else { - assert (0); - } - } - - // Get prolongation parameter name - char prolong_param_string[1000]; - bool have_prolong_param_string = false; - { - const int prolong_param_length = Util_TableGetString - (gp.tagstable, sizeof prolong_param_string, prolong_param_string, - "ProlongationParameter"); - if (prolong_param_length >= 0) { - have_prolong_param_string = true; - } else if (prolong_param_length == UTIL_ERROR_TABLE_NO_SUCH_KEY) { - // do nothing - } else { - assert (0); - } - } - - // Complain if both are given - if (have_prolong_string && have_prolong_param_string) { - char * const groupname = CCTK_GroupName (group); - CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, - "Group \"%s\" has both the tags \"Prolongation\" and \"ProlongationParameter\". This is not possible.", - groupname); - free (groupname); - } - - // Map the parameter name - if (have_prolong_param_string) { - char * thorn; - char * name; - ierr = CCTK_DecomposeName (prolong_param_string, þ, &name); - if (ierr < 0) { - char * const groupname = CCTK_GroupName (group); - CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, - "Group \"%s\" has the \"ProlongationParameter\" tag \"%s\". This is not a valid parameter name.", - groupname, prolong_param_string); - free (groupname); - } - int type; - char const * const * const value - = (static_cast - (CCTK_ParameterGet (name, thorn, &type))); - if (! value || ! *value) { - char * const groupname = CCTK_GroupName (group); - CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, - "Group \"%s\" has the \"ProlongationParameter\" tag \"%s\". This parameter does not exist.", - groupname, prolong_param_string); - free (groupname); - } - if (type != PARAMETER_KEYWORD && type != PARAMETER_STRING) { - char * const groupname = CCTK_GroupName (group); - CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, - "Group \"%s\" has the \"ProlongationParameter\" tag \"%s\". This parameter has the wrong type; it must be either KEYWORD or STRING.", - groupname, prolong_param_string); - free (groupname); - } - free (thorn); - free (name); - assert (strlen(*value) < sizeof prolong_string); - strcpy (prolong_string, *value); - have_prolong_string = true; - } - - // Select a default, if necessary - if (! have_prolong_string) { - if (can_transfer) { - // Use the default - if (gp.numtimelevels == 1) { - // Only one time level: do not prolongate - char * const groupname = CCTK_GroupName (group); - CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, - "Group \"%s\" has only one time level; therefore it will not be prolongated or restricted.", - groupname); - free (groupname); - return op_none; - } else { - // Several time levels: use the default - return op_Lagrange; - } - } else { - if (gp.grouptype == CCTK_GF) { - char * const groupname = CCTK_GroupName (group); - CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, - "Group \"%s\" has the variable type \"%s\" which cannot be prolongated or restricted.", - groupname, CCTK_VarTypeName(gp.vartype)); - free (groupname); - return op_none; - } else { - return op_error; - } - } - } - - // Select the prolongation method - assert (have_prolong_string); - if (CCTK_Equals(prolong_string, "none")) { - return op_none; - } else if (CCTK_Equals(prolong_string, "Lagrange")) { - return op_Lagrange; - } else if (CCTK_Equals(prolong_string, "TVD")) { - return op_TVD; - } else if (CCTK_Equals(prolong_string, "ENO")) { - return op_ENO; - } else { - char * const groupname = CCTK_GroupName (group); - CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, - "Group \"%s\" has the unknown prolongation method \"%s\".", - groupname, prolong_string); - free (groupname); - return op_error; - } - return op_error; - } - - - void* SetupGH (tFleshConfig* fc, int convLevel, cGH* cgh) { DECLARE_CCTK_PARAMETERS; - int ierr; - assert (cgh->cctk_dim == dim); // Not sure what to do with that assert (convLevel==0); dist::pseudoinit(); - - // Initialise current position - mglevel = -1; - reflevel = -1; - map = -1; - component = -1; - - - - Waypoint ("Setting up the grid hierarchy"); - - // Processor information - Output ("Carpet is running on %d processors", CCTK_nProcs(cgh)); - - // Multigrid information - basemglevel = convergence_level; - mglevels = num_convergence_levels; - mgfact = convergence_factor; - maxmglevelfact = ipow(mgfact, mglevels-1); - cgh->cctk_convfac = mgfact; + Checkpoint ("starting SetupGH..."); // Refinement information maxreflevels = max_refinement_levels; - reffact = refinement_factor; - maxreflevelfact = ipow(reffact, maxreflevels-1); + maxreflevelfact + = (int)floor(pow((double)refinement_factor, maxreflevels-1) + 0.5); + + // Ghost zones + vect lghosts, ughosts; + if (ghost_size == -1) { + lghosts = vect(ghost_size_x, ghost_size_y, ghost_size_z); + ughosts = vect(ghost_size_x, ghost_size_y, ghost_size_z); + } else { + lghosts = vect(ghost_size, ghost_size, ghost_size); + ughosts = vect(ghost_size, ghost_size, ghost_size); + } - // Map information - carpetGH.maps = maps = num_maps;; + // Grid size + const int stride = maxreflevelfact; + vect npoints; + if (global_nsize == -1) { + npoints = vect(global_nx, global_ny, global_nz); + } else { + npoints = vect(global_nsize, global_nsize, global_nsize); + } + const vect str(stride); + const vect lb(lghosts * str); + const vect ub = (npoints - ughosts - 1) * str; + const bbox baseext(lb, ub, str); - // Allocate space for groups - groupdata.resize(CCTK_NumGroups()); - arrdata.resize(CCTK_NumGroups()); + // Allocate grid hierarchy + hh = new gh(refinement_factor, vertex_centered, + multigrid_factor, vertex_centered, + baseext); - vhh.resize(maps); - vdd.resize(maps); - vtt.resize(maps); + // Allocate time hierarchy + tt = new th(hh, maxreflevelfact); - // Loop over maps - for (int m=0; m 0.001)) { - CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, - "The domain size for map %d scaled for convergence level %d with convergence factor %d is not integer", - m, basemglevel, convergence_factor); - } - - // Sanity check - // (if this fails, someone requested an insane amount of memory) - assert (all(npoints <= INT_MAX)); - { - int max = INT_MAX; - for (int d=0; d(refinement_factor, vertex_centered, - convergence_factor, vertex_centered, baseext); - - // Allocate data hierarchy - vdd.at(m) = new dh(*vhh.at(m), lghosts, ughosts, - prolongation_order_space, buffer_width); - - // Allocate time hierarchy - vtt.at(m) = new th(*vhh.at(m), 1.0); - - if (max_refinement_levels > 1) { - const int prolongation_stencil_size - = vdd.at(m)->prolongation_stencil_size(); - const int min_nghosts - = ((prolongation_stencil_size + refinement_factor - 1) - / (refinement_factor-1)); - if (any(min(lghosts,ughosts) < min_nghosts)) { - CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, - "There are not enough ghost zones for the desired spatial prolongation order on map %d. With Carpet::prolongation_order_space=%d, you need at least %d ghost zones.", - m, prolongation_order_space, min_nghosts); - } - - } - - - - // Set initial refinement structure - - vector bbs; - vector obs; - if (strcmp(base_extents, "") == 0) { - - // default: one grid component covering everything - bbs.push_back (vhh.at(m)->baseextent); - obs.push_back (bbvect(true)); - - } else { - - // explicit grid components - // TODO: invent something for the other convergence levels - istringstream ext_str(base_extents); - try { - ext_str >> bbs; - } catch (input_error) { - CCTK_WARN (0, "Could not parse parameter \"base_extents\""); - } - CCTK_VInfo (CCTK_THORNSTRING, "Using %d grid patches", bbs.size()); - cout << "grid-patches-are " << bbs << endl; - if (bbs.size()<=0) { - CCTK_WARN (0, "Cannot evolve with 0 grid patches"); - } - istringstream ob_str (base_outerbounds); - try { - ob_str >> obs; - } catch (input_error) { - CCTK_WARN (0, "Could not parse parameter \"base_outerbounds\""); - } - assert (obs.size() == bbs.size()); - - } - - // Distribute onto the processors - // (TODO: this should be done globally for all maps) - vector ps; - SplitRegions (cgh, bbs, obs, ps); - - // Create all multigrid levels - vector > bbss; - MakeMultigridBoxes (cgh, bbs, obs, bbss); - - // Only one refinement level - vector > > bbsss(1); - vector > obss(1); - vector > pss(1); - bbsss.at(0) = bbss; - obss.at(0) = obs; - pss.at(0) = ps; - - // Check the regions - CheckRegions (bbsss, obss, pss); - -#if 0 - // Do this later, because CactusBase/IO might not yet be initialised - // Write grid structure to file - OutputGridStructure (cgh, m, bbsss, obss, pss); -#endif - - // Recompose grid hierarchy - vhh.at(m)->recompose (bbsss, obss, pss, false); - - CCTK_INFO ("Grid structure (grid points):"); - const int rl = 0; - for (int c=0; ccomponents(rl); ++c) { - for (int ml=0; mlmglevels(rl,c); ++ml) { - const ivect lower = vhh.at(m)->extents.at(rl).at(c).at(ml).lower(); - const ivect upper = vhh.at(m)->extents.at(rl).at(c).at(ml).upper(); - const int convfact = ipow(mgfact, ml); - assert (all(lower % maxreflevelfact == 0)); - assert (all(upper % maxreflevelfact == 0)); - assert (all(((upper - lower) / maxreflevelfact) % convfact == 0)); - cout << " [" << ml << "][" << rl << "][" << m << "][" << c << "]" - << " exterior extent: " << lower / maxreflevelfact - << " : " << upper / maxreflevelfact - << " (" << (upper - lower) / maxreflevelfact / convfact + 1 - << ")" << endl; - } - } - - } // loop over maps + // Allocate data hierarchy + dd = new dh(*hh, lghosts, ughosts, prolongation_order_space); - reflevels = 1; - for (int m=0; mreflevels() == reflevels); + if (max_refinement_levels > 1) { + const int prolongation_stencil_size = dd->prolongation_stencil_size(); + const int min_nghosts + = ((prolongation_stencil_size + refinement_factor - 1) + / (refinement_factor-1)); + if (any(min(lghosts,ughosts) < min_nghosts)) { + CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, + "There are not enough ghost zones for the desired spatial prolongation order. With Carpet::prolongation_order_space=%d, you need at least %d ghost zones.", + prolongation_order_space, min_nghosts); + } } - + // Allocate space for groups + scdata.resize(CCTK_NumGroups()); + arrdata.resize(CCTK_NumGroups()); + gfdata.resize(CCTK_NumGroups()); // Allocate space for variables in group (but don't enable storage // yet) for (int group=0; group=1 || gp.dim<=dim); + + switch (gp.disttype) { + case CCTK_DISTRIB_CONSTANT: + case CCTK_DISTRIB_DEFAULT: break; - - case CCTK_ARRAY: { - assert (gp.dim>=1 || gp.dim<=dim); - const CCTK_INT * const * const sz = CCTK_GroupSizesI(group); - const CCTK_INT * const * const gsz = CCTK_GroupGhostsizesI(group); - for (int d=0; d sizes(1), ghostsizes(0); + for (int d=0; d alb(0), aub(stride), astr(stride); + for (int d=0; d arrext(alb, aub-astr, astr); + + arrdata[group].hh = new gh(refinement_factor, vertex_centered, + multigrid_factor, vertex_centered, + arrext); + + arrdata[group].tt = new th (arrdata[group].hh, maxreflevelfact); + + vect alghosts(0), aughosts(0); + for (int d=0; d(*arrdata[group].hh, alghosts, aughosts, + prolongation_order_space); + + if (max_refinement_levels > 1) { + const int prolongation_stencil_size + = arrdata[group].dd->prolongation_stencil_size(); + const int min_nghosts + = ((prolongation_stencil_size + refinement_factor - 2) + / (refinement_factor-1)); + if (any(min(alghosts,aughosts) < min_nghosts)) { + CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, + "There are not enough ghost zones for the desired spatial prolongation order for the grid function group \"%s\". With Carpet::prolongation_order_space=%d, you need at least %d ghost zones.", + CCTK_GroupName(group), + prolongation_order_space, min_nghosts); + } } - - ivect alghosts(0), aughosts(0); - for (int d=0; d= 0) { - int status; - - status = Util_TableGetIntArray - (gp.tagstable, gp.dim, &convpowers[0], "convergence_power"); - if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY) { - // keep default: independent of convergence level - } else if (status == 1) { - // a scalar was given - convpowers = convpowers[0]; - } else if (status == gp.dim) { - // do nothing - } else { - char * const groupname = CCTK_GroupName(group); - CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, - "The key \"convergence_power\" in the tags table of group \"%s\" is wrong", - groupname); - free (groupname); - } - assert (all (convpowers >= 0)); - - status = Util_TableGetIntArray - (gp.tagstable, gp.dim, &convoffsets[0], "convergence_offset"); - if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY) { - // keep default: offset is 0 - } else if (status == 1) { - // a scalar was given - convoffsets = convoffsets[0]; - } else if (status == gp.dim) { - // do nothing - } else { - char * const groupname = CCTK_GroupName(group); - CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, - "The key \"convergence_offset\" in the tags table of group \"%s\" is wrong", - groupname); - free (groupname); - } - - } // if there is a group tags table - - rvect real_sizes - = ((sizes - convoffsets) - / pow(rvect(convergence_factor), convpowers * basemglevel) - + convoffsets); - for (int d=gp.dim; d 0.001)) { - char * const groupname = CCTK_GroupName(group); - CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, - "The shape of group \"%s\" scaled for convergence level %d with convergence factor %d is not integer", - groupname, basemglevel, convergence_factor); - free (groupname); - } - - - - assert (gp.disttype==CCTK_DISTRIB_CONSTANT - || gp.disttype==CCTK_DISTRIB_DEFAULT); - - if (gp.disttype==CCTK_DISTRIB_CONSTANT) { - if (! all (ghostsizes == 0)) { - char * const groupname = CCTK_GroupName(group); - CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, - "The group \"%s\" has DISTRIB=constant, but its ghostsize is not 0", - groupname); - free (groupname); - } - assert (all (ghostsizes == 0)); - const int d = gp.dim==0 ? 0 : gp.dim-1; - sizes[d] = (sizes[d] - 2*ghostsizes[d]) * CCTK_nProcs(cgh) + 2*ghostsizes[d]; - assert (sizes[d] >= 0); - } - - - - const ivect alb(0); - const ivect aub(sizes-1); - const ivect astr(1); - const ibbox abaseext(alb, aub, astr); - - assert (all(convpowers == convpowers[0])); - const int amgfact1 = ipow(mgfact, convpowers[0]); - - arrdata.at(group).at(0).hh - = new gh(refinement_factor, vertex_centered, - amgfact1, vertex_centered, - abaseext); - - arrdata.at(group).at(0).dd - = new dh(*arrdata.at(group).at(0).hh, alghosts, aughosts, 0, 0); - - arrdata.at(group).at(0).tt - = new th(*arrdata.at(group).at(0).hh, 1.0); - - - // Set refinement structure for scalars and arrays - vector bbs; - vector obs; - bbs.push_back (abaseext); - obs.push_back (bbvect(true)); - - // Split it into components, one for each processor - vector ps; - if (gp.disttype==CCTK_DISTRIB_CONSTANT) { - SplitRegions_AlongDir (cgh, bbs, obs, ps, gp.dim==0 ? 0 : gp.dim-1); - } else { - SplitRegions_Automatic (cgh, bbs, obs, ps); - } - - // Create all multigrid levels - vector > bbss (bbs.size()); - ivect amgfact; - iivect aoffset; - for (int d=0; d > > bbsss(1); - vector > obss(1); - vector > pss(1); - bbsss.at(0) = bbss; - obss.at(0) = obs; - pss.at(0) = ps; + arrdata[group].data.resize(CCTK_NumVarsInGroupI(group)); + for (int var=0; var<(int)scdata[group].data.size(); ++var) { + arrdata[group].data[var] = 0; + } + break; + } - // Recompose for this map - char * const groupname = CCTK_GroupName (group); - assert (groupname); - Checkpoint ("Recomposing grid array group \"%s\"", groupname); - free (groupname); - arrdata.at(group).at(0).hh->recompose (bbsss, obss, pss, false); + case CCTK_GF: { + assert (CCTK_GroupDimI(group) == dim); + gfdata[group].data.resize(CCTK_NumVarsInGroupI(group)); + for (int var=0; var<(int)scdata[group].data.size(); ++var) { + gfdata[group].data[var] = 0; + } break; - } // case of scalar or array + } default: - assert (0); - } // switch on group type - - // Initialise group information - groupdata.at(group).info.dim = gp.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 = GetTransportOperator (cgh, group); - - // Initialise group variables - for (int m=0; m<(int)arrdata.at(group).size(); ++m) { - - arrdata.at(group).at(m).data.resize(CCTK_NumVarsInGroupI(group)); - for (int var=0; var<(int)arrdata.at(group).at(m).data.size(); ++var) { - arrdata.at(group).at(m).data.at(var) = 0; - } - + abort(); } - - } // for group - - - - // Allocate level times - leveltimes.resize (mglevels); - for (int ml=0; mlcctk_time = 0; - cgh->cctk_delta_time = 1.0; + // Initialise cgh for (int d=0; dcctk_origin_space[d] = 0.0; - cgh->cctk_delta_space[d] = 1.0; + cgh->cctk_nghostzones[d] = dd->lghosts[d]; } - mglevel = 0; + // Initialise current position reflevel = 0; - map = 0; - component = 0; - - leave_local_mode (cgh); - leave_singlemap_mode (cgh); - leave_level_mode (cgh); - leave_global_mode (cgh); - + mglevel = 0; + component = -1; + // Recompose grid hierarchy + Recompose (cgh); - // Some statistics - if (verbose || veryverbose) { - - int num_gf_groups = 0; - int num_gf_vars = 0; - vector 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; - } - - for (int group=0; groupGH[convlev]; - for (int rl=reflevels-1; rl>=0; --rl) { - BEGIN_REVERSE_MGLEVEL_LOOP(cgh) { - enter_level_mode (cgh, rl); - do_global_mode = reflevel==0; - do_meta_mode = do_global_mode && mglevel==mglevels-1; - - Checkpoint ("Shutdown at iteration %d time %g%s%s", - cgh->cctk_iteration, (double)cgh->cctk_time, - (do_global_mode ? " (global)" : ""), - (do_meta_mode ? " (meta)" : "")); - - // Terminate - Checkpoint ("Scheduling TERMINATE"); - CCTK_ScheduleTraverse ("CCTK_TERMINATE", cgh, CallFunction); - - leave_level_mode (cgh); - } END_REVERSE_MGLEVEL_LOOP; - } // for rl + // Terminate + BEGIN_REFLEVEL_LOOP(cgh) { + Checkpoint ("%*sScheduling TERMINATE", 2*reflevel, ""); + CCTK_ScheduleTraverse ("CCTK_TERMINATE", cgh, CallFunction); + } END_REFLEVEL_LOOP(cgh); - BEGIN_REVERSE_MGLEVEL_LOOP(cgh) { - do_global_mode = true; - do_meta_mode = mglevel==mglevels-1; - - // Shutdown - Checkpoint ("Scheduling SHUTDOWN"); + // Shutdown + BEGIN_REFLEVEL_LOOP(cgh) { + Checkpoint ("%*sScheduling SHUTDOWN", 2*reflevel, ""); CCTK_ScheduleTraverse ("CCTK_SHUTDOWN", cgh, CallFunction); - - } END_REVERSE_MGLEVEL_LOOP; + } END_REFLEVEL_LOOP(cgh); CCTK_PRINTSEPARATOR; printf ("Done.\n"); - // earlier checkpoint before finalising MPI - Waypoint ("Done with shutdown"); + // earlier checkpoint before calling finalising MPI + Checkpoint ("done with Shutdown."); dist::finalize(); diff --git a/Carpet/Carpet/src/Storage.cc b/Carpet/Carpet/src/Storage.cc index 16847f077..a741ff5e7 100644 --- a/Carpet/Carpet/src/Storage.cc +++ b/Carpet/Carpet/src/Storage.cc @@ -2,17 +2,13 @@ #include #include "cctk.h" -#include "cctk_Parameters.h" -#include "dh.hh" -#include "gf.hh" +#include "Carpet/CarpetLib/src/dh.hh" +#include "Carpet/CarpetLib/src/gf.hh" #include "carpet.hh" -extern "C" { - static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Storage.cc,v 1.36 2004/06/02 09:04:43 schnetter Exp $"; - CCTK_FILEVERSION(Carpet_Carpet_Storage_cc); -} +static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Storage.cc,v 1.1 2001/07/04 12:29:47 schnetter Exp $"; @@ -22,44 +18,31 @@ namespace Carpet { - int EnableGroupStorage (const cGH* cgh, const char* groupname) + int EnableGroupStorage (cGH* cgh, const char* groupname) { - DECLARE_CCTK_PARAMETERS; - - Checkpoint ("EnableGroupStorage \"%s\"", groupname); + Checkpoint ("%*sEnableGroupStorage %s", 2*reflevel, "", groupname); // TODO: Enabling storage for one refinement level has to enable // it for all other refinement levels as well. Disabling must // wait until all refinement levels have been disabled. + // TODO: Invent a mode "reflevel==-1" that is global, i. e. has + // effect for all refinement levels. This mode is used during + // INITIAL, and en-/disabling storage in it is also global. + const int group = CCTK_GroupIndex(groupname); assert (group>=0 && grouplocal_components(reflevel) == 1)); - } - if (CCTK_QueryGroupStorageI(cgh, group)) { // storage was enabled previously - const int n0 = CCTK_FirstVarIndexI(group); - assert (n0>=0); - const int num_tl = CCTK_NumTimeLevelsFromVarI(n0); - assert (num_tl>0); - return num_tl; + return 1; } // There is a difference between the Cactus time levels and the // Carpet time levels. If there are n time levels, then the // Cactus time levels are numbered 0 ... n-1, with the current - // time level being 0. In Carpet, the time levels are numbered - // -(n-1) ... 0, where the current time level is also 0. + // time level being n-1. In Carpet, the time levels are numbered + // -(n-1) ... 0, where the current time level is always 0. const int n0 = CCTK_FirstVarIndexI(group); assert (n0>=0); const int num_tl = CCTK_NumTimeLevelsFromVarI(n0); @@ -67,62 +50,78 @@ namespace Carpet { const int tmin = 1 - num_tl; const int tmax = 0; - if (grouptype == CCTK_GF) { - if (max_refinement_levels > 1) { - if (groupdata.at(group).transport_operator != op_none) { - if (num_tl <= prolongation_order_time) { - CCTK_VWarn (0, __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.", - CCTK_GroupName(group), - prolongation_order_time, prolongation_order_time+1); - } - } - } - } - - cGroup gp; - const int ierr = CCTK_GroupData (group, &gp); - assert (!ierr); - const int vectorlength = gp.vectorgroup ? gp.vectorlength : 1; - - assert (arrdata.at(group).at(0).data.size()==0 - || arrdata.at(group).at(0).data.at(0) == 0); - for (int m=0; m<(int)arrdata.at(group).size(); ++m) { - for (int var=0; var* vectorleader - = (gp.vectorgroup && vectorindex>0 - ? arrdata.at(group).at(m).data.at(var - vectorindex) - : NULL); - const int n = n0 + var; - switch (CCTK_VarTypeI(n)) { -#define TYPECASE(N,T) \ - case N: \ - arrdata.at(group).at(m).data.at(var) = new gf \ - (n, groupdata.at(group).transport_operator, \ - *arrdata.at(group).at(m).tt, *arrdata.at(group).at(m).dd, \ - tmin, tmax, prolongation_order_time, \ - vectorlength, vectorindex, (gf*)vectorleader); \ - break; + switch (CCTK_GroupTypeI(group)) { + + case CCTK_SCALAR: + assert (scdata[group].data.size()==0 + || scdata[group].data[0].size()==0 + || scdata[group].data[0][0].size()==0 + || scdata[group].data[0][0][0] == 0); + for (int var=0; var<(int)scdata[group].data.size(); ++var) { + for (int rl=0; rl<(int)scdata[group].data[var].size(); ++rl) { + for (int tl=0; tl<(int)scdata[group].data[var][rl].size(); ++tl) { + const int n = n0 + var; + switch (CCTK_VarTypeI(n)) { +#define TYPECASE(N,T) \ + case N: \ + scdata[group].data[var][rl][tl] = new T; \ + break; +#include "typecase" +#undef TYPECASE + default: + UnsupportedVarType(n); + } // switch + } // for + } // for + } // for + break; + + case CCTK_ARRAY: + assert (arrdata[group].data.size()==0 + || arrdata[group].data[0] == 0); + for (int var=0; var<(int)arrdata[group].data.size(); ++var) { + const int n = n0 + var; + switch (CCTK_VarTypeI(n)) { +#define TYPECASE(N,T) \ + case N: \ + arrdata[group].data[var] = new gf \ + (CCTK_VarName(n), *arrdata[group].tt, *arrdata[group].dd, \ + tmin, tmax); \ + break; #include "typecase" #undef TYPECASE - default: - UnsupportedVarType(n); - } // switch - - if (grouptype != CCTK_GF) { - for (int tl=0; tldata[n][tl] = ((*arrdata.at(group).at(m).data.at(var)) - (-tl, 0, c, 0)->storage()); - } - } - - } // for - } // for - -// PoisonGroup (cgh, group, alltimes); + default: + UnsupportedVarType(n); + } // switch + } // for + break; + + case CCTK_GF: + assert (gfdata[group].data.size()==0 + || gfdata[group].data[0] == 0); + for (int var=0; var<(int)gfdata[group].data.size(); ++var) { + const int n = n0 + var; + switch (CCTK_VarTypeI(n)) { +#define TYPECASE(N,T) \ + case N: \ + gfdata[group].data[var] = new gf \ + (CCTK_VarName(n), *tt, *(dh*)dd, tmin, tmax); \ + break; +#include "typecase" +#undef TYPECASE + default: + UnsupportedVarType(n); + } // switch + } // for + break; + + default: + abort(); + } + + // Reinitialise Cactus variables + set_component (cgh, component); + PoisonGroup (cgh, group, alltimes); // storage was not enabled previously return 0; @@ -130,122 +129,147 @@ namespace Carpet { - int DisableGroupStorage (const cGH* cgh, const char* groupname) + int DisableGroupStorage (cGH* cgh, const char* groupname) { - Checkpoint ("DisableGroupStorage \"%s\"", groupname); + Checkpoint ("%*sDisableGroupStorage %s", 2*reflevel, "", groupname); const int group = CCTK_GroupIndex(groupname); assert (group>=0 && group=0); - const int num_tl = CCTK_NumTimeLevelsFromVarI(n0); - assert (num_tl>0); - for (int m=0; m<(int)arrdata.at(group).size(); ++m) { - for (int var=0; var*)arrdata.at(group).at(m).data.at(var); \ - arrdata.at(group).at(m).data.at(var) = NULL; \ - break; + switch (CCTK_GroupTypeI(group)) { + + case CCTK_SCALAR: + if (scdata[group].data.size()==0 + || scdata[group].data[0].size()==0 + || scdata[group].data[0][0].size()==0 + || scdata[group].data[0][0][0] == 0) { + // group already has no storage + break; + } + for (int var=0; var<(int)scdata[group].data.size(); ++var) { + const int n = n0 + var; + for (int rl=0; rl<(int)scdata[group].data[var].size(); ++rl) { + for (int tl=0; tl<(int)scdata[group].data[var][rl].size(); ++tl) { + switch (CCTK_VarTypeI(n)) { +#define TYPECASE(N,T) \ + case N: \ + delete (T*)scdata[group].data[var][rl][tl]; \ + break; +#include "typecase" +#undef TYPECASE + default: + UnsupportedVarType(n); + } + scdata[group].data[var][rl][tl] = 0; + } + } + } + break; + + case CCTK_ARRAY: + if (arrdata[group].data.size()==0 || arrdata[group].data[0] == 0) { + // group already has no storage + break; + } + for (int var=0; var<(int)arrdata[group].data.size(); ++var) { + const int n = CCTK_FirstVarIndexI(group) + var; + switch (CCTK_VarTypeI(n)) { +#define TYPECASE(N,T) \ + case N: \ + delete (gf*)arrdata[group].data[var]; \ + break; +#include "typecase" +#undef TYPECASE + default: + UnsupportedVarType(n); + } // switch + arrdata[group].data[var] = 0; + } // for + break; + + case CCTK_GF: + if (gfdata[group].data.size()==0 + || gfdata[group].data[0] == 0) { + // group already has no storage + break; + } + for (int var=0; var<(int)gfdata[group].data.size(); ++var) { + const int n = CCTK_FirstVarIndexI(group) + var; + switch (CCTK_VarTypeI(n)) { +#define TYPECASE(N,T) \ + case N: \ + delete (gf*)gfdata[group].data[var]; \ + break; #include "typecase" #undef TYPECASE - default: - UnsupportedVarType(n); - } // switch - arrdata.at(group).at(m).data.at(var) = NULL; - - if (CCTK_GroupTypeI(group) != CCTK_GF) { - for (int tl=0; tldata[n][tl] = NULL; - } - } - - } // for - } // for - - // storage was enabled previously + default: + UnsupportedVarType(n); + } // switch + gfdata[group].data[var] = 0; + } // for + break; + + default: + abort(); + } + + // Reinitialise Cactus variables + set_component (cgh, component); + + // storage was not disabled previously return 1; } - int QueryGroupStorageB (const cGH* cgh, int group, const char* groupname) + int QueryGroupStorageB (cGH* cgh, int group, const char* groupname) { if (groupname) { group = CCTK_GroupIndex(groupname); } assert (group>=0 && group=0 && n=0 && group=0 && dir=0 && group -#include -#include -#include -#include -#include - -#include -#include - -#include - -#include "cctk.h" -#include "cctk_Parameters.h" -#include "cctki_GHExtensions.h" -#include "cctki_ScheduleBindings.h" -#include "cctki_WarnLevel.h" - -#include "Carpet/CarpetLib/src/bbox.hh" -#include "Carpet/CarpetLib/src/defs.hh" -#include "Carpet/CarpetLib/src/dh.hh" -#include "Carpet/CarpetLib/src/dist.hh" -#include "Carpet/CarpetLib/src/gf.hh" -#include "Carpet/CarpetLib/src/gh.hh" -#include "Carpet/CarpetLib/src/th.hh" -#include "Carpet/CarpetLib/src/vect.hh" - -#include "carpet.hh" - -static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Attic/carpet.cc,v 1.28 2001/07/02 13:22:09 schnetter Exp $"; - - - -namespace Carpet { - - using namespace std; - - static void Recompose (cGH* cgh); - static void CycleTimeLevels (cGH* cgh); - static void Restrict (cGH* cgh); - - enum checktimes { currenttime, - currenttimebutnotifonly, - allbutlasttime, - allbutcurrenttime, - alltimes }; - - static void Poison (cGH* cgh, checktimes where); - static void PoisonGroup (cGH* cgh, int group, checktimes where); - static void PoisonCheck (cGH* cgh, checktimes where); - - static void CalculateChecksums (cGH* cgh, checktimes where); - static void CheckChecksums (cGH* cgh, checktimes where); - - static int mintl (checktimes where, int num_tl); - static int maxtl (checktimes where, int num_tl); - - // Debugging output - static void Checkpoint (const char* fmt, ...); - - // Error output - static void UnsupportedVarType (int vindex); - - - - // Handle from CCTK_RegisterGHExtension - int GHExtension; - - // Maximum number of refinement levels - int maxreflevels; - - // Refinement factor on finest grid - int maxreflevelfact; - - // Current iteration per refinement level - vector iteration; - - // Current position on the grid hierarchy - int mglevel; - int reflevel; - int component; - - // refinement factor of current level: pow(refinement_factor, reflevel) - int reflevelfact; - - // Time step on base grid - CCTK_REAL base_delta_time; - - - - // Data for scalars - vector scdata; // [group] - - // Data for arrays - // TODO: have replicated arrays - vector arrdata; // [group] - - // The grid hierarchy - gh* hh; - th* tt; - dh* dd; - int gfsize[gfdim]; - - // Data for grid functions - vector gfdata; // [group] - - // Checksums - vector > > > checksums; // [n][rl][tl][c] - - - - int CarpetStartup() - { - CCTK_RegisterBanner ("AMR driver provided by Carpet"); - - GHExtension = CCTK_RegisterGHExtension("Carpet"); - CCTK_RegisterGHExtensionSetupGH (GHExtension, SetupGH); - - CCTK_OverloadInitialise (Initialise); - CCTK_OverloadEvolve (Evolve); - CCTK_OverloadShutdown (Shutdown); - - CCTK_OverloadSyncGroup (SyncGroup); - CCTK_OverloadEnableGroupStorage (EnableGroupStorage); - CCTK_OverloadDisableGroupStorage (DisableGroupStorage); - CCTK_OverloadEnableGroupComm (EnableGroupComm); - CCTK_OverloadDisableGroupComm (DisableGroupComm); - CCTK_OverloadBarrier (Barrier); - CCTK_OverloadExit (Exit); - CCTK_OverloadAbort (Abort); - CCTK_OverloadMyProc (MyProc); - CCTK_OverloadnProcs (nProcs); - CCTK_OverloadArrayGroupSizeB (ArrayGroupSizeB); - CCTK_OverloadQueryGroupStorageB (QueryGroupStorageB); - - return 0; - } - - - - void* SetupGH (tFleshConfig* fc, int convLevel, cGH* cgh) - { - DECLARE_CCTK_PARAMETERS; - - assert (maxdim>0 && gfdim>0 && gfdim <= maxdim); - assert (cgh->cctk_dim == gfdim); - - // Not sure what to do with that - assert (convLevel==0); - - dist::pseudoinit(); - Checkpoint ("starting SetupGH..."); - - // Refinement information - maxreflevels = max_refinement_levels; - maxreflevelfact - = (int)floor(pow((double)refinement_factor, maxreflevels-1) + 0.5); - - // Ghost zones - assert (gfdim == 3); - vect lghosts, ughosts; - if (ghost_size == -1) { - lghosts = vect(ghost_size_x, ghost_size_y, ghost_size_z); - ughosts = vect(ghost_size_x, ghost_size_y, ghost_size_z); - } else { - lghosts = vect(ghost_size, ghost_size, ghost_size); - ughosts = vect(ghost_size, ghost_size, ghost_size); - } - - // Grid size - const int stride = maxreflevelfact; - vect npoints; - if (global_nsize == -1) { - npoints = vect(global_nx, global_ny, global_nz); - } else { - npoints = vect(global_nsize, global_nsize, global_nsize); - } - - const vect str(stride); - const vect lb(lghosts * str); - const vect ub = (npoints - ughosts - 1) * str; - - const bbox baseext(lb, ub, str); - - // Allocate grid hierarchy - hh = new gh(refinement_factor, vertex_centered, - multigrid_factor, vertex_centered, - baseext); - - // Allocate time hierarchy - tt = new th(hh, maxreflevelfact); - - // Allocate data hierarchy - dd = new dh(*hh, lghosts, ughosts, prolongation_order_space); - - if (max_refinement_levels > 1) { - const int prolongation_stencil_size = dd->prolongation_stencil_size(); - const int min_nghosts - = ((prolongation_stencil_size + refinement_factor - 1) - / (refinement_factor-1)); - if (any(min(lghosts,ughosts) < min_nghosts)) { - CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, - "There are not enough ghost zones for the desired spatial prolongation order. With Carpet::prolongation_order_space=%d, you need at least %d ghost zones.", - prolongation_order_space, min_nghosts); - } - } - - // Allocate space for groups - scdata.resize(CCTK_NumGroups()); - arrdata.resize(CCTK_NumGroups()); - gfdata.resize(CCTK_NumGroups()); - - // Allocate space for variables in group (but don't enable storage - // yet) - for (int group=0; group=1 || gp.dim<=maxdim); - - // TODO -// assert (gp.disttype == CCTK_DISTRIB_CONSTANT); - assert (gp.disttype == CCTK_DISTRIB_DEFAULT); - - switch (gp.dim) { - -#define CODE \ - do { \ - vect alb, aub, astr; \ - for (int d=0; d arrext(alb, aub-astr, astr); \ - \ - arrdata[group].hh = new gh \ - (refinement_factor, vertex_centered, \ - multigrid_factor, vertex_centered, \ - arrext); \ - \ - arrdata[group].tt = new th (arrdata[group].hh, maxreflevelfact); \ - \ - vect alghosts, aughosts; \ - for (int d=0; d(*(gh*)arrdata[group].hh, \ - alghosts, aughosts, \ - prolongation_order_space); \ - \ - if (max_refinement_levels > 1) { \ - const int prolongation_stencil_size \ - = arrdata[group].dd->prolongation_stencil_size(); \ - const int min_nghosts \ - = ((prolongation_stencil_size + refinement_factor - 2) \ - / (refinement_factor-1)); \ - if (any(min(alghosts,aughosts) < min_nghosts)) { \ - CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, \ - "There are not enough ghost zones for the desired spatial prolongation order for the grid function group \"%s\". With Carpet::prolongation_order_space=%d, you need at least %d ghost zones.", \ - CCTK_GroupName(group), \ - prolongation_order_space, min_nghosts); \ - } \ - } \ - \ - } while (0) - - case 1: { - const int dim = 1; - CODE; - break; - } - - case 2: { - const int dim = 2; - CODE; - break; - } - - case 3: { - const int dim = 3; - CODE; - break; - } - - default: - abort(); - -#undef CODE - - } - - arrdata[group].data.resize(CCTK_NumVarsInGroupI(group)); - for (int var=0; var<(int)scdata[group].data.size(); ++var) { - arrdata[group].data[var] = 0; - } - break; - } - - case CCTK_GF: { - assert (CCTK_GroupDimI(group) == gfdim); - - gfdata[group].data.resize(CCTK_NumVarsInGroupI(group)); - for (int var=0; var<(int)scdata[group].data.size(); ++var) { - gfdata[group].data[var] = 0; - } - break; - } - - default: - abort(); - } - } - - // Initialise cgh - for (int d=0; dcctk_nghostzones[d] = dd->lghosts[d]; - } - - // Initialise current position - reflevel = 0; - mglevel = 0; - component = -1; - - // Recompose grid hierarchy - Recompose (cgh); - - // Initialise time step on coarse grid - base_delta_time = 0; - - // Current iteration - iteration.resize(maxreflevels, 0); - - // Set current position (this time for real) - set_reflevel (cgh, 0); - set_mglevel (cgh, 0); - set_component (cgh, -1); - - // Enable storage for all groups if desired - // XXX - if (true || enable_all_storage) { - for (int group=0; groupcctk_iteration = 0; - cgh->cctk_time = cctk_initial_time; - - // Enable storage and communtication - CCTKi_ScheduleGHInit (cgh); - - // Initialise stuff - CCTKi_InitGHExtensions (cgh); - - // Check parameters - Checkpoint ("PARAMCHECK"); - CCTK_ScheduleTraverse ("CCTK_PARAMCHECK", cgh, CallFunction); - CCTKi_FinaliseParamWarn(); - - BEGIN_REFLEVEL_LOOP(cgh) { - - // Checking - Poison (cgh, allbutlasttime); - - // Set up the grid - Checkpoint ("%*sScheduling BASEGRID", 2*reflevel, ""); - CCTK_ScheduleTraverse ("CCTK_BASEGRID", cgh, CallFunction); - if (reflevel==0) { - base_delta_time = cgh->cctk_delta_time; - } else { -// assert (abs(cgh->cctk_delta_time - base_delta_time / reflevelfactor) -// < 1e-6 * base_delta_time); - // This circumvents a bug in CactusBase/Time - cgh->cctk_delta_time = base_delta_time / reflevelfact; - } - - // Set up the initial data - Checkpoint ("%*sScheduling INITIAL", 2*reflevel, ""); - CCTK_ScheduleTraverse ("CCTK_INITIAL", cgh, CallFunction); - Checkpoint ("%*sScheduling POSTINITIAL", 2*reflevel, ""); - CCTK_ScheduleTraverse ("CCTK_POSTINITIAL", cgh, CallFunction); - - // Recover - Checkpoint ("%*sScheduling RECOVER_VARIABLES", 2*reflevel, ""); - CCTK_ScheduleTraverse ("CCTK_RECOVER_VARIABLES", cgh, CallFunction); - Checkpoint ("%*sScheduling CPINITIAL", 2*reflevel, ""); - CCTK_ScheduleTraverse ("CCTK_CPINITIAL", cgh, CallFunction); - - // Poststep - Checkpoint ("%*sScheduling POSTSTEP", 2*reflevel, ""); - CCTK_ScheduleTraverse ("CCTK_POSTSTEP", cgh, CallFunction); - - // Checking - PoisonCheck (cgh, allbutlasttime); - - } END_REFLEVEL_LOOP(cgh); - - BEGIN_REVERSE_REFLEVEL_LOOP(cgh) { - - // Restrict - Restrict (cgh); - - // Checking - CalculateChecksums (cgh, allbutlasttime); - - // Analysis - Checkpoint ("%*sScheduling ANALYSIS", 2*reflevel, ""); - CCTK_ScheduleTraverse ("CCTK_ANALYSIS", cgh, CallFunction); - - // Output - Checkpoint ("%*sOutputGH", 2*reflevel, ""); - CCTK_OutputGH (cgh); - - // Checking - CheckChecksums (cgh, allbutlasttime); - - } END_REVERSE_REFLEVEL_LOOP(cgh); - - Checkpoint ("done with Initialise."); - - return 0; - } - - - - int Evolve (tFleshConfig* fc) - { - DECLARE_CCTK_PARAMETERS; - - Checkpoint ("starting Evolve..."); - - const int convlev = 0; - cGH* cgh = fc->GH[convlev]; - - // Main loop - while (cgh->cctk_iteration < cctk_itlast - || (cctk_final_time >= cctk_initial_time - && cgh->cctk_time < cctk_final_time)) { - - // Advance time - ++cgh->cctk_iteration; - cgh->cctk_time += base_delta_time / maxreflevelfact; - - Checkpoint ("Evolving iteration %d...", cgh->cctk_iteration); - - BEGIN_REFLEVEL_LOOP(cgh) { - if ((cgh->cctk_iteration-1) % (maxreflevelfact/reflevelfact) == 0) { - - // Cycle time levels - CycleTimeLevels (cgh); - - // Advance level times - tt->advance_time (reflevel, mglevel); - for (int group=0; groupadvance_time (reflevel, mglevel); - break; - case CCTK_GF: - break; - default: - abort(); - } - } - - // Checking - CalculateChecksums (cgh, allbutcurrenttime); - Poison (cgh, currenttimebutnotifonly); - - // Evolve - Checkpoint ("%*sScheduling PRESTEP", 2*reflevel, ""); - CCTK_ScheduleTraverse ("CCTK_PRESTEP", cgh, CallFunction); - Checkpoint ("%*sScheduling EVOL", 2*reflevel, ""); - CCTK_ScheduleTraverse ("CCTK_EVOL", cgh, CallFunction); - Checkpoint ("%*sScheduling POSTSTEP", 2*reflevel, ""); - CCTK_ScheduleTraverse ("CCTK_POSTSTEP", cgh, CallFunction); - - // Checking - PoisonCheck (cgh, currenttimebutnotifonly); - - } - } END_REFLEVEL_LOOP(cgh); - - BEGIN_REVERSE_REFLEVEL_LOOP(cgh) { - if (cgh->cctk_iteration % (maxreflevelfact/reflevelfact) == 0) { - - // Restrict - Restrict (cgh); - - // Checking - CalculateChecksums (cgh, currenttime); - - // Checkpoint - Checkpoint ("%*sScheduling CHECKPOINT", 2*reflevel, ""); - CCTK_ScheduleTraverse ("CCTK_CHECKPOINT", cgh, CallFunction); - - // Analysis - Checkpoint ("%*sScheduling ANALYSIS", 2*reflevel, ""); - CCTK_ScheduleTraverse ("CCTK_ANALYSIS", cgh, CallFunction); - - // Output - Checkpoint ("%*sOutputGH", 2*reflevel, ""); - CCTK_OutputGH (cgh); - - // Checking - CheckChecksums (cgh, alltimes); - - } - } END_REVERSE_REFLEVEL_LOOP(cgh); - - } // main loop - - Checkpoint ("done with Evolve."); - - return 0; - } - - - - int Shutdown (tFleshConfig* fc) - { - DECLARE_CCTK_PARAMETERS; - - Checkpoint ("starting Shutdown..."); - - const int convlev = 0; - cGH* cgh = fc->GH[convlev]; - - // Terminate - BEGIN_REFLEVEL_LOOP(cgh) { - Checkpoint ("%*sScheduling TERMINATE", 2*reflevel, ""); - CCTK_ScheduleTraverse ("CCTK_TERMINATE", cgh, CallFunction); - } END_REFLEVEL_LOOP(cgh); - - // Shutdown - BEGIN_REFLEVEL_LOOP(cgh) { - Checkpoint ("%*sScheduling SHUTDOWN", 2*reflevel, ""); - CCTK_ScheduleTraverse ("CCTK_SHUTDOWN", cgh, CallFunction); - } END_REFLEVEL_LOOP(cgh); - - CCTK_PRINTSEPARATOR; - printf ("Done.\n"); - - // earlier checkpoint before calling finalising MPI - Checkpoint ("done with Shutdown."); - - dist::finalize(); - - return 0; - } - - - - int CallFunction (void* function, cFunctionData* attribute, void* data) - { - // Traverse one function on all components of one refinement level - // of one multigrid level - - assert (mglevel>=0); - assert (reflevel>=0); - -// Checkpoint ("%*sStarting CallFunction...", 2*reflevel, ""); - - cGH* cgh = (cGH*)data; - - if (attribute->global) { - // Global operation: call once per refinement level - - const int res = CCTK_CallFunction (function, attribute, data); - assert (res==0); - - } else { - // Local operation: call once per component - - BEGIN_COMPONENT_LOOP(cgh) { - // This requires that all processors have the same number of - // local components - if (hh->is_local(reflevel, component)) { - - const int res = CCTK_CallFunction (function, attribute, data); - assert (res==0); - - } // if is_local - - } END_COMPONENT_LOOP(cgh); - - } - -// Checkpoint ("%*sdone with CallFunction.", 2*reflevel, ""); - - // The return value indicates whether the grid functions have been - // synchronised. - // return 0: let the flesh do the synchronisation, if necessary - return 0; - } - - - - int SyncGroup (cGH* cgh, const char* groupname) - { - DECLARE_CCTK_PARAMETERS; - - if (hh->components(reflevel) > 1) assert (component == -1); - - Checkpoint ("%*sSyncGroup %s", 2*reflevel, "", groupname); - - const int group = CCTK_GroupIndex(groupname); - assert (group>=0 && group=0); - const int num_tl = CCTK_NumTimeLevelsFromVarI(n0); - assert (num_tl>0); - const int tl = 0; - - switch (CCTK_GroupTypeI(group)) { - - case CCTK_SCALAR: - // TODO: Check whether the local values are consistent - break; - - case CCTK_ARRAY: - assert (group<(int)arrdata.size()); - for (int var=0; var<(int)arrdata[group].data.size(); ++var) { - if (num_tl>1 && reflevel>0) { - for (int c=0; ccomponents(reflevel); ++c) { - arrdata[group].data[var]->ref_bnd_prolongate - (tl, reflevel, c, mglevel, - prolongation_order_space, prolongation_order_time); - } - } - for (int c=0; ccomponents(reflevel); ++c) { - arrdata[group].data[var]->sync (tl, reflevel, c, mglevel); - } - } - break; - - case CCTK_GF: - assert (group<(int)gfdata.size()); - for (int var=0; var<(int)gfdata[group].data.size(); ++var) { - if (num_tl>1 && reflevel>0) { - for (int c=0; ccomponents(reflevel); ++c) { - gfdata[group].data[var]->ref_bnd_prolongate - (tl, reflevel, c, mglevel, - prolongation_order_space, prolongation_order_time); - } - } - for (int c=0; ccomponents(reflevel); ++c) { - gfdata[group].data[var]->sync (tl, reflevel, c, mglevel); - } - } - break; - - default: - abort(); - } - - return 0; - } - - - - int EnableGroupStorage (cGH* cgh, const char* groupname) - { - Checkpoint ("%*sEnableGroupStorage %s", 2*reflevel, "", groupname); - - // TODO: Enabling storage for one refinement level has to enable - // it for all other refinement levels as well. Disabling must - // wait until all refinement levels have been disabled. - - // TODO: Invent a mode "reflevel==-1" that is global, i. e. has - // effect for all refinement levels. This mode is used during - // INITIAL, and en-/disabling storage in it is also global. - - const int group = CCTK_GroupIndex(groupname); - assert (group>=0 && group=0); - const int num_tl = CCTK_NumTimeLevelsFromVarI(n0); - assert (num_tl>0); - const int tmin = 1 - num_tl; - const int tmax = 0; - - switch (CCTK_GroupTypeI(group)) { - - case CCTK_SCALAR: - assert (scdata[group].data.size()==0 - || scdata[group].data[0].size()==0 - || scdata[group].data[0][0].size()==0 - || scdata[group].data[0][0][0] == 0); - for (int var=0; var<(int)scdata[group].data.size(); ++var) { - for (int rl=0; rl<(int)scdata[group].data[var].size(); ++rl) { - for (int tl=0; tl<(int)scdata[group].data[var][rl].size(); ++tl) { - const int n = n0 + var; - switch (CCTK_VarTypeI(n)) { -#define TYPECASE(N,T) \ - case N: \ - scdata[group].data[var][rl][tl] = new T; \ - break; -#include "typecase" -#undef TYPECASE - default: - UnsupportedVarType(n); - } - } - } - } - break; - - case CCTK_ARRAY: - assert (arrdata[group].data.size()==0 - || arrdata[group].data[0] == 0); - for (int var=0; var<(int)arrdata[group].data.size(); ++var) { - const int n = n0 + var; - switch (CCTK_GroupDimI(group)) { - case 1: { - const int dim=1; - switch (CCTK_VarTypeI(n)) { -#define TYPECASE(N,T) \ - case N: \ - arrdata[group].data[var] = new gf \ - (CCTK_VarName(n), *arrdata[group].tt, *(dh*)arrdata[group].dd, \ - tmin, tmax); \ - break; -#include "typecase" -#undef TYPECASE - default: - UnsupportedVarType(n); - } - break; - } - case 2: { - const int dim=2; - switch (CCTK_VarTypeI(n)) { -#define TYPECASE(N,T) \ - case N: \ - arrdata[group].data[var] = new gf \ - (CCTK_VarName(n), *arrdata[group].tt, *(dh*)arrdata[group].dd, \ - tmin, tmax); \ - break; -#include "typecase" -#undef TYPECASE - default: - UnsupportedVarType(n); - } - break; - } - case 3: { - const int dim=3; - switch (CCTK_VarTypeI(n)) { -#define TYPECASE(N,T) \ - case N: \ - arrdata[group].data[var] = new gf \ - (CCTK_VarName(n), *arrdata[group].tt, *(dh*)arrdata[group].dd, \ - tmin, tmax); \ - break; -#include "typecase" -#undef TYPECASE - default: - UnsupportedVarType(n); - } - break; - } - default: abort(); - } - } - break; - - case CCTK_GF: - assert (gfdata[group].data.size()==0 - || gfdata[group].data[0] == 0); - for (int var=0; var<(int)gfdata[group].data.size(); ++var) { - const int n = n0 + var; - switch (CCTK_GroupDimI(group)) { - case 1: { - const int dim=1; - switch (CCTK_VarTypeI(n)) { -#define TYPECASE(N,T) \ - case N: \ - gfdata[group].data[var] = new gf \ - (CCTK_VarName(n), *tt, *(dh*)dd, tmin, tmax); \ - break; -#include "typecase" -#undef TYPECASE - default: - UnsupportedVarType(n); - } - break; - } - case 2: { - const int dim=2; - switch (CCTK_VarTypeI(n)) { -#define TYPECASE(N,T) \ - case N: \ - gfdata[group].data[var] = new gf \ - (CCTK_VarName(n), *tt, *(dh*)dd, tmin, tmax); \ - break; -#include "typecase" -#undef TYPECASE - default: - UnsupportedVarType(n); - } - break; - } - case 3: { - const int dim=3; - switch (CCTK_VarTypeI(n)) { -#define TYPECASE(N,T) \ - case N: \ - gfdata[group].data[var] = new gf \ - (CCTK_VarName(n), *tt, *(dh*)dd, tmin, tmax); \ - break; -#include "typecase" -#undef TYPECASE - default: - UnsupportedVarType(n); - } - break; - } - default: abort(); - } - } - break; - - default: - abort(); - } - - // Reinitialise Cactus variables - set_component (cgh, component); - PoisonGroup (cgh, group, alltimes); - - // storage was not enabled previously - return 0; - } - - - - int DisableGroupStorage (cGH* cgh, const char* groupname) - { - Checkpoint ("%*sDisableGroupStorage %s", 2*reflevel, "", groupname); - - const int group = CCTK_GroupIndex(groupname); - assert (group>=0 && group*)arrdata[group].data[var]; \ - break; -#include "typecase" -#undef TYPECASE - default: - UnsupportedVarType(n); - } - break; - } - case 2: { - const int dim=2; - switch (CCTK_VarTypeI(n)) { -#define TYPECASE(N,T) \ - case N: \ - delete (gf*)arrdata[group].data[var]; \ - break; -#include "typecase" -#undef TYPECASE - default: - UnsupportedVarType(n); - } - break; - } - case 3: { - const int dim=3; - switch (CCTK_VarTypeI(n)) { -#define TYPECASE(N,T) \ - case N: \ - delete (gf*)arrdata[group].data[var]; \ - break; -#include "typecase" -#undef TYPECASE - default: - UnsupportedVarType(n); - } - break; - } - default: abort(); - } - arrdata[group].data[var] = 0; - } - break; - - case CCTK_GF: - if (gfdata[group].data.size()==0 - || gfdata[group].data[0] == 0) { - // group already has no storage - break; - } - for (int var=0; var<(int)gfdata[group].data.size(); ++var) { - const int n = CCTK_FirstVarIndexI(group) + var; - switch (CCTK_GroupDimI(group)) { - case 1: { - const int dim=1; - switch (CCTK_VarTypeI(n)) { -#define TYPECASE(N,T) \ - case N: \ - delete (gf*)gfdata[group].data[var]; \ - break; -#include "typecase" -#undef TYPECASE - default: - UnsupportedVarType(n); - } - break; - } - case 2: { - const int dim=2; - switch (CCTK_VarTypeI(n)) { -#define TYPECASE(N,T) \ - case N: \ - delete (gf*)gfdata[group].data[var]; \ - break; -#include "typecase" -#undef TYPECASE - default: - UnsupportedVarType(n); - } - break; - } - case 3: { - const int dim=3; - switch (CCTK_VarTypeI(n)) { -#define TYPECASE(N,T) \ - case N: \ - delete (gf*)gfdata[group].data[var]; \ - break; -#include "typecase" -#undef TYPECASE - default: - UnsupportedVarType(n); - } - break; - } - default: abort(); - } - gfdata[group].data[var] = 0; - } - break; - - default: - abort(); - } - - // Reinitialise Cactus variables - set_component (cgh, component); - - // storage was not disabled previously - return 1; - } - - - - int EnableGroupComm (cGH* cgh, const char* groupname) - { - // Communication is always enabled - return 0; - } - - int DisableGroupComm (cGH* cgh, const char* groupname) - { - // Communication is always enabled - return -1; - } - - - - static void Recompose (cGH* cgh) - { - DECLARE_CCTK_PARAMETERS; - - assert (component == -1); - Checkpoint ("%*sRecompose", 2*reflevel, ""); - - const int dim = gfdim; - const int nprocs = CCTK_nProcs(cgh); - const int reflevels = max_refinement_levels; // arbitrary value - const int mglevels = 1; // arbitrary value - vector > > bbss(reflevels); - // note: what this routine calls "ub" is "ub+str" elsewhere - vect rstr = hh->baseextent.stride(); - vect rlb = hh->baseextent.lower(); - vect rub = hh->baseextent.upper() + rstr; - for (int rl=0; rl0) { - // save old values - const vect oldrlb = rlb; - const vect oldrub = rub; - // calculate extent and centre - const vect rextent = rub - rlb + rstr; - const vect rcentre = rlb + (rextent / 2 / rstr) * rstr; - // calculate new extent - assert (all(rextent % hh->reffact == 0)); - const vect newrextent = rextent / hh->reffact; - // refined boxes have smaller stride - assert (all(rstr%hh->reffact == 0)); - rstr /= hh->reffact; - // refine (arbitrarily) around the center only - rlb = rcentre - (newrextent/2 / rstr) * rstr; - rub = rlb + newrextent - rstr; -// // refine (arbitrarily) around the lower boundary only -// rlb = rlb; -// rub = rlb + newrextent - rstr; - assert (all(rlb >= oldrlb && rub <= oldrub)); - } - vector > bbs(nprocs); - for (int c=0; c cstr = rstr; - vect clb = rlb; - vect cub = rub; - const int glonpz = (rub[dim-1] - rlb[dim-1]) / cstr[dim-1]; - const int locnpz = (glonpz + nprocs - 1) / nprocs; - const int zstep = locnpz * cstr[dim-1]; - clb[dim-1] = rlb[dim-1] + zstep * c; - cub[dim-1] = rlb[dim-1] + zstep * (c+1); - if (c == nprocs-1) cub[dim-1] = rub[dim-1]; - bbs[c] = bbox(clb, cub-cstr, cstr); - } - bbss[rl] = bbs; - } - vector > > > bbsss - = hh->make_multigrid_boxes(bbss, mglevels); - vector > pss(bbss.size()); - for (int rl=0; rl(bbss[rl].size()); - // make sure all processors have the same number of components - assert (bbss[rl].size() % nprocs == 0); - for (int c=0; c<(int)bbss[rl].size(); ++c) { - pss[rl][c] = c % nprocs; // distribute among processors - } - } - hh->recompose(bbsss, pss); - // TODO: recompose arrays too - } - - - - static void CycleTimeLevels (cGH* cgh) - { - Checkpoint ("%*sCycleTimeLevels", 2*reflevel, ""); - - for (int group=0; group0); - void* tmpdata = scdata[group].data[var][reflevel][0]; - for (int tl=1; tlcomponents(reflevel); ++c) { - arrdata[group].data[var]->cycle (reflevel, c, mglevel); - } - break; - } - case CCTK_GF: { - assert (group<(int)gfdata.size()); - assert (var<(int)gfdata[group].data.size()); - for (int c=0; ccomponents(reflevel); ++c) { - gfdata[group].data[var]->cycle (reflevel, c, mglevel); - } - break; - } - default: - abort(); - } - } - } - } - } - - - - static void Restrict (cGH* cgh) - { - assert (component == -1); - - Checkpoint ("%*sRestrict", 2*reflevel, ""); - - if (reflevel == hh->reflevels()-1) return; - - for (int group=0; groupcomponents(reflevel); ++c) { - arrdata[group].data[var]->ref_restrict - (tl, reflevel, c, mglevel); - } - for (int c=0; ccomponents(reflevel); ++c) { - arrdata[group].data[var]->sync (tl, reflevel, c, mglevel); - } - } - break; - - case CCTK_GF: - for (int var=0; var<(int)gfdata[group].data.size(); ++var) { - for (int c=0; ccomponents(reflevel); ++c) { - gfdata[group].data[var]->ref_restrict - (tl, reflevel, c, mglevel); - } - for (int c=0; ccomponents(reflevel); ++c) { - gfdata[group].data[var]->sync (tl, reflevel, c, mglevel); - } - } - break; - - default: - abort(); - } - - } // if group has storage - - } // loop over groups - } - - - - int Barrier (cGH* cgh) - { - MPI_Barrier (dist::comm); - return 0; - } - - - - int Exit (cGH* cgh, int retval) - { - CCTK_Barrier (cgh); - dist::finalize(); - exit (retval); - } - - int Abort (cGH* cgh, int retval) - { - MPI_Abort (dist::comm, retval); - abort (); - } - - - - int MyProc (cGH* cgh) - { - int rank; - MPI_Comm_rank (dist::comm, &rank); - return rank; - } - - int nProcs (cGH* cgh) - { - int size; - MPI_Comm_size (dist::comm, &size); - return size; - } - - - - const int* ArrayGroupSizeB (cGH* cgh, int dir, int group, - const char* groupname) - { - static const int zero = 0, one = 1; - - if (component == -1) { - // global routine - return &zero; - } - - if (groupname) { - group = CCTK_GroupIndex(groupname); - } - assert (group>=0 && group=0 && dir=0 && var=0 && group=0 && n=0 && group0); - - const int num_tl = CCTK_NumTimeLevelsFromVarI(n); - assert (num_tl>0); - const int min_tl = mintl(where, num_tl); - const int max_tl = maxtl(where, num_tl); - - for (int tl=min_tl; tl<=max_tl; ++tl) { - - switch (CCTK_GroupTypeFromVarI(n)) { - case CCTK_SCALAR: { - assert (group<(int)scdata.size()); - assert (var<(int)scdata[group].data.size()); - memset (cgh->data[n][tl], poison_value, sz); - break; - } - case CCTK_ARRAY: - case CCTK_GF: { - BEGIN_COMPONENT_LOOP(cgh) { - if (hh->is_local(reflevel,component)) { - const int dim = CCTK_GroupDimI(group); - int np = 1; - for (int d=0; ddata[n][tl], poison_value, np*sz); - } - } END_COMPONENT_LOOP(cgh); - break; - } - default: - abort(); - } - - } // for tl - - } // for var - } - - - - void PoisonCheck (cGH* cgh, const checktimes where) - { - DECLARE_CCTK_PARAMETERS; - - if (! check_for_poison) return; - - Checkpoint ("%*sPoisonCheck", 2*reflevel, ""); - - for (int group=0; group0); - const int min_tl = mintl(where, num_tl); - const int max_tl = maxtl(where, num_tl); - - for (int tl=min_tl; tl<=max_tl; ++tl) { - - switch (CCTK_GroupTypeFromVarI(n)) { - - case CCTK_SCALAR: { - bool poisoned=false; - switch (CCTK_VarTypeI(n)) { -#define TYPECASE(N,T) \ - case N: { \ - T worm; \ - memset (&worm, poison_value, sizeof(worm)); \ - poisoned = *(const T*)cgh->data[n][tl] == worm; \ - break; \ - } -#include "typecase" -#undef TYPECASE - default: - UnsupportedVarType(n); - } - if (poisoned) { - char* fullname = CCTK_FullName(n); - CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, - "The variable \"%s\" contains poison in timelevel %d", - fullname, tl); - free (fullname); - } - break; - } - - case CCTK_ARRAY: - case CCTK_GF: { - BEGIN_COMPONENT_LOOP(cgh) { - if (hh->is_local(reflevel,component)) { - const int dim = CCTK_GroupDimI(group); - int size[dim]; - for (int d=0; ddata[n][tl]; - int numpoison=0; - for (int k=0; k10) { - char* fullname = CCTK_FullName(n); - CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, - "The variable \"%s\" contains poison at %d locations in timelevel %d; not all locations were printed.", - fullname, numpoison, tl); - free (fullname); - } - } // if is local - } END_COMPONENT_LOOP(cgh); - break; - } - - default: - abort(); - } - - } // for tl - - } // for var - } // if has storage - } // for group - } - - - - void CalculateChecksums (cGH* cgh, const checktimes where) - { - DECLARE_CCTK_PARAMETERS; - - if (! checksum_timelevels) return; - - Checkpoint ("%*sCalculateChecksums", 2*reflevel, ""); - - checksums.resize(CCTK_NumVars()); - for (int group=0; group0); - const int min_tl = mintl(where, num_tl); - const int max_tl = maxtl(where, num_tl); - - checksums[n].resize(maxreflevels); - checksums[n][reflevel].resize(num_tl); - for (int tl=min_tl; tl<=max_tl; ++tl) { - switch (CCTK_GroupTypeFromVarI(n)) { - - case CCTK_SCALAR: { - checksums[n][reflevel][tl].resize(1); - const int c=0; - checksums[n][reflevel][tl][c].valid = false; - break; - } - - case CCTK_ARRAY: - case CCTK_GF: { - checksums[n][reflevel][tl].resize(hh->components(reflevel)); - BEGIN_COMPONENT_LOOP(cgh) { - checksums[n][reflevel][tl][component].valid = false; - } END_COMPONENT_LOOP(cgh); - break; - } - - default: - abort(); - } - } - } - } - - for (int group=0; group0); - const int num_tl = CCTK_NumTimeLevelsFromVarI(n); - assert (num_tl>0); - const int min_tl = mintl(where, num_tl); - const int max_tl = maxtl(where, num_tl); - - for (int tl=min_tl; tl<=max_tl; ++tl) { - switch (CCTK_GroupTypeFromVarI(n)) { - - case CCTK_SCALAR: { - const int c=0; - int chk = 0; - const void* data = cgh->data[n][tl]; - for (int i=0; iis_local(reflevel,component)) { - const int dim = CCTK_GroupDimI(group); - int np = 1; - for (int d=0; ddata[n][tl]; - int chk = 0; - for (int i=0; i0); - const int num_tl = CCTK_NumTimeLevelsFromVarI(n); - assert (num_tl>0); - const int min_tl = mintl(where, num_tl); - const int max_tl = maxtl(where, num_tl); - - assert ((int)checksums[n].size()==maxreflevels); - assert ((int)checksums[n][reflevel].size()==num_tl); - - for (int tl=min_tl; tl<=max_tl; ++tl) { - - bool unexpected_change = false; - - switch (CCTK_GroupTypeFromVarI(n)) { - - case CCTK_SCALAR: { - assert ((int)checksums[n][reflevel][tl].size()==1); - const int c=0; - if (checksums[n][reflevel][tl][c].valid) { - int chk = 0; - const void* data = cgh->data[n][tl]; - for (int i=0; icomponents(reflevel)); - BEGIN_COMPONENT_LOOP(cgh) { - if (checksums[n][reflevel][tl][component].valid) { - if (hh->is_local(reflevel,component)) { - const int dim = CCTK_GroupDimI(group); - int np = 1; - for (int d=0; ddata[n][tl]; - int chk = 0; - for (int i=0; i0); - 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; - case allbutlasttime: - // do include current time if there is only one time level - return 0; - case allbutcurrenttime: - return 1; - case alltimes: - return 0; - default: - abort(); - } - } - - int maxtl (const checktimes where, const int num_tl) - { - assert (num_tl>0); - switch (where) { - case currenttime: - return 0; - case currenttimebutnotifonly: - return 0; - case allbutlasttime: - return num_tl-2; - case allbutcurrenttime: - return num_tl-1; - case alltimes: - return num_tl-1; - default: - abort(); - } - } - - - - void Checkpoint (const char* fmt, ...) - { - DECLARE_CCTK_PARAMETERS; - if (verbose) { - va_list args; - char msg[1000]; - va_start (args, fmt); - vsnprintf (msg, sizeof(msg), fmt, args); - va_end (args); - CCTK_INFO (msg); - } - if (barriers) { - MPI_Barrier (dist::comm); - } - } - - - - void UnsupportedVarType (const int vindex) - { - assert (vindex>=0 && vindex=0 && rlreflevels()); - assert (component == -1); - - // Change - reflevel = rl; - const int dim = gfdim; - const bbox& base = hh->baseextent; - reflevelfact = (int)floor(pow((double)hh->reffact, reflevel)+0.5); - cgh->cctk_delta_time = base_delta_time / reflevelfact; - for (int d=0; dcctk_gsh[d] - = ((base.shape() / base.stride() - + dd->lghosts + dd->ughosts)[d] - 1) * reflevelfact + 1; - cgh->cctk_levfac[d] = reflevelfact; - } - } - - - - void set_mglevel (cGH* cgh, const int ml) - { - assert (ml==0); - assert (component==-1); - mglevel = ml; - cgh->cctk_convlevel = mglevel; - } - - - - void set_component (cGH* cgh, const int c) - { - assert (c==-1 || (c>=0 && ccomponents(reflevel))); - component = c; - - const int dim = gfdim; - - if (component == -1) { - // Global mode -- no component is active - - // Set Cactus parameters to pseudo values - for (int d=0; dcctk_lsh[d] = 0xdeadbeef; - cgh->cctk_bbox[2*d ] = 0xdeadbeef; - cgh->cctk_bbox[2*d+1] = 0xdeadbeef; - cgh->cctk_lbnd[d] = 0xdeadbeef; - cgh->cctk_ubnd[d] = 0xdeadbeef; - for (int stg=0; stgcctk_lssh[CCTK_LSSH_IDX(stg,d)] = 0xdeadbeef; - } - } - - // Set local grid function and array sizes - for (int d=0; d=0); - const int var = n - CCTK_FirstVarIndexI(group); - assert (var>=0); - - if (CCTK_QueryGroupStorageI(cgh, group)) { - // Group has storage - - switch (CCTK_GroupTypeFromVarI(n)) { - - case CCTK_SCALAR: - // Scalar variables can be accessed - assert (group<(int)scdata.size()); - assert (var<(int)scdata[group].data.size()); - assert (reflevel<(int)scdata[group].data[var].size()); - assert (tl<(int)scdata[group].data[var][reflevel].size()); - cgh->data[n][tl] = scdata[group].data[var][reflevel][tl]; - break; - - case CCTK_ARRAY: - case CCTK_GF: - // Arrays and grid functions cannot be accessed - cgh->data[n][tl] = 0; - break; - - default: - abort(); - } - - } else { - // Group has no storage - - cgh->data[n][tl] = 0; - - } - - } // for tl - } // for n - - } else { - // Local mode -- a component is active - - // Set Cactus parameters - for (int d=0; d& ext - = dd->boxes[reflevel][component][mglevel].exterior; - cgh->cctk_lsh[d] = (ext.shape() / ext.stride())[d]; - cgh->cctk_lbnd[d] = (ext.lower() / ext.stride())[d]; - cgh->cctk_ubnd[d] = (ext.upper() / ext.stride())[d]; - assert (cgh->cctk_lsh[d]>=0 && cgh->cctk_lsh[d]<=cgh->cctk_gsh[d]); - assert (cgh->cctk_lbnd[d]>=0 && cgh->cctk_ubnd[d]cctk_gsh[d]); - assert (cgh->cctk_lbnd[d]<=cgh->cctk_ubnd[d]+1); - // No outer boundaries on the finer grids - cgh->cctk_bbox[2*d ] - = reflevel==0 && cgh->cctk_lbnd[d] == 0; - cgh->cctk_bbox[2*d+1] - = reflevel==0 && cgh->cctk_ubnd[d] == cgh->cctk_gsh[d]-1; -#if 0 - // Do allow outer boundaries on the finer grids (but this is - // generally inconsistent -- c. f. periodicity) - const bbox& base = hh->baseextent; - cgh->cctk_bbox[2*d ] = (ext.lower() < base.lower())[d]; - cgh->cctk_bbox[2*d+1] = (ext.upper() > base.upper())[d]; -#endif - for (int stg=0; stgcctk_lssh[CCTK_LSSH_IDX(stg,d)] = cgh->cctk_lsh[d]; - } - } - - // Set local grid function and array sizes - const bbox ext = - dd->boxes[reflevel][component][mglevel].exterior; - for (int d=0; d ext = \ - ((dh*)arrdata[group].dd)-> \ - boxes[reflevel][component][mglevel].exterior; \ - for (int d=0; d=0); - const int var = n - CCTK_FirstVarIndexI(group); - assert (var>=0); - const int num_tl = CCTK_NumTimeLevelsFromVarI(n); - assert (num_tl>0); - - for (int tl=0; tldata[n][tl] = scdata[group].data[var][reflevel][tl]; - break; - - case CCTK_ARRAY: - assert (group<(int)arrdata.size()); - assert (var<(int)arrdata[group].data.size()); - cgh->data[n][tl] - = ((*arrdata[group].data[var]) - (-tl, reflevel, component, mglevel)->storage()); - break; - - case CCTK_GF: - assert (group<(int)gfdata.size()); - assert (var<(int)gfdata[group].data.size()); - cgh->data[n][tl] - = ((*gfdata[group].data[var]) - (-tl, reflevel, component, mglevel)->storage()); - break; - - default: - abort(); - } - if (CCTK_GroupTypeFromVarI(n)==CCTK_SCALAR - || hh->is_local(reflevel,component)) { - assert (cgh->data[n][tl]); - } else { - assert (! cgh->data[n][tl]); - } - - } else { - // Group has no storage - - cgh->data[n][tl] = 0; - - } // if ! has storage - - } // for tl - } // for n - - } - } - - - - MPI_Comm CarpetMPICommunicator () - { - return dist::comm; - } - -} // namespace Carpet diff --git a/Carpet/Carpet/src/carpet.hh b/Carpet/Carpet/src/carpet.hh index a221b300c..f3005fb49 100644 --- a/Carpet/Carpet/src/carpet.hh +++ b/Carpet/Carpet/src/carpet.hh @@ -1,13 +1,17 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/carpet.hh,v 1.10 2001/07/02 13:22:10 schnetter Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/carpet.hh,v 1.11 2001/07/04 12:29:48 schnetter Exp $ + +// It is assumed that the number of components of all arrays is equal +// to the number of components of the grid functions, and that their +// distribution onto the processors is the same, and that all +// processors own the same number of components. + +// TODO: treat scalars as arrays with a dimension of zero #include #include "cctk.h" #include "cctk_Schedule.h" -#include "Carpet/CarpetLib/src/dgdh.hh" -#include "Carpet/CarpetLib/src/dggf.hh" -#include "Carpet/CarpetLib/src/dggh.hh" #include "Carpet/CarpetLib/src/dh.hh" #include "Carpet/CarpetLib/src/ggf.hh" #include "Carpet/CarpetLib/src/gh.hh" @@ -17,8 +21,7 @@ namespace Carpet { - const int maxdim = 3; - const int gfdim = 3; + const int dim = 3; @@ -55,25 +58,24 @@ namespace Carpet { // Data for arrays struct arrdesc { - int dim; - dimgeneric_gh* hh; + gh* hh; th* tt; - dimgeneric_dh* dd; - vector data; // [var] - int size[maxdim]; + dh* dd; + vector* > data; // [var] + int size[dim]; }; extern vector arrdata; // [group] // Data for grid functions // The grid hierarchy - extern gh* hh; + extern gh* hh; extern th* tt; - extern dh* dd; - extern int gfsize[gfdim]; + extern dh* dd; + extern int gfsize[dim]; struct gfdesc { - vector data; // [var] + vector*> data; // [var] }; extern vector gfdata; // [group] @@ -187,4 +189,34 @@ namespace Carpet { MPI_Comm CarpetMPICommunicator(); } + + + // These are private; don't use these from the outside + + void Recompose (cGH* cgh); + void CycleTimeLevels (cGH* cgh); + void Restrict (cGH* cgh); + + enum checktimes { currenttime, + currenttimebutnotifonly, + allbutlasttime, + allbutcurrenttime, + alltimes }; + + int mintl (checktimes where, int num_tl); + int maxtl (checktimes where, int num_tl); + + void Poison (cGH* cgh, checktimes where); + void PoisonGroup (cGH* cgh, int group, checktimes where); + void PoisonCheck (cGH* cgh, checktimes where); + + void CalculateChecksums (cGH* cgh, checktimes where); + void CheckChecksums (cGH* cgh, checktimes where); + + // Debugging output + void Checkpoint (const char* fmt, ...); + + // Error output + void UnsupportedVarType (int vindex); + } // namespace Carpet diff --git a/Carpet/Carpet/src/helpers.cc b/Carpet/Carpet/src/helpers.cc index 70e892b4d..c7af53a78 100644 --- a/Carpet/Carpet/src/helpers.cc +++ b/Carpet/Carpet/src/helpers.cc @@ -1,208 +1,123 @@ #include -#include -#include #include -#include #include #include "cctk.h" -#include "cctk_FortranString.h" #include "cctk_Parameters.h" -#include "defs.hh" -#include "dist.hh" -#include "ggf.hh" +#include "Carpet/CarpetLib/src/dist.hh" #include "carpet.hh" -extern "C" { - static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/helpers.cc,v 1.49 2004/07/08 12:43:34 tradke Exp $"; - CCTK_FILEVERSION(Carpet_Carpet_helpers_cc); -} +static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/helpers.cc,v 1.1 2001/07/04 12:29:48 schnetter Exp $"; namespace Carpet { - + using namespace std; - - - - // Get Carpet's GH extension - CarpetGH const * GetCarpetGH (const cGH * const cgh) + + + + int Barrier (cGH* cgh) { - assert (cgh); - return &carpetGH; - } - - - - // Enable or disable prolongating - CCTK_INT CarpetEnableProlongating (const CCTK_INT flag) - { - assert (flag==0 || flag==1); - do_prolongate = flag; - if (do_prolongate) { - Checkpoint ("Prolongating enabled"); - } else { - Checkpoint ("Prolongating disabled"); - } - return 0; - } - - - - // Communication - - int Barrier (const cGH* cgh) - { - void *dummy = &dummy; - dummy = &cgh; - MPI_Barrier (dist::comm); return 0; } - - - + + + int Exit (cGH* cgh, int retval) { CCTK_Barrier (cgh); dist::finalize(); exit (retval); - return -999; } - + int Abort (cGH* cgh, int retval) { - void *dummy = &dummy; - dummy = &cgh; - - MPI_Comm comm = dist::comm; - if (comm == MPI_COMM_NULL) - { - comm = MPI_COMM_WORLD; - } - MPI_Abort (comm, retval); + MPI_Abort (dist::comm, retval); abort (); - return -999; } - - - - int MyProc (const cGH* cgh) + + + + int MyProc (cGH* cgh) { - // if there is no cgh yet, assume nothing has been initialised - // yet, and don't use dist::comm int rank; - MPI_Comm_rank (cgh ? dist::comm : MPI_COMM_WORLD, &rank); + MPI_Comm_rank (dist::comm, &rank); return rank; } - - int nProcs (const cGH* cgh) + + int nProcs (cGH* cgh) { - // if there is no cgh yet, assume nothing has been initialised - // yet, and don't use dist::comm int size; - MPI_Comm_size (cgh ? dist::comm : MPI_COMM_WORLD, &size); + MPI_Comm_size (dist::comm, &size); return size; } - - - - MPI_Comm CarpetMPIComm () + + + + const int* ArrayGroupSizeB (cGH* cgh, int dir, int group, + const char* groupname) { - return dist::comm; - } - - - - // Datatypes - - MPI_Datatype CarpetMPIDatatype (const int vartype) - { - switch (vartype) { -#define TYPECASE(N,T) \ - case N: { \ - T dummy; \ - return dist::datatype(dummy); \ + static const int zero = 0, one = 1; + + if (component == -1) { + // global routine + return &zero; } -#include "typecase" -#undef TYPECASE - default: - CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, - "Carpet does not support the variable type %d.", vartype); + + if (groupname) { + group = CCTK_GroupIndex(groupname); } - // notreached - return MPI_CHAR; - } - - MPI_Datatype CarpetSimpleMPIDatatype (const int vartype) - { - switch (vartype) { -#ifdef CARPET_COMPLEX - case CCTK_VARIABLE_COMPLEX: - return CarpetMPIDatatype (CCTK_VARIABLE_REAL); -#endif -#ifdef CARPET_COMPLEX8 -# ifdef CCTK_REAL4 - case CCTK_VARIABLE_COMPLEX8: - return CarpetMPIDatatype (CCTK_VARIABLE_REAL4); -# endif -#endif -#ifdef CARPET_COMPLEX16 -# ifdef CCTK_REAL8 - case CCTK_VARIABLE_COMPLEX16: - return CarpetMPIDatatype (CCTK_VARIABLE_REAL8); -# endif -#endif -#ifdef CARPET_COMPLEX32 -# ifdef CCTK_REAL16 - case CCTK_VARIABLE_COMPLEX32: - return CarpetMPIDatatype (CCTK_VARIABLE_REAL16); -# endif -#endif - default: - return CarpetMPIDatatype (vartype); + assert (group>=0 && group=0 && dir=0 && var0); @@ -211,9 +126,7 @@ namespace Carpet { return 0; case currenttimebutnotifonly: // don't include current time if there is only one time level - return num_tl>1 ? 0 : 1; - case previoustime: - return 1; + return num_tl>1 ? 0: 1; case allbutlasttime: // do include current time if there is only one time level return 0; @@ -222,11 +135,10 @@ namespace Carpet { case alltimes: return 0; default: - assert (0); + abort(); } - return -999; } - + int maxtl (const checktimes where, const int num_tl) { assert (num_tl>0); @@ -235,8 +147,6 @@ namespace Carpet { return 0; case currenttimebutnotifonly: return 0; - case previoustime: - return num_tl>1 ? 1 : 0; case allbutlasttime: return num_tl-2; case allbutcurrenttime: @@ -244,76 +154,20 @@ namespace Carpet { case alltimes: return num_tl-1; default: - assert (0); + abort(); } - return -999; } - - - - // Diagnostic output - - static void prepend_id (char* const msg, size_t const msglen) - { - if (mglevel!=-1) { - snprintf (msg+strlen(msg), msglen-strlen(msg), "[%d]", mglevel); - if (reflevel!=-1) { - snprintf (msg+strlen(msg), msglen-strlen(msg), "[%d]", reflevel); - if (map!=-1) { - snprintf (msg+strlen(msg), msglen-strlen(msg), "[%d]", map); - if (component!=-1) { - snprintf (msg+strlen(msg), msglen-strlen(msg), "[%d]", component); - } - } - } - snprintf (msg+strlen(msg), msglen-strlen(msg), " "); - } - } - - void Output (const char* fmt, ...) - { - DECLARE_CCTK_PARAMETERS; - va_list args; - char msg[1000]; - snprintf (msg, sizeof msg, "%s", ""); - prepend_id (msg + strlen(msg), sizeof msg - strlen(msg)); - va_start (args, fmt); - vsnprintf (msg + strlen(msg), sizeof msg - strlen(msg), fmt, args); - va_end (args); - CCTK_INFO (msg); - if (barriers) { - MPI_Barrier (dist::comm); - } - } - - void Waypoint (const char* fmt, ...) - { - DECLARE_CCTK_PARAMETERS; - if (verbose || veryverbose) { - va_list args; - char msg[1000]; - snprintf (msg, sizeof msg, "%s", ""); - prepend_id (msg + strlen(msg), sizeof msg - strlen(msg)); - va_start (args, fmt); - vsnprintf (msg + strlen(msg), sizeof msg - strlen(msg), fmt, args); - va_end (args); - CCTK_INFO (msg); - } - if (barriers) { - MPI_Barrier (dist::comm); - } - } - + + + void Checkpoint (const char* fmt, ...) { DECLARE_CCTK_PARAMETERS; - if (veryverbose) { + if (verbose) { va_list args; char msg[1000]; - snprintf (msg, sizeof msg, "%s", ""); - prepend_id (msg + strlen(msg), sizeof msg - strlen(msg)); va_start (args, fmt); - vsnprintf (msg + strlen(msg), sizeof msg - strlen(msg), fmt, args); + vsnprintf (msg, sizeof(msg), fmt, args); va_end (args); CCTK_INFO (msg); } @@ -321,9 +175,9 @@ namespace Carpet { MPI_Barrier (dist::comm); } } - - - + + + void UnsupportedVarType (const int vindex) { assert (vindex>=0 && vindex=0 && rlreflevels()); + assert (component == -1); + + // Change + reflevel = rl; + const bbox& base = hh->baseextent; + reflevelfact = (int)floor(pow((double)hh->reffact, reflevel)+0.5); + cgh->cctk_delta_time = base_delta_time / reflevelfact; + for (int d=0; dcctk_gsh[d] + = ((base.shape() / base.stride() + + dd->lghosts + dd->ughosts)[d] - 1) * reflevelfact + 1; + cgh->cctk_levfac[d] = reflevelfact; + } + } + + + + void set_mglevel (cGH* cgh, const int ml) + { + assert (ml==0); + assert (component==-1); + mglevel = ml; + cgh->cctk_convlevel = mglevel; + } + + + + void set_component (cGH* cgh, const int c) + { + assert (c==-1 || (c>=0 && ccomponents(reflevel))); + component = c; + + if (component == -1) { + // Global mode -- no component is active + + // Set Cactus parameters to pseudo values + for (int d=0; dcctk_lsh[d] = 0xdeadbeef; + cgh->cctk_bbox[2*d ] = 0xdeadbeef; + cgh->cctk_bbox[2*d+1] = 0xdeadbeef; + cgh->cctk_lbnd[d] = 0xdeadbeef; + cgh->cctk_ubnd[d] = 0xdeadbeef; + for (int stg=0; stgcctk_lssh[CCTK_LSSH_IDX(stg,d)] = 0xdeadbeef; + } + } + + // Set local grid function and array sizes + for (int d=0; d=0); + const int var = n - CCTK_FirstVarIndexI(group); + assert (var>=0); + + if (CCTK_QueryGroupStorageI(cgh, group)) { + // Group has storage + + switch (CCTK_GroupTypeFromVarI(n)) { + + case CCTK_SCALAR: + // Scalar variables can be accessed + assert (group<(int)scdata.size()); + assert (var<(int)scdata[group].data.size()); + assert (reflevel<(int)scdata[group].data[var].size()); + assert (tl<(int)scdata[group].data[var][reflevel].size()); + cgh->data[n][tl] = scdata[group].data[var][reflevel][tl]; + break; + + case CCTK_ARRAY: + case CCTK_GF: + // Arrays and grid functions cannot be accessed + cgh->data[n][tl] = 0; + break; + + default: + abort(); + } + + } else { + // Group has no storage + + cgh->data[n][tl] = 0; + + } + + } // for tl + } // for n + + } else { + // Local mode -- a component is active + + // Set Cactus parameters + for (int d=0; d& ext + = dd->boxes[reflevel][component][mglevel].exterior; + cgh->cctk_lsh[d] = (ext.shape() / ext.stride())[d]; + cgh->cctk_lbnd[d] = (ext.lower() / ext.stride())[d]; + cgh->cctk_ubnd[d] = (ext.upper() / ext.stride())[d]; + assert (cgh->cctk_lsh[d]>=0 && cgh->cctk_lsh[d]<=cgh->cctk_gsh[d]); + assert (cgh->cctk_lbnd[d]>=0 && cgh->cctk_ubnd[d]cctk_gsh[d]); + assert (cgh->cctk_lbnd[d]<=cgh->cctk_ubnd[d]+1); + // No outer boundaries on the finer grids + cgh->cctk_bbox[2*d ] + = reflevel==0 && cgh->cctk_lbnd[d] == 0; + cgh->cctk_bbox[2*d+1] + = reflevel==0 && cgh->cctk_ubnd[d] == cgh->cctk_gsh[d]-1; +#if 0 + // Do allow outer boundaries on the finer grids (but this is + // generally inconsistent -- c. f. periodicity) + const bbox& base = hh->baseextent; + cgh->cctk_bbox[2*d ] = (ext.lower() < base.lower())[d]; + cgh->cctk_bbox[2*d+1] = (ext.upper() > base.upper())[d]; +#endif + for (int stg=0; stgcctk_lssh[CCTK_LSSH_IDX(stg,d)] = cgh->cctk_lsh[d]; + } + } + + // Set local grid function and array sizes + const bbox ext = + dd->boxes[reflevel][component][mglevel].exterior; + for (int d=0; dboxes.size()); + assert (component < (int)arrdata[group].dd->boxes[reflevel].size()); + assert (mglevel + < (int)arrdata[group].dd->boxes[reflevel][component].size()); + const bbox ext = + arrdata[group].dd->boxes[reflevel][component][mglevel].exterior; + for (int d=0; d=0); + const int var = n - CCTK_FirstVarIndexI(group); + assert (var>=0); + const int num_tl = CCTK_NumTimeLevelsFromVarI(n); + assert (num_tl>0); + + for (int tl=0; tldata[n][tl] = scdata[group].data[var][reflevel][tl]; + break; + + case CCTK_ARRAY: + assert (group<(int)arrdata.size()); + assert (var<(int)arrdata[group].data.size()); + cgh->data[n][tl] + = ((*arrdata[group].data[var]) + (-tl, reflevel, component, mglevel)->storage()); + break; + + case CCTK_GF: + assert (group<(int)gfdata.size()); + assert (var<(int)gfdata[group].data.size()); + cgh->data[n][tl] + = ((*gfdata[group].data[var]) + (-tl, reflevel, component, mglevel)->storage()); + break; + + default: + abort(); + } + if (CCTK_GroupTypeFromVarI(n)==CCTK_SCALAR + || hh->is_local(reflevel,component)) { + assert (cgh->data[n][tl]); + } else { + assert (! cgh->data[n][tl]); + } + + } else { + // Group has no storage + + cgh->data[n][tl] = 0; + + } // if ! has storage + + } // for tl + } // for n + + } + } + } // namespace Carpet diff --git a/Carpet/Carpet/src/make.code.defn b/Carpet/Carpet/src/make.code.defn index f07a00e65..0c1146197 100644 --- a/Carpet/Carpet/src/make.code.defn +++ b/Carpet/Carpet/src/make.code.defn @@ -1,8 +1,22 @@ # Main make.code.defn file for thorn Carpet -*-Makefile-*- -# $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/make.code.defn,v 1.1 2001/03/01 13:40:10 eschnett Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/make.code.defn,v 1.2 2001/07/04 12:29:48 schnetter Exp $ # Source files in this directory -SRCS = carpet.cc +SRCS = CallFunction.cc \ + CarpetStartup.cc \ + Checksum.cc \ + Comm.cc \ + Cycle.cc \ + Evolve.cc \ + Initialise.cc \ + Poison.cc \ + Recompose.cc \ + Restrict.cc \ + SetupGH.cc \ + Shutdown.cc \ + Storage.cc \ + helpers.cc \ + variables.cc # Subdirectories containing source files SUBDIRS = diff --git a/Carpet/Carpet/src/variables.cc b/Carpet/Carpet/src/variables.cc index f23c8d3f6..f884950f1 100644 --- a/Carpet/Carpet/src/variables.cc +++ b/Carpet/Carpet/src/variables.cc @@ -1,10 +1,13 @@ +#include +#include "Carpet/CarpetLib/src/dh.hh" +#include "Carpet/CarpetLib/src/gh.hh" +#include "Carpet/CarpetLib/src/th.hh" + +#include "carpet.hh" + +static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/variables.cc,v 1.1 2001/07/04 12:29:48 schnetter Exp $"; -#include "variables.hh" -extern "C" { - static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/variables.cc,v 1.22 2004/05/21 18:16:23 schnetter Exp $"; - CCTK_FILEVERSION(Carpet_Carpet_variables_cc); -} namespace Carpet { @@ -18,81 +21,42 @@ namespace Carpet { // Maximum number of refinement levels int maxreflevels; - // Refinement levels - int reflevels; - - // Refinement factor - int reffact; - // Refinement factor on finest grid int maxreflevelfact; - // Base multigrid level - int basemglevel; - - // Multigrid levels - int mglevels; - - // Multigrid factor - int mgfact; - - // Multigrid factor on coarsest grid - int maxmglevelfact; - - // Maps - int maps; - - + // Current iteration per refinement level + vector iteration; // Current position on the grid hierarchy - int reflevel; int mglevel; - int map; + int reflevel; int component; - // refinement factor of current level: ipow(refinement_factor, reflevel) + // refinement factor of current level: pow(refinement_factor, reflevel) int reflevelfact; - // multigrid factor of current level: ipow(multigrid_factor, mglevel) - int mglevelfact; + // Time step on base grid + CCTK_REAL base_delta_time; - // Carpet's GH - CarpetGH carpetGH; - - - - // Times and spaces on the refinement levels - CCTK_REAL global_time; - vector > leveltimes; // [mglevel][reflevel] - CCTK_REAL delta_time; - - vector > origin_space; // [mglevel] - vect delta_space; - - - - // Is this the time for a global mode call? - bool do_meta_mode; - bool do_global_mode; - - // Is prolongation enabled? - bool do_prolongate; + // Data for scalars + vector scdata; // [group] + // Data for arrays + // TODO: have replicated arrays + vector arrdata; // [group] + // The grid hierarchy + gh* hh; + th* tt; + dh* dd; + int gfsize[dim]; // Data for grid functions + vector gfdata; // [group] - // The grid hierarchy - vector*> vhh; // [map] - vector*> vdd; // [map] - vector*> vtt; // [map] - - // Data for the groups - vector groupdata; // [group] - - // Data for everything - vector > arrdata; // [group][map] + // Checksums + vector > > > checksums; // [n][rl][tl][c] } // namespace Carpet diff --git a/Carpet/CarpetIOASCII/src/ioascii.cc b/Carpet/CarpetIOASCII/src/ioascii.cc index 674bf9361..348c81630 100644 --- a/Carpet/CarpetIOASCII/src/ioascii.cc +++ b/Carpet/CarpetIOASCII/src/ioascii.cc @@ -24,7 +24,7 @@ #include "ioascii.hh" -static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOASCII/src/ioascii.cc,v 1.13 2001/07/02 13:22:10 schnetter Exp $"; +static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOASCII/src/ioascii.cc,v 1.14 2001/07/04 12:29:49 schnetter Exp $"; @@ -253,229 +253,109 @@ int CarpetIOASCII } } - switch (CCTK_GroupDimI(group)) { - - case 1: { - const int dim=1; - - // Find the output offset - vect offset(0); - switch (outdim) { - case 1: - // The offset doesn't matter in this case + assert (outdim <= CCTK_GroupDimI(group)); + + // Find the output offset + vect offset(0); + switch (outdim) { + case 1: + switch (dirs[0]) { + case 0: + offset[1] = GetGridOffset (cgh, 2, + "out%dD_xline_yi", "out_xline_yi", + "out%dD_xline_y", "out_xline_y", + out_xline_y); + offset[2] = GetGridOffset (cgh, 3, + "out%dD_xline_zi", "out_xline_zi", + "out%dD_xline_z", "out_xline_z", + out_xline_z); break; - default: - abort(); - } - - // Traverse all components on this refinement and - // multigrid level - BEGIN_COMPONENT_LOOP(cgh) { - - generic_gf* ff = 0; - - switch (CCTK_GroupTypeI(group)) { - - case CCTK_ARRAY: - assert (var < (int)arrdata[group].data.size()); - ff = (generic_gf*)arrdata[group].data[var]; - break; - - case CCTK_GF: - assert (var < (int)gfdata[group].data.size()); - ff = (generic_gf*)gfdata[group].data[var]; - break; - - default: - abort(); - } - - const generic_data* const data - = (*ff) (tl, reflevel, component, mglevel); - const bbox ext = data->extent(); - const vect offset1 = offset * ext.stride(); - - data->write_ascii (filename, cgh->cctk_iteration, offset1, dirs, - tl, reflevel, component, mglevel); - - } END_COMPONENT_LOOP(cgh); - break; - } - - case 2: { - const int dim=2; - - // Find the output offset - vect offset(0); - switch (outdim) { case 1: - switch (dirs[0]) { - case 0: - offset[1] = GetGridOffset (cgh, 2, - "out%dD_xline_yi", "out_xline_yi", - "out%dD_xline_y", "out_xline_y", - out_xline_y); - break; - case 1: - offset[0] = GetGridOffset (cgh, 1, - "out%dD_yline_xi", "out_yline_xi", - "out%dD_yline_x", "out_yline_x", - out_yline_x); - break; - default: - abort(); - } + offset[0] = GetGridOffset (cgh, 1, + "out%dD_yline_xi", "out_yline_xi", + "out%dD_yline_x", "out_yline_x", + out_yline_x); + offset[2] = GetGridOffset (cgh, 3, + "out%dD_yline_zi", "out_yline_zi", + "out%dD_yline_z", "out_yline_z", + out_yline_z); break; case 2: - // The offset doesn't matter in this case + offset[0] = GetGridOffset (cgh, 1, + "out%dD_zline_xi", "out_zline_xi", + "out%dD_zline_x", "out_zline_x", + out_zline_x); + offset[1] = GetGridOffset (cgh, 2, + "out%dD_zline_yi", "out_zline_yi", + "out%dD_zline_y", "out_zline_y", + out_zline_y); break; default: abort(); } - - // Traverse all components on this refinement and - // multigrid level - BEGIN_COMPONENT_LOOP(cgh) { - - generic_gf* ff = 0; - - switch (CCTK_GroupTypeI(group)) { - - case CCTK_ARRAY: - assert (var < (int)arrdata[group].data.size()); - ff = (generic_gf*)arrdata[group].data[var]; - break; - - case CCTK_GF: - assert (var < (int)gfdata[group].data.size()); - ff = (generic_gf*)gfdata[group].data[var]; - break; - - default: - abort(); - } - - const generic_data* const data - = (*ff) (tl, reflevel, component, mglevel); - const bbox ext = data->extent(); - const vect offset1 = offset * ext.stride(); - - data->write_ascii (filename, cgh->cctk_iteration, offset1, dirs, - tl, reflevel, component, mglevel); - - } END_COMPONENT_LOOP(cgh); break; + case 2: + if (dirs[0]==0 && dirs[1]==1) { + offset[2] = GetGridOffset + (cgh, 3, + "out%dD_xyplane_zi", "out_xyplane_zi", + "out%dD_xyplane_z", "out_xyplane_z", + out_xyplane_z); + } else if (dirs[0]==0 && dirs[1]==2) { + offset[1] = GetGridOffset + (cgh, 2, + "out%dD_xzplane_yi", "out_xzplane_yi", + "out%dD_xzplane_y", "out_xzplane_y", + out_xzplane_y); + } else if (dirs[0]==1 && dirs[1]==2) { + offset[0] = GetGridOffset + (cgh, 1, + "out%dD_yzplane_xi", "out_yzplane_xi", + "out%dD_yzplane_x", "out_yzplane_x", + out_yzplane_x); + } else { + abort(); + } + break; + case 3: + // The offset doesn't matter in this case + break; + default: + abort(); } + + // Traverse all components on this refinement and + // multigrid level + BEGIN_COMPONENT_LOOP(cgh) { - case 3: { - const int dim=3; + generic_gf* ff = 0; - // Find the output offset - vect offset(0); - switch (outdim) { - case 1: - switch (dirs[0]) { - case 0: - offset[1] = GetGridOffset (cgh, 2, - "out%dD_xline_yi", "out_xline_yi", - "out%dD_xline_y", "out_xline_y", - out_xline_y); - offset[2] = GetGridOffset (cgh, 3, - "out%dD_xline_zi", "out_xline_zi", - "out%dD_xline_z", "out_xline_z", - out_xline_z); - break; - case 1: - offset[0] = GetGridOffset (cgh, 1, - "out%dD_yline_xi", "out_yline_xi", - "out%dD_yline_x", "out_yline_x", - out_yline_x); - offset[2] = GetGridOffset (cgh, 3, - "out%dD_yline_zi", "out_yline_zi", - "out%dD_yline_z", "out_yline_z", - out_yline_z); - break; - case 2: - offset[0] = GetGridOffset (cgh, 1, - "out%dD_zline_xi", "out_zline_xi", - "out%dD_zline_x", "out_zline_x", - out_zline_x); - offset[1] = GetGridOffset (cgh, 2, - "out%dD_zline_yi", "out_zline_yi", - "out%dD_zline_y", "out_zline_y", - out_zline_y); - break; - default: - abort(); - } - break; - case 2: - if (dirs[0]==0 && dirs[1]==1) { - offset[2] = GetGridOffset - (cgh, 3, - "out%dD_xyplane_zi", "out_xyplane_zi", - "out%dD_xyplane_z", "out_xyplane_z", - out_xyplane_z); - } else if (dirs[0]==0 && dirs[1]==2) { - offset[1] = GetGridOffset - (cgh, 2, - "out%dD_xzplane_yi", "out_xzplane_yi", - "out%dD_xzplane_y", "out_xzplane_y", - out_xzplane_y); - } else if (dirs[0]==1 && dirs[1]==2) { - offset[0] = GetGridOffset - (cgh, 1, - "out%dD_yzplane_xi", "out_yzplane_xi", - "out%dD_yzplane_x", "out_yzplane_x", - out_yzplane_x); - } else { - abort(); - } + switch (CCTK_GroupTypeI(group)) { + + case CCTK_ARRAY: + assert (var < (int)arrdata[group].data.size()); + ff = (generic_gf*)arrdata[group].data[var]; break; - case 3: - // The offset doesn't matter in this case + + case CCTK_GF: + assert (var < (int)gfdata[group].data.size()); + ff = (generic_gf*)gfdata[group].data[var]; break; + default: abort(); } - // Traverse all components on this refinement and - // multigrid level - BEGIN_COMPONENT_LOOP(cgh) { - - generic_gf* ff = 0; - - switch (CCTK_GroupTypeI(group)) { - - case CCTK_ARRAY: - assert (var < (int)arrdata[group].data.size()); - ff = (generic_gf*)arrdata[group].data[var]; - break; - - case CCTK_GF: - assert (var < (int)gfdata[group].data.size()); - ff = (generic_gf*)gfdata[group].data[var]; - break; - - default: - abort(); - } - - const generic_data* const data - = (*ff) (tl, reflevel, component, mglevel); - const bbox ext = data->extent(); - const vect offset1 = offset * ext.stride(); - - data->write_ascii (filename, cgh->cctk_iteration, offset1, dirs, - tl, reflevel, component, mglevel); - - } END_COMPONENT_LOOP(cgh); - break; - } + const generic_data* const data + = (*ff) (tl, reflevel, component, mglevel); + const bbox ext = data->extent(); + const vect offset1 = offset * ext.stride(); - default: abort(); + data->write_ascii (filename, cgh->cctk_iteration, offset1, dirs, + tl, reflevel, component, mglevel); - } // switch (CCTK_GroupDimI(group)) + } END_COMPONENT_LOOP(cgh); + break; // Append EOL after every complete set of components if (CCTK_MyProc(cgh)==0) { diff --git a/Carpet/CarpetLib/src/bboxset.cc b/Carpet/CarpetLib/src/bboxset.cc index 835a46e96..5713f2d35 100644 --- a/Carpet/CarpetLib/src/bboxset.cc +++ b/Carpet/CarpetLib/src/bboxset.cc @@ -5,7 +5,7 @@ copyright : (C) 2000 by Erik Schnetter email : schnetter@astro.psu.edu - $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bboxset.cc,v 1.7 2001/07/02 13:22:12 schnetter Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bboxset.cc,v 1.8 2001/07/04 12:29:51 schnetter Exp $ ***************************************************************************/ @@ -303,8 +303,6 @@ void bboxset::output (ostream& os) const { #if defined(TMPL_EXPLICIT) -template class bboxset; -template class bboxset; template class bboxset; #endif diff --git a/Carpet/CarpetLib/src/data.cc b/Carpet/CarpetLib/src/data.cc index d95f0dd6f..d21eb0bdd 100644 --- a/Carpet/CarpetLib/src/data.cc +++ b/Carpet/CarpetLib/src/data.cc @@ -5,7 +5,7 @@ copyright : (C) 2000 by Erik Schnetter email : schnetter@astro.psu.edu - $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.cc,v 1.12 2001/07/02 13:22:12 schnetter Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.cc,v 1.13 2001/07/04 12:29:51 schnetter Exp $ ***************************************************************************/ @@ -586,8 +586,6 @@ ostream& data::output (ostream& os) const { #if defined(TMPL_EXPLICIT) #define INSTANTIATE(T) \ -template class data; \ -template class data; \ template class data; #include "instantiate" diff --git a/Carpet/CarpetLib/src/defs.cc b/Carpet/CarpetLib/src/defs.cc index d2fed76b5..ac868e5a9 100644 --- a/Carpet/CarpetLib/src/defs.cc +++ b/Carpet/CarpetLib/src/defs.cc @@ -5,7 +5,7 @@ copyright : (C) 2000 by Erik Schnetter email : schnetter@astro.psu.edu - $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/defs.cc,v 1.6 2001/07/02 13:22:12 schnetter Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/defs.cc,v 1.7 2001/07/04 12:29:51 schnetter Exp $ ***************************************************************************/ @@ -76,20 +76,6 @@ ostream& output (ostream& os, const vector& v) { #include "bbox.hh" #include "bboxset.hh" -template ostream& output (ostream& os, const list >& l); -template ostream& output (ostream& os, const set >& s); -template ostream& output (ostream& os, const set >& s); -template ostream& output (ostream& os, const vector > >& v); -template ostream& output (ostream& os, const vector > >& v); -template ostream& output (ostream& os, const vector > > >& v); - -template ostream& output (ostream& os, const list >& l); -template ostream& output (ostream& os, const set >& s); -template ostream& output (ostream& os, const set >& s); -template ostream& output (ostream& os, const vector > >& v); -template ostream& output (ostream& os, const vector > >& v); -template ostream& output (ostream& os, const vector > > >& v); - template ostream& output (ostream& os, const list >& l); template ostream& output (ostream& os, const set >& s); template ostream& output (ostream& os, const set >& s); diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc index a22240ea8..f1512be42 100644 --- a/Carpet/CarpetLib/src/dh.cc +++ b/Carpet/CarpetLib/src/dh.cc @@ -6,7 +6,7 @@ copyright : (C) 2000 by Erik Schnetter email : schnetter@astro.psu.edu - $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.cc,v 1.15 2001/07/02 13:22:13 schnetter Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.cc,v 1.16 2001/07/04 12:29:51 schnetter Exp $ ***************************************************************************/ @@ -360,7 +360,5 @@ void dh::output (ostream& os) const { #if defined(TMPL_EXPLICIT) -template class dh<1>; -template class dh<2>; template class dh<3>; #endif diff --git a/Carpet/CarpetLib/src/gdata.cc b/Carpet/CarpetLib/src/gdata.cc index de80207c0..347bfb52c 100644 --- a/Carpet/CarpetLib/src/gdata.cc +++ b/Carpet/CarpetLib/src/gdata.cc @@ -5,7 +5,7 @@ copyright : (C) 2000 by Erik Schnetter email : schnetter@astro.psu.edu - $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.cc,v 1.13 2001/07/02 13:22:13 schnetter Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.cc,v 1.14 2001/07/04 12:29:52 schnetter Exp $ ***************************************************************************/ @@ -216,40 +216,6 @@ void generic_data::write_ascii (const string name, const int time, #if defined(TMPL_EXPLICIT) -template class generic_data<1>; - -template void generic_data<1> -::write_ascii (const string name, const int time, - const vect& org, const vect& dirs, - const int tl, const int rl, const int c, const int ml) const; -template void generic_data<1> -::write_ascii (const string name, const int time, - const vect& org, const vect& dirs, - const int tl, const int rl, const int c, const int ml) const; -template void generic_data<1> -::write_ascii (const string name, const int time, - const vect& org, const vect& dirs, - const int tl, const int rl, const int c, const int ml) const; - - - -template class generic_data<2>; - -template void generic_data<2> -::write_ascii (const string name, const int time, - const vect& org, const vect& dirs, - const int tl, const int rl, const int c, const int ml) const; -template void generic_data<2> -::write_ascii (const string name, const int time, - const vect& org, const vect& dirs, - const int tl, const int rl, const int c, const int ml) const; -template void generic_data<2> -::write_ascii (const string name, const int time, - const vect& org, const vect& dirs, - const int tl, const int rl, const int c, const int ml) const; - - - template class generic_data<3>; template void generic_data<3> diff --git a/Carpet/CarpetLib/src/gf.cc b/Carpet/CarpetLib/src/gf.cc index 42743ff29..30eb4a5b4 100644 --- a/Carpet/CarpetLib/src/gf.cc +++ b/Carpet/CarpetLib/src/gf.cc @@ -6,7 +6,7 @@ copyright : (C) 2000 by Erik Schnetter email : schnetter@astro.psu.edu - $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gf.cc,v 1.6 2001/07/02 13:22:14 schnetter Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gf.cc,v 1.7 2001/07/04 12:29:52 schnetter Exp $ ***************************************************************************/ @@ -80,8 +80,6 @@ ostream& gf::output (ostream& os) const { #if defined(TMPL_EXPLICIT) #define INSTANTIATE(T) \ -template class gf; \ -template class gf; \ template class gf; #include "instantiate" diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc index 4afe6983b..a81bb9bb2 100644 --- a/Carpet/CarpetLib/src/ggf.cc +++ b/Carpet/CarpetLib/src/ggf.cc @@ -6,7 +6,7 @@ copyright : (C) 2000 by Erik Schnetter email : schnetter@astro.psu.edu - $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.cc,v 1.11 2001/07/02 13:22:14 schnetter Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.cc,v 1.12 2001/07/04 12:29:52 schnetter Exp $ ***************************************************************************/ @@ -465,7 +465,5 @@ void generic_gf::ref_prolongate (int tl, int rl, int c, int ml, #if defined(TMPL_EXPLICIT) -template class generic_gf<1>; -template class generic_gf<2>; template class generic_gf<3>; #endif diff --git a/Carpet/CarpetLib/src/gh.cc b/Carpet/CarpetLib/src/gh.cc index 52f225e61..6cb6fcdda 100644 --- a/Carpet/CarpetLib/src/gh.cc +++ b/Carpet/CarpetLib/src/gh.cc @@ -7,7 +7,7 @@ copyright : (C) 2000 by Erik Schnetter email : schnetter@astro.psu.edu - $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.cc,v 1.8 2001/07/02 13:22:14 schnetter Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.cc,v 1.9 2001/07/04 12:29:52 schnetter Exp $ ***************************************************************************/ @@ -236,8 +236,6 @@ ostream& gh::output (ostream& os) const { #if defined(TMPL_EXPLICIT) -template class gh<1>; -template class gh<2>; template class gh<3>; #endif diff --git a/Carpet/CarpetSlab/src/carpetslab.cc b/Carpet/CarpetSlab/src/carpetslab.cc index 41cac072a..0bc745d3b 100644 --- a/Carpet/CarpetSlab/src/carpetslab.cc +++ b/Carpet/CarpetSlab/src/carpetslab.cc @@ -1,4 +1,4 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/Attic/carpetslab.cc,v 1.9 2001/07/02 13:22:15 schnetter Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/Attic/carpetslab.cc,v 1.10 2001/07/04 12:29:52 schnetter Exp $ #include #include @@ -18,7 +18,7 @@ #include "carpetslab.hh" -static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/Attic/carpetslab.cc,v 1.9 2001/07/02 13:22:15 schnetter Exp $"; +static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/Attic/carpetslab.cc,v 1.10 2001/07/04 12:29:52 schnetter Exp $"; @@ -66,13 +66,11 @@ namespace CarpetSlab { // Get info about group cGroup gp; CCTK_GroupData (group, &gp); - assert (gp.dim<=maxdim); + assert (gp.dim<=dim); assert (CCTK_QueryGroupStorageI(cgh, group)); const int typesize = CCTK_VarTypeSize(gp.vartype); assert (typesize>0); - const int dim = gp.dim; - // Check dimension assert (hdim>=0 && hdim<=gp.dim); @@ -107,9 +105,9 @@ namespace CarpetSlab { // } // Get insider information about variable - const dimgeneric_gh* myhh; - const dimgeneric_dh* mydd; - const dimgeneric_gf* myff; + const gh* myhh; + const dh* mydd; + const generic_gf* myff; switch (gp.grouptype) { case CCTK_SCALAR: abort(); @@ -155,129 +153,89 @@ namespace CarpetSlab { if (hh->components(reflevel) > 0) { - switch (CCTK_GroupDimI(group)) { - -#define CODE \ - do { \ - \ - /* convert types */ \ - const gh* myhh1 = dynamic_cast*>(myhh); \ - const dh* mydd1 = dynamic_cast*>(mydd); \ - const generic_gf* myff1 \ - = dynamic_cast*>(myff); \ - \ - /* Only temporarily */ \ - component = 0; \ - \ - /* Get sample data */ \ - const generic_data* mydata; \ - mydata = (*myff1)(tl, reflevel, component, mglevel); \ - \ - /* Stride of data in memory */ \ - const vect str = mydata->extent().stride(); \ - \ - /* Stride of collected data */ \ - vect hstr = str; \ - for (int dd=0; dd hlb; \ - for (int d=0; d hub = hlb; \ - for (int dd=0; dd hextent (hlb, hub, hstr); \ - assert (hextent.num_points() == totalsize); \ - \ - /* Create collector data object */ \ - void* myhdata = rank==collect_proc ? hdata : 0; \ - generic_data* const alldata = mydata->make_typed(); \ - alldata->allocate (hextent, collect_proc, myhdata); \ - \ - /* Done with the temporary stuff */ \ - mydata = 0; \ - component = -1; \ - \ - /* Loop over all components, copying data from them */ \ - assert (component == -1); \ - for (component=0; \ - componentcomponents(reflevel); \ - ++component) { \ - \ - /* Get data object */ \ - mydata = (*myff1)(tl, reflevel, component, mglevel); \ - \ - /* Calculate overlapping extents */ \ - const bboxset myextents \ - = ((mydd1->boxes[reflevel][component][mglevel].sync_not \ - | mydd1->boxes[reflevel][component][mglevel].interior) \ - & hextent); \ - \ - /* Loop over overlapping extents */ \ - for (bboxset::const_iterator \ - ext_iter = myextents.begin(); \ - ext_iter != myextents.end(); \ - ++ext_iter) { \ - \ - /* Copy data */ \ - alldata->copy_from (mydata, *ext_iter); \ - \ - } \ - \ - } /* Loop over components */ \ - component = -1; \ - \ - /* Copy result to all processors */ \ - if (dest_proc == -1) { \ - for (int proc=0; proc* const tmpdata = mydata->make_typed(); \ - tmpdata->allocate (alldata->extent(), proc, myhdata); \ - tmpdata->copy_from (alldata, alldata->extent()); \ - delete tmpdata; \ - \ - } \ - } \ - } /* Copy result */ \ - \ - delete alldata; \ - \ - } while(0) - - case 1: { - const int dim=1; - CODE; - break; + // Only temporarily + component = 0; + + // Get sample data + const generic_data* mydata; + mydata = (*myff)(tl, reflevel, component, mglevel); + + // Stride of data in memory + const vect str = mydata->extent().stride(); + + // Stride of collected data + vect hstr = str; + for (int dd=0; dd hlb(0); + for (int d=0; d hub = hlb; + for (int dd=0; dd hextent (hlb, hub, hstr); + assert (hextent.num_points() == totalsize); + + // Create collector data object + void* myhdata = rank==collect_proc ? hdata : 0; + generic_data* const alldata = mydata->make_typed(); + alldata->allocate (hextent, collect_proc, myhdata); + + // Done with the temporary stuff + mydata = 0; + component = -1; + + // Loop over all components, copying data from them + assert (component == -1); + for (component=0; componentcomponents(reflevel); ++component) { - default: abort(); + // Get data object + mydata = (*myff)(tl, reflevel, component, mglevel); -#undef CODE + // Calculate overlapping extents + const bboxset myextents + = ((mydd->boxes[reflevel][component][mglevel].sync_not + | mydd->boxes[reflevel][component][mglevel].interior) + & hextent); - } // switch + // Loop over overlapping extents + for (bboxset::const_iterator ext_iter = myextents.begin(); + ext_iter != myextents.end(); + ++ext_iter) { + + // Copy data + alldata->copy_from (mydata, *ext_iter); + + } + + } // Loop over components + component = -1; + + // Copy result to all processors + if (dest_proc == -1) { + for (int proc=0; proc* const tmpdata = mydata->make_typed(); + tmpdata->allocate (alldata->extent(), proc, myhdata); + tmpdata->copy_from (alldata, alldata->extent()); + delete tmpdata; + + } + } + } // Copy result + + delete alldata; } // if components>0 @@ -302,8 +260,8 @@ namespace CarpetSlab { void** const hdata, int hsize [/*hdim*/]) { - const int dim = CCTK_GroupDimFromVarI(vindex); - assert (dim>=1 && dim<=maxdim); + const int gpdim = CCTK_GroupDimFromVarI(vindex); + assert (gpdim>=1 && gpdim<=dim); // Check some arguments assert (hdim>=0 && hdim<=dim); diff --git a/Carpet/CarpetTest/interface.ccl b/Carpet/CarpetTest/interface.ccl index 55cd937c5..eb5323639 100644 --- a/Carpet/CarpetTest/interface.ccl +++ b/Carpet/CarpetTest/interface.ccl @@ -1,5 +1,5 @@ # Interface definition for thorn CarpetTest -# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetTest/interface.ccl,v 1.2 2001/07/09 09:00:26 schnetter Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetTest/interface.ccl,v 1.1 2001/07/04 12:29:53 schnetter Exp $ implements: CarpetTest @@ -22,8 +22,3 @@ CCTK_REAL arrg1 type=ARRAY dim=1 size=i8 { arr1 } "1D array" - -CCTK_REAL scg type=SCALAR -{ - sc -} "scalar" diff --git a/Carpet/CarpetTest/schedule.ccl b/Carpet/CarpetTest/schedule.ccl index 1ade016b9..49b161fba 100644 --- a/Carpet/CarpetTest/schedule.ccl +++ b/Carpet/CarpetTest/schedule.ccl @@ -1,14 +1,14 @@ # Schedule definitions for thorn CarpetTest -# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetTest/schedule.ccl,v 1.3 2003/11/05 16:18:39 schnetter Exp $ - -STORAGE: gfg arrg1 arrg2 arrg3 +# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetTest/schedule.ccl,v 1.1 2001/07/04 12:29:54 schnetter Exp $ schedule carpettest_check_sizes AT initial { LANG: C + STORAGE: gfg arrg1 arrg2 arrg3 } "Check grid function and grid array sizes" schedule carpettest_check_arguments AT initial { LANG: Fortran + STORAGE: gfg arrg1 arrg2 arrg3 } "Check grid function and grid array arguments" diff --git a/Carpet/CarpetTest/src/carpettest_check_arguments.F77 b/Carpet/CarpetTest/src/carpettest_check_arguments.F77 index 297e38977..7c4fab201 100644 --- a/Carpet/CarpetTest/src/carpettest_check_arguments.F77 +++ b/Carpet/CarpetTest/src/carpettest_check_arguments.F77 @@ -1,44 +1,16 @@ -c -*-Fortran-*- -c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetTest/src/carpettest_check_arguments.F77,v 1.5 2004/01/25 14:57:31 schnetter Exp $ +c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetTest/src/carpettest_check_arguments.F77,v 1.1 2001/07/04 12:29:54 schnetter Exp $ #include "cctk.h" #include "cctk_Arguments.h" -#include "cctk_Functions.h" #include "cctk_Parameters.h" subroutine carpettest_check_arguments (CCTK_ARGUMENTS) implicit none DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_FUNCTIONS DECLARE_CCTK_PARAMETERS - integer i,j,k - print *, "Xgfg ", Xgfg0, Xgfg1, Xgfg2 - print *, "Xarrg3 ", Xarrg30, Xarrg31, Xarrg32 - print *, "Xarrg2 ", Xarrg20, Xarrg21 - print *, "Xarrg1 ", Xarrg10 - print *, "Xscg" + print *, "gfg ", gfg0, gfg1, gfg2 + print *, "arrg3 ", arrg30, arrg31, arrg32 + print *, "arrg2 ", arrg20, arrg21 + print *, "arrg1 ", arrg10 print * - do k=1,Xgfg2 - do j=1,Xgfg1 - do i=1,Xgfg0 - gf(i,j,k) = i*10000 + j*100 + k - end do - end do - end do - do k=1,Xarrg32 - do j=1,Xarrg31 - do i=1,Xarrg30 - arr3(i,j,k) = i*10000 + j*100 + k - end do - end do - end do - do j=1,Xarrg21 - do i=1,Xarrg20 - arr2(i,j) = i*100 + j - end do - end do - do i=1,Xarrg10 - arr1(i) = i - end do - sc = 42 end diff --git a/Carpet/CarpetTest/src/carpettest_check_sizes.c b/Carpet/CarpetTest/src/carpettest_check_sizes.c index c4e17baf7..79ca08929 100644 --- a/Carpet/CarpetTest/src/carpettest_check_sizes.c +++ b/Carpet/CarpetTest/src/carpettest_check_sizes.c @@ -1,4 +1,4 @@ -/* $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetTest/src/carpettest_check_sizes.c,v 1.2 2001/07/09 09:00:27 schnetter Exp $ */ +/* $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetTest/src/carpettest_check_sizes.c,v 1.1 2001/07/04 12:29:54 schnetter Exp $ */ #include #include @@ -10,98 +10,125 @@ void carpettest_check_sizes (CCTK_ARGUMENTS); -static void print_scalar (const char *name, int sc); -static void print_scalar_descr (const char *name, int sc, const char *descr); -static void print_array (const char *name, int dim, const int *arr); - -static const char *grouptype_string (int grouptype); -static const char *disttype_string (int disttype); - void carpettest_check_sizes (CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS + char msg[1000]; int group; + int var; int dim; - cGroup data; - cGroupDynamicData dyndata; + CCTK_INT **sizes, **ghostsizes; + int d; + int gsh[3], lsh[3], lbnd[3], ubnd[3]; + const int *cip; + int size[3]; - dim = cctk_dim; - print_scalar ("cctk_dim", cctk_dim); - print_array ("cctk_gsh", dim, cctk_gsh); - print_array ("cctk_lsh", dim, cctk_lsh); - print_array ("cctk_lbnd", dim, cctk_lbnd); - print_array ("cctk_ubnd", dim, cctk_ubnd); - print_array ("cctk_bbox", 2*dim, cctk_bbox); - print_array ("cctk_nghostzones", dim, cctk_nghostzones); + sprintf (msg, "gsh: %d %d %d\n", cctk_gsh[0], cctk_gsh[1], cctk_gsh[2]); + printf (msg); + sprintf (msg, "lsh: %d %d %d\n", cctk_lsh[0], cctk_lsh[1], cctk_lsh[2]); + printf (msg); + sprintf (msg, "lbnd: %d %d %d\n", cctk_lbnd[0], cctk_lbnd[1], cctk_lbnd[2]); + printf (msg); + sprintf (msg, "ubnd: %d %d %d\n", cctk_ubnd[0], cctk_ubnd[1], cctk_ubnd[2]); + printf (msg); printf ("\n"); +#if 0 for (group=0; groupsetType (IObase::Float64); - const int dim = CCTK_GroupDimI(group); + const int gpdim = CCTK_GroupDimI(group); // Set coordinate information - CCTK_REAL lower[maxdim], upper[maxdim]; - double origin[maxdim], delta[maxdim], timestep; + CCTK_REAL lower[dim], upper[dim]; + double origin[dim], delta[dim], timestep; for (int d=0; dsetTopLevelParameters - (dim, origin, delta, timestep, maxreflevels); + (gpdim, origin, delta, timestep, maxreflevels); // Set refinement information int interlevel_timerefinement; - int interlevel_spacerefinement[maxdim]; - int initial_gridplacementrefinement[maxdim]; + int interlevel_spacerefinement[dim]; + int initial_gridplacementrefinement[dim]; interlevel_timerefinement = hh->reffact; for (int d=0; dreffact; @@ -261,77 +263,46 @@ namespace CarpetIOFlexIO { // level BEGIN_COMPONENT_LOOP(cgh) { - switch (CCTK_GroupDimI(group)) { - -#define CODE \ - do { \ - const generic_gf* ff = 0; \ - \ - switch (CCTK_GroupTypeI(group)) { \ - \ - case CCTK_ARRAY: \ - assert (var < (int)arrdata[group].data.size()); \ - ff = (generic_gf*)arrdata[group].data[var]; \ - break; \ - \ - case CCTK_GF: \ - assert (var < (int)gfdata[group].data.size()); \ - ff = (generic_gf*)gfdata[group].data[var]; \ - break; \ - \ - default: \ - abort(); \ - } \ - \ - const generic_data* const data \ - = (*ff) (tl, reflevel, component, mglevel); \ - \ - /* Make temporary copy on processor 0 */ \ - const bbox ext = data->extent(); \ - generic_data* const tmp = data->make_typed (); \ - tmp->allocate (ext, 0); \ - tmp->copy_from (data, ext); \ - \ - /* Write data */ \ - if (CCTK_MyProc(cgh)==0) { \ - int origin[dim], dims[dim]; \ - for (int d=0; dwrite (origin, dims, (void*)tmp->storage()); \ - } \ - \ - /* Delete temporary copy */ \ - delete tmp; \ - \ - } while (0) + const generic_gf* ff = 0; + + switch (CCTK_GroupTypeI(group)) { - case 1: { - const int dim=1; - CODE; + case CCTK_ARRAY: + assert (var < (int)arrdata[group].data.size()); + ff = (generic_gf*)arrdata[group].data[var]; break; - } - case 2: { - const int dim=2; - CODE; + case CCTK_GF: + assert (var < (int)gfdata[group].data.size()); + ff = (generic_gf*)gfdata[group].data[var]; break; - } - - case 3: { - const int dim=3; - CODE; - break; - } default: abort(); - -#undef CODE - } + const generic_data* const data + = (*ff) (tl, reflevel, component, mglevel); + + /* Make temporary copy on processor 0 */ + const bbox ext = data->extent(); + generic_data* const tmp = data->make_typed (); + tmp->allocate (ext, 0); + tmp->copy_from (data, ext); + + /* Write data */ + if (CCTK_MyProc(cgh)==0) { + int origin[dim], dims[dim]; + for (int d=0; dwrite (origin, dims, (void*)tmp->storage()); + } + + /* Delete temporary copy */ + delete tmp; + } END_COMPONENT_LOOP(cgh); // Close the file @@ -481,14 +452,16 @@ namespace CarpetIOFlexIO { } extension = GetStringParameter ("in3D_extension", extension); - char* const filename = (char*)alloca(strlen(myindir)+strlen(alias)+strlen(extension)+100); + char* const filename = (char*)alloca(strlen(myindir) + +strlen(alias) + +strlen(extension)+100); sprintf (filename, "%s/%s.%s", myindir, alias, extension); IObase* reader = 0; AmrGridReader* amrreader = 0; int ndatasets = -1; - const int dim = CCTK_GroupDimI(group); + const int gpdim = CCTK_GroupDimI(group); // Read the file only on the root processor if (CCTK_MyProc(cgh)==0) { @@ -511,11 +484,11 @@ namespace CarpetIOFlexIO { // Read information about dataset IObase::DataType numbertype; int rank; - int dims[maxdim]; + int dims[dim]; reader->readInfo (numbertype, rank, dims); // Check rank - assert (rank==dim); + assert (rank==gpdim); // Check datatype // TODO: Check datatype correctly @@ -539,8 +512,8 @@ namespace CarpetIOFlexIO { // Read grid AmrGrid* amrgrid = 0; - int amr_origin[maxdim]; - int amr_dims[maxdim]; + int amr_origin[dim]; + int amr_dims[dim]; if (CCTK_MyProc(cgh)==0) { @@ -554,15 +527,19 @@ namespace CarpetIOFlexIO { IObase::DataType atype; int alength; if (reader->readAttributeInfo("iorigin", atype, alength) < 0) { - for (int d=0; diorigin[d] = 0; } } - for (int d=0; diorigin[d]; amr_dims[d] = amrgrid->dims[d]; } + for (int d=gpdim; d* ff = 0; \ - \ - switch (CCTK_GroupTypeI(group)) { \ - \ - case CCTK_ARRAY: \ - assert (var < (int)arrdata[group].data.size()); \ - ff = (generic_gf*)arrdata[group].data[var]; \ - break; \ - \ - case CCTK_GF: \ - assert (var < (int)gfdata[group].data.size()); \ - ff = (generic_gf*)gfdata[group].data[var]; \ - break; \ - \ - default: \ - abort(); \ - } \ - \ - generic_data* const data \ - = (*ff) (tl, reflevel, component, mglevel); \ - \ - /* Create temporary data storage on processor 0 */ \ - const vect str = vect(reflevelfact); \ - const vect lb = vect(amr_origin) * str; \ - const vect ub \ - = lb + (vect(amr_dims) - 1) * str; \ - const bbox ext(lb,ub,str); \ - generic_data* const tmp = data->make_typed (); \ - \ - if (CCTK_MyProc(cgh)==0) { \ - tmp->allocate (ext, 0, amrgrid->data); \ - } else { \ - tmp->allocate (ext, 0, 0); \ - } \ - \ - /* Copy into grid function */ \ - data->copy_from (tmp, ext); \ - \ - /* Delete temporary copy */ \ - delete tmp; \ - \ - } while (0) - - case 1: { - const int dim=1; - CODE; - break; - } + generic_gf* ff = 0; + + switch (CCTK_GroupTypeI(group)) { - case 2: { - const int dim=2; - CODE; + case CCTK_ARRAY: + assert (var < (int)arrdata[group].data.size()); + ff = (generic_gf*)arrdata[group].data[var]; break; - } - case 3: { - const int dim=3; - CODE; + case CCTK_GF: + assert (var < (int)gfdata[group].data.size()); + ff = (generic_gf*)gfdata[group].data[var]; break; - } default: abort(); - -#undef CODE - } + generic_data* const data + = (*ff) (tl, reflevel, component, mglevel); + + /* Create temporary data storage on processor 0 */ + const vect str = vect(reflevelfact); + const vect lb = vect(amr_origin) * str; + const vect ub + = lb + (vect(amr_dims) - 1) * str; + const bbox ext(lb,ub,str); + generic_data* const tmp = data->make_typed (); + + if (CCTK_MyProc(cgh)==0) { + tmp->allocate (ext, 0, amrgrid->data); + } else { + tmp->allocate (ext, 0, 0); + } + + /* Copy into grid function */ + data->copy_from (tmp, ext); + + /* Delete temporary copy */ + delete tmp; + } END_COMPONENT_LOOP(cgh); if (CCTK_MyProc(cgh)==0) { -- cgit v1.2.3