diff options
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); |