diff options
| author | Erik Schnetter <schnetter@cct.lsu.edu> | 2010-03-18 15:10:52 -0700 |
|---|---|---|
| committer | Barry Wardell <barry.wardell@gmail.com> | 2011-12-14 16:45:33 +0000 |
| commit | 0707cfeb82f565afec604f39202dafcbef889db0 (patch) | |
| tree | 820f1ac0bc04d5f9827e4e00bf1efccb6d76727e | |
| parent | 6b2d5314e79a5547f2fdd1dc4898bdeef74da9f2 (diff) | |
Re-organise time level handling
Store the current Cactus time (and not a fake Carpet time) in the th
"time hiearchy". This removes the now redundant "leveltimes" data
structure in Carpet.
Add past time levels to th, so that it can store the time for past
time levels instead of assuming the time step size is constant. This
allows changing the time step size during evolution.
Share the time hierarchy between all maps, instead of having one time
hierarchy per map.
Simplify the time level cycling and time stepping code used during
evolution.
Improve structure of the code that loops over time levels for certain
schedule bins. Introduce a new Carpet variable "timelevel", similar
to "reflevel".
This also makes it possible to avoid time interpolation for the past
time levels during regridding. The past time levels of the fine grid
then remain aligned (in time) with the past time levels of the coarse
grid. This is controlled by a new parameter
"time_interpolation_during_regridding", which defaults to "yes" for
backwards compatibility.
Simplify the three time level initialisation. Instead of initialising
all three time levels by taking altogether three time steps (forwards
and backwards), initialise only one past time level by taking one time
step backwards. The remaining time level is initialised during the
first time step of the evolution, which begins by cycling time levels,
which drops the non-initialised last time level anyway.
Update Carpet and the mode handling correspondingly.
Update the CarpetIOHDF5 checkpoint format correspondingly.
Update CarpetInterp, CarpetReduce, and CarpetRegrid2 correspondingly.
Update CarpetJacobi and CarpetMG correspondingly.
29 files changed, 651 insertions, 745 deletions
diff --git a/Carpet/Carpet/param.ccl b/Carpet/Carpet/param.ccl index 1cbabcfc0..d08cb98aa 100644 --- a/Carpet/Carpet/param.ccl +++ b/Carpet/Carpet/param.ccl @@ -494,6 +494,10 @@ BOOLEAN regrid_in_level_mode "Regrid in level mode (instead of singlemap mode), { } "yes" +BOOLEAN time_interpolation_during_regridding "Interpolate finer levels in time during regridding" STEERABLE=recover +{ +} "yes" + BOOLEAN init_3_timelevels "Set up 3 timelevels of initial data" STEERABLE=always diff --git a/Carpet/Carpet/src/Comm.cc b/Carpet/Carpet/src/Comm.cc index 64f2f0480..c01a59fa3 100644 --- a/Carpet/Carpet/src/Comm.cc +++ b/Carpet/Carpet/src/Comm.cc @@ -22,7 +22,6 @@ namespace Carpet { static void ProlongateGroupBoundaries (const cGH* cctkGH, - CCTK_REAL initial_time, const vector<int>& groups); @@ -137,7 +136,7 @@ namespace Carpet { if (local_do_prolongate) { static Timer timer ("Evolve::Prolongate"); timer.start(); - ProlongateGroupBoundaries (cctkGH, cctk_initial_time, goodgroups); + ProlongateGroupBoundaries (cctkGH, goodgroups); timer.stop(); } @@ -156,11 +155,9 @@ namespace Carpet { // Prolongate the boundaries of all CCTK_GF groups in the given set static void ProlongateGroupBoundaries (const cGH* cctkGH, - const CCTK_REAL initial_time, const vector<int>& groups) { DECLARE_CCTK_PARAMETERS; - const int tl = 0; if (reflevel == 0) return; @@ -169,8 +166,7 @@ namespace Carpet { assert (groups.size() > 0); // use the current time here (which may be modified by the user) - const CCTK_REAL time = - (cctkGH->cctk_time - initial_time) / delta_time; + const CCTK_REAL time = cctkGH->cctk_time; for (comm_state state; not state.done(); state.step()) { for (int group = 0; group < (int)groups.size(); ++group) { @@ -180,11 +176,16 @@ namespace Carpet { continue; } assert (reflevel>=0 and reflevel<reflevels); + const int active_tl = + groupdata.AT(g).activetimelevels.AT(mglevel).AT(reflevel); + assert (active_tl>=0); + const int tl = active_tl > 1 ? timelevel : 0; for (int m = 0; m < (int)arrdata.AT(g).size(); ++m) { for (int v = 0; v < (int)arrdata.AT(g).AT(m).data.size(); ++v) { ggf *const gv = arrdata.AT(g).AT(m).data.AT(v); - gv->ref_bnd_prolongate_all (state, tl, reflevel, mglevel, time); + gv->ref_bnd_prolongate_all + (state, tl, reflevel, mglevel, time); } } } @@ -196,7 +197,6 @@ namespace Carpet { void SyncGroups (const cGH* cctkGH, const vector<int>& groups) { DECLARE_CCTK_PARAMETERS; - const int tl = 0; Checkpoint ("SyncGroups"); @@ -208,6 +208,10 @@ namespace Carpet { const int grouptype = CCTK_GroupTypeI (g); const int ml = grouptype == CCTK_GF ? mglevel : 0; const int rl = grouptype == CCTK_GF ? reflevel : 0; + const int active_tl = + groupdata.AT(g).activetimelevels.AT(mglevel).AT(reflevel); + assert (active_tl>=0); + const int tl = active_tl > 1 ? timelevel : 0; for (int m = 0; m < (int)arrdata.AT(g).size(); ++m) { for (int v = 0; v < (int)arrdata.AT(g).AT(m).data.size(); ++v) { arrdesc& array = arrdata.AT(g).AT(m); diff --git a/Carpet/Carpet/src/Cycle.cc b/Carpet/Carpet/src/Cycle.cc index acf4d19db..582a966e5 100644 --- a/Carpet/Carpet/src/Cycle.cc +++ b/Carpet/Carpet/src/Cycle.cc @@ -17,17 +17,21 @@ namespace Carpet { - void CycleTimeLevels (const cGH* cgh) + void CycleTimeLevels (cGH* const cctkGH) { DECLARE_CCTK_PARAMETERS; Checkpoint ("CycleTimeLevels"); assert (is_level_mode()); + assert (timelevel == 0); + tt->advance_time (mglevel, reflevel); + cctkGH->cctk_time = tt->get_time (mglevel, reflevel, timelevel); + for (int group=0; group<CCTK_NumGroups(); ++group) { - if (CCTK_QueryGroupStorageI(cgh, group)) { + if (CCTK_QueryGroupStorageI(cctkGH, group)) { - int const activetimelevels = CCTK_ActiveTimeLevelsGI (cgh, group); + int const activetimelevels = CCTK_ActiveTimeLevelsGI (cctkGH, group); if (activetimelevels > 1) { if (activetimelevels < prolongation_order_time+1) { char * const groupname = CCTK_GroupName (group); @@ -62,7 +66,7 @@ namespace Carpet { { int const varindex = firstvarindex + var; for (int tl=0; tl<numtimelevels; ++tl) { - cgh->data[varindex][tl] + cctkGH->data[varindex][tl] = (tl < groupdata.AT(group).info.activetimelevels ? ((*arrdata.AT(group).AT(0).data.AT(var)) (tl, 0, 0, 0)->storage()) @@ -82,17 +86,91 @@ namespace Carpet { - void FlipTimeLevels (const cGH* cgh) + void UncycleTimeLevels (cGH* const cctkGH) + { + DECLARE_CCTK_PARAMETERS; + + Checkpoint ("UncycleTimeLevels"); + assert (is_level_mode()); + + assert (timelevel == 0); + tt->retreat_time (mglevel, reflevel); + cctkGH->cctk_time = tt->get_time (mglevel, reflevel, timelevel); + + for (int group=0; group<CCTK_NumGroups(); ++group) { + if (CCTK_QueryGroupStorageI(cctkGH, group)) { + + int const activetimelevels = CCTK_ActiveTimeLevelsGI (cctkGH, group); + if (activetimelevels > 1) { + if (activetimelevels < prolongation_order_time+1) { + char * const groupname = CCTK_GroupName (group); + CCTK_VWarn (CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING, + "Group \"%s\" has %d only active time levels. Groups with more than one active time level need at least %d active time levels for prolongation_order_time=%d", + groupname, + activetimelevels, int(prolongation_order_time+1), + int(prolongation_order_time)); + free (groupname); + } + } + + switch (CCTK_GroupTypeI(group)) { + + case CCTK_GF: + assert (reflevel>=0 and reflevel<reflevels); + for (int m=0; m<(int)arrdata.AT(group).size(); ++m) { + for (int var=0; var<CCTK_NumVarsInGroupI(group); ++var) { + arrdata.AT(group).AT(m).data.AT(var)-> + uncycle_all (reflevel, mglevel); + } + } + break; + + case CCTK_SCALAR: + case CCTK_ARRAY: + if (do_global_mode) { + int const numtimelevels = CCTK_NumTimeLevelsI (group); + int const firstvarindex = CCTK_FirstVarIndexI (group); + for (int var=0; var<CCTK_NumVarsInGroupI(group); ++var) { + arrdata.AT(group).AT(0).data.AT(var)->uncycle_all (0, mglevel); + { + int const varindex = firstvarindex + var; + for (int tl=0; tl<numtimelevels; ++tl) { + cctkGH->data[varindex][tl] + = (tl < groupdata.AT(group).info.activetimelevels + ? ((*arrdata.AT(group).AT(0).data.AT(var)) + (tl, 0, 0, 0)->storage()) + : NULL); + } + } + } + } + break; + + default: + assert (0); + } // switch grouptype + } // if storage + } // for group + } + + + + void FlipTimeLevels (cGH* const cctkGH) { DECLARE_CCTK_PARAMETERS; Checkpoint ("FlipTimeLevels"); assert (is_level_mode()); + assert (timelevel == 0); + tt->flip_timelevels (mglevel, reflevel); + cctkGH->cctk_time = tt->get_time (mglevel, reflevel, timelevel); + cctkGH->cctk_delta_time *= -1; + for (int group=0; group<CCTK_NumGroups(); ++group) { - if (CCTK_QueryGroupStorageI(cgh, group)) { + if (CCTK_QueryGroupStorageI(cctkGH, group)) { - int const activetimelevels = CCTK_ActiveTimeLevelsGI (cgh, group); + int const activetimelevels = CCTK_ActiveTimeLevelsGI (cctkGH, group); if (activetimelevels > 1) { if (activetimelevels < prolongation_order_time+1) { char * const groupname = CCTK_GroupName (group); @@ -127,7 +205,7 @@ namespace Carpet { { int const varindex = firstvarindex + var; for (int tl=0; tl<numtimelevels; ++tl) { - cgh->data[varindex][tl] + cctkGH->data[varindex][tl] = (tl < groupdata.AT(group).info.activetimelevels ? ((*arrdata.AT(group).AT(0).data.AT(var)) (tl, 0, 0, 0)->storage()) @@ -147,13 +225,13 @@ namespace Carpet { - void FillTimeLevels (const cGH* cgh) + void FillTimeLevels (const cGH* const cctkGH) { Checkpoint ("FillTimeLevels"); assert (is_level_mode()); for (int group=0; group<CCTK_NumGroups(); ++group) { - if (CCTK_QueryGroupStorageI(cgh, group)) { + if (CCTK_QueryGroupStorageI(cctkGH, group)) { switch (CCTK_GroupTypeI(group)) { case CCTK_GF: diff --git a/Carpet/Carpet/src/Evolve.cc b/Carpet/Carpet/src/Evolve.cc index baadd19bd..3ed5ae848 100644 --- a/Carpet/Carpet/src/Evolve.cc +++ b/Carpet/Carpet/src/Evolve.cc @@ -95,7 +95,9 @@ namespace Carpet { int const do_every = ipow (mgfact, ml) * (maxtimereflevelfact / timereffacts.AT(rl)); if (cctkGH->cctk_iteration % do_every == 0) { - assert (abs (leveltimes.AT(ml).AT(rl) - global_time) <= + // assert (abs (leveltimes.AT(ml).AT(rl) - global_time) <= + // eps * global_time); + assert (abs (tt->get_time(ml,rl,0) - global_time) <= eps * global_time); } } @@ -296,57 +298,19 @@ namespace Carpet { do_global_mode = do_late_global_mode; do_meta_mode = do_late_meta_mode; - Waypoint ("Postregrid at iteration %d time %g%s%s", - cctkGH->cctk_iteration, (double)cctkGH->cctk_time, - (do_global_mode ? " (global)" : ""), - (do_meta_mode ? " (meta)" : "")); - - int const num_tl = prolongation_order_time+1; - - bool const old_do_allow_past_timelevels = - do_allow_past_timelevels; - do_allow_past_timelevels = false; - - // Rewind times - for (int m=0; m<maps; ++m) { - CCTK_REAL const old_delta = - vtt.AT(m)->get_delta (reflevel, mglevel); - vtt.AT(m)->set_delta (reflevel, mglevel, - old_delta); - } - FlipTimeLevels (cctkGH); - for (int tl=0; tl<num_tl; ++tl) { - for (int m=0; m<maps; ++m) { - vtt.AT(m)->advance_time (reflevel, mglevel); - } - CycleTimeLevels (cctkGH); - } - for (int m=0; m<maps; ++m) { - CCTK_REAL const old_delta = - vtt.AT(m)->get_delta (reflevel, mglevel); - vtt.AT(m)->set_delta (reflevel, mglevel, - old_delta); - } - FlipTimeLevels (cctkGH); - CCTK_REAL const old_cctk_time = cctkGH->cctk_time; - cctkGH->cctk_time -= - num_tl * (cctkGH->cctk_delta_time / cctkGH->cctk_timefac); - - for (int tl=0; tl<num_tl; ++tl) { + BEGIN_TIMELEVEL_LOOP (cctkGH) { - // Advance times - for (int m=0; m<maps; ++m) { - vtt.AT(m)->advance_time (reflevel, mglevel); - } - CycleTimeLevels (cctkGH); - cctkGH->cctk_time += - cctkGH->cctk_delta_time / cctkGH->cctk_timefac; + Waypoint ("Postregrid at iteration %d time %g timelevel %d%s%s", + cctkGH->cctk_iteration, + (double)cctkGH->cctk_time, + timelevel, + (do_global_mode ? " (global)" : ""), + (do_meta_mode ? " (meta)" : "")); // Postregrid ScheduleTraverse (where, "CCTK_POSTREGRID", cctkGH); - } // for tl - cctkGH->cctk_time = old_cctk_time; - - do_allow_past_timelevels = old_do_allow_past_timelevels; + } END_TIMELEVEL_LOOP; EndTimingLevel (cctkGH); } LEAVE_LEVEL_MODE; @@ -413,43 +377,14 @@ namespace Carpet { } // Advance times - cctkGH->cctk_time - = (global_time - - delta_time / maxtimereflevelfact - + delta_time * mglevelfact / timereflevelfact); - CCTK_REAL const carpet_time = cctkGH->cctk_time / delta_time; - for (int m=0; m<maps; ++m) { - vtt.AT(m)->advance_time (reflevel, mglevel); - if (not adaptive_stepsize) { -#if 0 - // We must not perform this check, since the - // relative accuracy of incrementally adding to the - // current time cannot be good enough. Just setting - // the time (see below) is fine. - CCTK_REAL const eps = 1.0e-12; - static_assert (abs(0.1) > 0, - "Function CarpetLib::abs has wrong signature"); - CCTK_REAL const level_time = - vtt.AT(m)->get_time (reflevel, mglevel); - if (not (abs (level_time - carpet_time) <= - eps * max (carpet_time, 1.0))) { - int const oldprecision = cerr.precision(); - cerr.precision (17); - cerr << "ml: " << ml << endl - << "rl: " << rl << endl - << "m: " << m << endl - << "level_time: " << level_time << endl - << "carpet_time: " << carpet_time << endl - << "(level_time - carpet_time): " << (level_time - carpet_time) << endl; - cerr.precision (oldprecision); - } - assert (abs (level_time - carpet_time) <= - eps * max (carpet_time, 1.0)); -#endif - vtt.AT(m)->set_time (reflevel, mglevel, carpet_time); - } - } CycleTimeLevels (cctkGH); + if (not adaptive_stepsize) { + cctkGH->cctk_time + = (global_time + - delta_time / maxtimereflevelfact + + delta_time * mglevelfact / timereflevelfact); + tt->set_time (mglevel, reflevel, timelevel, cctkGH->cctk_time); + } Waypoint ("Evolution I at iteration %d time %g%s%s%s", cctkGH->cctk_iteration, (double)cctkGH->cctk_time, @@ -633,7 +568,7 @@ namespace Carpet { streamsize const oldprecision = cout.precision(); cout.precision (17); cout << " global_time: " << global_time << endl - << " leveltimes: " << leveltimes << endl + // << " leveltimes: " << leveltimes << endl << " delta_time: " << delta_time << endl; cout.precision (oldprecision); } diff --git a/Carpet/Carpet/src/Initialise.cc b/Carpet/Carpet/src/Initialise.cc index f81c0d368..1ff189dd7 100644 --- a/Carpet/Carpet/src/Initialise.cc +++ b/Carpet/Carpet/src/Initialise.cc @@ -64,8 +64,14 @@ namespace Carpet { global_time = cctk_initial_time; delta_time = 1.0; for (int ml = 0; ml < mglevel; ++ ml) { - assert (leveltimes.AT(ml).size() == 1); - leveltimes.AT(ml).AT(0) = global_time; + // assert (leveltimes.AT(ml).size() == 1); + // leveltimes.AT(ml).AT(0) = global_time; + for (int rl = 0; rl < reflevels; ++ rl) { + CCTK_REAL const dt = delta_time / timereffacts.AT(rl); + for (int tl = 0; tl < tt->timelevels; ++ tl) { + tt->set_time (ml, rl, tl, global_time - tl * dt); + } + } } cctkGH->cctk_iteration = 0; @@ -99,6 +105,10 @@ namespace Carpet { CallRecoverVariables (cctkGH); CallPostRecoverVariables (cctkGH); + // TODO: We should probably restrict here. + // CallRestrict (cctkGH); + // TODO: Should we also execute another bin after this? + // CallPostRestrictRecover (cctkGH); print_internal_data (); } else { @@ -263,47 +273,14 @@ namespace Carpet { do_global_mode = do_early_global_mode; // on first iteration, coarsest grid do_meta_mode = do_early_meta_mode; // on first iteration, coarsest grid - Waypoint ("Recovering II at iteration %d time %g%s%s", - cctkGH->cctk_iteration, (double)cctkGH->cctk_time, - (do_global_mode ? " (global)" : ""), - (do_meta_mode ? " (meta)" : "")); - - int const num_tl = prolongation_order_time+1; - - bool const old_do_allow_past_timelevels = do_allow_past_timelevels; - do_allow_past_timelevels = false; - - // Rewind times - for (int m=0; m<maps; ++m) { - CCTK_REAL const old_delta = - vtt.AT(m)->get_delta (reflevel, mglevel); - vtt.AT(m)->set_delta (reflevel, mglevel, - old_delta); - } - FlipTimeLevels (cctkGH); - for (int tl=0; tl<num_tl; ++tl) { - for (int m=0; m<maps; ++m) { - vtt.AT(m)->advance_time (reflevel, mglevel); - } - CycleTimeLevels (cctkGH); - } - for (int m=0; m<maps; ++m) { - CCTK_REAL const old_delta = - vtt.AT(m)->get_delta (reflevel, mglevel); - vtt.AT(m)->set_delta (reflevel, mglevel, - old_delta); - } - FlipTimeLevels (cctkGH); - CCTK_REAL const old_cctk_time = cctkGH->cctk_time; - cctkGH->cctk_time -= - num_tl * (cctkGH->cctk_delta_time / cctkGH->cctk_timefac); - - for (int tl=0; tl<num_tl; ++tl) { + BEGIN_TIMELEVEL_LOOP (cctkGH) { - // Advance times - for (int m=0; m<maps; ++m) { - vtt.AT(m)->advance_time (reflevel, mglevel); - } - CycleTimeLevels (cctkGH); - cctkGH->cctk_time += cctkGH->cctk_delta_time / cctkGH->cctk_timefac; + Waypoint ("Recovering II at iteration %d time %g timelevel %d%s%s", + cctkGH->cctk_iteration, + (double)cctkGH->cctk_time, + timelevel, + (do_global_mode ? " (global)" : ""), + (do_meta_mode ? " (meta)" : "")); // Post recover variables ScheduleTraverse (where, "CCTK_POST_RECOVER_VARIABLES", cctkGH); @@ -311,10 +288,7 @@ namespace Carpet { // Checking PoisonCheck (cctkGH, currenttime); - } // for tl - cctkGH->cctk_time = old_cctk_time; - - do_allow_past_timelevels = old_do_allow_past_timelevels; + } END_TIMELEVEL_LOOP; CheckChecksums (cctkGH, allbutcurrenttime); @@ -384,61 +358,36 @@ namespace Carpet { // Set up the grids ScheduleTraverse (where, "CCTK_BASEGRID", cctkGH); - int const num_tl = - init_each_timelevel ? prolongation_order_time+1 : 1; - - bool const old_do_allow_past_timelevels = do_allow_past_timelevels; - do_allow_past_timelevels = - not CCTK_EQUALS (initial_data_setup_method, "init_single_level"); - - // Rewind times - for (int m=0; m<maps; ++m) { - CCTK_REAL const old_delta = - vtt.AT(m)->get_delta (reflevel, mglevel); - vtt.AT(m)->set_delta (reflevel, mglevel, - old_delta); - } - FlipTimeLevels (cctkGH); - for (int tl=0; tl<num_tl; ++tl) { - for (int m=0; m<maps; ++m) { - vtt.AT(m)->advance_time (reflevel, mglevel); - } - CycleTimeLevels (cctkGH); - } - for (int m=0; m<maps; ++m) { - CCTK_REAL const old_delta = - vtt.AT(m)->get_delta (reflevel, mglevel); - vtt.AT(m)->set_delta (reflevel, mglevel, - old_delta); - } - FlipTimeLevels (cctkGH); - CCTK_REAL const old_cctk_time = cctkGH->cctk_time; - cctkGH->cctk_time -= - num_tl * (cctkGH->cctk_delta_time / cctkGH->cctk_timefac); - - for (int tl=0; tl<num_tl; ++tl) { + if (init_each_timelevel) { - // Advance times - for (int m=0; m<maps; ++m) { - vtt.AT(m)->advance_time (reflevel, mglevel); - } - CycleTimeLevels (cctkGH); - cctkGH->cctk_time += cctkGH->cctk_delta_time / cctkGH->cctk_timefac; + BEGIN_TIMELEVEL_LOOP (cctkGH) { + + // Set up the initial data + ScheduleTraverse (where, "CCTK_INITIAL", cctkGH); + ScheduleTraverse (where, "CCTK_POSTINITIAL", cctkGH); + + } END_TIMELEVEL_LOOP; + + } else { // not init_each_timelevel + + assert (do_allow_past_timelevels); + do_allow_past_timelevels = + not CCTK_EQUALS (initial_data_setup_method, "init_single_level"); // Set up the initial data ScheduleTraverse (where, "CCTK_INITIAL", cctkGH); ScheduleTraverse (where, "CCTK_POSTINITIAL", cctkGH); - if (init_fill_timelevels) { - assert (tl==0); - FillTimeLevels (cctkGH); - } - - // Checking - PoisonCheck (cctkGH, currenttime); + do_allow_past_timelevels = true; - } // for tl - cctkGH->cctk_time = old_cctk_time; + } // not init_each_timelevel - do_allow_past_timelevels = old_do_allow_past_timelevels; + if (init_fill_timelevels) { + FillTimeLevels (cctkGH); + } + + // Checking + PoisonCheck (cctkGH, currenttime); if (regrid_during_initialisation and mglevel==0) { // Regrid after initialising each level @@ -875,57 +824,19 @@ namespace Carpet { do_global_mode = do_late_global_mode; do_meta_mode = do_late_meta_mode; - Waypoint ("Postregrid at iteration %d time %g%s%s", - cctkGH->cctk_iteration, (double)cctkGH->cctk_time, - (do_global_mode ? " (global)" : ""), - (do_meta_mode ? " (meta)" : "")); - - int const num_tl = prolongation_order_time + 1; - - bool const old_do_allow_past_timelevels = - do_allow_past_timelevels; - do_allow_past_timelevels = false; - - // Rewind times - for (int m=0; m<maps; ++m) { - CCTK_REAL const old_delta = - vtt.AT(m)->get_delta (reflevel, mglevel); - vtt.AT(m)->set_delta (reflevel, mglevel, - old_delta); - } - FlipTimeLevels (cctkGH); - for (int tl=0; tl<num_tl; ++tl) { - for (int m=0; m<maps; ++m) { - vtt.AT(m)->advance_time (reflevel, mglevel); - } - CycleTimeLevels (cctkGH); - } - for (int m=0; m<maps; ++m) { - CCTK_REAL const old_delta = - vtt.AT(m)->get_delta (reflevel, mglevel); - vtt.AT(m)->set_delta (reflevel, mglevel, - old_delta); - } - FlipTimeLevels (cctkGH); - CCTK_REAL const old_cctk_time = cctkGH->cctk_time; - cctkGH->cctk_time -= - num_tl * (cctkGH->cctk_delta_time / cctkGH->cctk_timefac); + BEGIN_TIMELEVEL_LOOP (cctkGH) { - for (int tl=0; tl<num_tl; ++tl) { - - // Advance times - for (int m=0; m<maps; ++m) { - vtt.AT(m)->advance_time (reflevel, mglevel); - } - CycleTimeLevels (cctkGH); - cctkGH->cctk_time += - cctkGH->cctk_delta_time / cctkGH->cctk_timefac; + Waypoint ("Postregrid at iteration %d time %g timelevel %d%s%s", + cctkGH->cctk_iteration, + (double)cctkGH->cctk_time, + timelevel, + (do_global_mode ? " (global)" : ""), + (do_meta_mode ? " (meta)" : "")); // Postregrid ScheduleTraverse (where, "CCTK_POSTREGRID", cctkGH); - } // for tl - cctkGH->cctk_time = old_cctk_time; - - do_allow_past_timelevels = old_do_allow_past_timelevels; + } END_TIMELEVEL_LOOP; EndTimingLevel (cctkGH); } LEAVE_LEVEL_MODE; @@ -1079,53 +990,26 @@ namespace Carpet { (do_global_mode ? " (global)" : ""), (do_meta_mode ? " (meta)" : "")); - int const num_tl = - init_each_timelevel ? prolongation_order_time+1 : 1; - - bool const old_do_allow_past_timelevels = - do_allow_past_timelevels; - do_allow_past_timelevels = false; - - // Rewind times - for (int m=0; m<maps; ++m) { - CCTK_REAL const old_delta = - vtt.AT(m)->get_delta (reflevel, mglevel); - vtt.AT(m)->set_delta (reflevel, mglevel, - old_delta); - } - FlipTimeLevels (cctkGH); - for (int tl=0; tl<num_tl; ++tl) { - for (int m=0; m<maps; ++m) { - vtt.AT(m)->advance_time (reflevel, mglevel); - } - CycleTimeLevels (cctkGH); - } - for (int m=0; m<maps; ++m) { - CCTK_REAL const old_delta = - vtt.AT(m)->get_delta (reflevel, mglevel); - vtt.AT(m)->set_delta (reflevel, mglevel, - old_delta); - } - FlipTimeLevels (cctkGH); - CCTK_REAL const old_cctk_time = cctkGH->cctk_time; - cctkGH->cctk_time -= - num_tl * (cctkGH->cctk_delta_time / cctkGH->cctk_timefac); - - for (int tl=0; tl<num_tl; ++tl) { + if (init_each_timelevel) { - // Advance times - for (int m=0; m<maps; ++m) { - vtt.AT(m)->advance_time (reflevel, mglevel); - } - CycleTimeLevels (cctkGH); - cctkGH->cctk_time += - cctkGH->cctk_delta_time / cctkGH->cctk_timefac; + BEGIN_TIMELEVEL_LOOP (cctkGH) { + + // Postregrid + ScheduleTraverse (where, "CCTK_POSTREGRIDINITIAL", cctkGH); + + } END_TIMELEVEL_LOOP; + + } else { // not init_each_timelevel + + assert (do_allow_past_timelevels); + do_allow_past_timelevels = false; // Postregrid ScheduleTraverse (where, "CCTK_POSTREGRIDINITIAL", cctkGH); - } // for tl - cctkGH->cctk_time = old_cctk_time; - - do_allow_past_timelevels = old_do_allow_past_timelevels; + do_allow_past_timelevels = true; + + } // not init_each_timelevel EndTimingLevel (cctkGH); } LEAVE_LEVEL_MODE; @@ -1156,274 +1040,107 @@ namespace Carpet { - // Use Scott Hawley's algorithm to get two extra timelevels of data - - static void initialise_3tl_advance_time (cGH * const cctkGH); - static void initialise_3tl_evolve_Ia (cGH * const cctkGH); static void initialise_3tl_flip_timelevels (cGH * const cctkGH); - static void initialise_3tl_evolve_Ib (cGH * const cctkGH); - static void initialise_3tl_evolve_IIb (cGH * const cctkGH); - static void initialise_3tl_advance_time_2 (cGH * const cctkGH); - static void initialise_3tl_evolve_Ic (cGH * const cctkGH); - static void initialise_3tl_reset_time (cGH * const cctkGH); + static void initialise_3tl_evolve (cGH * const cctkGH); + static void initialise_3tl_recycle (cGH * const cctkGH); + + void Initialise3tl (cGH * const cctkGH) { DECLARE_CCTK_PARAMETERS; - Waypoint ("Initialising three timelevels"); + Waypoint ("Initialising three timelevels:"); - // TODO: ensure that there are 3 timelevels - if (not (prolongation_order_time == 2)) { - CCTK_WARN (CCTK_WARN_ABORT, - "The 3 timelevel initialisation scheme (init_3_timelevels=yes) requires 3 timelevels and a temporal prolongation order of 2 (prolongation_order_time=2)"); - } - assert (prolongation_order_time == 2); + initialise_3tl_flip_timelevels (cctkGH); + initialise_3tl_evolve (cctkGH); + // TODO: May want to restrict here if possible (i.e. if the time + // refinement factor is one) + initialise_3tl_recycle (cctkGH); + initialise_3tl_flip_timelevels (cctkGH); - for (int rl=0; rl<int(timereffacts.size()); ++rl) { - if (not (timereffacts.AT(rl) == ipow (2, rl))) { - CCTK_WARN (CCTK_WARN_ABORT, - "The 3 timelevel initialisation scheme (init_3_timelevels=yes) requires temporal refinement factors of 2 for all refinement levels (time_refinement_factors[rl]=pow(2,rl))"); - } - assert (timereffacts.AT(rl) == ipow (2, rl)); - } + Waypoint ("Finished initialising three timelevels"); + } + + void + initialise_3tl_flip_timelevels (cGH * const cctkGH) + { + Waypoint ("Initialise3TL::Flip"); + + delta_time *= -1; BEGIN_MGLEVEL_LOOP(cctkGH) { - BEGIN_REFLEVEL_LOOP(cctkGH) { - BeginTimingLevel (cctkGH); + BEGIN_REFLEVEL_LOOP (cctkGH) { - do_early_global_mode = reflevel==0; - do_late_global_mode = reflevel==reflevels-1; - do_early_meta_mode = do_early_global_mode and mglevel==mglevels-1; - do_late_meta_mode = do_late_global_mode and mglevel==0; - do_global_mode = do_early_global_mode; - do_meta_mode = do_early_meta_mode; + FlipTimeLevels (cctkGH); - initialise_3tl_advance_time (cctkGH); - initialise_3tl_evolve_Ia (cctkGH); - initialise_3tl_flip_timelevels (cctkGH); - initialise_3tl_evolve_Ib (cctkGH); - initialise_3tl_flip_timelevels (cctkGH); - - EndTimingLevel (cctkGH); } END_REFLEVEL_LOOP; } END_MGLEVEL_LOOP; - - Waypoint ("Hourglass structure in place"); - - initialise_3tl_flip_timelevels (cctkGH); + } + + void + initialise_3tl_evolve (cGH * const cctkGH) + { + char const * const where = "Initialise3TL::Evolve"; + static Timer timer (where); + timer.start(); BEGIN_MGLEVEL_LOOP(cctkGH) { - BEGIN_REVERSE_REFLEVEL_LOOP(cctkGH) { + BEGIN_REFLEVEL_LOOP(cctkGH) { BeginTimingLevel (cctkGH); - do_early_global_mode = reflevel==reflevels-1; - do_late_global_mode = reflevel==0; + do_early_global_mode = reflevel==0; + do_late_global_mode = reflevel==reflevels-1; do_early_meta_mode = do_early_global_mode and mglevel==mglevels-1; do_late_meta_mode = do_late_global_mode and mglevel==0; do_global_mode = do_early_global_mode; do_meta_mode = do_early_meta_mode; - initialise_3tl_evolve_IIb (cctkGH); - initialise_3tl_advance_time_2 (cctkGH); - initialise_3tl_evolve_Ic (cctkGH); + Waypoint ("Initialisation 3TL evolution", + (do_global_mode ? " (global)" : ""), + (do_meta_mode ? " (meta)" : "")); - EndTimingLevel (cctkGH); - } END_REVERSE_REFLEVEL_LOOP; - } END_MGLEVEL_LOOP; - - initialise_3tl_flip_timelevels (cctkGH); - - BEGIN_MGLEVEL_LOOP(cctkGH) { - BEGIN_REVERSE_REFLEVEL_LOOP(cctkGH) { - BeginTimingLevel (cctkGH); + CycleTimeLevels (cctkGH); - do_early_global_mode = reflevel==reflevels-1; - do_late_global_mode = reflevel==0; - do_early_meta_mode = do_early_global_mode and mglevel==mglevels-1; - do_late_meta_mode = do_late_global_mode and mglevel==0; - do_global_mode = do_early_global_mode; - do_meta_mode = do_early_meta_mode; + CalculateChecksums (cctkGH, allbutcurrenttime); + Poison (cctkGH, currenttimebutnotifonly); + + // Evolve + ScheduleTraverse (where, "CCTK_PRESTEP", cctkGH); + ScheduleTraverse (where, "CCTK_EVOL", cctkGH); - initialise_3tl_reset_time (cctkGH); + PoisonCheck (cctkGH, currenttime); EndTimingLevel (cctkGH); - } END_REVERSE_REFLEVEL_LOOP; + } END_REFLEVEL_LOOP; } END_MGLEVEL_LOOP; - Waypoint ("Finished initialising three timelevels"); - } - - - - void - initialise_3tl_advance_time (cGH * const cctkGH) - { - Waypoint ("Advancing time"); - - cctkGH->cctk_time - = global_time + delta_time * mglevelfact / timereflevelfact; - for (int m=0; m<maps; ++m) { - vtt.AT(m)->advance_time (reflevel, mglevel); - } - - CycleTimeLevels (cctkGH); - } - - void - initialise_3tl_evolve_Ia (cGH * const cctkGH) - { - char const * const where = "Initialise3TL::EvolveIa"; - static Timer timer (where); - timer.start(); - - Waypoint ("Initialisation 3TL evolution I (a) (forwards) at iteration" - " %d time %g%s%s", - cctkGH->cctk_iteration, (double)cctkGH->cctk_time, - (do_global_mode ? " (global)" : ""), - (do_meta_mode ? " (meta)" : "")); - - CalculateChecksums (cctkGH, allbutcurrenttime); - Poison (cctkGH, currenttimebutnotifonly); - - // Evolve forward - ScheduleTraverse (where, "CCTK_PRESTEP", cctkGH); - ScheduleTraverse (where, "CCTK_EVOL", cctkGH); - - PoisonCheck (cctkGH, currenttime); - - timer.stop(); - } - - void - initialise_3tl_flip_timelevels (cGH * const cctkGH) - { - Waypoint ("Flipping timelevels"); - - BEGIN_META_MODE(cctkGH) { - - delta_time *= -1; - - BEGIN_MGLEVEL_LOOP(cctkGH) { - BEGIN_REFLEVEL_LOOP (cctkGH) { - - cctkGH->cctk_time - = global_time + delta_time * mglevelfact / timereflevelfact; - - FlipTimeLevels (cctkGH); - - } END_REFLEVEL_LOOP; - } END_MGLEVEL_LOOP; - } END_META_MODE; - } - - void - initialise_3tl_evolve_Ib (cGH * const cctkGH) - { - char const * const where = "Initialise3TL::EvolveIb"; - static Timer timer (where); - timer.start(); - - Waypoint ("Initialisation 3TL evolution I (b) (backwards) at iteration" - " %d time %g%s%s", - cctkGH->cctk_iteration, (double)cctkGH->cctk_time, - (do_global_mode ? " (global)" : ""), - (do_meta_mode ? " (meta)" : "")); - - // Checking - CalculateChecksums (cctkGH, allbutcurrenttime); - Poison (cctkGH, currenttimebutnotifonly); - - // Evolve backward - ScheduleTraverse (where, "CCTK_PRESTEP", cctkGH); - ScheduleTraverse (where, "CCTK_EVOL", cctkGH); - - // Checking - PoisonCheck (cctkGH, alltimes); - timer.stop(); } - // Evolve backwards one more timestep - // Starting with the finest level and proceeding to the coarsest void - initialise_3tl_evolve_IIb (cGH * const cctkGH) + initialise_3tl_recycle (cGH * const cctkGH) { - char const * const where = "Initialise3TL::EvolveIIb"; + char const * const where = "Initialise3TL::Recycle"; static Timer timer (where); timer.start(); - Waypoint ("Initialisation 3TL evolution II (b) (backwards) at iteration" - " %d time %g%s%s", - cctkGH->cctk_iteration, (double)cctkGH->cctk_time, - (do_global_mode ? " (global)" : ""), - (do_meta_mode ? " (meta)" : "")); - - Restrict (cctkGH); - - if (reflevel < reflevels-1) { - ScheduleTraverse (where, "CCTK_POSTRESTRICT", cctkGH); - } - - ScheduleTraverse (where, "CCTK_POSTSTEP", cctkGH); - - PoisonCheck (cctkGH, alltimes); - - timer.stop(); - } - - void - initialise_3tl_advance_time_2 (cGH * const cctkGH) - { - Waypoint ("Advancing time"); - - cctkGH->cctk_time - = global_time + 2 * delta_time * mglevelfact / timereflevelfact; - for (int m=0; m<maps; ++m) { - vtt.AT(m)->advance_time (reflevel, mglevel); - } - - CycleTimeLevels (cctkGH); - } - - void - initialise_3tl_evolve_Ic (cGH * const cctkGH) - { - char const * const where = "Initialise3TL::EvolveIc"; - static Timer timer (where); - timer.start(); - - Waypoint ("Initialisation 3TL evolution I (c) (backwards) at iteration" - " %d time %g%s%s", - cctkGH->cctk_iteration, (double)cctkGH->cctk_time, - (do_global_mode ? " (global)" : ""), - (do_meta_mode ? " (meta)" : "")); - - CalculateChecksums (cctkGH, allbutcurrenttime); - Poison (cctkGH, currenttimebutnotifonly); - - // Evolve backward - ScheduleTraverse (where, "CCTK_PRESTEP", cctkGH); - ScheduleTraverse (where, "CCTK_EVOL", cctkGH); - ScheduleTraverse (where, "CCTK_POSTSTEP", cctkGH); - - PoisonCheck (cctkGH, alltimes); + BEGIN_MGLEVEL_LOOP(cctkGH) { + BEGIN_REFLEVEL_LOOP(cctkGH) { + BeginTimingLevel (cctkGH); + + Waypoint ("Initialisation 3TL recycling"); + + UncycleTimeLevels (cctkGH); + + EndTimingLevel (cctkGH); + } END_REFLEVEL_LOOP; + } END_MGLEVEL_LOOP; timer.stop(); } - void - initialise_3tl_reset_time (cGH * const cctkGH) - { - Waypoint ("Resetting time"); - - cctkGH->cctk_time = global_time; - for (int m=0; m<maps; ++m) { - vtt.AT(m)->set_time (reflevel, mglevel, 0); - } - } - ////////////////////////////////////////////////////////////////////////////// @@ -1442,7 +1159,7 @@ namespace Carpet { streamsize const oldprecision = cout.precision(); cout.precision (17); cout << " global_time: " << global_time << endl - << " leveltimes: " << leveltimes << endl + // << " leveltimes: " << leveltimes << endl << " delta_time: " << delta_time << endl; cout.precision (oldprecision); } diff --git a/Carpet/Carpet/src/Recompose.cc b/Carpet/Carpet/src/Recompose.cc index 598d48a14..94b1f90f7 100644 --- a/Carpet/Carpet/src/Recompose.cc +++ b/Carpet/Carpet/src/Recompose.cc @@ -356,10 +356,10 @@ namespace Carpet { } // Set new level times - for (int ml=0; ml<mglevels; ++ml) { - leveltimes.AT(ml).resize - (reflevels, leveltimes.AT(ml).AT(oldreflevels-1)); - } + // for (int ml=0; ml<mglevels; ++ml) { + // leveltimes.AT(ml).resize + // (reflevels, leveltimes.AT(ml).AT(oldreflevels-1)); + // } ++ regridding_epoch; diff --git a/Carpet/Carpet/src/Restrict.cc b/Carpet/Carpet/src/Restrict.cc index bad055641..8ede1ebcd 100644 --- a/Carpet/Carpet/src/Restrict.cc +++ b/Carpet/Carpet/src/Restrict.cc @@ -62,28 +62,16 @@ namespace Carpet { static void RestrictGroups (const cGH* cctkGH, const vector<int>& groups) { DECLARE_CCTK_PARAMETERS; - const int tl = 0; - for (comm_state state; not state.done(); state.step()) { for (int group = 0; group < (int)groups.size(); ++group) { const int g = groups.AT(group); + const int active_tl = CCTK_ActiveTimeLevelsGI (cctkGH, g); + assert (active_tl>=0); + const int tl = active_tl > 1 ? timelevel : 0; for (int m=0; m<(int)arrdata.AT(g).size(); ++m) { - - // use background time here (which may not be modified - // by the user) - const CCTK_REAL time = vtt.AT(m)->time (tl, reflevel, mglevel); - - const CCTK_REAL time1 = vtt.AT(m)->time (0, reflevel, mglevel); - const CCTK_REAL time2 = - (cctkGH->cctk_time - cctk_initial_time) / delta_time; - const CCTK_REAL time0 = - abs(time1) + abs(time2) + abs(cctkGH->cctk_delta_time); - const CCTK_REAL eps = 1.0e-12; - assert (abs(time1 - time2) <= eps * time0); - for (int v = 0; v < (int)arrdata.AT(g).AT(m).data.size(); ++v) { ggf *const gv = arrdata.AT(g).AT(m).data.AT(v); - gv->ref_restrict_all (state, tl, reflevel, mglevel, time); + gv->ref_restrict_all (state, tl, reflevel, mglevel); } } } // loop over groups diff --git a/Carpet/Carpet/src/SetupGH.cc b/Carpet/Carpet/src/SetupGH.cc index af5f13422..4961cec3b 100644 --- a/Carpet/Carpet/src/SetupGH.cc +++ b/Carpet/Carpet/src/SetupGH.cc @@ -54,8 +54,7 @@ namespace Carpet { allocate_data_hierarchy (cGH const * cctkGH, int m); static void - allocate_time_hierarchy (cGH const * cctkGH, - int m); + allocate_time_hierarchy (cGH const * cctkGH); static void setup_grid_hierarchy (cGH const * cctkGH); static void @@ -272,6 +271,9 @@ namespace Carpet { cctkGH->cctk_mode = CCTK_MODE_META; #endif + timelevels = prolongation_order_time + 1; + timelevel = 0; + // Say hello Waypoint ("Setting up the grid hierarchy"); @@ -448,9 +450,10 @@ namespace Carpet { // Allocate hierarchies allocate_grid_hierarchy (cctkGH, m); allocate_data_hierarchy (cctkGH, m); - allocate_time_hierarchy (cctkGH, m); } // for m + + allocate_time_hierarchy (cctkGH); } @@ -601,15 +604,14 @@ namespace Carpet { void - allocate_time_hierarchy (cGH const * const cctkGH, - int const m) + allocate_time_hierarchy (cGH const * const cctkGH) { DECLARE_CCTK_PARAMETERS; - vtt.resize (maps); - vtt.AT(m) = new th (* vhh.AT(m), - timereffacts, - 1.0); + // We are using the gh of the first map. This works because all + // maps have the same number of refinement levels. + tt = new th (* vhh.AT(0), timelevels, timereffacts, + time_interpolation_during_regridding); } @@ -774,7 +776,7 @@ namespace Carpet { for (int m=0; m<maps; ++m) { arrdata.AT(group).AT(m).hh = vhh.AT(m); arrdata.AT(group).AT(m).dd = vdd.AT(m); - arrdata.AT(group).AT(m).tt = vtt.AT(m); + arrdata.AT(group).AT(m).tt = tt; } break; @@ -865,9 +867,9 @@ namespace Carpet { new dh (*arrdata.AT(group).AT(m).hh, ghosts, buffers, my_prolongation_orders_space); - CCTK_REAL const basedelta = 1.0; arrdata.AT(group).AT(m).tt = - new th (*arrdata.AT(group).AT(m).hh, grouptimereffacts, basedelta); + new th (*arrdata.AT(group).AT(m).hh, timelevels, grouptimereffacts, + time_interpolation_during_regridding); } @@ -1053,11 +1055,11 @@ namespace Carpet { void set_state (cGH * const cctkGH) { - // Allocate level times - leveltimes.resize (mglevels); - for (int ml=0; ml<mglevels; ++ml) { - leveltimes.AT(ml).resize (1); - } + // // Allocate level times + // leveltimes.resize (mglevels); + // for (int ml=0; ml<mglevels; ++ml) { + // leveltimes.AT(ml).resize (1); + // } // Allocate orgin and spacings origin_space.resize (maps); diff --git a/Carpet/Carpet/src/Timing.cc b/Carpet/Carpet/src/Timing.cc index 51ecc45e8..43c230588 100644 --- a/Carpet/Carpet/src/Timing.cc +++ b/Carpet/Carpet/src/Timing.cc @@ -585,9 +585,9 @@ namespace Carpet { << " Grid hierarchy:" << eol; for (int m = 0; m < Carpet::maps; ++ m) { cout << " gh[" << m << "]: " << PRINTMEM(*vhh.AT(m)) << eol - << " dh[" << m << "]: " << PRINTMEM(*vdd.AT(m)) << eol - << " th[" << m << "]: " << PRINTMEM(*vtt.AT(m)) << eol; + << " dh[" << m << "]: " << PRINTMEM(*vdd.AT(m)) << eol; } + cout << " th: " << PRINTMEM(*tt) << eol; #if 0 for (int g = 0; g < (int)arrdata.size(); ++ g) { if (CCTK_GroupTypeI(g) != CCTK_GF) { diff --git a/Carpet/Carpet/src/carpet.hh b/Carpet/Carpet/src/carpet.hh index bbcb4296a..3ca8535fe 100644 --- a/Carpet/Carpet/src/carpet.hh +++ b/Carpet/Carpet/src/carpet.hh @@ -38,8 +38,9 @@ namespace Carpet { // Other functions bool Regrid (cGH const * cctkGH, bool force_recompose, bool do_init); - void CycleTimeLevels (const cGH* cgh); - void FlipTimeLevels (const cGH* cgh); + void CycleTimeLevels (cGH* cgh); + void UncycleTimeLevels (cGH* cgh); + void FlipTimeLevels (cGH* cgh); void FillTimeLevels (const cGH* cgh); void SyncGroups (const cGH* cgh, const vector<int>& groups); int SyncProlongateGroups (const cGH* cgh, const vector<int>& groups); diff --git a/Carpet/Carpet/src/modes.cc b/Carpet/Carpet/src/modes.cc index 68db7beaa..038815393 100644 --- a/Carpet/Carpet/src/modes.cc +++ b/Carpet/Carpet/src/modes.cc @@ -330,12 +330,17 @@ namespace Carpet { } // Set current time - assert (mglevel>=0 and mglevel<(int)leveltimes.size()); - assert (reflevel>=0 and reflevel<(int)leveltimes.AT(mglevel).size()); + // assert (mglevel>=0 and mglevel<(int)leveltimes.size()); + // assert (reflevel>=0 and reflevel<(int)leveltimes.AT(mglevel).size()); + // if (not adaptive_stepsize) { + // cctkGH->cctk_time = leveltimes.AT(mglevel).AT(reflevel); + // } else { + // leveltimes.AT(mglevel).AT(reflevel) = cctkGH->cctk_time; + // } if (not adaptive_stepsize) { - cctkGH->cctk_time = leveltimes.AT(mglevel).AT(reflevel); + cctkGH->cctk_time = tt->get_time(mglevel, reflevel, timelevel); } else { - leveltimes.AT(mglevel).AT(reflevel) = cctkGH->cctk_time; + tt->set_time (mglevel, reflevel, timelevel, cctkGH->cctk_time); } assert (is_level_mode()); @@ -354,9 +359,15 @@ namespace Carpet { CCTK_INT const deadbeef = get_deadbeef(); // Save and unset current time - assert (mglevel>=0 and mglevel<(int)leveltimes.size()); - assert (reflevel>=0 and reflevel<(int)leveltimes.AT(mglevel).size()); - leveltimes.AT(mglevel).AT(reflevel) = cctkGH->cctk_time; + // assert (mglevel>=0 and mglevel<(int)leveltimes.size()); + // assert (reflevel>=0 and reflevel<(int)leveltimes.AT(mglevel).size()); + // leveltimes.AT(mglevel).AT(reflevel) = cctkGH->cctk_time; + // if (not adaptive_stepsize) { + // cctkGH->cctk_time = global_time; + // } else { + // global_time = cctkGH->cctk_time; + // } + tt->set_time (mglevel, reflevel, timelevel, cctkGH->cctk_time); if (not adaptive_stepsize) { cctkGH->cctk_time = global_time; } else { @@ -606,10 +617,30 @@ namespace Carpet { assert (firstvar>=0); const int max_tl = CCTK_MaxTimeLevelsGI (group); assert (max_tl>=0); - const int active_tl = CCTK_ActiveTimeLevelsGI (cctkGH, group); + const int active_tl = + groupdata.AT(group).activetimelevels.AT(mglevel).AT(reflevel); assert (active_tl>=0 and active_tl<=max_tl); - const int available_tl = - do_allow_past_timelevels ? active_tl : min (1, active_tl); + int available_tl; + int tl_offset; + if (active_tl == 0) { + // gropu has no storage + available_tl = active_tl; + tl_offset = 0; + } else if (do_allow_past_timelevels) { + // regular case; all timelevels are accessible + available_tl = active_tl; + tl_offset = 0; + } else { + // only one timelevel is accessible + available_tl = 1; + if (timelevel < active_tl) { + // timelevel "timelevel" exists + tl_offset = timelevel; + } else { + // timelevel "timelevel" does not exist + tl_offset = active_tl - 1; + } + } // assert (vhh.AT(map)->is_local(reflevel,component)); @@ -620,7 +651,7 @@ namespace Carpet { for (int tl=0; tl<max_tl; ++tl) { if (ff and tl<available_tl) { gdata * const data = - (*ff) (tl, reflevel, local_component, mglevel); + (*ff) (tl_offset+tl, reflevel, local_component, mglevel); assert (data); cctkGH->data[firstvar+var][tl] = data->storage(); } else { diff --git a/Carpet/Carpet/src/modes.hh b/Carpet/Carpet/src/modes.hh index c2631d952..868c54172 100644 --- a/Carpet/Carpet/src/modes.hh +++ b/Carpet/Carpet/src/modes.hh @@ -247,6 +247,25 @@ namespace Carpet { assert (local_component_loop_); \ local_component_loop_ = false; \ } while (false) + +#define BEGIN_TIMELEVEL_LOOP(cctkGH) \ + do { \ + bool timelevel_loop_ = true; \ + assert (do_allow_past_timelevels); \ + do_allow_past_timelevels = false; \ + assert (timelevel == 0); \ + for (timelevel = 0; timelevel < timelevels; ++ timelevel) { \ + cctkGH->cctk_time = tt->get_time (mglevel, reflevel, timelevel); \ + { +#define END_TIMELEVEL_LOOP \ + } \ + } \ + assert (timelevel_loop_); \ + timelevel_loop_ = false; \ + timelevel = 0; \ + cctkGH->cctk_time = tt->get_time (mglevel, reflevel, timelevel); \ + do_allow_past_timelevels = true; \ + } while (false) diff --git a/Carpet/Carpet/src/variables.cc b/Carpet/Carpet/src/variables.cc index 0198a6018..c9d82cb16 100644 --- a/Carpet/Carpet/src/variables.cc +++ b/Carpet/Carpet/src/variables.cc @@ -42,6 +42,9 @@ namespace Carpet { // Maps int maps; + // Timelevels + int timelevels; + // Current position on the grid hierarchy @@ -51,6 +54,7 @@ namespace Carpet { int map; int component; int local_component; + int timelevel; // Current refinement factors int timereflevelfact; @@ -68,7 +72,7 @@ namespace Carpet { // Times and spaces on the refinement levels CCTK_REAL global_time; - vector<vector<CCTK_REAL> > leveltimes; // [mglevel][reflevel] + // vector<vector<CCTK_REAL> > leveltimes; // [mglevel][reflevel] CCTK_REAL delta_time; vector<vector<vect<CCTK_REAL,dim> > > origin_space; // [map][mglevel] @@ -109,7 +113,7 @@ namespace Carpet { // The grid hierarchy vector<gh*> vhh; // [map] vector<dh*> vdd; // [map] - vector<th*> vtt; // [map] + th* tt; int regridding_epoch; // Data for the groups @@ -122,6 +126,6 @@ namespace Carpet { // MPI Communicators MPI_Comm comm_universe = MPI_COMM_NULL; - MPI_Comm comm_world = MPI_COMM_NULL; + MPI_Comm comm_world = MPI_COMM_NULL; } // namespace Carpet diff --git a/Carpet/Carpet/src/variables.hh b/Carpet/Carpet/src/variables.hh index 2746ec295..61ac33cc2 100644 --- a/Carpet/Carpet/src/variables.hh +++ b/Carpet/Carpet/src/variables.hh @@ -72,6 +72,9 @@ namespace Carpet { // Maps extern int maps; + // Timelevels + extern int timelevels; + // Current position on the grid hierarchy @@ -81,6 +84,7 @@ namespace Carpet { extern int map; extern int component; extern int local_component; // -1 for non-local + extern int timelevel; // Current refinement factors extern int timereflevelfact; @@ -98,7 +102,7 @@ namespace Carpet { // Times and spaces on the refinement levels extern CCTK_REAL global_time; - extern vector<vector<CCTK_REAL> > leveltimes; // [mglevel][reflevel] + // extern vector<vector<CCTK_REAL> > leveltimes; // [mglevel][reflevel] extern CCTK_REAL delta_time; extern vector<vector<vect<CCTK_REAL,dim> > > origin_space; // [map][mglevel] @@ -146,7 +150,7 @@ namespace Carpet { // The grid hierarchy extern vector<gh*> vhh; // [map] extern vector<dh*> vdd; // [map] - extern vector<th*> vtt; // [map] + extern th* tt; extern int regridding_epoch; // Data for the groups diff --git a/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.cc b/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.cc index dcb2d97e4..944fbaa89 100644 --- a/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.cc +++ b/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.cc @@ -1091,14 +1091,14 @@ int WriteMetadata (const cGH * const cctkGH, int const nioprocs, "Cactus version", CCTK_FullVersion()); // all times on all refinement levels - error_count += WriteAttribute (group, - "numberofmgtimes", mglevels); - for (int i = 0; i < mglevels; i++) { - char name[100]; - snprintf (name, sizeof (name), "mgleveltimes %d", i); - error_count += WriteAttribute - (group, name, &leveltimes.at(i).front(), leveltimes.at(i).size()); - } + // error_count += WriteAttribute (group, + // "numberofmgtimes", mglevels); + // for (int i = 0; i < mglevels; i++) { + // char name[100]; + // snprintf (name, sizeof (name), "mgleveltimes %d", i); + // error_count += WriteAttribute + // (group, name, &leveltimes.at(i).front(), leveltimes.at(i).size()); + // } // unique configuration identifier if (CCTK_IsFunctionAliased ("UniqueConfigID")) { @@ -1196,21 +1196,23 @@ int WriteMetadata (const cGH * const cctkGH, int const nioprocs, if (called_from_checkpoint or not CCTK_Equals (out_save_parameters, "no")) { vector <vector <vector <region_t> > > grid_superstructure (maps); vector <vector <vector <region_t> > > grid_structure (maps); - vector <vector <vector <CCTK_REAL> > > grid_times (maps); vector <vector <i2vect> > grid_ghosts (maps); vector <vector <i2vect> > grid_buffers (maps); vector <vector <int> > grid_prolongation_orders (maps); for (int m = 0; m < maps; ++ m) { grid_superstructure.at(m) = vhh.at(m)->superregions; grid_structure.at(m) = vhh.at(m)->regions.at(0); - grid_times.at(m).resize(mglevels); grid_ghosts.at(m) = vdd.at(m)->ghost_widths; grid_buffers.at(m) = vdd.at(m)->buffer_widths; grid_prolongation_orders.at(m) = vdd.at(m)->prolongation_orders_space; - for (int ml = 0; ml < mglevels; ++ ml) { - grid_times.at(m).at(ml).resize(vhh.at(m)->reflevels()); - for (int rl = 0; rl < vhh.at(m)->reflevels(); ++ rl) { - grid_times.at(m).at(ml).at(rl) = vtt.at(m)->get_time(rl, ml); + } + vector <vector <vector <CCTK_REAL> > > grid_times (mglevels); + for (int ml = 0; ml < mglevels; ++ ml) { + grid_times.at(ml).resize(vhh.at(0)->reflevels()); + for (int rl = 0; rl < vhh.at(0)->reflevels(); ++ rl) { + grid_times.at(ml).at(rl).resize(tt->timelevels); + for (int tl = 0; tl < tt->timelevels; ++ tl) { + grid_times.at(ml).at(rl).at(tl) = tt->get_time(ml, rl, tl); } } } @@ -1224,7 +1226,7 @@ int WriteMetadata (const cGH * const cctkGH, int const nioprocs, // only into one of the checkpoint files gs_buf << "grid_structure:" << grid_structure << ","; gs_buf << "grid_times:" << grid_times << ","; - gs_buf << "grid_leveltimes:" << leveltimes << ","; + // gs_buf << "grid_leveltimes:" << leveltimes << ","; gs_buf << "grid_ghosts:" << grid_ghosts << ","; gs_buf << "grid_buffers:" << grid_buffers << ","; gs_buf << "grid_prolongation_orders:" << grid_prolongation_orders << "."; diff --git a/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.hh b/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.hh index 81ae3d943..3a30ba48f 100644 --- a/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.hh +++ b/Carpet/CarpetIOHDF5/src/CarpetIOHDF5.hh @@ -17,7 +17,7 @@ // some macros for HDF5 group names #define METADATA_GROUP "Parameters and Global Attributes" #define ALL_PARAMETERS "All Parameters" -#define GRID_STRUCTURE "Grid Structure v3" +#define GRID_STRUCTURE "Grid Structure v4" // atomic HDF5 datatypes for the generic CCTK datatypes // (the one for CCTK_COMPLEX is created at startup as a compound HDF5 datatype) diff --git a/Carpet/CarpetIOHDF5/src/Input.cc b/Carpet/CarpetIOHDF5/src/Input.cc index 340651cb0..336655025 100644 --- a/Carpet/CarpetIOHDF5/src/Input.cc +++ b/Carpet/CarpetIOHDF5/src/Input.cc @@ -59,12 +59,12 @@ typedef struct { int main_loop_index; double global_time; double delta_time; - vector<double> mgleveltimes; // [num_mglevels*num_reflevels] + // vector<double> mgleveltimes; // [num_mglevels*num_reflevels] vector<vector<vector<region_t> > > grid_superstructure; // [map][reflevel][component] vector<vector<vector<region_t> > > grid_structure; // [map][reflevel][component] - vector<vector<vector<CCTK_REAL> > > grid_times; // [map][mglevel][reflevel] - vector<vector<CCTK_REAL> > leveltimes; // [mglevel][reflevel] + vector<vector<vector<CCTK_REAL> > > grid_times; // [mglevel][reflevel][timelevel] + // vector<vector<CCTK_REAL> > leveltimes; // [mglevel][reflevel] vector <vector <i2vect> > grid_ghosts; // [map] vector <vector <i2vect> > grid_buffers; // [map] vector <vector <int> > grid_prolongation_orders; // [map] @@ -193,23 +193,23 @@ void CarpetIOHDF5_RecoverGridStructure (CCTK_ARGUMENTS) } // if regrid_in_level_mode for (int m = 0; m < maps; ++ m) { - // Regrid RegridMap (cctkGH, m, superregsss.at(m), regssss.at(m), false); - - // Set time hierarchy correctly after RegridMap created it - for (int ml = 0; ml < mglevels; ++ ml) { - for (size_t rl = 0; rl < fileset.grid_times.at(m).at(mglevel).size(); ++ rl) { - vtt.at(m)->set_time (rl, ml, fileset.grid_times.at(m).at(ml).at(rl)); + } // for m + + // Set time hierarchy correctly after RegridMap created it + for (int ml = 0; ml < vhh.at(0)->mglevels(); ++ ml) { + for (int rl = 0; rl < vhh.at(0)->reflevels(); ++ rl) { + for (int tl = 0; tl < tt->timelevels; ++ tl) { + tt->set_time (ml, rl, tl, fileset.grid_times.at(ml).at(rl).at(tl)); } } - - } // for m + } PostRegrid (cctkGH); // Set level times correctly after PostRegrid created them - leveltimes = fileset.leveltimes; + // leveltimes = fileset.leveltimes; for (int rl = 0; rl < reflevels; ++ rl) { Recompose (cctkGH, rl, false); @@ -363,13 +363,14 @@ int Recover (cGH* cctkGH, const char *basefilename, int called_from) cctkGH->cctk_iteration = fileset->cctk_iteration; int const idx = mglevel*fileset->num_reflevels + reflevel; - cctkGH->cctk_time = fileset->mgleveltimes.at(idx); + // cctkGH->cctk_time = fileset->mgleveltimes.at(idx); + cctkGH->cctk_time = global_time; if (use_grid_structure_from_checkpoint) { // recover the grid structure only once - static bool is_first = true; - if (is_first) { - is_first = false; + static bool have_grid_structure = false; + if (not have_grid_structure) { + have_grid_structure = true; CarpetIOHDF5_RecoverGridStructure (cctkGH); } } @@ -886,9 +887,9 @@ static void ReadMetadata (fileset_t& fileset, hid_t file) HDF5_ERROR (attr = H5Aopen_name (metadata, "carpet_reflevels")); HDF5_ERROR (H5Aread (attr, H5T_NATIVE_INT, &fileset.num_reflevels)); HDF5_ERROR (H5Aclose (attr)); - HDF5_ERROR (attr = H5Aopen_name (metadata, "numberofmgtimes")); - HDF5_ERROR (H5Aread (attr, H5T_NATIVE_INT, &fileset.num_mglevels)); - HDF5_ERROR (H5Aclose (attr)); +// HDF5_ERROR (attr = H5Aopen_name (metadata, "numberofmgtimes")); +// HDF5_ERROR (H5Aread (attr, H5T_NATIVE_INT, &fileset.num_mglevels)); +// HDF5_ERROR (H5Aclose (attr)); HDF5_ERROR (attr = H5Aopen_name (metadata, "GH$iteration")); HDF5_ERROR (H5Aread (attr, H5T_NATIVE_INT, &fileset.cctk_iteration)); HDF5_ERROR (H5Aclose (attr)); @@ -935,12 +936,12 @@ static void ReadMetadata (fileset_t& fileset, hid_t file) skipws (gs_buf); consume (gs_buf, ","); - skipws (gs_buf); - consume (gs_buf, "grid_leveltimes:"); - skipws (gs_buf); - gs_buf >> fileset.leveltimes; - skipws (gs_buf); - consume (gs_buf, ","); + // skipws (gs_buf); + // consume (gs_buf, "grid_leveltimes:"); + // skipws (gs_buf); + // gs_buf >> fileset.leveltimes; + // skipws (gs_buf); + // consume (gs_buf, ","); skipws (gs_buf); consume (gs_buf, "grid_ghosts:"); @@ -964,19 +965,19 @@ static void ReadMetadata (fileset_t& fileset, hid_t file) consume (gs_buf, "."); } - fileset.mgleveltimes.resize (fileset.num_mglevels * fileset.num_reflevels); - for (int i = 0; i < fileset.num_mglevels; i++) { - char buffer[32]; +// fileset.mgleveltimes.resize (fileset.num_mglevels * fileset.num_reflevels); +// for (int i = 0; i < fileset.num_mglevels; i++) { +// char buffer[32]; - snprintf (buffer, sizeof (buffer), "mgleveltimes %d", i); - HDF5_ERROR (attr = H5Aopen_name (metadata, buffer)); - HDF5_ERROR (dataspace = H5Aget_space (attr)); - assert (H5Sget_simple_extent_npoints (dataspace) == fileset.num_reflevels); - HDF5_ERROR (H5Sclose (dataspace)); - HDF5_ERROR (H5Aread (attr, H5T_NATIVE_DOUBLE, - &fileset.mgleveltimes[i * fileset.num_reflevels])); - HDF5_ERROR (H5Aclose (attr)); - } +// snprintf (buffer, sizeof (buffer), "mgleveltimes %d", i); +// HDF5_ERROR (attr = H5Aopen_name (metadata, buffer)); +// HDF5_ERROR (dataspace = H5Aget_space (attr)); +// assert (H5Sget_simple_extent_npoints (dataspace) == fileset.num_reflevels); +// HDF5_ERROR (H5Sclose (dataspace)); +// HDF5_ERROR (H5Aread (attr, H5T_NATIVE_DOUBLE, +// &fileset.mgleveltimes[i * fileset.num_reflevels])); +// HDF5_ERROR (H5Aclose (attr)); +// } if (is_old_fashioned_file) { HDF5_ERROR (H5Dclose (metadata)); @@ -1245,9 +1246,12 @@ static int ReadVar (const cGH* const cctkGH, } END_LOCAL_COMPONENT_LOOP; if (in_recovery) { - data.tt->set_time (reflevel, mglevel, - ((cctkGH->cctk_time - cctk_initial_time) - / (delta_time * mglevelfact)) ); + if (group.grouptype != CCTK_GF) { + assert (data.tt->timelevels == 1); + data.tt->set_time (mglevel, reflevel, 0, + ((cctkGH->cctk_time - cctk_initial_time) + / (delta_time * mglevelfact)) ); + } } } END_LOCAL_MAP_LOOP; diff --git a/Carpet/CarpetInterp/param.ccl b/Carpet/CarpetInterp/param.ccl index 090ef1bcb..f73f4cc61 100644 --- a/Carpet/CarpetInterp/param.ccl +++ b/Carpet/CarpetInterp/param.ccl @@ -21,7 +21,3 @@ BOOLEAN tree_search "Use a tree search to find the source processor" STEERABLE=a BOOLEAN check_tree_search "Cross-check the result of the tree search" STEERABLE=always { } "no" - -SHARES: Cactus - -USES CCTK_REAL cctk_initial_time diff --git a/Carpet/CarpetInterp/src/interp.cc b/Carpet/CarpetInterp/src/interp.cc index 936b7d5d6..418408110 100644 --- a/Carpet/CarpetInterp/src/interp.cc +++ b/Carpet/CarpetInterp/src/interp.cc @@ -888,9 +888,8 @@ namespace CarpetInterp { assert (partype == PARAMETER_INTEGER); prolongation_order_time = *(CCTK_INT const*) parptr; - current_time = - (cctkGH->cctk_time - cctk_initial_time) / cctkGH->cctk_delta_time; - delta_time = cctkGH->cctk_delta_time; + current_time = cctkGH->cctk_time; + delta_time = cctkGH->cctk_delta_time; iret = Util_TableGetIntArray (param_table_handle, N_output_arrays, &time_deriv_order.front(), "time_deriv_order"); @@ -956,8 +955,10 @@ namespace CarpetInterp { ivect const ipos = ivect(floor((pos - lower) / (delta * rvect(fact)) + rhalf)) * fact; - ivect const & stride = hh->baseextent(ml,rl).stride(); - assert (all (ipos % stride == 0)); + ivect const & istride = hh->baseextent(ml,rl).stride(); + if (hh->refcent == cell_centered) { + assert (all (istride % 2 == 0)); + } gh::cregs const & regs = hh->regions.AT(ml).AT(rl); @@ -1011,14 +1012,16 @@ namespace CarpetInterp { // Try finer levels first for (rl = maxrl-1; rl >= minrl; --rl) { + ivect const & istride = hh->baseextent(ml,rl).stride(); + if (hh->refcent == cell_centered) { + assert (all (istride % 2 == 0)); + } + ivect const fact = maxspacereflevelfact / spacereffacts.AT(rl) * ipow(mgfact, ml); ivect const ipos = ivect(floor((pos - lower) / (delta * rvect(fact)) + rhalf)) * fact; - ivect const & stride = hh->baseextent(ml,rl).stride(); - assert (all (ipos % stride == 0)); - gh::cregs const & regs = hh->superregions.AT(rl); // Search all superregions linearly. Each superregion @@ -1415,7 +1418,7 @@ namespace CarpetInterp { : vector<CCTK_REAL> (num_timelevels_ ) { for (int tl=0; tl<num_timelevels_; ++tl) { - at(tl) = vtt.AT(Carpet::map)->time (tl, reflevel, mglevel); + at(tl) = tt->get_time (mglevel, reflevel, tl); } } diff --git a/Carpet/CarpetLib/src/defs.cc b/Carpet/CarpetLib/src/defs.cc index 3bfdcd585..a2b366ea2 100644 --- a/Carpet/CarpetLib/src/defs.cc +++ b/Carpet/CarpetLib/src/defs.cc @@ -354,6 +354,7 @@ template size_t memoryof (vector<vector<ibbox> > const & v); template size_t memoryof (vector<vector<dh::dboxes> > const & v); template size_t memoryof (vector<vector<dh::fast_dboxes> > const & v); template size_t memoryof (vector<vector<region_t> > const & v); +template size_t memoryof (vector<vector<vector<CCTK_REAL> > > const & v); template size_t memoryof (vector<vector<vector<dh::dboxes> > > const & v); template size_t memoryof (vector<vector<vector<dh::fast_dboxes> > > const & v); template size_t memoryof (vector<vector<vector<region_t> > > const & v); @@ -394,6 +395,7 @@ template ostream& output (ostream& os, const vector<CCTK_REAL>& v); template ostream& output (ostream& os, const vector<ibbox>& v); template ostream& output (ostream& os, const vector<rbbox>& v); template ostream& output (ostream& os, const vector<ivect>& v); +template ostream& output (ostream& os, const vector<rvect>& v); template ostream& output (ostream& os, const vector<bbvect>& v); template ostream& output (ostream& os, const vector<i2vect>& v); template ostream& output (ostream& os, const vector<dh::dboxes> & v); diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc index 7d947fae5..8a7c643da 100644 --- a/Carpet/CarpetLib/src/ggf.cc +++ b/Carpet/CarpetLib/src/ggf.cc @@ -38,7 +38,7 @@ ggf::ggf (const int varindex_, const operator_type transport_operator_, vectorlength(vectorlength_), vectorindex(vectorindex_), vectorleader(vectorleader_) { - assert (&t.h == &d.h); + // assert (&t.h == &d.h); assert (vectorlength >= 1); assert (vectorindex >= 0 and vectorindex < vectorlength); @@ -176,18 +176,6 @@ void ggf::recompose_fill (comm_state & state, int const rl, assert (d.fast_boxes.AT(ml).AT(rl).do_init); - vector <int> tls; - if (do_prolongate and rl > 0 and - transport_operator != op_none and transport_operator != op_sync and - transport_operator != op_restrict) - { - int const numtl = timelevels (ml, rl); - tls.resize (numtl); - for (int tl = 0; tl < numtl; ++ tl) { - tls.AT(tl) = tl; - } - } - // Initialise from the same level of the old hierarchy, where // possible if (rl < (int)oldstorage.AT(ml).size()) { @@ -207,12 +195,18 @@ void ggf::recompose_fill (comm_state & state, int const rl, if (transport_operator != op_none and transport_operator != op_sync and transport_operator != op_restrict) { + int const numtl = timelevels (ml, rl); + vector <int> tls (numtl); + for (int tl = 0; tl < numtl; ++ tl) { + tls.AT(tl) = tl; + } + for (int tl = 0; tl < timelevels (ml, rl); ++tl) { transfer_from_all (state, tl, rl, ml, & dh::fast_dboxes::fast_old2new_ref_prol_sendrecv, tls, rl - 1, ml, - t.time (tl, rl, ml)); + t.get_time (ml, rl, tl)); } // for tl } // if transport_operator } // if rl @@ -277,6 +271,22 @@ void ggf::cycle_all (int const rl, int const ml) { } } +// Uncycle the time levels by rotating the data sets +void ggf::uncycle_all (int const rl, int const ml) { + assert (rl>=0 and rl<h.reflevels()); + assert (ml>=0 and ml<h.mglevels()); + int const ntl = timelevels(ml,rl); + assert (ntl > 0); + for (int lc=0; lc<(int)storage.AT(ml).AT(rl).size(); ++lc) { + fdata & fdatas = storage.AT(ml).AT(rl).AT(lc); + gdata * const tmpdata = fdatas.AT(0); + for (int tl=0; tl<ntl-1; ++tl) { + fdatas.AT(tl) = fdatas.AT(tl+1); + } + fdatas.AT(ntl-1) = tmpdata; + } +} + // Flip the time levels by exchanging the data sets void ggf::flip_all (int const rl, int const ml) { assert (rl>=0 and rl<h.reflevels()); @@ -285,9 +295,9 @@ void ggf::flip_all (int const rl, int const ml) { assert (ntl > 0); for (int lc=0; lc<(int)storage.AT(ml).AT(rl).size(); ++lc) { fdata & fdatas = storage.AT(ml).AT(rl).AT(lc); - for (int tl=0; tl<ntl/2; ++tl) { - const int tl1 = tl; - const int tl2 = ntl-1 - tl; + for (int tl=1; tl<(ntl+1)/2; ++tl) { + const int tl1 = tl; + const int tl2 = ntl - tl; assert (tl1 < tl2); gdata * const tmpdata = fdatas.AT(tl1); fdatas.AT(tl1) = fdatas.AT(tl2); @@ -383,15 +393,15 @@ ref_bnd_prolongate_all (comm_state & state, void ggf:: mg_restrict_all (comm_state & state, - int const tl, int const rl, int const ml, + int const tl, int const rl, int const ml, CCTK_REAL const time) { static Timer timer ("mg_restrict_all"); timer.start (); // Require same times static_assert (abs(0.1) > 0, "Function CarpetLib::abs has wrong signature"); - assert (abs(t.get_time(rl,ml) - t.get_time(rl,ml-1)) - <= 1.0e-8 * (1.0 + abs(t.get_time(rl,ml)))); + assert (abs(t.get_time(ml,rl,0) - t.get_time(ml-1,rl,0)) + <= 1.0e-8 * (1.0 + abs(t.get_time(ml,rl,0)))); vector<int> const tl2s(1,tl); transfer_from_all (state, tl ,rl,ml, @@ -414,8 +424,8 @@ mg_prolongate_all (comm_state & state, timer.start (); // Require same times static_assert (abs(0.1) > 0, "Function CarpetLib::abs has wrong signature"); - assert (abs(t.get_time(rl,ml) - t.get_time(rl,ml+1)) - <= 1.0e-8 * (1.0 + abs(t.get_time(rl,ml)))); + assert (abs(t.get_time(ml,rl,0) - t.get_time(ml+1,rl,0)) + <= 1.0e-8 * (1.0 + abs(t.get_time(ml,rl,0)))); vector<int> const tl2s(1,tl); transfer_from_all (state, tl ,rl,ml, @@ -431,22 +441,19 @@ mg_prolongate_all (comm_state & state, void ggf:: ref_restrict_all (comm_state & state, - int const tl, int const rl, int const ml, - CCTK_REAL const time) + int const tl, int const rl, int const ml) { // Require same times static_assert (abs(0.1) > 0, "Function CarpetLib::abs has wrong signature"); - assert (abs(t.get_time(rl,ml) - t.get_time(rl+1,ml)) - <= 1.0e-8 * (1.0 + abs(t.get_time(rl,ml)))); + assert (abs(t.get_time(ml,rl,tl) - t.get_time(ml,rl+1,tl)) <= + 1.0e-8 * (1.0 + abs(t.get_time(ml,rl,tl)))); if (transport_operator == op_none or transport_operator == op_sync) return; static Timer timer ("ref_restrict_all"); timer.start (); - vector<int> const tl2s(1,tl); transfer_from_all (state, - tl ,rl ,ml, + tl,rl ,ml, & dh::fast_dboxes::fast_ref_rest_sendrecv, - tl2s,rl+1,ml, - time); + tl,rl+1,ml); timer.stop (0); } @@ -522,7 +529,7 @@ transfer_from_all (comm_state & state, for (size_t j=0; j<i; ++j) { assert (tl2s.AT(i) != tl2s.AT(j)); } - times.AT(i) = t.time(tl2s.AT(i),rl2,ml2); + times.AT(i) = t.get_time(ml2,rl2,tl2s.AT(i)); } // Interpolation orders diff --git a/Carpet/CarpetLib/src/ggf.hh b/Carpet/CarpetLib/src/ggf.hh index 43c5c2e93..5538f3941 100644 --- a/Carpet/CarpetLib/src/ggf.hh +++ b/Carpet/CarpetLib/src/ggf.hh @@ -105,6 +105,9 @@ public: // Cycle the time levels by rotating the data sets void cycle_all (int rl, int ml); + + // Uncycle the time levels by rotating the data sets + void uncycle_all (int rl, int ml); // Flip the time levels by exchanging the data sets void flip_all (int rl, int ml); @@ -135,7 +138,7 @@ public: // Restrict a refinement level void ref_restrict_all (comm_state& state, - int tl, int rl, int ml, CCTK_REAL time); + int tl, int rl, int ml); // Prolongate a refinement level void ref_prolongate_all (comm_state& state, @@ -171,7 +174,7 @@ protected: { vector <int> tl2s(1); tl2s.AT(0) = tl2; - CCTK_REAL const time = t.time (tl2,rl2,ml2); + CCTK_REAL const time = t.get_time (ml2,rl2,tl2); transfer_from_all (state, tl1, rl1, ml1, sendrecvs, diff --git a/Carpet/CarpetLib/src/gh.cc b/Carpet/CarpetLib/src/gh.cc index 21cb399a6..05c1b55ea 100644 --- a/Carpet/CarpetLib/src/gh.cc +++ b/Carpet/CarpetLib/src/gh.cc @@ -549,7 +549,7 @@ gh:: output (ostream & os) const { - os << "gh:" + os << "gh:{" << "reffacts=" << reffacts << ",refcentering=" << refcent << "," << "mgfactor=" << mgfact << ",mgcentering=" << mgcent << "," << "baseextents=" << baseextents << "," @@ -565,6 +565,6 @@ output (ostream & os) os << *d; } } - os << "}"; + os << "}}"; return os; } diff --git a/Carpet/CarpetLib/src/th.cc b/Carpet/CarpetLib/src/th.cc index dd441e455..bed473245 100644 --- a/Carpet/CarpetLib/src/th.cc +++ b/Carpet/CarpetLib/src/th.cc @@ -19,8 +19,10 @@ list<th*> th::allth; // Constructors -th::th (gh& h_, const vector<int> & reffacts_, const CCTK_REAL basedelta) - : h(h_), reffacts(reffacts_), delta(basedelta) +th::th (gh& h_, int const timelevels_, vector<int> const& reffacts_, + bool const time_interpolation_during_regridding_) + : h(h_), timelevels(timelevels_), reffacts(reffacts_), + time_interpolation_during_regridding (time_interpolation_during_regridding_) { assert (reffacts.size() >= 1); assert (reffacts.front() == 1); @@ -42,6 +44,9 @@ th::~th () // Modifiers void th::regrid () { + CCTK_REAL const basetime = 0.0; + CCTK_REAL const basedelta = 1.0; + const int old_mglevels = times.size(); times.resize(h.mglevels()); deltas.resize(h.mglevels()); @@ -50,17 +55,64 @@ void th::regrid () times.AT(ml).resize(h.reflevels()); deltas.AT(ml).resize(h.reflevels()); for (int rl=0; rl<h.reflevels(); ++rl) { + if (ml==0) { + deltas.AT(ml).AT(rl) = basedelta / reffacts.AT(rl); + } else { + deltas.AT(ml).AT(rl) = deltas.AT(ml-1).AT(rl) * h.mgfact; + } if (old_mglevels==0) { - times.AT(ml).AT(rl) = 0; + times.AT(ml).AT(rl).resize(timelevels); + for (int tl=0; tl<timelevels; ++tl) { + times.AT(ml).AT(rl).AT(tl) = basetime - tl * deltas.AT(ml).AT(rl); + } } else if (rl < old_reflevels) { // do nothing + } else if (rl == 0) { + assert (0); + times.AT(ml).AT(rl).resize(timelevels); + for (int tl=0; tl<timelevels; ++tl) { + times.AT(ml).AT(rl).AT(tl) = basetime - tl * deltas.AT(ml).AT(rl); + } } else { - times.AT(ml).AT(rl) = times.AT(ml).AT(rl-1); + if (time_interpolation_during_regridding) { +#warning "We probably don't want to do this, but it is nice for compatibility" + times.AT(ml).AT(rl).resize(timelevels); + for (int tl=0; tl<timelevels; ++tl) { + // linear interpolation between the two surrounding coarse + // grid times + assert (reffacts.AT(rl) % reffacts.AT(rl-1) == 0); + int const rf = reffacts.AT(rl) / reffacts.AT(rl-1); + int const ctl = tl / rf; + int const mtl = tl % rf; + if (mtl == 0) { + times.AT(ml).AT(rl).AT(tl) = times.AT(ml).AT(rl-1).AT(ctl); + } else { + assert (ctl+1<timelevels); + CCTK_REAL const alpha = (CCTK_REAL)mtl / rf; + assert (alpha>0 and alpha<1); + times.AT(ml).AT(rl).AT(tl) = + (1-alpha) * times.AT(ml).AT(rl-1).AT(ctl ) + + ( alpha) * times.AT(ml).AT(rl-1).AT(ctl+1); + } + } + } else { + times.AT(ml).AT(rl) = times.AT(ml).AT(rl-1); + } } - if (ml==0) { - deltas.AT(ml).AT(rl) = delta / reffacts.AT(rl); - } else { - deltas.AT(ml).AT(rl) = deltas.AT(ml-1).AT(rl) * h.mgfact; + } + } + for (int ml=0; ml<h.mglevels(); ++ml) { + for (int rl=0; rl<h.reflevels(); ++rl) { + for (int tl=1; tl<timelevels; ++tl) { + assert (times.AT(ml).AT(rl).AT(tl) < times.AT(ml).AT(rl).AT(tl-1)); + } + } + } + for (int ml=0; ml<h.mglevels(); ++ml) { + for (int rl=0; rl<h.reflevels(); ++rl) { + // assert (isfinite(deltas.AT(ml).AT(rl))); + for (int tl=0; tl<timelevels; ++tl) { + // assert (isfinite(times.AT(ml).AT(rl).AT(tl))); } } } @@ -72,6 +124,37 @@ void th::regrid_free () +void th::advance_time (int const ml, int const rl) +{ + for (int tl=timelevels-1; tl>0; --tl) { + set_time (ml,rl,tl, get_time(ml,rl,tl-1)); + } + set_time (ml,rl,0, get_time(ml,rl,0) + get_delta(ml,rl)); +} + +void th::retreat_time (int const ml, int const rl) +{ + CCTK_REAL const t = get_time(ml,rl,0); + for (int tl=0; tl<timelevels-1; ++tl) { + set_time (ml,rl,tl, get_time(ml,rl,tl+1)); + } + set_time (ml,rl,timelevels-1, t); +} + +void th::flip_timelevels (int const ml, int const rl) +{ + for (int tl=1; tl<(timelevels+1)/2; ++tl) { + int const tl2 = timelevels - tl; + CCTK_REAL const t = get_time(ml,rl,tl ); + CCTK_REAL const t2 = get_time(ml,rl,tl2); + set_time (ml,rl,tl ,t2); + set_time (ml,rl,tl2,t ); + } + set_delta (ml,rl, -get_delta(ml,rl)); +} + + + // Memory usage size_t th:: @@ -80,7 +163,6 @@ memory () { return memoryof (reffacts) + - memoryof (delta) + memoryof (times) + memoryof (deltas); } @@ -100,19 +182,33 @@ allmemory () +// Input +istream& th::input (istream& is) +{ + skipws (is); + consume (is, "th:{"); + consume (is, "timelevels="); + is >> timelevels; + consume (is, ","); + consume (is, "reffacts="); + is >> reffacts; + consume (is, ","); + consume (is, "times="); + is >> times; + consume (is, ","); + consume (is, "deltas="); + is >> deltas; + consume (is, "}"); + return is; +} + // Output -void th::output (ostream& os) const +ostream& th::output (ostream& os) const { - os << "th:" - << "times={"; - const char * sep = ""; - for (int ml=0; ml<h.mglevels(); ++ml) { - for (int rl=0; rl<h.reflevels(); ++rl) { - os << sep - << ml << ":" << rl << ":" - << times.AT(ml).AT(rl) << "(" << deltas.AT(ml).AT(rl) << ")"; - sep = ","; - } - } - os << "}"; + os << "th:{" + << "timelevels=" << timelevels << "," + << "reffacts=" << reffacts << "," + << "times=" << times << "," + << "deltas=" << deltas << "}"; + return os; } diff --git a/Carpet/CarpetLib/src/th.hh b/Carpet/CarpetLib/src/th.hh index f1a7b72b7..c1453cf4d 100644 --- a/Carpet/CarpetLib/src/th.hh +++ b/Carpet/CarpetLib/src/th.hh @@ -2,10 +2,11 @@ #define TH_HH #include <cassert> +#include <cmath> #include <iostream> #include <vector> -#include "cctk.h" +#include <cctk.h> #include "defs.hh" #include "gh.hh" @@ -17,6 +18,8 @@ using namespace std; // Forward declaration class th; +// Input +istream& operator>> (istream& is, th& t); // Output ostream& operator<< (ostream& os, const th& t); @@ -34,18 +37,22 @@ public: // should be readonly gh& h; // hierarchy gh::th_handle gh_handle; + int timelevels; // const + private: - const vector<int> reffacts; - - CCTK_REAL delta; // time step - vector<vector<CCTK_REAL> > times; // current times [ml][rl] - vector<vector<CCTK_REAL> > deltas; // time steps [ml][rl] + vector<int> reffacts; // const + + bool const time_interpolation_during_regridding; + + vector<vector<vector<CCTK_REAL> > > times; // current times [ml][rl][tl] + vector<vector<CCTK_REAL> > deltas; // time steps [ml][rl] public: // Constructors - th (gh& h, const vector<int> & reffacts, const CCTK_REAL basedelta); + th (gh& h, int timelevels, vector<int> const& reffacts, + bool time_interpolation_during_regridding); // Destructors ~th (); @@ -55,50 +62,52 @@ public: void regrid_free (); // Time management - CCTK_REAL get_time (const int rl, const int ml) const CCTK_ATTRIBUTE_PURE + void set_time (int const ml, int const rl, int const tl, CCTK_REAL const& t) { - assert (rl>=0 and rl<h.reflevels()); assert (ml>=0 and ml<h.mglevels()); - return times.AT(ml).AT(rl); - } - - void set_time (const int rl, const int ml, const CCTK_REAL t) - { assert (rl>=0 and rl<h.reflevels()); - assert (ml>=0 and ml<h.mglevels()); - times.AT(ml).AT(rl) = t; - } - - void advance_time (const int rl, const int ml) - { - set_time(rl,ml, get_time(rl,ml) + get_delta(rl,ml)); + assert (tl>=0 and tl<timelevels); + // assert (isfinite(t)); + times.AT(ml).AT(rl).AT(tl) = t; } - CCTK_REAL get_delta (const int rl, const int ml) const CCTK_ATTRIBUTE_PURE + CCTK_REAL get_time (int const ml, int const rl, int const tl) + const CCTK_ATTRIBUTE_PURE { - assert (rl>=0 and rl<h.reflevels()); assert (ml>=0 and ml<h.mglevels()); - return deltas.AT(ml).AT(rl); + assert (rl>=0 and rl<h.reflevels()); + assert (tl>=0 and tl<timelevels); + CCTK_REAL const t = times.AT(ml).AT(rl).AT(tl); + // assert (isfinite(t)); + return t; } - void set_delta (const int rl, const int ml, const CCTK_REAL dt) + void set_delta (int const ml, int const rl, CCTK_REAL const& dt) { - assert (rl>=0 and rl<h.reflevels()); assert (ml>=0 and ml<h.mglevels()); + assert (rl>=0 and rl<h.reflevels()); + // assert (isfinite(dt)); deltas.AT(ml).AT(rl) = dt; } - CCTK_REAL time (const int tl, const int rl, const int ml) const CCTK_ATTRIBUTE_PURE + CCTK_REAL get_delta (int const ml, int const rl) const CCTK_ATTRIBUTE_PURE { - assert (rl>=0 and rl<h.reflevels()); assert (ml>=0 and ml<h.mglevels()); - return get_time(rl, ml) - tl * get_delta(rl, ml); + assert (rl>=0 and rl<h.reflevels()); + CCTK_REAL const dt = deltas.AT(ml).AT(rl); + // assert (isfinite(dt)); + return dt; } + void advance_time (int const ml, int const rl); + void retreat_time (int const ml, int const rl); + void flip_timelevels (int const ml, int const rl); + // Output size_t memory () const CCTK_ATTRIBUTE_PURE; static size_t allmemory () CCTK_ATTRIBUTE_PURE; - void output (ostream& os) const; + istream& input (istream& is); + ostream& output (ostream& os) const; }; @@ -108,9 +117,11 @@ inline size_t memoryof (th const & t) { return t.memory (); } +inline istream& operator>> (istream& is, th& t) { + return t.input(is); +} inline ostream& operator<< (ostream& os, const th& t) { - t.output(os); - return os; + return t.output(os); } diff --git a/Carpet/CarpetReduce/src/reduce.cc b/Carpet/CarpetReduce/src/reduce.cc index cebda71db..577fe2478 100644 --- a/Carpet/CarpetReduce/src/reduce.cc +++ b/Carpet/CarpetReduce/src/reduce.cc @@ -1275,7 +1275,7 @@ namespace CarpetReduce { CCTK_REAL const time = current_time; vector<CCTK_REAL> times(num_tl); for (int tl=0; tl<num_tl; ++tl) { - times.at(tl) = vtt.at(0)->time (tl, reflevel, mglevel); + times.at(tl) = tt->get_time (mglevel, reflevel, tl); } // Calculate interpolation weights diff --git a/Carpet/CarpetRegrid2/src/regrid.cc b/Carpet/CarpetRegrid2/src/regrid.cc index 360b6e0c4..c5ea20a80 100644 --- a/Carpet/CarpetRegrid2/src/regrid.cc +++ b/Carpet/CarpetRegrid2/src/regrid.cc @@ -344,10 +344,8 @@ namespace CarpetRegrid2 { if (freeze_unaligned_parent_levels) { // Assume that non-existing levels are always aligned if (min_rl < reflevels) { - CCTK_REAL const mytime = - vtt.at(Carpet::map)->time (0, min_rl, mglevel); - CCTK_REAL const parenttime = - vtt.at(Carpet::map)->time (0, min_rl - 1, mglevel); + CCTK_REAL const mytime = tt->get_time (mglevel, min_rl, 0); + CCTK_REAL const parenttime = tt->get_time (mglevel, min_rl - 1, 0); CCTK_REAL const eps = 1.0e-12; in_sync = abs (mytime - parenttime) <= eps * abs (delta_time); diff --git a/CarpetDev/CarpetJacobi/src/Jacobi.cc b/CarpetDev/CarpetJacobi/src/Jacobi.cc index 0f1149918..b6969b6a5 100644 --- a/CarpetDev/CarpetJacobi/src/Jacobi.cc +++ b/CarpetDev/CarpetJacobi/src/Jacobi.cc @@ -23,7 +23,7 @@ namespace Carpet { // TODO: fix this - void CycleTimeLevels (const cGH* cctkGH); + void CycleTimeLevels (cGH* cctkGH); void Restrict (const cGH* cctkGH); }; @@ -202,12 +202,12 @@ namespace CarpetJacobi { * const_cast<CCTK_REAL *> (& cctkGH->cctk_delta_time) = delta_time; BEGIN_REFLEVEL_LOOP(cctkGH) { * const_cast<CCTK_REAL *> (& cctkGH->cctk_time) = global_time; - for (int m=0; m<maps; ++m) { - vtt[m]->set_time (reflevel, mglevel, 0); - vtt[m]->set_delta - (reflevel, mglevel, - 1.0 / ipow (maxval (spacereflevelfact), reflevelpower)); + for (int tl=0; tl<timelevels; ++tl) { + tt->set_time (mglevel, reflevel, tl, global_time - tl * delta_time); } + tt->set_delta + (mglevel, reflevel, + 1.0 / ipow (maxval (spacereflevelfact), reflevelpower)); } END_REFLEVEL_LOOP; @@ -252,13 +252,11 @@ namespace CarpetJacobi { if (reflevel <= solve_level && (iter-istep) % do_every == 0) { // Advance time - for (int m=0; m<maps; ++m) { - vtt[m]->advance_time (reflevel, mglevel); - } + tt->advance_time (mglevel, reflevel); * const_cast<CCTK_REAL *> (& cctkGH->cctk_time) = (1.0 * (iter - istep + do_every) / ipow (maxval (maxspacereflevelfact), reflevelpower)); - CycleTimeLevels (cctkGH); + CycleTimeLevels (const_cast<cGH*>(cctkGH)); if (DEBUG) cout << "CJ residual iter " << iter << " reflevel " << reflevel << " time " << cctkGH->cctk_time << flush << endl; // Advance time levels @@ -569,10 +567,10 @@ namespace CarpetJacobi { * const_cast<CCTK_REAL *> (& cctkGH->cctk_delta_time) = delta_time; BEGIN_REFLEVEL_LOOP(cctkGH) { * const_cast<CCTK_REAL *> (& cctkGH->cctk_time) = global_time; - for (int m=0; m<maps; ++m) { - vtt[m]->set_time (reflevel, mglevel, 0); - vtt[m]->set_delta (reflevel, mglevel, 1.0 / timereflevelfact); + for (int tl=0; tl<timelevels; ++tl) { + tt->set_time (mglevel, reflevel, tl, global_time - tl * delta_time); } + tt->set_delta (mglevel, reflevel, 1.0 / timereflevelfact); } END_REFLEVEL_LOOP; } END_GLOBAL_MODE; diff --git a/CarpetDev/CarpetMG/src/mg.cc b/CarpetDev/CarpetMG/src/mg.cc index d76a10212..07e889e02 100644 --- a/CarpetDev/CarpetMG/src/mg.cc +++ b/CarpetDev/CarpetMG/src/mg.cc @@ -1015,7 +1015,6 @@ namespace CarpetMG { int const tl = 0; for (comm_state state; !state.done(); state.step()) { for (int m=0; m<maps; ++m) { - const CCTK_REAL time = vtt.at(m)->time (tl, reflevel, mglevel); for (int n=0; n<var.size(); ++n) { int const vi = var.at(n); assert (vi >= 0); @@ -1024,7 +1023,7 @@ namespace CarpetMG { int const v0 = CCTK_FirstVarIndexI (gi); assert (v0 >= 0); arrdata.at(gi).at(m).data.at(vi-v0)->ref_restrict_all - (state, tl, reflevel, mglevel, time); + (state, tl, reflevel, mglevel); } // for n } // for m } // for state @@ -1040,9 +1039,9 @@ namespace CarpetMG { { assert (reflevel > 0); int const tl = 0; + CCTK_REAL const time = tt->get_time (mglevel, reflevel, tl); for (comm_state state; !state.done(); state.step()) { for (int m=0; m<maps; ++m) { - const CCTK_REAL time = vtt.at(m)->time (tl, reflevel, mglevel); for (int n=0; n<var.size(); ++n) { int const vi = var.at(n); assert (vi >= 0); |