diff options
53 files changed, 2244 insertions, 663 deletions
diff --git a/Carpet/Carpet/param.ccl b/Carpet/Carpet/param.ccl index 6742612e2..4c226ed5b 100644 --- a/Carpet/Carpet/param.ccl +++ b/Carpet/Carpet/param.ccl @@ -1,5 +1,5 @@ # Parameter definitions for thorn Carpet -# $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/param.ccl,v 1.33 2003/08/15 09:34:36 schnetter Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/param.ccl,v 1.34 2003/11/05 16:18:37 schnetter Exp $ shares: Cactus @@ -194,6 +194,7 @@ KEYWORD processor_topology "How to determine the processor topology" { "manual" :: "Specified by processor_topology_*" "along-z" :: "Split the region along the z direction only" + "along-dir" :: "Split the region along one direction only" "automatic" :: "Choose the topology automatically" } "automatic" @@ -212,6 +213,11 @@ CCTK_INT processor_topology_3d_z "Number of processors in z-direction" 1:* :: "must be positive" } 1 +CCTK_INT split_direction "Direction in which the domain should be split" +{ + 0:* :: "0 for x, 1 for y, 2 for z, etc." +} 2 + STRING grid_structure_filename "File name to output grid structure to (empty = no output)" diff --git a/Carpet/Carpet/schedule.ccl b/Carpet/Carpet/schedule.ccl index 228dced2b..33aab5abc 100644 --- a/Carpet/Carpet/schedule.ccl +++ b/Carpet/Carpet/schedule.ccl @@ -1,5 +1,5 @@ # Schedule definitions for thorn Carpet -# $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/schedule.ccl,v 1.5 2003/07/20 21:03:43 schnetter Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/schedule.ccl,v 1.6 2003/11/05 16:18:37 schnetter Exp $ schedule CarpetStartup at STARTUP as Driver_Startup { @@ -8,6 +8,6 @@ schedule CarpetStartup at STARTUP as Driver_Startup -schedule group PostRestrict +schedule group POSTRESTRICT { } "Apply boundary conditions after restricting" diff --git a/Carpet/Carpet/src/Comm.cc b/Carpet/Carpet/src/Comm.cc index 99f52b5f9..e110bcbc1 100644 --- a/Carpet/Carpet/src/Comm.cc +++ b/Carpet/Carpet/src/Comm.cc @@ -10,7 +10,7 @@ #include "carpet.hh" extern "C" { - static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Comm.cc,v 1.22 2003/08/10 21:59:51 schnetter Exp $"; + static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Comm.cc,v 1.23 2003/11/05 16:18:37 schnetter Exp $"; CCTK_FILEVERSION(Carpet_Carpet_Comm_cc); } @@ -65,37 +65,49 @@ namespace Carpet { assert (num_tl>0); const int tl = 0; - assert (group<(int)arrdata.size()); - for (int var=0; var<(int)arrdata[group].data.size(); ++var) { + // Prolongate the boundaries + if (reflevel>0) { if (grouptype == CCTK_GF) { if (do_prolongate) { if (arrdata[group].do_transfer) { - if (reflevel>0) { - // use the current time here (which may be modified by the - // user) - const CCTK_REAL time = (cgh->cctk_time - cctk_initial_time) / delta_time; + for (comm_state<dim> state; !state.done(); state.step()) { + assert (group<(int)arrdata.size()); + for (int var=0; var<(int)arrdata[group].data.size(); ++var) { + // use the current time here (which may be modified by + // the user) + const CCTK_REAL time = (cgh->cctk_time - cctk_initial_time) / delta_time; #if 0 - const CCTK_REAL time1 = tt->time (tl, reflevel, mglevel); - const CCTK_REAL time2 = (cgh->cctk_time - cctk_initial_time) / delta_time; - assert (fabs((time1 - time2) / (fabs(time1) + fabs(time2) + fabs(cgh->cctk_delta_time))) < 1e-12); + const CCTK_REAL time1 = tt->time (tl, reflevel, mglevel); + const CCTK_REAL time2 = (cgh->cctk_time - cctk_initial_time) / delta_time; + assert (fabs((time1 - time2) / (fabs(time1) + fabs(time2) + fabs(cgh->cctk_delta_time))) < 1e-12); #endif - for (int c=0; c<arrdata[group].hh->components(reflevel); ++c) { - arrdata[group].data[var]->ref_bnd_prolongate - (tl, reflevel, c, mglevel, time); - } - } + for (int c=0; c<arrdata[group].hh->components(reflevel); ++c) { + arrdata[group].data[var]->ref_bnd_prolongate + (state, tl, reflevel, c, mglevel, time); + } + } // for var + } // for state } else { Checkpoint ("%*s(no prolongating for group %s)", 2*reflevel, "", groupname); - } + } // if ! do_transfer } // if do_prolongate - for (int c=0; c<arrdata[group].hh->components(reflevel); ++c) { - arrdata[group].data[var]->sync (tl, reflevel, c, mglevel); - } - } else { // grouptype != CCTK_GF - arrdata[group].data[var]->sync (0, 0, 0, 0); - } - } + } // if grouptype == CCTK_GF + } // if reflevel>0 + + // Sync + for (comm_state<dim> state; !state.done(); state.step()) { + assert (group<(int)arrdata.size()); + for (int var=0; var<(int)arrdata[group].data.size(); ++var) { + if (grouptype == CCTK_GF) { + for (int c=0; c<arrdata[group].hh->components(reflevel); ++c) { + arrdata[group].data[var]->sync (state, tl, reflevel, c, mglevel); + } + } else { // grouptype != CCTK_GF + arrdata[group].data[var]->sync (state, 0, 0, 0, 0); + } // grouptype != CCTK_GF + } // for var + } // for state return 0; } diff --git a/Carpet/Carpet/src/Cycle.cc b/Carpet/Carpet/src/Cycle.cc index 170fe6dc9..932aab288 100644 --- a/Carpet/Carpet/src/Cycle.cc +++ b/Carpet/Carpet/src/Cycle.cc @@ -9,7 +9,7 @@ #include "carpet.hh" extern "C" { - static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Cycle.cc,v 1.14 2003/06/18 18:28:07 schnetter Exp $"; + static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Cycle.cc,v 1.15 2003/11/05 16:18:37 schnetter Exp $"; CCTK_FILEVERSION(Carpet_Carpet_Cycle_cc); } @@ -63,7 +63,7 @@ namespace Carpet { assert (group<(int)arrdata.size()); assert (var<(int)arrdata[group].data.size()); - for (int rl=0; rl<hh->reflevels(); ++rl) { + for (int rl=0; rl<arrdata[group].hh->reflevels(); ++rl) { for (int c=0; c<arrdata[group].hh->components(rl); ++c) { arrdata[group].data[var]->flip (rl, c, mglevel); } diff --git a/Carpet/Carpet/src/Evolve.cc b/Carpet/Carpet/src/Evolve.cc index 672d0b160..91308d2b2 100644 --- a/Carpet/Carpet/src/Evolve.cc +++ b/Carpet/Carpet/src/Evolve.cc @@ -5,10 +5,23 @@ #include "cctk_Parameters.h" #include "cctk_Termination.h" -#ifdef HAVE_TIME_GETTIMEOFDAY -# ifdef HAVE_SYS_TIME_H +// IRIX wants this before <time.h> +#if HAVE_SYS_TYPES_H +# include <sys/types.h> +#endif + +#if TIME_WITH_SYS_TIME +# include <sys/time.h> +# include <time.h> +#else +# if HAVE_SYS_TIME_H # include <sys/time.h> +# elif HAVE_TIME_H +# include <time.h> # endif +#endif + +#if HAVE_UNISTD_H # include <unistd.h> #endif @@ -18,7 +31,7 @@ #include "carpet.hh" extern "C" { - static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Evolve.cc,v 1.28 2003/08/10 21:59:51 schnetter Exp $"; + static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Evolve.cc,v 1.29 2003/11/05 16:18:37 schnetter Exp $"; CCTK_FILEVERSION(Carpet_Carpet_Evolve_cc); } @@ -134,7 +147,7 @@ namespace Carpet { if ((cgh->cctk_iteration-1) % do_every == 0) { BEGIN_MGLEVEL_LOOP(cgh) { - const int do_every = mglevelfact * maxreflevelfact/reflevelfact; + const int do_every = mglevelfact * (maxreflevelfact/reflevelfact); if ((cgh->cctk_iteration-1) % do_every == 0) { do_global_mode @@ -175,7 +188,7 @@ namespace Carpet { // Regrid Waypoint ("%*sRegrid", 2*reflevel, ""); - Regrid (cgh); + Regrid (cgh, reflevel+1, true); } } END_REFLEVEL_LOOP; @@ -187,7 +200,7 @@ namespace Carpet { if (cgh->cctk_iteration % do_every == 0) { BEGIN_MGLEVEL_LOOP(cgh) { - const int do_every = mglevelfact * maxreflevelfact/reflevelfact; + const int do_every = mglevelfact * (maxreflevelfact/reflevelfact); if (cgh->cctk_iteration % do_every == 0) { do_global_mode @@ -201,8 +214,8 @@ namespace Carpet { do_global_mode ? " (global time)" : ""); Restrict (cgh); - Waypoint ("%*sScheduling PostRestrict", 2*reflevel, ""); - CCTK_ScheduleTraverse ("PostRestrict", cgh, CallFunction); + Waypoint ("%*sScheduling POSTRESTRICT", 2*reflevel, ""); + CCTK_ScheduleTraverse ("CCTK_POSTRESTRICT", cgh, CallFunction); // Checking CalculateChecksums (cgh, currenttime); @@ -220,7 +233,7 @@ namespace Carpet { if (cgh->cctk_iteration % do_every == 0) { BEGIN_MGLEVEL_LOOP(cgh) { - const int do_every = mglevelfact * maxreflevelfact/reflevelfact; + const int do_every = mglevelfact * (maxreflevelfact/reflevelfact); if (cgh->cctk_iteration % do_every == 0) { do_global_mode diff --git a/Carpet/Carpet/src/Initialise.cc b/Carpet/Carpet/src/Initialise.cc index e0468bb28..231bbf5d1 100644 --- a/Carpet/Carpet/src/Initialise.cc +++ b/Carpet/Carpet/src/Initialise.cc @@ -12,7 +12,7 @@ #include "carpet.hh" extern "C" { - static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Initialise.cc,v 1.33 2003/09/02 13:11:16 tradke Exp $"; + static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Initialise.cc,v 1.34 2003/11/05 16:18:37 schnetter Exp $"; CCTK_FILEVERSION(Carpet_Carpet_Initialise_cc); } @@ -47,12 +47,12 @@ namespace Carpet { CCTKi_InitGHExtensions (cgh); // Register coordinates - Waypoint ("CCTK_WRAGH"); + Waypoint ("Scheduling CCTK_WRAGH"); CCTK_ScheduleTraverse ("CCTK_WRAGH", cgh, CallFunction); // Check parameters Waypoint ("Current time is %g", cgh->cctk_time); - Waypoint ("PARAMCHECK"); + Waypoint ("Scheduling PARAMCHECK"); CCTK_ScheduleTraverse ("CCTK_PARAMCHECK", cgh, CallFunction); CCTKi_FinaliseParamWarn(); @@ -153,7 +153,7 @@ namespace Carpet { // Regrid Waypoint ("%*sRegrid", 2*reflevel, ""); - Regrid (cgh, reflevel); + Regrid (cgh, reflevel+1, prolongate_initial_data); BEGIN_MGLEVEL_LOOP(cgh) { @@ -294,8 +294,8 @@ namespace Carpet { cout << "3TL rl=" << reflevel << " restricting" << endl; Restrict (cgh); - Waypoint ("%*sScheduling PostRestrict", 2*reflevel, ""); - CCTK_ScheduleTraverse ("PostRestrict", cgh, CallFunction); + Waypoint ("%*sScheduling POSTRESTRICT", 2*reflevel, ""); + CCTK_ScheduleTraverse ("CCTK_POSTRESTRICT", cgh, CallFunction); // Checking CalculateChecksums (cgh, currenttime); @@ -370,8 +370,8 @@ namespace Carpet { // Restrict Restrict (cgh); - Waypoint ("%*sScheduling PostRestrict", 2*reflevel, ""); - CCTK_ScheduleTraverse ("PostRestrict", cgh, CallFunction); + Waypoint ("%*sScheduling POSTRESTRICT", 2*reflevel, ""); + CCTK_ScheduleTraverse ("CCTK_POSTRESTRICT", cgh, CallFunction); } END_MGLEVEL_LOOP; } END_REVERSE_REFLEVEL_LOOP; diff --git a/Carpet/Carpet/src/Recompose.cc b/Carpet/Carpet/src/Recompose.cc index 23642b8b3..9ec0714af 100644 --- a/Carpet/Carpet/src/Recompose.cc +++ b/Carpet/Carpet/src/Recompose.cc @@ -26,10 +26,12 @@ #include "carpet.hh" extern "C" { - static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Recompose.cc,v 1.46 2003/10/16 17:00:04 schnetter Exp $"; + static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Recompose.cc,v 1.47 2003/11/05 16:18:37 schnetter Exp $"; CCTK_FILEVERSION(Carpet_Carpet_Recompose_cc); } +#define DEBUG false // false or true + namespace Carpet { @@ -70,9 +72,16 @@ namespace Carpet { - void SplitRegions_AlongZ (const cGH* cgh, - vector<ibbox>& bbs, - vector<bbvect>& obs); + void SplitRegions (const cGH* cgh, + vector<ibbox>& bbs, + vector<bbvect>& obs); + void SplitRegions_AlongZ (const cGH* cgh, + vector<ibbox>& bbs, + vector<bbvect>& obs); + void SplitRegions_AlongDir (const cGH* cgh, + vector<ibbox>& bbs, + vector<bbvect>& obs, + const int dir); static void SplitRegions_Automatic_Recursively (bvect const & dims, int const nprocs, dvect const dshape, @@ -80,12 +89,12 @@ namespace Carpet { bbvect const & ob, vector<ibbox> & bbs, vector<bbvect> & obs); - static void SplitRegions_Automatic (const cGH* cgh, - vector<ibbox>& bbs, - vector<bbvect>& obs); - static void SplitRegions_AsSpecified (const cGH* cgh, - vector<ibbox>& bbs, - vector<bbvect>& obs); + void SplitRegions_Automatic (const cGH* cgh, + vector<ibbox>& bbs, + vector<bbvect>& obs); + static void SplitRegions_AsSpecified (const cGH* cgh, + vector<ibbox>& bbs, + vector<bbvect>& obs); static void MakeProcessors_RoundRobin (const cGH* cgh, const gh<dim>::rexts& bbss, @@ -149,7 +158,8 @@ namespace Carpet { - void Regrid (const cGH* cgh, const int initialise_upto) + void Regrid (const cGH* cgh, + const int initialise_from, const bool do_prolongate) { assert (mglevel == -1); assert (component == -1); @@ -170,7 +180,7 @@ namespace Carpet { int do_recompose = (*regrid_routine) (cgh, bbsss, obss, pss); assert (do_recompose >= 0); if (do_recompose == 0) return; - Recompose (cgh, bbsss, obss, pss, initialise_upto); + Recompose (cgh, bbsss, obss, pss, initialise_from, do_prolongate); } @@ -179,7 +189,8 @@ namespace Carpet { const gh<dim>::rexts& bbsss, const gh<dim>::rbnds& obss, const gh<dim>::rprocs& pss, - const int initialise_upto) + const int initialise_from, + const bool do_prolongate) { assert (mglevel == -1); assert (component == -1); @@ -191,7 +202,7 @@ namespace Carpet { OutputGridStructure (cgh, bbsss, obss, pss); // Recompose - hh->recompose (bbsss, obss, pss, initialise_upto); + hh->recompose (bbsss, obss, pss, initialise_from, do_prolongate); Output (cgh, hh); } @@ -274,7 +285,7 @@ namespace Carpet { } } - hh->recompose(bbsss, obss, pss); + hh->recompose(bbsss, obss, pss, 0, true); } @@ -386,6 +397,8 @@ namespace Carpet { if (CCTK_EQUALS (processor_topology, "along-z")) { SplitRegions_AlongZ (cgh, bbs, obs); + } else if (CCTK_EQUALS (processor_topology, "along-dir")) { + SplitRegions_AlongDir (cgh, bbs, obs, split_direction); } else if (CCTK_EQUALS (processor_topology, "automatic")) { SplitRegions_Automatic (cgh, bbs, obs); } else if (CCTK_EQUALS (processor_topology, "manual")) { @@ -400,6 +413,15 @@ namespace Carpet { void SplitRegions_AlongZ (const cGH* cgh, vector<ibbox>& bbs, vector<bbvect>& obs) { + SplitRegions_AlongDir (cgh, bbs, obs, 2); + } + + + + void SplitRegions_AlongDir (const cGH* cgh, vector<ibbox>& bbs, + vector<bbvect>& obs, + const int dir) + { // Something to do? if (bbs.size() == 0) return; @@ -409,6 +431,8 @@ namespace Carpet { assert (bbs.size() == 1); + assert (dir>=0 && dir<dim); + const ivect rstr = bbs[0].stride(); const ivect rlb = bbs[0].lower(); const ivect rub = bbs[0].upper() + rstr; @@ -420,19 +444,19 @@ namespace Carpet { ivect cstr = rstr; ivect clb = rlb; ivect cub = rub; - const int glonpz = (rub[dim-1] - rlb[dim-1]) / cstr[dim-1]; + const int glonpz = (rub[dir] - rlb[dir]) / cstr[dir]; const int locnpz = (glonpz + nprocs - 1) / nprocs; - const int zstep = locnpz * cstr[dim-1]; - clb[dim-1] = rlb[dim-1] + zstep * c; - cub[dim-1] = rlb[dim-1] + zstep * (c+1); - if (clb[dim-1] > rub[dim-1]) clb[dim-1] = rub[dim-1]; - if (cub[dim-1] > rub[dim-1]) cub[dim-1] = rub[dim-1]; - assert (clb[dim-1] <= cub[dim-1]); - assert (cub[dim-1] <= rub[dim-1]); + const int zstep = locnpz * cstr[dir]; + clb[dir] = rlb[dir] + zstep * c; + cub[dir] = rlb[dir] + zstep * (c+1); + if (clb[dir] > rub[dir]) clb[dir] = rub[dir]; + if (cub[dir] > rub[dir]) cub[dir] = rub[dir]; + assert (clb[dir] <= cub[dir]); + assert (cub[dir] <= rub[dir]); bbs[c] = ibbox(clb, cub-cstr, cstr); obs[c] = obnd; - if (c>0) obs[c][dim-1][0] = false; - if (c<nprocs-1) obs[c][dim-1][1] = false; + if (c>0) obs[c][dir][0] = false; + if (c<nprocs-1) obs[c][dir][1] = false; } } @@ -446,11 +470,13 @@ namespace Carpet { vector<ibbox> & bbs, vector<bbvect> & obs) { + if (DEBUG) cout << "SRAR enter" << endl; // check preconditions assert (nprocs >= 1); // are we done? if (all(dims)) { + if (DEBUG) cout << "SRAR bottom" << endl; // check precondition assert (nprocs == 1); @@ -460,6 +486,7 @@ namespace Carpet { obs.assign (1, ob); // return + if (DEBUG) cout << "SRAR exit" << endl; return; } @@ -479,7 +506,37 @@ namespace Carpet { } } assert (mydim>=0 && mydim<dim); - assert (mysize>0); + assert (mysize>=0); + if (DEBUG) cout << "SRAR mydim " << mydim << endl; + if (DEBUG) cout << "SRAR mysize " << mysize << endl; + + if (mysize == 0) { + // the bbox is empty + if (DEBUG) cout << "SRAR empty" << endl; + + // create the bboxes + bbs.clear(); + obs.clear(); + bbs.reserve(nprocs); + obs.reserve(nprocs); + + // create a new bbox + assert (bb.empty()); + bbvect const newob (vect<bool,2>(false)); + ibbox const newbb (bb); + if (DEBUG) cout << "SRAR " << mydim << " newbb " << newbb << endl; + if (DEBUG) cout << "SRAR " << mydim << " newob " << newob << endl; + + // store + bbs.insert (bbs.end(), nprocs, newbb); + obs.insert (obs.end(), nprocs, newob); + + // check postconditions + assert (bbs.size() == nprocs); + assert (obs.size() == nprocs); + if (DEBUG) cout << "SRAR exit" << endl; + return; + } // mark this direction as done assert (! dims[mydim]); @@ -488,6 +545,8 @@ namespace Carpet { // choose a number of slices for this direction int const nslices = min(nprocs, (int)floor(mysize * pow(nprocs/allsizes, 1.0/alldims) + 0.5)); assert (nslices <= nprocs); + if (DEBUG) cout << "SRAR " << mydim << " nprocs " << nprocs << endl; + if (DEBUG) cout << "SRAR " << mydim << " nslices " << nslices << endl; // split the remaining processors vector<int> mynprocs(nslices); @@ -501,6 +560,7 @@ namespace Carpet { sum_mynprocs += mynprocs[n]; } assert (sum_mynprocs == nprocs); + if (DEBUG) cout << "SRAR " << mydim << " mynprocs " << mynprocs << endl; // split the region vector<int> myslice(nslices); @@ -512,19 +572,23 @@ namespace Carpet { } else { myslice[n] = (int)floor(1.0 * slice_left * mynprocs[n] / nprocs_left + 0.5); } + assert (myslice[n] > 0); slice_left -= myslice[n]; nprocs_left -= mynprocs[n]; } assert (slice_left == 0); assert (nprocs_left == 0); + if (DEBUG) cout << "SRAR " << mydim << " myslice " << myslice << endl; // create the bboxes and recurse + if (DEBUG) cout << "SRAR " << mydim << ": create bboxes" << endl; bbs.clear(); obs.clear(); bbs.reserve(nprocs); obs.reserve(nprocs); ivect last_up; for (int n=0; n<nslices; ++n) { + if (DEBUG) cout << "SRAR " << mydim << " n " << n << endl; // create a new bbox ivect lo = bb.lower(); @@ -541,12 +605,16 @@ namespace Carpet { last_up = up; } ibbox newbb(lo, up, str); + if (DEBUG) cout << "SRAR " << mydim << " newbb " << newbb << endl; + if (DEBUG) cout << "SRAR " << mydim << " newob " << newob << endl; // recurse vector<ibbox> newbbs; vector<bbvect> newobs; SplitRegions_Automatic_Recursively (newdims, mynprocs[n], dshape, newbb, newob, newbbs, newobs); + if (DEBUG) cout << "SRAR " << mydim << " newbbs " << newbbs << endl; + if (DEBUG) cout << "SRAR " << mydim << " newobs " << newobs << endl; // store assert (newbbs.size() == mynprocs[n]); @@ -558,21 +626,28 @@ namespace Carpet { // check postconditions assert (bbs.size() == nprocs); assert (obs.size() == nprocs); + if (DEBUG) cout << "SRAR exit" << endl; } void SplitRegions_Automatic (const cGH* cgh, vector<ibbox>& bbs, vector<bbvect>& obs) { + if (DEBUG) cout << "SRA enter" << endl; // Something to do? if (bbs.size() == 0) return; const int nprocs = CCTK_nProcs(cgh); + if (DEBUG) cout << "SRA nprocs " << nprocs << endl; // if (nprocs==1) return; // split the processors + // nslices: number of disjoint bboxes int const nslices = bbs.size(); + if (DEBUG) cout << "SRA nslices " << nslices << endl; + // ncomps: number of components per processor int const ncomps = (nslices + nprocs - 1) / nprocs; + if (DEBUG) cout << "SRA ncomps " << ncomps << endl; assert (ncomps > 0); vector<int> mysize(nslices); for (int c=0; c<nslices; ++c) { @@ -580,32 +655,41 @@ namespace Carpet { } vector<int> mynprocs(nslices); { + if (DEBUG) cout << "SRA: distributing processors to slices" << endl; int ncomps_left = nprocs * ncomps; for (int c=0; c<nslices; ++c) { mynprocs[c] = 1; -- ncomps_left; } while (ncomps_left > 0) { + if (DEBUG) cout << "SRA ncomps_left " << ncomps_left << endl; int maxc = -1; - double maxratio = 0; + double maxratio = -1; for (int c=0; c<nslices; ++c) { double const ratio = (double)mysize[c] / mynprocs[c]; if (ratio > maxratio) { maxc=c; maxratio=ratio; } } assert (maxc>=0 && maxc<nslices); ++ mynprocs[maxc]; + if (DEBUG) cout << "SRA maxc " << maxc << endl; + if (DEBUG) cout << "SRA mynprocs[maxc] " << mynprocs[maxc] << endl; -- ncomps_left; } assert (ncomps_left == 0); } + if (DEBUG) cout << "SRA mynprocs " << mynprocs << endl; vector<ibbox> allbbs; vector<bbvect> allobs; + if (DEBUG) cout << "SRA: splitting regions" << endl; for (int c=0; c<nslices; ++c) { const ibbox bb = bbs[c]; const bbvect ob = obs[c]; + if (DEBUG) cout << "SRA c " << c << endl; + if (DEBUG) cout << "SRA bb " << bb << endl; + if (DEBUG) cout << "SRA ob " << ob << endl; const ivect rstr = bb.stride(); const ivect rlb = bb.lower(); @@ -613,20 +697,28 @@ namespace Carpet { // calculate real shape factors dvect dshape; - for (int d=0; d<dim; ++d) { - dshape[d] = (double)(rub[d]-rlb[d]) / (rub[0]-rlb[0]); + if (any(rub == rlb)) { + // the bbox is empty + dshape = 0.0; + } else { + for (int d=0; d<dim; ++d) { + dshape[d] = (double)(rub[d]-rlb[d]) / (rub[0]-rlb[0]); + } + const double dfact = pow(nprocs / prod(dshape), 1.0/dim); + dshape *= dfact; + assert (abs(prod(dshape) - nprocs) < 1e-6); } - const double dfact = pow(nprocs / prod(dshape), 1.0/dim); - dshape *= dfact; - assert (abs(prod(dshape) - nprocs) < 1e-6); + if (DEBUG) cout << "SRA shapes " << dshape << endl; bvect const dims = false; vector<ibbox> thebbs; vector<bbvect> theobs; - + SplitRegions_Automatic_Recursively (dims, mynprocs[c], dshape, bb, ob, thebbs, theobs); + if (DEBUG) cout << "SRA thebbs " << thebbs << endl; + if (DEBUG) cout << "SRA theobs " << theobs << endl; allbbs.insert(allbbs.end(), thebbs.begin(), thebbs.end()); allobs.insert(allobs.end(), theobs.begin(), theobs.end()); @@ -635,6 +727,8 @@ namespace Carpet { bbs = allbbs; obs = allobs; + + if (DEBUG) cout << "SRA exit" << endl; } @@ -649,8 +743,6 @@ namespace Carpet { const int nprocs = CCTK_nProcs(cgh); - if (nprocs==1) return; - assert (bbs.size() == 1); const ivect rstr = bbs[0].stride(); @@ -668,6 +760,8 @@ namespace Carpet { } assert (prod(nprocs_dir) == nprocs); + if (nprocs==1) return; + bbs.resize(nprocs); obs.resize(nprocs); const ivect cstr = rstr; @@ -688,11 +782,11 @@ namespace Carpet { // cub = min (cub, rub); for (int d=0; d<dim; ++d) { if (ipos[d]<rem[d]) { - clb[d] += step[d] * ipos[d]; - cub[d] += step[d] * (ipos[d]+1); + clb[d] += cstr[d] * ipos[d]; + cub[d] += cstr[d] * (ipos[d]+1); } else { - clb[d] += step[d] * rem[d]; - cub[d] += step[d] * rem[d]; + clb[d] += cstr[d] * rem[d]; + cub[d] += cstr[d] * rem[d]; } } assert (all (clb >= 0)); diff --git a/Carpet/Carpet/src/Restrict.cc b/Carpet/Carpet/src/Restrict.cc index e81a01771..06d87ef1a 100644 --- a/Carpet/Carpet/src/Restrict.cc +++ b/Carpet/Carpet/src/Restrict.cc @@ -10,7 +10,7 @@ #include "carpet.hh" extern "C" { - static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Restrict.cc,v 1.21 2003/08/10 21:59:51 schnetter Exp $"; + static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Restrict.cc,v 1.22 2003/11/05 16:18:37 schnetter Exp $"; CCTK_FILEVERSION(Carpet_Carpet_Restrict_cc); } @@ -24,47 +24,76 @@ namespace Carpet { void Restrict (const cGH* cgh) { + assert (reflevel != -1); + assert (mglevel != -1); assert (component == -1); Checkpoint ("%*sRestrict", 2*reflevel, ""); - // Loop over grid functions with storage - for (int group=0; group<CCTK_NumGroups(); ++group) { - if (CCTK_GroupTypeI(group) == CCTK_GF - && CCTK_QueryGroupStorageI(cgh, group)) { - if (arrdata[group].do_transfer) { - for (int var=0; var<(int)arrdata[group].data.size(); ++var) { - - const int tl = 0; - - // use background time here (which may not be modified by - // the user) - const CCTK_REAL time = tt->time (tl, reflevel, mglevel); - const CCTK_REAL time1 = tt->time (0, reflevel, mglevel); - const CCTK_REAL time2 = cgh->cctk_time / cgh->cctk_delta_time; - assert (fabs((time1 - time2) / (fabs(time1) + fabs(time2) + fabs(cgh->cctk_delta_time))) < 1e-12); - - if (reflevel < hh->reflevels()-1) { - - for (int c=0; c<hh->components(reflevel); ++c) { - arrdata[group].data[var]->ref_restrict - (tl, reflevel, c, mglevel, time); - } - for (int c=0; c<arrdata[group].hh->components(reflevel); ++c) { - arrdata[group].data[var]->sync (tl, reflevel, c, mglevel); - } - - } // if not finest refinement level - - } // loop over variables - } else { - char * groupname = CCTK_GroupName(group); - Checkpoint ("%*s(no restricting for group %s)", - 2*reflevel, "", groupname); - free (groupname); - } - } // if group has storage - } // loop over groups + // Restrict + if (reflevel < hh->reflevels()-1) { + for (comm_state<dim> state; !state.done(); state.step()) { + for (int group=0; group<CCTK_NumGroups(); ++group) { + if (CCTK_GroupTypeI(group) == CCTK_GF) { + if (CCTK_QueryGroupStorageI(cgh, group)) { + if (arrdata[group].do_transfer) { + for (int var=0; var<(int)arrdata[group].data.size(); ++var) { + + const int tl = 0; + + // use background time here (which may not be + // modified by the user) + const CCTK_REAL time = tt->time (tl, reflevel, mglevel); + const CCTK_REAL time1 = tt->time (0, reflevel, mglevel); + const CCTK_REAL time2 = cgh->cctk_time / cgh->cctk_delta_time; + assert (fabs((time1 - time2) / (fabs(time1) + fabs(time2) + fabs(cgh->cctk_delta_time))) < 1e-12); + + + for (int c=0; c<hh->components(reflevel); ++c) { + arrdata[group].data[var]->ref_restrict + (state, tl, reflevel, c, mglevel, time); + } + + } // loop over variables + } else { + if (state.thestate==state_recv) { + char * groupname = CCTK_GroupName(group); + Checkpoint ("%*s(no restricting for group %s)", + 2*reflevel, "", groupname); + free (groupname); + } + } // if ! do_transfer + } // if group has storage + } // if grouptype == CCTK_GF + } // loop over groups + } // for state + } // if not finest refinement level + + + + // Sync + if (reflevel < hh->reflevels()-1) { + for (comm_state<dim> state; !state.done(); state.step()) { + for (int group=0; group<CCTK_NumGroups(); ++group) { + if (CCTK_GroupTypeI(group) == CCTK_GF) { + if (CCTK_QueryGroupStorageI(cgh, group)) { + if (arrdata[group].do_transfer) { + for (int var=0; var<(int)arrdata[group].data.size(); ++var) { + + const int tl = 0; + + for (int c=0; c<arrdata[group].hh->components(reflevel); ++c) { + arrdata[group].data[var]->sync + (state, tl, reflevel, c, mglevel); + } + + } // loop over variables + } // if do_transfer + } // if group has storage + } // if grouptype == CCTK_GF + } // loop over groups + } // for state + } // if not finest refinement level } diff --git a/Carpet/Carpet/src/SetupGH.cc b/Carpet/Carpet/src/SetupGH.cc index 89d73343e..676ba350a 100644 --- a/Carpet/Carpet/src/SetupGH.cc +++ b/Carpet/Carpet/src/SetupGH.cc @@ -20,7 +20,7 @@ #include "carpet.hh" extern "C" { - static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/SetupGH.cc,v 1.53 2003/09/19 16:08:37 schnetter Exp $"; + static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/SetupGH.cc,v 1.54 2003/11/05 16:18:37 schnetter Exp $"; CCTK_FILEVERSION(Carpet_Carpet_SetupGH_cc); } @@ -245,7 +245,7 @@ namespace Carpet { case CCTK_SCALAR: // treat scalars as DIM=0, DISTRIB=const arrays - arrdata[group].info.dim = 0; + assert (arrdata[group].info.dim == 0); disttype = CCTK_DISTRIB_CONSTANT; break; @@ -270,7 +270,12 @@ namespace Carpet { || disttype==CCTK_DISTRIB_DEFAULT); if (disttype==CCTK_DISTRIB_CONSTANT) { - sizes[dim-1] = (sizes[dim-1] - 2*ghostsizes[dim-1]) * CCTK_nProcs(cgh) + 2*ghostsizes[dim-1]; + for (int d=0; d<dim; ++d) { + ghostsizes[d] = 0; + } + const int d = gp.dim==0 ? 0 : gp.dim-1; + sizes[d] = (sizes[d] - 2*ghostsizes[d]) * CCTK_nProcs(cgh) + 2*ghostsizes[d]; + assert (sizes[d] >= 0); } const vect<int,dim> alb(0); @@ -302,7 +307,11 @@ namespace Carpet { obs.push_back (vect<vect<bool,2>,dim>(vect<bool,2>(true))); // Split it into components, one for each processor - SplitRegions_AlongZ (cgh, bbs, obs); + if (disttype==CCTK_DISTRIB_CONSTANT) { + SplitRegions_AlongDir (cgh, bbs, obs, gp.dim==0 ? 0 : gp.dim-1); + } else { + SplitRegions_Automatic (cgh, bbs, obs); + } // For all refinement levels (but there is only one) vector<vector<bbox<int,dim> > > bbss(1); @@ -323,7 +332,7 @@ namespace Carpet { assert (groupname); Checkpoint ("Recomposing grid array group %s", groupname); free (groupname); - arrdata[group].hh->recompose (bbsss, obss, pss); + arrdata[group].hh->recompose (bbsss, obss, pss, 1, false); break; } @@ -383,14 +392,24 @@ namespace Carpet { for (int d=0; d<dim; ++d) { ((int*)arrdata[group].info.gsh )[d] = (base.shape() / base.stride())[d]; ((int*)arrdata[group].info.lsh )[d] = (ext.shape() / ext.stride())[d]; - ((int*)arrdata[group].info.lbnd)[d] = (ext.lower() / ext.stride())[d]; - ((int*)arrdata[group].info.ubnd)[d] = (ext.upper() / ext.stride())[d]; - for (int f=0; f<2; ++f) { - ((int*)arrdata[group].info.bbox)[2*d+f] = obnds[d][f]; - } + ((int*)arrdata[group].info.lbnd)[d] = ((ext.lower() - baseext.lower()) / ext.stride())[d]; + ((int*)arrdata[group].info.ubnd)[d] = ((ext.upper() - baseext.lower()) / ext.stride())[d]; + ((int*)arrdata[group].info.bbox)[2*d ] = obnds[d][0]; + ((int*)arrdata[group].info.bbox)[2*d+1] = obnds[d][1]; assert (arrdata[group].info.lsh[d]>=0); assert (arrdata[group].info.lsh[d]<=arrdata[group].info.gsh[d]); + if (d>=arrdata[group].info.dim) + if (! (arrdata[group].info.lsh[d]==1)) { + cout << "group " << group << " " << CCTK_GroupName(group) << " " + << "dim " << arrdata[group].info.dim << " " + << "d " << d << endl; + cout << "gsh " << arrdata[group].info.gsh[0] << " " << arrdata[group].info.gsh[1] << " " << arrdata[group].info.gsh[2] << endl; + cout << "lsh " << arrdata[group].info.lsh[0] << " " << arrdata[group].info.lsh[1] << " " << arrdata[group].info.lsh[2] << endl; + cout << "lbnd " << arrdata[group].info.lbnd[0] << " " << arrdata[group].info.lbnd[1] << " " << arrdata[group].info.lbnd[2] << endl; + } + if (d>=arrdata[group].info.dim) + assert (arrdata[group].info.lsh[d]==1); assert (arrdata[group].info.lbnd[d]>=0); assert (arrdata[group].info.lbnd[d]<=arrdata[group].info.ubnd[d]+1); assert (arrdata[group].info.ubnd[d]<arrdata[group].info.gsh[d]); @@ -403,7 +422,10 @@ namespace Carpet { if (numvars>0) { const int firstvar = CCTK_FirstVarIndexI (group); assert (firstvar>=0); + const int num_tl = CCTK_NumTimeLevelsFromVarI (firstvar); + // We don't know when to cycle arrays or scalars + assert (num_tl==1); assert (rl>=0 && rl<(int)arrdata[group].dd->boxes.size()); assert (c>=0 && c<(int)arrdata[group].dd->boxes[rl].size()); @@ -471,7 +493,7 @@ namespace Carpet { // Recompose grid hierarchy Checkpoint ("Recomposing grid functions"); - Recompose (cgh, bbsss, obss, pss); + Recompose (cgh, bbsss, obss, pss, maxreflevels, false); diff --git a/Carpet/Carpet/src/carpet.hh b/Carpet/Carpet/src/carpet.hh index cb13a6281..9d50f0eb0 100644 --- a/Carpet/Carpet/src/carpet.hh +++ b/Carpet/Carpet/src/carpet.hh @@ -1,4 +1,4 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/carpet.hh,v 1.25 2003/08/10 21:59:51 schnetter Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/carpet.hh,v 1.26 2003/11/05 16:18:37 schnetter Exp $ #ifndef CARPET_HH #define CARPET_HH @@ -11,7 +11,8 @@ namespace Carpet { - void Regrid (const cGH* cgh, const int initialise_upto = -1); + void Regrid (const cGH* cgh, + const int initialise_from, const bool do_prolongate); void CycleTimeLevels (const cGH* cgh); void FlipTimeLevels (const cGH* cgh); void Restrict (const cGH* cgh); @@ -20,7 +21,8 @@ namespace Carpet { const gh<dim>::rexts& bbsss, const gh<dim>::rbnds& obss, const gh<dim>::rprocs& pss, - const int initialise_upto = -1); + const int initialise_from, + const bool do_prolongate); enum checktimes { currenttime, currenttimebutnotifonly, diff --git a/Carpet/Carpet/src/carpet_public.hh b/Carpet/Carpet/src/carpet_public.hh index 61e4fbe28..64ce2703b 100644 --- a/Carpet/Carpet/src/carpet_public.hh +++ b/Carpet/Carpet/src/carpet_public.hh @@ -1,4 +1,4 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/carpet_public.hh,v 1.38 2003/10/13 11:44:51 schnetter Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/carpet_public.hh,v 1.39 2003/11/05 16:18:37 schnetter Exp $ // It is assumed that the number of components of all arrays is equal // to the number of components of the grid functions, and that their @@ -149,6 +149,11 @@ namespace Carpet { vector<vect<vect<bool,2>,dim> >& obs); void SplitRegions_AlongZ (const cGH* cgh, vector<bbox<int,dim> >& bbs, vector<vect<vect<bool,2>,dim> >& obs); + void SplitRegions_AlongDir (const cGH* cgh, vector<bbox<int,dim> >& bbs, + vector<vect<vect<bool,2>,dim> >& obs, + const int dir); + void SplitRegions_Automatic (const cGH* cgh, vector<bbox<int,dim> >& bbs, + vector<vect<vect<bool,2>,dim> >& obs); void MakeProcessors (const cGH* cgh, const gh<dim>::rexts& bbsss, gh<dim>::rprocs& pss); diff --git a/Carpet/Carpet/src/helpers.cc b/Carpet/Carpet/src/helpers.cc index 028ed3c6b..937f745a0 100644 --- a/Carpet/Carpet/src/helpers.cc +++ b/Carpet/Carpet/src/helpers.cc @@ -15,7 +15,7 @@ #include "carpet.hh" extern "C" { - static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/helpers.cc,v 1.41 2003/08/10 21:59:51 schnetter Exp $"; + static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/helpers.cc,v 1.42 2003/11/05 16:18:37 schnetter Exp $"; CCTK_FILEVERSION(Carpet_Carpet_helpers_cc); } @@ -266,8 +266,8 @@ namespace Carpet { assert (mglevelfact==1); cgh->cctk_timefac = reflevelfact / mglevelfact; for (int d=0; d<dim; ++d) { - assert (baseext.lower()[d] * reflevelfact % maxreflevelfact == 0); - cgh->cctk_levoff[d] = baseext.lower()[d] * reflevelfact / maxreflevelfact; + assert (baseext.lower()[d] % (maxreflevelfact/reflevelfact) == 0); + cgh->cctk_levoff[d] = baseext.lower()[d] / (maxreflevelfact/reflevelfact); cgh->cctk_levoffdenom[d] = 1; } @@ -384,7 +384,7 @@ namespace Carpet { for (int d=0; d<dim; ++d) { - ((int*)arrdata[group].info.lsh)[d] = (ext.shape() / ext.stride())[d]; + ((int*)arrdata[group].info.lsh )[d] = (ext.shape() / ext.stride())[d]; ((int*)arrdata[group].info.lbnd)[d] = ((ext.lower() - baseext.lower()) / ext.stride())[d]; ((int*)arrdata[group].info.ubnd)[d] = ((ext.upper() - baseext.lower()) / ext.stride())[d]; ((int*)arrdata[group].info.bbox)[2*d ] = obnds[d][0]; @@ -392,6 +392,8 @@ namespace Carpet { assert (arrdata[group].info.lsh[d]>=0); assert (arrdata[group].info.lsh[d]<=arrdata[group].info.gsh[d]); + if (d>=arrdata[group].info.dim) + assert (arrdata[group].info.lsh[d]==1); assert (arrdata[group].info.lbnd[d]>=0); assert (arrdata[group].info.lbnd[d]<=arrdata[group].info.ubnd[d]+1); assert (arrdata[group].info.ubnd[d]<arrdata[group].info.gsh[d]); @@ -458,10 +460,10 @@ namespace Carpet { } extern "C" void CCTK_FCALL CCTK_FNAME(CallScheduleGroup) - (int * const ierr, cGH * const cgh, ONE_FORTSTRING_ARG) + (int * const ierr, cGH * const * const cgh, ONE_FORTSTRING_ARG) { ONE_FORTSTRING_CREATE (group); - *ierr = CallScheduleGroup (cgh, group); + *ierr = CallScheduleGroup (*cgh, group); free (group); } diff --git a/Carpet/CarpetIOASCII/schedule.ccl b/Carpet/CarpetIOASCII/schedule.ccl index aff6826a2..3c7f55c59 100644 --- a/Carpet/CarpetIOASCII/schedule.ccl +++ b/Carpet/CarpetIOASCII/schedule.ccl @@ -1,5 +1,5 @@ # Schedule definitions for thorn CarpetIOASCII -# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOASCII/schedule.ccl,v 1.3 2001/03/22 18:42:05 eschnett Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOASCII/schedule.ccl,v 1.4 2003/11/05 16:18:37 schnetter Exp $ schedule CarpetIOASCIIStartup at STARTUP after IOUtil_Startup { diff --git a/Carpet/CarpetIOASCII/src/ioascii.cc b/Carpet/CarpetIOASCII/src/ioascii.cc index 896fcce42..eabe10ef7 100644 --- a/Carpet/CarpetIOASCII/src/ioascii.cc +++ b/Carpet/CarpetIOASCII/src/ioascii.cc @@ -30,7 +30,7 @@ #include "ioascii.hh" extern "C" { - static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOASCII/src/ioascii.cc,v 1.53 2003/10/14 16:39:16 schnetter Exp $"; + static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOASCII/src/ioascii.cc,v 1.54 2003/11/05 16:18:37 schnetter Exp $"; CCTK_FILEVERSION(Carpet_CarpetIOASCII_ioascii_cc); } @@ -829,7 +829,9 @@ namespace CarpetIOASCII { gdata<D>* const tmp = gfdata->make_typed(vi); tmp->allocate(gfdata->extent(), 0); - tmp->copy_from (gfdata, gfdata->extent()); + for (comm_state<dim> state; !state.done(); state.step()) { + tmp->copy_from (state, gfdata, gfdata->extent()); + } WriteASCII (os, tmp, gfext, vi, time, org, dirs, tl, rl, c, ml, coord_time, coord_lower, coord_upper); delete tmp; diff --git a/Carpet/CarpetIOASCII/src/ioascii.h b/Carpet/CarpetIOASCII/src/ioascii.h index 5e677843c..d2e5c549b 100644 --- a/Carpet/CarpetIOASCII/src/ioascii.h +++ b/Carpet/CarpetIOASCII/src/ioascii.h @@ -1,4 +1,4 @@ -/* $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOASCII/src/ioascii.h,v 1.4 2002/09/01 14:52:25 schnetter Exp $ */ +/* $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOASCII/src/ioascii.h,v 1.5 2003/11/05 16:18:37 schnetter Exp $ */ #ifndef CARPETIOASCII_H #define CARPETIOASCII_H diff --git a/Carpet/CarpetInterp/src/interp.cc b/Carpet/CarpetInterp/src/interp.cc index 75ae5675d..b018768dc 100644 --- a/Carpet/CarpetInterp/src/interp.cc +++ b/Carpet/CarpetInterp/src/interp.cc @@ -1,4 +1,4 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetInterp/src/interp.cc,v 1.13 2003/10/14 16:39:16 schnetter Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetInterp/src/interp.cc,v 1.14 2003/11/05 16:18:38 schnetter Exp $ #include <assert.h> #include <math.h> @@ -12,6 +12,8 @@ #include "bbox.hh" #include "data.hh" +#include "defs.hh" +#include "ggf.hh" #include "vect.hh" #include "carpet.hh" @@ -19,7 +21,7 @@ #include "interp.hh" extern "C" { - static char const * const rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetInterp/src/interp.cc,v 1.13 2003/10/14 16:39:16 schnetter Exp $"; + static char const * const rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetInterp/src/interp.cc,v 1.14 2003/11/05 16:18:38 schnetter Exp $"; CCTK_FILEVERSION(Carpet_CarpetInterp_interp_cc); } @@ -112,9 +114,10 @@ namespace CarpetInterp { delta[d] = (upper[d] - lower[d]) / (hh->baseextent.shape()[d] - hh->baseextent.stride()[d]); } + assert (N_interp_points >= 0); assert (interp_coords); for (int d=0; d<dim; ++d) { - assert (interp_coords[d]); + assert (N_interp_points==0 || interp_coords[d]); } @@ -135,6 +138,7 @@ namespace CarpetInterp { +#if 0 // Assert that all refinement levels have the same time // TODO: interpolate in time { @@ -152,13 +156,21 @@ namespace CarpetInterp { "Cannot interpolate in global mode at iteration %d (time %g), because not all refinement levels exist at this time. Interpolation in time is not yet implemented.", cgh->cctk_iteration, (double)cgh->cctk_time); } else { - CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, - "Cannot interpolate in level mode at iteration %d (time %g), because refinement level %d does not exist at this time. Interpolation in time is not yet implemented.", - cgh->cctk_iteration, (double)cgh->cctk_time, reflevel); + // called in level mode at a time where the level does not + // exist + CCTK_WARN (0, "internal error"); } assert (0); } } +#endif + + // Find the time interpolation order + int partype; + void const * const parptr = CCTK_ParameterGet ("prolongation_order_time", "Carpet", &partype); + assert (parptr); + assert (partype == PARAMETER_INTEGER); + int const prolongation_order_time = * (CCTK_INT const *) parptr; @@ -181,8 +193,10 @@ namespace CarpetInterp { home[n] = -1; for (int rl=maxrl-1; rl>=minrl; --rl) { + const int fact = maxreflevelfact / ipow(reffact,rl); CCTK_REAL (* const rfloor) (CCTK_REAL const) = floor; - ivect const ipos = ivect(map(rfloor, (pos - lower) / delta + 0.5)); + ivect const ipos = ivect(map(rfloor, (pos - lower) / (delta * fact) + 0.5)) * fact; + assert (all(ipos % hh->bases[rl][ml].stride() == 0)); // TODO: use something faster than a linear search for (int c=0; c<hh->components(rl); ++c) { @@ -251,10 +265,13 @@ namespace CarpetInterp { } // Transfer coordinate patches - for (int p=0; p<nprocs; ++p) { - for (int rl=minrl; rl<maxrl; ++rl) { - for (int c=0; c<hh->components(rl); ++c) { - allcoords[ind_prc(p,rl,c)].change_processor (hh->processors[rl][c]); + for (comm_state<dim> state; !state.done(); state.step()) { + for (int p=0; p<nprocs; ++p) { + for (int rl=minrl; rl<maxrl; ++rl) { + for (int c=0; c<hh->components(rl); ++c) { + allcoords[ind_prc(p,rl,c)].change_processor + (state, hh->processors[rl][c]); + } } } } @@ -298,6 +315,15 @@ namespace CarpetInterp { if (reflevel>=minrl && reflevel<maxrl) { BEGIN_MGLEVEL_LOOP(cgh) { + // Number of necessary time levels + int const tl = 0; + CCTK_REAL const time1 = tt->time (tl, reflevel, mglevel); + CCTK_REAL const time2 = cgh->cctk_time / cgh->cctk_delta_time; + bool const need_time_interp + = fabs((time1 - time2) / (fabs(time1) + fabs(time2) + fabs(cgh->cctk_delta_time))) > 1e-12; + int const num_tl + = need_time_interp ? prolongation_order_time + 1 : 1; + BEGIN_LOCAL_COMPONENT_LOOP(cgh, CCTK_GF) { // Find out about the local geometry @@ -313,13 +339,13 @@ namespace CarpetInterp { // Find out about the grid functions vector<CCTK_INT> input_array_type_codes(N_input_arrays); - vector<const void *> input_arrays(N_input_arrays); + vector<const void *> input_arrays(N_input_arrays * num_tl); for (int n=0; n<N_input_arrays; ++n) { if (input_array_variable_indices[n] == -1) { // Ignore this entry - input_array_type_codes[n] = -1; - input_arrays[n] = 0; + input_array_type_codes[tl+n*num_tl] = -1; + input_arrays[tl+n*num_tl] = 0; } else { @@ -337,8 +363,12 @@ namespace CarpetInterp { assert (group.disttype == CCTK_DISTRIB_DEFAULT); assert (group.stagtype == 0); // not staggered - input_array_type_codes[n] = group.vartype; - input_arrays[n] = CCTK_VarDataPtrI (cgh, 0, vi); + assert (group.numtimelevels >= num_tl); + + for (int tl=0; tl<num_tl; ++tl) { + input_array_type_codes[tl+n*num_tl] = group.vartype; + input_arrays[tl+n*num_tl] = CCTK_VarDataPtrI (cgh, tl, vi); + } } } // for input arrays @@ -353,30 +383,85 @@ namespace CarpetInterp { assert (allhomecnts[ind_prc(p,reflevel,component)] == alloutputs[ind_prc(p,reflevel,component)].shape()[0]); + int const npoints = allhomecnts[ind_prc(p,reflevel,component)]; + // Do the processor-local interpolation vector<const void *> tmp_interp_coords (dim); for (int d=0; d<dim; ++d) { tmp_interp_coords[d] = &allcoords[ind_prc(p,reflevel,component)][ivect(0,d,0)]; } - vector<void *> tmp_output_arrays (N_output_arrays); - for (int m=0; m<N_output_arrays; ++m) { - assert (output_array_type_codes[m] == CCTK_VARIABLE_REAL); - tmp_output_arrays[m] - = &alloutputs[ind_prc(p,reflevel,component)][ivect(0,m,0)]; + vector<void *> tmp_output_arrays (N_output_arrays * num_tl); + if (need_time_interp) { + for (int m=0; m<N_output_arrays; ++m) { + assert (output_array_type_codes[m] == CCTK_VARIABLE_REAL); + for (int tl=0; tl<num_tl; ++tl) { + tmp_output_arrays[tl+m*num_tl] = new CCTK_REAL [npoints]; + } + } + } else { + for (int m=0; m<N_output_arrays; ++m) { + assert (output_array_type_codes[m] == CCTK_VARIABLE_REAL); + tmp_output_arrays[m] + = &alloutputs[ind_prc(p,reflevel,component)][ivect(0,m,0)]; + } } ierr = CCTK_InterpLocalUniform (N_dims, local_interp_handle, param_table_handle, &coord_origin[0], &coord_delta[0], - allhomecnts[ind_prc(p,reflevel,component)], + npoints, interp_coords_type_code, &tmp_interp_coords[0], - N_input_arrays, &lsh[0], + N_input_arrays * num_tl, &lsh[0], &input_array_type_codes[0], &input_arrays[0], - N_output_arrays, + N_output_arrays * num_tl, output_array_type_codes, &tmp_output_arrays[0]); assert (!ierr); + // Interpolate in time, if necessary + if (need_time_interp) { + CCTK_REAL const time = cgh->cctk_time / cgh->cctk_delta_time; + CCTK_REAL times[num_tl]; + for (int tl=0; tl<num_tl; ++tl) { + times[tl] = tt->time (tl, reflevel, mglevel); + } + CCTK_REAL tfacs[num_tl]; + switch (num_tl) { + case 1: + // no interpolation + assert (fabs((time - times[0]) / fabs(time + times[0] + cgh->cctk_delta_time)) < 1e-12); + tfacs[0] = 1.0; + break; + case 2: + // linear (2-point) interpolation + tfacs[0] = (time - times[1]) / (times[0] - times[1]); + tfacs[1] = (time - times[0]) / (times[1] - times[0]); + break; + case 3: + // quadratic (3-point) interpolation + tfacs[0] = (time - times[1]) * (time - times[2]) / ((times[0] - times[1]) * (times[0] - times[2])); + tfacs[1] = (time - times[0]) * (time - times[2]) / ((times[1] - times[0]) * (times[1] - times[2])); + tfacs[2] = (time - times[0]) * (time - times[1]) / ((times[2] - times[0]) * (times[2] - times[1])); + break; + default: + assert (0); + } + for (int m=0; m<N_output_arrays; ++m) { + assert (output_array_type_codes[m] == CCTK_VARIABLE_REAL); + for (int k=0; k<npoints; ++k) { + CCTK_REAL & dest = alloutputs[ind_prc(p,reflevel,component)][ivect(k,m,0)]; + dest = 0; + for (int tl=0; tl<num_tl; ++tl) { + CCTK_REAL const src = ((CCTK_REAL const *)tmp_output_arrays[tl+m*num_tl])[k]; + dest += tfacs[tl] * src; + } + } + for (int tl=0; tl<num_tl; ++tl) { + delete [] (CCTK_REAL *)tmp_output_arrays[tl+m*num_tl]; + } + } + } + } // for processors } END_LOCAL_COMPONENT_LOOP; @@ -396,14 +481,18 @@ namespace CarpetInterp { // Transfer output patches back - for (int p=0; p<nprocs; ++p) { - for (int rl=0; rl<minrl; ++rl) { - for (int c=0; c<hh->components(rl); ++c) { - alloutputs[ind_prc(p,rl,c)].change_processor (p); + for (comm_state<dim> state; !state.done(); state.step()) { + for (int p=0; p<nprocs; ++p) { + for (int rl=minrl; rl<maxrl; ++rl) { + for (int c=0; c<hh->components(rl); ++c) { + alloutputs[ind_prc(p,rl,c)].change_processor (state, p); + } } } } + + // Read out local output patches { vector<int> tmpcnts ((maxrl-minrl) * maxncomps); @@ -412,7 +501,7 @@ namespace CarpetInterp { int const c = home[n]; for (int m=0; m<N_output_arrays; ++m) { assert (interp_coords_type_code == CCTK_VARIABLE_REAL); - assert (dim==3); + assert (alloutputs[ind_prc(myproc,rl,c)].owns_storage()); static_cast<CCTK_REAL *>(output_arrays[m])[n] = alloutputs[ind_prc(myproc,rl,c)][ivect(tmpcnts[ind_rc(rl,c)],m,0)]; } diff --git a/Carpet/CarpetLib/param.ccl b/Carpet/CarpetLib/param.ccl index 4316c4dd7..7731ef890 100644 --- a/Carpet/CarpetLib/param.ccl +++ b/Carpet/CarpetLib/param.ccl @@ -1,5 +1,5 @@ # Parameter definitions for thorn CarpetLib -# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/param.ccl,v 1.5 2003/07/20 21:03:43 schnetter Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/param.ccl,v 1.6 2003/11/05 16:18:38 schnetter Exp $ private: @@ -7,6 +7,10 @@ BOOLEAN verbose "Print info to the screen" { } "no" +BOOLEAN check_array_accesses "Check all array accesses in Fortran" +{ +} "no" + BOOLEAN barriers "Insert barriers at strategic places for debugging purposes (slows down execution)" { } "no" diff --git a/Carpet/CarpetLib/src/copy_3d_int4.F77 b/Carpet/CarpetLib/src/copy_3d_int4.F77 index b440a32b2..4cd070277 100644 --- a/Carpet/CarpetLib/src/copy_3d_int4.F77 +++ b/Carpet/CarpetLib/src/copy_3d_int4.F77 @@ -1,5 +1,5 @@ c -*-Fortran-*- -c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/copy_3d_int4.F77,v 1.2 2004/03/11 12:03:09 schnetter Exp $ +c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/copy_3d_int4.F77,v 1.1 2003/11/05 16:18:39 schnetter Exp $ #include "cctk.h" #include "cctk_Parameters.h" @@ -54,11 +54,6 @@ c bbox(:,3) is stride c This could be handled, but is likely to point to an error elsewhere call CCTK_WARN (0, "Internal error: region extent is empty") end if - if (mod(srcbbox(d,2) - srcbbox(d,1), srcbbox(d,3)).ne.0 - $ .or. mod(dstbbox(d,2) - dstbbox(d,1), dstbbox(d,3)).ne.0 - $ .or. mod(regbbox(d,2) - regbbox(d,1), regbbox(d,3)).ne.0) then - call CCTK_WARN (0, "Internal error: array extents are not integer multiples of the strides") - end if if (regbbox(d,1).lt.srcbbox(d,1) $ .or. regbbox(d,1).lt.dstbbox(d,1) $ .or. regbbox(d,2).gt.srcbbox(d,2) diff --git a/Carpet/CarpetLib/src/copy_3d_real8.F77 b/Carpet/CarpetLib/src/copy_3d_real8.F77 index d115418b8..43507cdd3 100644 --- a/Carpet/CarpetLib/src/copy_3d_real8.F77 +++ b/Carpet/CarpetLib/src/copy_3d_real8.F77 @@ -1,18 +1,8 @@ c -*-Fortran-*- -c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/copy_3d_real8.F77,v 1.6 2003/02/25 22:57:00 schnetter Exp $ +c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/copy_3d_real8.F77,v 1.7 2003/11/05 16:18:39 schnetter Exp $ #include "cctk.h" - - - -#define CHKIDX(i,j,k, imax,jmax,kmax, where) \ - if ((i).lt.1 .or. (i).gt.(imax) \ - .or. (j).lt.1 .or. (j).gt.(jmax) \ - .or. (k).lt.1 .or. (k).gt.(kmax)) then &&\ - write (msg, '(a, " array index out of bounds: shape is (",i4,",",i4,",",i4,"), index is (",i4,",",i4,",",i4,")")') \ - (where), (imax), (jmax), (kmax), (i), (j), (k) &&\ - call CCTK_WARN (0, msg(1:len_trim(msg))) &&\ - end if +#include "cctk_Parameters.h" @@ -23,6 +13,8 @@ c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/copy_3d_real8 implicit none + DECLARE_CCTK_PARAMETERS + integer srciext, srcjext, srckext CCTK_REAL8 src(srciext,srcjext,srckext) integer dstiext, dstjext, dstkext @@ -96,16 +88,19 @@ c This could be handled, but is likely to point to an error elsewhere c Loop over region - do k = 0, regkext-1 - do j = 0, regjext-1 - do i = 0, regiext-1 + do k = 1, regkext + do j = 1, regjext + do i = 1, regiext + + if (check_array_accesses.ne.0) then + call checkindex (srcioff+i, srcjoff+j+1, srckoff+k+1, 1,1,1, + $ "source") + call checkindex (dstioff+i, dstjoff+j+1, dstkoff+k+1, 1,1,1, + $ "destination") + end if - CHKIDX (srcioff+i+1, srcjoff+j+1, srckoff+k+1, \ - srciext,srcjext,srckext, "source") - CHKIDX (dstioff+i+1, dstjoff+j+1, dstkoff+k+1, \ - dstiext,dstjext,dstkext, "destination") - dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1) - $ = src (srcioff+i+1, srcjoff+j+1, srckoff+k+1) + dst (dstioff+i, dstjoff+j, dstkoff+k) + $ = src (srcioff+i, srcjoff+j, srckoff+k) end do end do diff --git a/Carpet/CarpetLib/src/data.cc b/Carpet/CarpetLib/src/data.cc index d2177364a..46ff5df79 100644 --- a/Carpet/CarpetLib/src/data.cc +++ b/Carpet/CarpetLib/src/data.cc @@ -1,4 +1,4 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.cc,v 1.32 2003/10/17 08:14:41 cvs_anon Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.cc,v 1.33 2003/11/05 16:18:39 schnetter Exp $ #include <assert.h> #include <limits.h> @@ -17,24 +17,36 @@ #include "data.hh" -#include "util_ErrorCodes.h" -#include "util_Table.h" - using namespace std; +// Hand out the next MPI tag +static int nexttag () +{ + static int last = 100; + ++last; + if (last > 30000) last = 100; + return last; +} + + + // Constructors template<class T, int D> data<T,D>::data (const int varindex_) : gdata<D>(varindex_), - _storage(0) + _storage(0), + comm_active(false), + tag(nexttag()) { } template<class T, int D> data<T,D>::data (const int varindex_, const ibbox& extent_, const int proc_) : gdata<D>(varindex_), - _storage(0) + _storage(0), + comm_active(false), + tag(nexttag()) { allocate(extent_, proc_); } @@ -102,7 +114,32 @@ void data<T,D>::transfer_from (gdata<D>* gsrc) { // Processor management template<class T, int D> -void data<T,D>::change_processor (const int newproc, void* const mem) { +void data<T,D>::change_processor (comm_state<D>& state, + const int newproc, void* const mem) +{ + switch (state.thestate) { + case state_recv: + change_processor_recv (newproc, mem); + break; + case state_send: + change_processor_send (newproc, mem); + break; + case state_wait: + change_processor_wait (newproc, mem); + break; + default: + assert(0); + } +} + + + +template<class T, int D> +void data<T,D>::change_processor_recv (const int newproc, void* const mem) +{ + assert (!comm_active); + comm_active = true; + if (newproc == this->_proc) { assert (!mem); return; @@ -121,21 +158,87 @@ void data<T,D>::change_processor (const int newproc, void* const mem) { } else { _storage = (T*)mem; } - + T dummy; - MPI_Status status; - MPI_Recv (_storage, this->_size, dist::datatype(dummy), this->_proc, - dist::tag, dist::comm, &status); + MPI_Irecv (_storage, this->_size, dist::datatype(dummy), this->_proc, + this->tag, dist::comm, &request); + + } else if (rank == this->_proc) { + // copy to other processor + + } else { + assert (!mem); + assert (!_storage); + } + } +} + + +template<class T, int D> +void data<T,D>::change_processor_send (const int newproc, void* const mem) +{ + assert (comm_active); + + if (newproc == this->_proc) { + assert (!mem); + return; + } + + if (this->_has_storage) { + int rank; + MPI_Comm_rank (dist::comm, &rank); + if (rank == newproc) { + // copy from other processor + } else if (rank == this->_proc) { // copy to other processor assert (!mem); assert (_storage); + T dummy; - MPI_Send (_storage, this->_size, dist::datatype(dummy), newproc, - dist::tag, dist::comm); + MPI_Isend (_storage, this->_size, dist::datatype(dummy), newproc, + this->tag, dist::comm, &request); + + } else { + assert (!mem); + assert (!_storage); + } + } +} + + +template<class T, int D> +void data<T,D>::change_processor_wait (const int newproc, void* const mem) +{ + assert (comm_active); + comm_active = false; + + if (newproc == this->_proc) { + assert (!mem); + return; + } + + if (this->_has_storage) { + int rank; + MPI_Comm_rank (dist::comm, &rank); + if (rank == newproc) { + // copy from other processor + + MPI_Status status; + MPI_Wait (&request, &status); + + } else if (rank == this->_proc) { + // copy to other processor + + assert (!mem); + assert (_storage); + + MPI_Status status; + MPI_Wait (&request, &status); + if (this->_owns_storage) { delete [] _storage; } @@ -146,7 +249,7 @@ void data<T,D>::change_processor (const int newproc, void* const mem) { assert (!_storage); } } - + this->_proc = newproc; } @@ -169,6 +272,13 @@ void data<T,D> && (box.lower()-src->extent().lower())%box.stride() == 0)); assert (this->proc() == src->proc()); + + const int groupindex = CCTK_GroupIndexFromVarI(varindex); + const int group_tags_table = CCTK_GroupTagsTableI(groupindex); + assert (group_tags_table >= 0); + + // Disallow this. + assert (0); int rank; MPI_Comm_rank (dist::comm, &rank); @@ -212,15 +322,27 @@ void data<T,D> MPI_Comm_rank (dist::comm, &rank); assert (rank == this->proc()); - T Tdummy; - CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING, - "There is no interpolator available for variable type %s, dimension %d, spatial interpolation order %d, temporal interpolation order %d. The interpolation will not be done.", - typestring(Tdummy), D, order_space, order_time); + assert (varindex >= 0); + const int groupindex = CCTK_GroupIndexFromVarI (varindex); + assert (groupindex >= 0); + char* groupname = CCTK_GroupName(groupindex); + CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, + "There is no interpolator available for the group \"%s\" with variable type %s, dimension %d, spatial interpolation order %d, temporal interpolation order %d.", + groupname, D, order_space, order_time); + ::free (groupname); } extern "C" { + void CCTK_FCALL CCTK_FNAME(copy_3d_int4) + (const CCTK_INT4* src, + const int& srciext, const int& srcjext, const int& srckext, + CCTK_INT4* dst, + const int& dstiext, const int& dstjext, const int& dstkext, + const int srcbbox[3][3], + const int dstbbox[3][3], + const int regbbox[3][3]); void CCTK_FCALL CCTK_FNAME(copy_3d_real8) (const CCTK_REAL8* src, const int& srciext, const int& srcjext, const int& srckext, @@ -232,6 +354,65 @@ extern "C" { } template<> +void data<CCTK_INT4,3> +::copy_from_innerloop (const gdata<3>* gsrc, const ibbox& box) +{ + const data* src = (const data*)gsrc; + assert (has_storage() && src->has_storage()); + assert (all(box.lower()>=extent().lower() + && box.lower()>=src->extent().lower())); + assert (all(box.upper()<=extent().upper() + && box.upper()<=src->extent().upper())); + assert (all(box.stride()==extent().stride() + && box.stride()==src->extent().stride())); + assert (all((box.lower()-extent().lower())%box.stride() == 0 + && (box.lower()-src->extent().lower())%box.stride() == 0)); + + assert (proc() == src->proc()); + + int rank; + MPI_Comm_rank (dist::comm, &rank); + assert (rank == proc()); + + const ibbox& sext = src->extent(); + const ibbox& dext = extent(); + + int srcshp[3], dstshp[3]; + int srcbbox[3][3], dstbbox[3][3], regbbox[3][3]; + + for (int d=0; d<3; ++d) { + srcshp[d] = (sext.shape() / sext.stride())[d]; + dstshp[d] = (dext.shape() / dext.stride())[d]; + + srcbbox[0][d] = sext.lower()[d]; + srcbbox[1][d] = sext.upper()[d]; + srcbbox[2][d] = sext.stride()[d]; + + dstbbox[0][d] = dext.lower()[d]; + dstbbox[1][d] = dext.upper()[d]; + dstbbox[2][d] = dext.stride()[d]; + + regbbox[0][d] = box.lower()[d]; + regbbox[1][d] = box.upper()[d]; + regbbox[2][d] = box.stride()[d]; + } + + assert (all(dext.stride() == box.stride())); + if (all(sext.stride() == dext.stride())) { + CCTK_FNAME(copy_3d_int4) ((const CCTK_INT4*)src->storage(), + srcshp[0], srcshp[1], srcshp[2], + (CCTK_INT4*)storage(), + dstshp[0], dstshp[1], dstshp[2], + srcbbox, + dstbbox, + regbbox); + + } else { + assert (0); + } +} + +template<> void data<CCTK_REAL8,3> ::copy_from_innerloop (const gdata<3>* gsrc, const ibbox& box) { @@ -302,6 +483,14 @@ extern "C" { const int srcbbox[3][3], const int dstbbox[3][3], const int regbbox[3][3]); + void CCTK_FCALL CCTK_FNAME(restrict_3d_real8_rf2) + (const CCTK_REAL8* src, + const int& srciext, const int& srcjext, const int& srckext, + CCTK_REAL8* dst, + const int& dstiext, const int& dstjext, const int& dstkext, + const int srcbbox[3][3], + const int dstbbox[3][3], + const int regbbox[3][3]); @@ -313,6 +502,14 @@ extern "C" { const int srcbbox[3][3], const int dstbbox[3][3], const int regbbox[3][3]); + void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_rf2) + (const CCTK_REAL8* src, + const int& srciext, const int& srcjext, const int& srckext, + CCTK_REAL8* dst, + const int& dstiext, const int& dstjext, const int& dstkext, + const int srcbbox[3][3], + const int dstbbox[3][3], + const int regbbox[3][3]); void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_o3) (const CCTK_REAL8* src, const int& srciext, const int& srcjext, const int& srckext, @@ -321,6 +518,14 @@ extern "C" { const int srcbbox[3][3], const int dstbbox[3][3], const int regbbox[3][3]); + void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_o3_rf2) + (const CCTK_REAL8* src, + const int& srciext, const int& srcjext, const int& srckext, + CCTK_REAL8* dst, + const int& dstiext, const int& dstjext, const int& dstkext, + const int srcbbox[3][3], + const int dstbbox[3][3], + const int regbbox[3][3]); void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_minmod) (const CCTK_REAL8* src, const int& srciext, const int& srcjext, const int& srckext, @@ -385,6 +590,16 @@ extern "C" { const int srcbbox[3][3], const int dstbbox[3][3], const int regbbox[3][3]); + void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_3tl_rf2) + (const CCTK_REAL8* src1, const CCTK_REAL8& t1, + const CCTK_REAL8* src2, const CCTK_REAL8& t2, + const CCTK_REAL8* src3, const CCTK_REAL8& t3, + const int& srciext, const int& srcjext, const int& srckext, + CCTK_REAL8* dst, const CCTK_REAL8& t, + const int& dstiext, const int& dstjext, const int& dstkext, + const int srcbbox[3][3], + const int dstbbox[3][3], + const int regbbox[3][3]); void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_3tl_o3) (const CCTK_REAL8* src1, const CCTK_REAL8& t1, const CCTK_REAL8* src2, const CCTK_REAL8& t2, @@ -395,6 +610,16 @@ extern "C" { const int srcbbox[3][3], const int dstbbox[3][3], const int regbbox[3][3]); + void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_3tl_o3_rf2) + (const CCTK_REAL8* src1, const CCTK_REAL8& t1, + const CCTK_REAL8* src2, const CCTK_REAL8& t2, + const CCTK_REAL8* src3, const CCTK_REAL8& t3, + const int& srciext, const int& srcjext, const int& srckext, + CCTK_REAL8* dst, const CCTK_REAL8& t, + const int& dstiext, const int& dstjext, const int& dstkext, + const int srcbbox[3][3], + const int dstbbox[3][3], + const int regbbox[3][3]); void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_3tl_minmod) (const CCTK_REAL8* src1, const CCTK_REAL8& t1, const CCTK_REAL8* src2, const CCTK_REAL8& t2, @@ -415,7 +640,6 @@ extern "C" { const int srcbbox[3][3], const int dstbbox[3][3], const int regbbox[3][3]); - } template<> @@ -468,230 +692,151 @@ void data<CCTK_REAL8,3> regbbox[1][d] = box.upper()[d]; regbbox[2][d] = box.stride()[d]; } - - const int varindex = (gsrcs[0])->var_index(); - const int groupindex = CCTK_GroupIndexFromVarI(varindex); - const int group_tags_table = CCTK_GroupTagsTableI(groupindex); - assert(group_tags_table >= 0); - - int typecode = -1; - int keysize = -1; - if (! Util_TableQueryValueInfo(group_tags_table, - &typecode, - &keysize, - "Prolongation")) - { - CCTK_VWarn(4, __LINE__, __FILE__, CCTK_THORNSTRING, - "Tags table for group %s does not contain 'Prolongation'" - " tag", CCTK_GroupName(groupindex)); - } - else - { - assert(typecode == CCTK_VARIABLE_CHAR); - } - -#define PROLONG_NONE 0 -#define PROLONG_LAGRANGE 1 -#define PROLONG_TVD 2 - - int prolong_method; - if (keysize < 0) { // No key in table - default to Lagrange. - prolong_method = 1; - } - else { - char prolong_string[keysize+10]; - const int error = Util_TableGetString(group_tags_table, - keysize+10, - prolong_string, - "Prolongation"); - if (error < 0) { - CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, - "Error code %d getting prolongation method" - " from tags table for group %s.", - error, CCTK_GroupName(groupindex)); - } - if (CCTK_Equals(prolong_string, "None")){ - prolong_method = PROLONG_NONE; - } - else if (CCTK_Equals(prolong_string, "Lagrange")){ - prolong_method = PROLONG_LAGRANGE; - } - else if (CCTK_Equals(prolong_string, "TVD")){ - prolong_method = PROLONG_TVD; - } - else { - CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING, - "Prolongation method %s for group %s not recognized.\n" - "Is this is a typo? Known values are:\n" - "None\nLagrange\nTVD", - prolong_method, CCTK_GroupName(groupindex)); - } - - } - -// CCTK_VInfo(CCTK_THORNSTRING, -// "Variable %d is %s with method %d", -// varindex, CCTK_VarName(varindex), prolong_method); assert (all(dext.stride() == box.stride())); if (all(sext.stride() < dext.stride())) { - assert (srcs.size() == 1); - CCTK_FNAME(restrict_3d_real8) - ((const CCTK_REAL8*)srcs[0]->storage(), - srcshp[0], srcshp[1], srcshp[2], - (CCTK_REAL8*)storage(), - dstshp[0], dstshp[1], dstshp[2], - srcbbox, dstbbox, regbbox); + switch (transport_operator) { + + case op_Lagrange: + case op_TVD: + assert (srcs.size() == 1); + if (all (dext.stride() == sext.stride() * 2)) { + CCTK_FNAME(restrict_3d_real8_rf2) + ((const CCTK_REAL8*)srcs[0]->storage(), + srcshp[0], srcshp[1], srcshp[2], + (CCTK_REAL8*)storage(), + dstshp[0], dstshp[1], dstshp[2], + srcbbox, dstbbox, regbbox); + } else { + CCTK_FNAME(restrict_3d_real8) + ((const CCTK_REAL8*)srcs[0]->storage(), + srcshp[0], srcshp[1], srcshp[2], + (CCTK_REAL8*)storage(), + dstshp[0], dstshp[1], dstshp[2], + srcbbox, dstbbox, regbbox); + } + break; + + default: + assert (0); + + } // switch (prolong_method) } else if (all(sext.stride() > dext.stride())) { - - switch (prolong_method) { - - case PROLONG_NONE: // PROLONG_NONE: NOTHING - break; - case PROLONG_LAGRANGE: // PROLONG_LAGRANGE: DEFAULT - switch (order_time) { - - case 0: - assert (srcs.size()>=1); - switch (order_space) { - case 0: - case 1: - CCTK_FNAME(prolongate_3d_real8) - ((const CCTK_REAL8*)srcs[0]->storage(), - srcshp[0], srcshp[1], srcshp[2], - (CCTK_REAL8*)storage(), - dstshp[0], dstshp[1], dstshp[2], - srcbbox, dstbbox, regbbox); - break; - case 2: - case 3: - CCTK_FNAME(prolongate_3d_real8_o3) - ((const CCTK_REAL8*)srcs[0]->storage(), - srcshp[0], srcshp[1], srcshp[2], - (CCTK_REAL8*)storage(), - dstshp[0], dstshp[1], dstshp[2], - srcbbox, dstbbox, regbbox); - break; - case 4: - case 5: - CCTK_FNAME(prolongate_3d_real8_o5) - ((const CCTK_REAL8*)srcs[0]->storage(), - srcshp[0], srcshp[1], srcshp[2], - (CCTK_REAL8*)storage(), - dstshp[0], dstshp[1], dstshp[2], - srcbbox, dstbbox, regbbox); - break; - default: - assert (0); - } - break; - - case 1: - assert (srcs.size()>=2); - switch (order_space) { - case 0: - case 1: - CCTK_FNAME(prolongate_3d_real8_2tl) - ((const CCTK_REAL8*)srcs[0]->storage(), times[0], - (const CCTK_REAL8*)srcs[1]->storage(), times[1], - srcshp[0], srcshp[1], srcshp[2], - (CCTK_REAL8*)storage(), time, - dstshp[0], dstshp[1], dstshp[2], - srcbbox, dstbbox, regbbox); - break; - case 2: - case 3: - CCTK_FNAME(prolongate_3d_real8_2tl_o3) - ((const CCTK_REAL8*)srcs[0]->storage(), times[0], - (const CCTK_REAL8*)srcs[1]->storage(), times[1], - srcshp[0], srcshp[1], srcshp[2], - (CCTK_REAL8*)storage(), time, - dstshp[0], dstshp[1], dstshp[2], - srcbbox, dstbbox, regbbox); - break; - case 4: - case 5: - CCTK_FNAME(prolongate_3d_real8_2tl_o5) - ((const CCTK_REAL8*)srcs[0]->storage(), times[0], - (const CCTK_REAL8*)srcs[1]->storage(), times[1], - srcshp[0], srcshp[1], srcshp[2], - (CCTK_REAL8*)storage(), time, - dstshp[0], dstshp[1], dstshp[2], - srcbbox, dstbbox, regbbox); - break; - default: - assert (0); - } - break; - - case 2: - assert (srcs.size()>=3); - switch (order_space) { - case 0: - case 1: - CCTK_FNAME(prolongate_3d_real8_3tl) - ((const CCTK_REAL8*)srcs[0]->storage(), times[0], - (const CCTK_REAL8*)srcs[1]->storage(), times[1], - (const CCTK_REAL8*)srcs[2]->storage(), times[2], - srcshp[0], srcshp[1], srcshp[2], - (CCTK_REAL8*)storage(), time, - dstshp[0], dstshp[1], dstshp[2], - srcbbox, dstbbox, regbbox); - break; - case 2: - case 3: - CCTK_FNAME(prolongate_3d_real8_3tl_o3) - ((const CCTK_REAL8*)srcs[0]->storage(), times[0], - (const CCTK_REAL8*)srcs[1]->storage(), times[1], - (const CCTK_REAL8*)srcs[2]->storage(), times[2], - srcshp[0], srcshp[1], srcshp[2], - (CCTK_REAL8*)storage(), time, - dstshp[0], dstshp[1], dstshp[2], - srcbbox, dstbbox, regbbox); - break; - case 4: - case 5: - CCTK_FNAME(prolongate_3d_real8_3tl_o5) - ((const CCTK_REAL8*)srcs[0]->storage(), times[0], - (const CCTK_REAL8*)srcs[1]->storage(), times[1], - (const CCTK_REAL8*)srcs[2]->storage(), times[2], - srcshp[0], srcshp[1], srcshp[2], - (CCTK_REAL8*)storage(), time, - dstshp[0], dstshp[1], dstshp[2], - srcbbox, dstbbox, regbbox); - break; - default: - assert (0); - } - break; - - default: - assert (0); - } // switch (order_time) - break; - case PROLONG_TVD: // PROLONG_TVD - switch (order_time) { - case 0: - CCTK_FNAME(prolongate_3d_real8_minmod) + + switch (transport_operator) { + + case op_Lagrange: + switch (order_time) { + + case 0: + assert (srcs.size()>=1); + switch (order_space) { + case 0: + case 1: + if (all (sext.stride() == dext.stride() * 2)) { + CCTK_FNAME(prolongate_3d_real8_rf2) + ((const CCTK_REAL8*)srcs[0]->storage(), + srcshp[0], srcshp[1], srcshp[2], + (CCTK_REAL8*)storage(), + dstshp[0], dstshp[1], dstshp[2], + srcbbox, dstbbox, regbbox); + } else { + CCTK_FNAME(prolongate_3d_real8) + ((const CCTK_REAL8*)srcs[0]->storage(), + srcshp[0], srcshp[1], srcshp[2], + (CCTK_REAL8*)storage(), + dstshp[0], dstshp[1], dstshp[2], + srcbbox, dstbbox, regbbox); + } + break; + case 2: + case 3: + if (all (sext.stride() == dext.stride() * 2)) { + CCTK_FNAME(prolongate_3d_real8_o3_rf2) ((const CCTK_REAL8*)srcs[0]->storage(), srcshp[0], srcshp[1], srcshp[2], (CCTK_REAL8*)storage(), dstshp[0], dstshp[1], dstshp[2], srcbbox, dstbbox, regbbox); - break; - case 1: - CCTK_FNAME(prolongate_3d_real8_2tl_minmod) + } else { + CCTK_FNAME(prolongate_3d_real8_o3) + ((const CCTK_REAL8*)srcs[0]->storage(), + srcshp[0], srcshp[1], srcshp[2], + (CCTK_REAL8*)storage(), + dstshp[0], dstshp[1], dstshp[2], + srcbbox, dstbbox, regbbox); + } + break; + case 4: + case 5: + CCTK_FNAME(prolongate_3d_real8_o5) + ((const CCTK_REAL8*)srcs[0]->storage(), + srcshp[0], srcshp[1], srcshp[2], + (CCTK_REAL8*)storage(), + dstshp[0], dstshp[1], dstshp[2], + srcbbox, dstbbox, regbbox); + break; + default: + assert (0); + } + break; + + case 1: + assert (srcs.size()>=2); + switch (order_space) { + case 0: + case 1: + CCTK_FNAME(prolongate_3d_real8_2tl) + ((const CCTK_REAL8*)srcs[0]->storage(), times[0], + (const CCTK_REAL8*)srcs[1]->storage(), times[1], + srcshp[0], srcshp[1], srcshp[2], + (CCTK_REAL8*)storage(), time, + dstshp[0], dstshp[1], dstshp[2], + srcbbox, dstbbox, regbbox); + break; + case 2: + case 3: + CCTK_FNAME(prolongate_3d_real8_2tl_o3) + ((const CCTK_REAL8*)srcs[0]->storage(), times[0], + (const CCTK_REAL8*)srcs[1]->storage(), times[1], + srcshp[0], srcshp[1], srcshp[2], + (CCTK_REAL8*)storage(), time, + dstshp[0], dstshp[1], dstshp[2], + srcbbox, dstbbox, regbbox); + break; + case 4: + case 5: + CCTK_FNAME(prolongate_3d_real8_2tl_o5) + ((const CCTK_REAL8*)srcs[0]->storage(), times[0], + (const CCTK_REAL8*)srcs[1]->storage(), times[1], + srcshp[0], srcshp[1], srcshp[2], + (CCTK_REAL8*)storage(), time, + dstshp[0], dstshp[1], dstshp[2], + srcbbox, dstbbox, regbbox); + break; + default: + assert (0); + } + break; + + case 2: + assert (srcs.size()>=3); + switch (order_space) { + case 0: + case 1: + if (all (sext.stride() == dext.stride() * 2)) { + CCTK_FNAME(prolongate_3d_real8_3tl_rf2) ((const CCTK_REAL8*)srcs[0]->storage(), times[0], (const CCTK_REAL8*)srcs[1]->storage(), times[1], + (const CCTK_REAL8*)srcs[2]->storage(), times[2], srcshp[0], srcshp[1], srcshp[2], (CCTK_REAL8*)storage(), time, dstshp[0], dstshp[1], dstshp[2], srcbbox, dstbbox, regbbox); - break; - case 2: - CCTK_FNAME(prolongate_3d_real8_3tl_minmod) + } else { + CCTK_FNAME(prolongate_3d_real8_3tl) ((const CCTK_REAL8*)srcs[0]->storage(), times[0], (const CCTK_REAL8*)srcs[1]->storage(), times[1], (const CCTK_REAL8*)srcs[2]->storage(), times[2], @@ -699,16 +844,88 @@ void data<CCTK_REAL8,3> (CCTK_REAL8*)storage(), time, dstshp[0], dstshp[1], dstshp[2], srcbbox, dstbbox, regbbox); - break; - default: - assert (0); + } + break; + case 2: + case 3: + if (all (sext.stride() == dext.stride() * 2)) { + CCTK_FNAME(prolongate_3d_real8_3tl_o3_rf2) + ((const CCTK_REAL8*)srcs[0]->storage(), times[0], + (const CCTK_REAL8*)srcs[1]->storage(), times[1], + (const CCTK_REAL8*)srcs[2]->storage(), times[2], + srcshp[0], srcshp[1], srcshp[2], + (CCTK_REAL8*)storage(), time, + dstshp[0], dstshp[1], dstshp[2], + srcbbox, dstbbox, regbbox); + } else { + CCTK_FNAME(prolongate_3d_real8_3tl_o3) + ((const CCTK_REAL8*)srcs[0]->storage(), times[0], + (const CCTK_REAL8*)srcs[1]->storage(), times[1], + (const CCTK_REAL8*)srcs[2]->storage(), times[2], + srcshp[0], srcshp[1], srcshp[2], + (CCTK_REAL8*)storage(), time, + dstshp[0], dstshp[1], dstshp[2], + srcbbox, dstbbox, regbbox); + } + break; + case 4: + case 5: + CCTK_FNAME(prolongate_3d_real8_3tl_o5) + ((const CCTK_REAL8*)srcs[0]->storage(), times[0], + (const CCTK_REAL8*)srcs[1]->storage(), times[1], + (const CCTK_REAL8*)srcs[2]->storage(), times[2], + srcshp[0], srcshp[1], srcshp[2], + (CCTK_REAL8*)storage(), time, + dstshp[0], dstshp[1], dstshp[2], + srcbbox, dstbbox, regbbox); + break; + default: + assert (0); } break; + default: - assert(0); + assert (0); + } // switch (order_time) + break; + + case op_TVD: + switch (order_time) { + case 0: + CCTK_FNAME(prolongate_3d_real8_minmod) + ((const CCTK_REAL8*)srcs[0]->storage(), + srcshp[0], srcshp[1], srcshp[2], + (CCTK_REAL8*)storage(), + dstshp[0], dstshp[1], dstshp[2], + srcbbox, dstbbox, regbbox); + break; + case 1: + CCTK_FNAME(prolongate_3d_real8_2tl_minmod) + ((const CCTK_REAL8*)srcs[0]->storage(), times[0], + (const CCTK_REAL8*)srcs[1]->storage(), times[1], + srcshp[0], srcshp[1], srcshp[2], + (CCTK_REAL8*)storage(), time, + dstshp[0], dstshp[1], dstshp[2], + srcbbox, dstbbox, regbbox); + case 2: + CCTK_FNAME(prolongate_3d_real8_3tl_minmod) + ((const CCTK_REAL8*)srcs[0]->storage(), times[0], + (const CCTK_REAL8*)srcs[1]->storage(), times[1], + (const CCTK_REAL8*)srcs[2]->storage(), times[2], + srcshp[0], srcshp[1], srcshp[2], + (CCTK_REAL8*)storage(), time, + dstshp[0], dstshp[1], dstshp[2], + srcbbox, dstbbox, regbbox); + break; + default: + assert (0); + } + break; + + default: + assert(0); } // switch (prolong_method) - } else { assert (0); } diff --git a/Carpet/CarpetLib/src/data.hh b/Carpet/CarpetLib/src/data.hh index f126c6c47..bb8e5c497 100644 --- a/Carpet/CarpetLib/src/data.hh +++ b/Carpet/CarpetLib/src/data.hh @@ -1,4 +1,4 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.hh,v 1.14 2003/10/14 16:39:16 schnetter Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.hh,v 1.15 2003/11/05 16:18:39 schnetter Exp $ #ifndef DATA_HH #define DATA_HH @@ -30,7 +30,12 @@ class data: public gdata<D> { // Fields T* _storage; // the data (if located on this processor) - + + bool comm_active; + MPI_Request request; + + int tag; // MPI tag for this object + public: // Constructors @@ -50,7 +55,13 @@ public: virtual void transfer_from (gdata<D>* gsrc); // Processor management - virtual void change_processor (const int newproc, void* const mem=0); + virtual void change_processor (comm_state<D>& state, + const int newproc, void* const mem=0); + private: + virtual void change_processor_recv (const int newproc, void* const mem=0); + virtual void change_processor_send (const int newproc, void* const mem=0); + virtual void change_processor_wait (const int newproc, void* const mem=0); + public: // Accessors virtual const void* storage () const { diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc index 83d43932b..4cdc1ae28 100644 --- a/Carpet/CarpetLib/src/dh.cc +++ b/Carpet/CarpetLib/src/dh.cc @@ -1,4 +1,4 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.cc,v 1.43 2003/10/16 17:00:04 schnetter Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.cc,v 1.44 2003/11/05 16:18:39 schnetter Exp $ #include <assert.h> @@ -31,7 +31,7 @@ dh<D>::dh (gh<D>& h, assert (buffer_width>=0); h.add(this); CHECKPOINT; - recompose(false); + recompose (0, true); } // Destructors @@ -51,7 +51,7 @@ int dh<D>::prolongation_stencil_size () const { // Modifiers template<int D> -void dh<D>::recompose (const int initialise_upto) { +void dh<D>::recompose (const int initialise_from, const bool do_prolongate) { DECLARE_CCTK_PARAMETERS; CHECKPOINT; @@ -562,7 +562,7 @@ void dh<D>::recompose (const int initialise_upto) { for (typename list<ggf<D>*>::iterator f=gfs.begin(); f!=gfs.end(); ++f) { - (*f)->recompose(initialise_upto); + (*f)->recompose (initialise_from, do_prolongate); } } diff --git a/Carpet/CarpetLib/src/dh.hh b/Carpet/CarpetLib/src/dh.hh index 4a23d5cba..0f9cbaaa5 100644 --- a/Carpet/CarpetLib/src/dh.hh +++ b/Carpet/CarpetLib/src/dh.hh @@ -1,4 +1,4 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.hh,v 1.15 2003/07/20 21:03:43 schnetter Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.hh,v 1.16 2003/11/05 16:18:39 schnetter Exp $ #ifndef DH_HH #define DH_HH @@ -115,7 +115,7 @@ public: int prolongation_stencil_size () const; // Modifiers - void recompose (const int initialise_upto); + void recompose (const int initialise_from, const bool do_prolongate); // Grid function management void add (ggf<D>* f); diff --git a/Carpet/CarpetLib/src/dist.hh b/Carpet/CarpetLib/src/dist.hh index 10d4a3cfc..cb3d422ca 100644 --- a/Carpet/CarpetLib/src/dist.hh +++ b/Carpet/CarpetLib/src/dist.hh @@ -1,4 +1,4 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dist.hh,v 1.8 2003/01/03 15:49:36 schnetter Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dist.hh,v 1.9 2003/11/05 16:18:39 schnetter Exp $ #ifndef DIST_HH #define DIST_HH @@ -19,8 +19,6 @@ using namespace std; namespace dist { - const int tag = 1; - extern MPI_Comm comm; extern MPI_Datatype mpi_complex_float; diff --git a/Carpet/CarpetLib/src/gdata.cc b/Carpet/CarpetLib/src/gdata.cc index e438a3147..4236b4b41 100644 --- a/Carpet/CarpetLib/src/gdata.cc +++ b/Carpet/CarpetLib/src/gdata.cc @@ -1,11 +1,15 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.cc,v 1.23 2003/10/14 16:39:16 schnetter Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.cc,v 1.24 2003/11/05 16:18:39 schnetter Exp $ #include <assert.h> +#include <stdlib.h> #include <iostream> #include "cctk.h" +#include "util_ErrorCodes.h" +#include "util_Table.h" + #include "bbox.hh" #include "defs.hh" #include "dist.hh" @@ -17,11 +21,43 @@ using namespace std; +// Communication state control +template<int D> +comm_state<D>::comm_state () + : thestate(state_recv), + current(0) +{ +} + +template<int D> +void comm_state<D>::step () +{ + assert (thestate!=state_done); + assert (current==tmps.size()); + thestate=astate(size_t(thestate)+1); + current=0; +} + +template<int D> +bool comm_state<D>::done () +{ + return thestate==state_done; +} + +template<int D> +comm_state<D>::~comm_state () +{ + assert (thestate==state_recv || thestate==state_done); +} + + + // Constructors template<int D> gdata<D>::gdata (const int varindex_) : varindex(varindex_), - _has_storage(false) + _has_storage(false), + transport_operator (find_transport_operator(varindex_)) { } // Destructors @@ -30,9 +66,93 @@ gdata<D>::~gdata () { } +// Transport operator types +template<int D> +gdata<D>::operator_type gdata<D>::find_transport_operator (const int varindex) { + const operator_type default_operator = op_Lagrange; + if (varindex == -1) return op_error; + assert (varindex >= 0); + const int groupindex = CCTK_GroupIndexFromVarI (varindex); + assert (groupindex >= 0); + const int group_tags_table = CCTK_GroupTagsTableI(groupindex); + assert (group_tags_table >= 0); + char prolong_string[100]; + const int ierr = Util_TableGetString + (group_tags_table, sizeof prolong_string, prolong_string, "Prolongation"); + if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY) { + char* groupname = CCTK_GroupName(groupindex); + CCTK_VWarn (4, __LINE__, __FILE__, CCTK_THORNSTRING, + "Tags table for group \"%s\" does not contain \"Prolongation\" tag", groupname); + ::free (groupname); + return default_operator; + } + assert (ierr >= 0); + if (CCTK_Equals(prolong_string, "None")) { + return op_none; + } else if (CCTK_Equals(prolong_string, "Lagrange")) { + return op_Lagrange; + } else if (CCTK_Equals(prolong_string, "TVD")) { + return op_TVD; + } else { + assert (0); + } + return op_error; +} + + + // Data manipulators template<int D> -void gdata<D>::copy_from (const gdata* src, const ibbox& box) +void gdata<D>::copy_from (comm_state<D>& state, + const gdata* src, const ibbox& box) +{ + switch (state.thestate) { + case state_recv: + copy_from_recv (state, src, box); + break; + case state_send: + copy_from_send (state, src, box); + break; + case state_wait: + copy_from_wait (state, src, box); + break; + default: + assert(0); + } +} + + + +template<int D> +void gdata<D>::copy_from_nocomm (const gdata* src, const ibbox& box) +{ + assert (has_storage() && src->has_storage()); + assert (all(box.lower()>=extent().lower() + && box.lower()>=src->extent().lower())); + assert (all(box.upper()<=extent().upper() + && box.upper()<=src->extent().upper())); + assert (all(box.stride()==extent().stride() + && box.stride()==src->extent().stride())); + assert (all((box.lower()-extent().lower())%box.stride() == 0 + && (box.lower()-src->extent().lower())%box.stride() == 0)); + + if (box.empty()) return; + + assert (proc() == src->proc()); + + // copy on same processor + int rank; + MPI_Comm_rank (dist::comm, &rank); + if (rank == proc()) { + copy_from_innerloop (src, box); + } +} + + + +template<int D> +void gdata<D>::copy_from_recv (comm_state<D>& state, + const gdata* src, const ibbox& box) { assert (has_storage() && src->has_storage()); assert (all(box.lower()>=extent().lower() @@ -49,20 +169,81 @@ void gdata<D>::copy_from (const gdata* src, const ibbox& box) if (proc() == src->proc()) { // copy on same processor - int rank; - MPI_Comm_rank (dist::comm, &rank); - if (rank == proc()) { - copy_from_innerloop (src, box); - } - } else { // copy to different processor - gdata* const tmp = make_typed(varindex); + gdata<D>* const tmp = make_typed(varindex); + // TODO: is this efficient? + state.tmps.push_back (tmp); + ++state.current; tmp->allocate (box, src->proc()); - tmp->copy_from (src, box); - tmp->change_processor (proc()); - copy_from (tmp, box); + tmp->change_processor_recv (proc()); + + } +} + + + +template<int D> +void gdata<D>::copy_from_send (comm_state<D>& state, + const gdata* src, const ibbox& box) +{ + assert (has_storage() && src->has_storage()); + assert (all(box.lower()>=extent().lower() + && box.lower()>=src->extent().lower())); + assert (all(box.upper()<=extent().upper() + && box.upper()<=src->extent().upper())); + assert (all(box.stride()==extent().stride() + && box.stride()==src->extent().stride())); + assert (all((box.lower()-extent().lower())%box.stride() == 0 + && (box.lower()-src->extent().lower())%box.stride() == 0)); + + if (box.empty()) return; + + if (proc() == src->proc()) { + // copy on same processor + + copy_from_nocomm (src, box); + + } else { + + // copy to different processor + gdata<D>* const tmp = state.tmps.at(state.current++); + assert (tmp); + tmp->copy_from_nocomm (src, box); + tmp->change_processor_send (proc()); + + } +} + + + +template<int D> +void gdata<D>::copy_from_wait (comm_state<D>& state, + const gdata* src, const ibbox& box) +{ + assert (has_storage() && src->has_storage()); + assert (all(box.lower()>=extent().lower() + && box.lower()>=src->extent().lower())); + assert (all(box.upper()<=extent().upper() + && box.upper()<=src->extent().upper())); + assert (all(box.stride()==extent().stride() + && box.stride()==src->extent().stride())); + assert (all((box.lower()-extent().lower())%box.stride() == 0 + && (box.lower()-src->extent().lower())%box.stride() == 0)); + + if (box.empty()) return; + + if (proc() == src->proc()) { + // copy on same processor + + } else { + + // copy to different processor + gdata<D>* const tmp = state.tmps.at(state.current++); + assert (tmp); + tmp->change_processor_wait (proc()); + copy_from_nocomm (tmp, box); delete tmp; } @@ -72,11 +253,79 @@ void gdata<D>::copy_from (const gdata* src, const ibbox& box) template<int D> void gdata<D> -::interpolate_from (const vector<const gdata*> srcs, - const vector<CCTK_REAL> times, - const ibbox& box, const CCTK_REAL time, - const int order_space, - const int order_time) +::interpolate_from (comm_state<D>& state, + const vector<const gdata*> srcs, + const vector<CCTK_REAL> times, + const ibbox& box, const CCTK_REAL time, + const int order_space, + const int order_time) +{ + assert (transport_operator != op_error); + if (transport_operator == op_none) return; + switch (state.thestate) { + case state_recv: + interpolate_from_recv (state, srcs, times, box, time, order_space, order_time); + break; + case state_send: + interpolate_from_send (state, srcs, times, box, time, order_space, order_time); + break; + case state_wait: + interpolate_from_wait (state, srcs, times, box, time, order_space, order_time); + break; + default: + assert(0); + } +} + + + +template<int D> +void gdata<D> +::interpolate_from_nocomm (const vector<const gdata*> srcs, + const vector<CCTK_REAL> times, + const ibbox& box, const CCTK_REAL time, + const int order_space, + const int order_time) +{ + assert (has_storage()); + assert (all(box.lower()>=extent().lower())); + assert (all(box.upper()<=extent().upper())); + assert (all(box.stride()==extent().stride())); + assert (all((box.lower()-extent().lower())%box.stride() == 0)); + assert (srcs.size() == times.size() && srcs.size()>0); + for (int t=0; t<(int)srcs.size(); ++t) { + assert (srcs[t]->has_storage()); + assert (all(box.lower()>=srcs[t]->extent().lower())); + assert (all(box.upper()<=srcs[t]->extent().upper())); + } + + assert (! box.empty()); + if (box.empty()) return; + + assert (proc() == srcs[0]->proc()); + + assert (transport_operator != op_error); + assert (transport_operator != op_none); + + // interpolate on same processor + int rank; + MPI_Comm_rank (dist::comm, &rank); + if (rank == proc()) { + interpolate_from_innerloop + (srcs, times, box, time, order_space, order_time); + } +} + + + +template<int D> +void gdata<D> +::interpolate_from_recv (comm_state<D>& state, + const vector<const gdata*> srcs, + const vector<CCTK_REAL> times, + const ibbox& box, const CCTK_REAL time, + const int order_space, + const int order_time) { assert (has_storage()); assert (all(box.lower()>=extent().lower())); @@ -96,21 +345,97 @@ void gdata<D> if (proc() == srcs[0]->proc()) { // interpolate on same processor - int rank; - MPI_Comm_rank (dist::comm, &rank); - if (rank == proc()) { - interpolate_from_innerloop - (srcs, times, box, time, order_space, order_time); - } - } else { // interpolate from other processor - gdata* const tmp = make_typed(varindex); + gdata<D>* const tmp = make_typed(varindex); + // TODO: is this efficient? + state.tmps.push_back (tmp); + ++state.current; tmp->allocate (box, srcs[0]->proc()); - tmp->interpolate_from (srcs, times, box, time, order_space, order_time); - tmp->change_processor (proc()); - copy_from (tmp, box); + tmp->change_processor_recv (proc()); + + } +} + + + +template<int D> +void gdata<D> +::interpolate_from_send (comm_state<D>& state, + const vector<const gdata*> srcs, + const vector<CCTK_REAL> times, + const ibbox& box, const CCTK_REAL time, + const int order_space, + const int order_time) +{ + assert (has_storage()); + assert (all(box.lower()>=extent().lower())); + assert (all(box.upper()<=extent().upper())); + assert (all(box.stride()==extent().stride())); + assert (all((box.lower()-extent().lower())%box.stride() == 0)); + assert (srcs.size() == times.size() && srcs.size()>0); + for (int t=0; t<(int)srcs.size(); ++t) { + assert (srcs[t]->has_storage()); + assert (all(box.lower()>=srcs[t]->extent().lower())); + assert (all(box.upper()<=srcs[t]->extent().upper())); + } + + assert (! box.empty()); + if (box.empty()) return; + + if (proc() == srcs[0]->proc()) { + // interpolate on same processor + + interpolate_from_nocomm (srcs, times, box, time, order_space, order_time); + + } else { + // interpolate from other processor + + gdata<D>* const tmp = state.tmps.at(state.current++); + assert (tmp); + tmp->interpolate_from_nocomm (srcs, times, box, time, order_space, order_time); + tmp->change_processor_send (proc()); + + } +} + + + +template<int D> +void gdata<D> +::interpolate_from_wait (comm_state<D>& state, + const vector<const gdata*> srcs, + const vector<CCTK_REAL> times, + const ibbox& box, const CCTK_REAL time, + const int order_space, + const int order_time) +{ + assert (has_storage()); + assert (all(box.lower()>=extent().lower())); + assert (all(box.upper()<=extent().upper())); + assert (all(box.stride()==extent().stride())); + assert (all((box.lower()-extent().lower())%box.stride() == 0)); + assert (srcs.size() == times.size() && srcs.size()>0); + for (int t=0; t<(int)srcs.size(); ++t) { + assert (srcs[t]->has_storage()); + assert (all(box.lower()>=srcs[t]->extent().lower())); + assert (all(box.upper()<=srcs[t]->extent().upper())); + } + + assert (! box.empty()); + if (box.empty()) return; + + if (proc() == srcs[0]->proc()) { + // interpolate on same processor + + } else { + // interpolate from other processor + + gdata<D>* const tmp = state.tmps.at(state.current++); + assert (tmp); + tmp->change_processor_wait (proc()); + copy_from_nocomm (tmp, box); delete tmp; } @@ -118,4 +443,5 @@ void gdata<D> +template class comm_state<3>; template class gdata<3>; diff --git a/Carpet/CarpetLib/src/gdata.hh b/Carpet/CarpetLib/src/gdata.hh index edc08c8e7..808a8eb51 100644 --- a/Carpet/CarpetLib/src/gdata.hh +++ b/Carpet/CarpetLib/src/gdata.hh @@ -1,4 +1,4 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.hh,v 1.19 2003/10/15 07:14:01 hawke Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.hh,v 1.20 2003/11/05 16:18:39 schnetter Exp $ #ifndef GDATA_HH #define GDATA_HH @@ -8,6 +8,7 @@ #include <iostream> #include <string> +#include <vector> #include "cctk.h" @@ -19,6 +20,29 @@ using namespace std; + + +template<int D> +class gdata; + + + +// State information for communications +enum astate { state_recv, state_send, state_wait, state_done }; + +template<int D> +struct comm_state { + astate thestate; + comm_state (); + void step (); + bool done (); + ~comm_state (); + + vector<gdata<D>*> tmps; + size_t current; +}; + + // A generic data storage without type information template<int D> @@ -57,7 +81,13 @@ public: virtual gdata<D>* make_typed (const int varindex) const = 0; // Processor management - virtual void change_processor (const int newproc, void* const mem=0) = 0; + virtual void change_processor (comm_state<D>& state, + const int newproc, void* const mem=0) = 0; + protected: + virtual void change_processor_recv (const int newproc, void* const mem=0) = 0; + virtual void change_processor_send (const int newproc, void* const mem=0) = 0; + virtual void change_processor_wait (const int newproc, void* const mem=0) = 0; + public: // Storage management virtual void transfer_from (gdata<D>* src) = 0; @@ -68,10 +98,6 @@ public: // Accessors - int var_index () const { - return varindex; - } - bool has_storage () const { return _has_storage; } @@ -117,14 +143,59 @@ public: assert (all(ind>=0 && ind<=shape())); return dot(ind, stride()); } - + + // Transport operator types + protected: + enum operator_type { op_error, op_none, op_Lagrange, op_TVD }; + // readonly + operator_type transport_operator; + private: + static operator_type find_transport_operator (const int varindex_); + // Data manipulators - void copy_from (const gdata* src, const ibbox& box); - void interpolate_from (const vector<const gdata*> srcs, - const vector<CCTK_REAL> times, - const ibbox& box, const CCTK_REAL time, - const int order_space, - const int order_time); + public: + void copy_from (comm_state<D>& state, const gdata* src, const ibbox& box); + private: + void copy_from_nocomm (const gdata* src, const ibbox& box); + void copy_from_recv (comm_state<D>& state, + const gdata* src, const ibbox& box); + void copy_from_send (comm_state<D>& state, + const gdata* src, const ibbox& box); + void copy_from_wait (comm_state<D>& state, + const gdata* src, const ibbox& box); + public: + void interpolate_from (comm_state<D>& state, + const vector<const gdata*> srcs, + const vector<CCTK_REAL> times, + const ibbox& box, const CCTK_REAL time, + const int order_space, + const int order_time); + private: + void interpolate_from_nocomm (const vector<const gdata*> srcs, + const vector<CCTK_REAL> times, + const ibbox& box, const CCTK_REAL time, + const int order_space, + const int order_time); + void interpolate_from_recv (comm_state<D>& state, + const vector<const gdata*> srcs, + const vector<CCTK_REAL> times, + const ibbox& box, const CCTK_REAL time, + const int order_space, + const int order_time); + void interpolate_from_send (comm_state<D>& state, + const vector<const gdata*> srcs, + const vector<CCTK_REAL> times, + const ibbox& box, const CCTK_REAL time, + const int order_space, + const int order_time); + void interpolate_from_wait (comm_state<D>& state, + const vector<const gdata*> srcs, + const vector<CCTK_REAL> times, + const ibbox& box, const CCTK_REAL time, + const int order_space, + const int order_time); + public: + protected: virtual void copy_from_innerloop (const gdata* src, const ibbox& box) = 0; diff --git a/Carpet/CarpetLib/src/gf.cc b/Carpet/CarpetLib/src/gf.cc index 8de652672..4a913d199 100644 --- a/Carpet/CarpetLib/src/gf.cc +++ b/Carpet/CarpetLib/src/gf.cc @@ -1,4 +1,4 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gf.cc,v 1.13 2003/10/14 16:39:16 schnetter Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gf.cc,v 1.14 2003/11/05 16:18:39 schnetter Exp $ #include <assert.h> @@ -20,7 +20,7 @@ gf<T,D>::gf (const int varindex, th<D>& t, dh<D>& d, : ggf<D>(varindex, t, d, tmin, tmax, prolongation_order_time) { // VGF - this->recompose(); + this->recompose (0, true); } // Destructors diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc index a096ce0fb..6786e5eb7 100644 --- a/Carpet/CarpetLib/src/ggf.cc +++ b/Carpet/CarpetLib/src/ggf.cc @@ -1,4 +1,4 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.cc,v 1.28 2003/10/14 16:39:16 schnetter Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.cc,v 1.29 2003/11/05 16:18:39 schnetter Exp $ #include <assert.h> #include <stdlib.h> @@ -51,7 +51,9 @@ bool ggf<D>::operator== (const ggf<D>& f) const { // Modifiers // VGF template<int D> -void ggf<D>::recompose (const int initialise_upto) { +void ggf<D>::recompose (const int initialise_from, const bool do_prolongate) { + + // TODO: restructure storage only when needed // Retain storage that might be needed fdata oldstorage = storage; @@ -60,19 +62,19 @@ void ggf<D>::recompose (const int initialise_upto) { storage.resize(tmax-tmin+1); for (int tl=tmin; tl<=tmax; ++tl) { storage[tl-tmin].resize(h.reflevels()); - for (int rl=0; rl<h.reflevels(); ++rl) { + for (int rl=initialise_from; rl<h.reflevels(); ++rl) { storage[tl-tmin][rl].resize(h.components(rl)); for (int c=0; c<h.components(rl); ++c) { storage[tl-tmin][rl][c].resize(h.mglevels(rl,c)); for (int ml=0; ml<h.mglevels(rl,c); ++ml) { - storage[tl-tmin][rl][c][ml] = 0; + storage[tl-tmin][rl][c][ml] = 0; } // for ml } // for c } // for rl } // for tl // Initialise the new storage - for (int rl=0; rl<h.reflevels(); ++rl) { + for (int rl=initialise_from; rl<h.reflevels(); ++rl) { for (int tl=tmin; tl<=tmax; ++tl) { for (int c=0; c<h.components(rl); ++c) { for (int ml=0; ml<h.mglevels(rl,c); ++ml) { @@ -84,30 +86,32 @@ void ggf<D>::recompose (const int initialise_upto) { storage[tl-tmin][rl][c][ml]->allocate (d.boxes[rl][c][ml].exterior, h.proc(rl,c)); - // Initialise only if desired - if (initialise_upto >= 0 && rl <= initialise_upto) { - + if (do_prolongate) { // Initialise from coarser level, if possible // TODO: init only un-copied regions if (rl>0) { - const CCTK_REAL time = t.time(tl,rl,ml); - ref_prolongate (tl,rl,c,ml,time); + for (comm_state<D> state; !state.done(); state.step()) { + const CCTK_REAL time = t.time(tl,rl,ml); + ref_prolongate (state,tl,rl,c,ml,time); + } } // if rl - - // Copy from old storage, if possible - if (rl<(int)oldstorage[tl-tmin].size()) { + } // if do_prolongate + + // Copy from old storage, if possible + // todo: copy only from interior regions? + if (rl<(int)oldstorage[tl-tmin].size()) { + for (comm_state<D> state; !state.done(); state.step()) { for (int cc=0; cc<(int)oldstorage[tl-tmin][rl].size(); ++cc) { if (ml<(int)oldstorage[tl-tmin][rl][cc].size()) { const ibbox ovlp = (d.boxes[rl][c][ml].exterior & oldstorage[tl-tmin][rl][cc][ml]->extent()); storage[tl-tmin][rl][c][ml]->copy_from - (oldstorage[tl-tmin][rl][cc][ml], ovlp); + (state, oldstorage[tl-tmin][rl][cc][ml], ovlp); } // if ml } // for cc - } // if rl - - } // if initialise + } // for step + } // if rl } // for ml } // for c @@ -117,7 +121,7 @@ void ggf<D>::recompose (const int initialise_upto) { // Delete old storage for (int tl=tmin; tl<=tmax; ++tl) { - for (int rl=0; rl<(int)oldstorage[tl-tmin].size(); ++rl) { + for (int rl=initialise_from; rl<(int)oldstorage[tl-tmin].size(); ++rl) { for (int c=0; c<(int)oldstorage[tl-tmin][rl].size(); ++c) { for (int ml=0; ml<(int)oldstorage[tl-tmin][rl][c].size(); ++ml) { delete oldstorage[tl-tmin][rl][c][ml]; @@ -126,30 +130,29 @@ void ggf<D>::recompose (const int initialise_upto) { } // for rl } // for tl - for (int rl=0; rl<h.reflevels(); ++rl) { - for (int tl=tmin; tl<=tmax; ++tl) { + for (int tl=tmin; tl<=tmax; ++tl) { + for (int rl=initialise_from; rl<h.reflevels(); ++rl) { // Set boundaries for (int c=0; c<h.components(rl); ++c) { for (int ml=0; ml<h.mglevels(rl,c); ++ml) { - // Initialise only if desired - if (initialise_upto >= 0 && rl <= initialise_upto) { - - sync (tl,rl,c,ml); - // TODO: assert that reflevel 0 boundaries are copied - if (rl>0) { + for (comm_state<D> state; !state.done(); state.step()) { + sync (state,tl,rl,c,ml); + } + // TODO: assert that reflevel 0 boundaries are copied + if (rl>0) { + for (comm_state<D> state; !state.done(); state.step()) { const CCTK_REAL time = t.time(tl,rl,ml); - ref_bnd_prolongate (tl,rl,c,ml,time); - } // if rl - - } // if initialise - + ref_bnd_prolongate (state,tl,rl,c,ml,time); + } + } // if rl + } // for ml } // for c - } // for tl - } // for rl + } // for rl + } // for tl } // Cycle the time levels by rotating the data sets @@ -181,6 +184,7 @@ void ggf<D>::flip (int rl, int c, int ml) { } } +#if 0 // Copy data from current time level to all previous levels template<int D> void ggf<D>::copytoprevs (int rl, int c, int ml) { @@ -192,6 +196,7 @@ void ggf<D>::copytoprevs (int rl, int c, int ml) { (storage[tmax-tmin][rl][c][ml], storage[tmax-tmin][rl][c][ml]->extent()); } } +#endif @@ -201,7 +206,8 @@ void ggf<D>::copytoprevs (int rl, int c, int ml) { // Copy a region template<int D> -void ggf<D>::copycat (int tl1, int rl1, int c1, int ml1, +void ggf<D>::copycat (comm_state<D>& state, + int tl1, int rl1, int c1, int ml1, const ibbox dh<D>::dboxes::* recv_list, int tl2, int rl2, int ml2, const ibbox dh<D>::dboxes::* send_list) @@ -220,12 +226,13 @@ void ggf<D>::copycat (int tl1, int rl1, int c1, int ml1, // copy the content assert (recv==send); storage[tl1-tmin][rl1][c1][ml1]->copy_from - (storage[tl2-tmin][rl2][c2][ml2], recv); + (state, storage[tl2-tmin][rl2][c2][ml2], recv); } // Copy regions template<int D> -void ggf<D>::copycat (int tl1, int rl1, int c1, int ml1, +void ggf<D>::copycat (comm_state<D>& state, + int tl1, int rl1, int c1, int ml1, const iblist dh<D>::dboxes::* recv_list, int tl2, int rl2, int ml2, const iblist dh<D>::dboxes::* send_list) @@ -247,13 +254,14 @@ void ggf<D>::copycat (int tl1, int rl1, int c1, int ml1, // (use the send boxes for communication) // copy the content storage[tl1-tmin][rl1][c1][ml1]->copy_from - (storage[tl2-tmin][rl2][c2][ml2], *r); + (state, storage[tl2-tmin][rl2][c2][ml2], *r); } } // Copy regions template<int D> -void ggf<D>::copycat (int tl1, int rl1, int c1, int ml1, +void ggf<D>::copycat (comm_state<D>& state, + int tl1, int rl1, int c1, int ml1, const iblistvect dh<D>::dboxes::* recv_listvect, int tl2, int rl2, int ml2, const iblistvect dh<D>::dboxes::* send_listvect) @@ -276,14 +284,15 @@ void ggf<D>::copycat (int tl1, int rl1, int c1, int ml1, // (use the send boxes for communication) // copy the content storage[tl1-tmin][rl1][c1][ml1]->copy_from - (storage[tl2-tmin][rl2][c2][ml2], *r); + (state, storage[tl2-tmin][rl2][c2][ml2], *r); } } } // Interpolate a region template<int D> -void ggf<D>::intercat (int tl1, int rl1, int c1, int ml1, +void ggf<D>::intercat (comm_state<D>& state, + int tl1, int rl1, int c1, int ml1, const ibbox dh<D>::dboxes::* recv_list, const vector<int> tl2s, int rl2, int ml2, const ibbox dh<D>::dboxes::* send_list, @@ -317,13 +326,14 @@ void ggf<D>::intercat (int tl1, int rl1, int c1, int ml1, // interpolate the content assert (recv==send); storage[tl1-tmin][rl1][c1][ml1]->interpolate_from - (gsrcs, times, recv, time, + (state, gsrcs, times, recv, time, d.prolongation_order_space, prolongation_order_time); } // Interpolate regions template<int D> -void ggf<D>::intercat (int tl1, int rl1, int c1, int ml1, +void ggf<D>::intercat (comm_state<D>& state, + int tl1, int rl1, int c1, int ml1, const iblist dh<D>::dboxes::* recv_list, const vector<int> tl2s, int rl2, int ml2, const iblist dh<D>::dboxes::* send_list, @@ -360,14 +370,15 @@ void ggf<D>::intercat (int tl1, int rl1, int c1, int ml1, // (use the send boxes for communication) // interpolate the content storage[tl1-tmin][rl1][c1][ml1]->interpolate_from - (gsrcs, times, *r, time, + (state, gsrcs, times, *r, time, d.prolongation_order_space, prolongation_order_time); } } // Interpolate regions template<int D> -void ggf<D>::intercat (int tl1, int rl1, int c1, int ml1, +void ggf<D>::intercat (comm_state<D>& state, + int tl1, int rl1, int c1, int ml1, const iblistvect dh<D>::dboxes::* recv_listvect, const vector<int> tl2s, int rl2, int ml2, const iblistvect dh<D>::dboxes::* send_listvect, @@ -405,7 +416,7 @@ void ggf<D>::intercat (int tl1, int rl1, int c1, int ml1, // (use the send boxes for communication) // interpolate the content storage[tl1-tmin][rl1][c1][ml1]->interpolate_from - (gsrcs, times, *r, time, + (state, gsrcs, times, *r, time, d.prolongation_order_space, prolongation_order_time); } } @@ -415,25 +426,28 @@ void ggf<D>::intercat (int tl1, int rl1, int c1, int ml1, // Copy a component from the next time level template<int D> -void ggf<D>::copy (int tl, int rl, int c, int ml) +void ggf<D>::copy (comm_state<D>& state, int tl, int rl, int c, int ml) { // Copy - copycat (tl ,rl,c,ml, &dh<D>::dboxes::exterior, + copycat (state, + tl ,rl,c,ml, &dh<D>::dboxes::exterior, tl+1,rl, ml, &dh<D>::dboxes::exterior); } // Synchronise the boundaries a component template<int D> -void ggf<D>::sync (int tl, int rl, int c, int ml) +void ggf<D>::sync (comm_state<D>& state, int tl, int rl, int c, int ml) { // Copy - copycat (tl,rl,c,ml, &dh<D>::dboxes::recv_sync, + copycat (state, + tl,rl,c,ml, &dh<D>::dboxes::recv_sync, tl,rl, ml, &dh<D>::dboxes::send_sync); } // Prolongate the boundaries of a component template<int D> -void ggf<D>::ref_bnd_prolongate (int tl, int rl, int c, int ml, +void ggf<D>::ref_bnd_prolongate (comm_state<D>& state, + int tl, int rl, int c, int ml, CCTK_REAL time) { // Interpolate @@ -443,53 +457,61 @@ void ggf<D>::ref_bnd_prolongate (int tl, int rl, int c, int ml, assert (tmax-tmin+1 >= prolongation_order_time+1); tl2s.resize(prolongation_order_time+1); for (int i=0; i<=prolongation_order_time; ++i) tl2s[i] = tmax - i; - intercat (tl ,rl ,c,ml, &dh<D>::dboxes::recv_ref_bnd_coarse, + intercat (state, + tl ,rl ,c,ml, &dh<D>::dboxes::recv_ref_bnd_coarse, tl2s,rl-1, ml, &dh<D>::dboxes::send_ref_bnd_fine, time); } // Restrict a multigrid level template<int D> -void ggf<D>::mg_restrict (int tl, int rl, int c, int ml, +void ggf<D>::mg_restrict (comm_state<D>& state, + int tl, int rl, int c, int ml, CCTK_REAL time) { // Require same times assert (t.get_time(rl,ml) == t.get_time(rl,ml-1)); const vector<int> tl2s(1,tl); - intercat (tl ,rl,c,ml, &dh<D>::dboxes::recv_mg_coarse, + intercat (state, + tl ,rl,c,ml, &dh<D>::dboxes::recv_mg_coarse, tl2s,rl, ml-1, &dh<D>::dboxes::send_mg_fine, time); } // Prolongate a multigrid level template<int D> -void ggf<D>::mg_prolongate (int tl, int rl, int c, int ml, +void ggf<D>::mg_prolongate (comm_state<D>& state, + int tl, int rl, int c, int ml, CCTK_REAL time) { // Require same times assert (t.get_time(rl,ml) == t.get_time(rl,ml+1)); const vector<int> tl2s(1,tl); - intercat (tl ,rl,c,ml, &dh<D>::dboxes::recv_mg_coarse, + intercat (state, + tl ,rl,c,ml, &dh<D>::dboxes::recv_mg_coarse, tl2s,rl, ml+1, &dh<D>::dboxes::send_mg_fine, time); } // Restrict a refinement level template<int D> -void ggf<D>::ref_restrict (int tl, int rl, int c, int ml, +void ggf<D>::ref_restrict (comm_state<D>& state, + int tl, int rl, int c, int ml, CCTK_REAL time) { // Require same times assert (t.get_time(rl,ml) == t.get_time(rl+1,ml)); const vector<int> tl2s(1,tl); - intercat (tl ,rl ,c,ml, &dh<D>::dboxes::recv_ref_fine, + intercat (state, + tl ,rl ,c,ml, &dh<D>::dboxes::recv_ref_fine, tl2s,rl+1, ml, &dh<D>::dboxes::send_ref_coarse, time); } // Prolongate a refinement level template<int D> -void ggf<D>::ref_prolongate (int tl, int rl, int c, int ml, +void ggf<D>::ref_prolongate (comm_state<D>& state, + int tl, int rl, int c, int ml, CCTK_REAL time) { assert (rl>=1); @@ -498,7 +520,8 @@ void ggf<D>::ref_prolongate (int tl, int rl, int c, int ml, assert (tmax-tmin+1 >= prolongation_order_time+1); tl2s.resize(prolongation_order_time+1); for (int i=0; i<=prolongation_order_time; ++i) tl2s[i] = tmax - i; - intercat (tl ,rl ,c,ml, &dh<D>::dboxes::recv_ref_coarse, + intercat (state, + tl ,rl ,c,ml, &dh<D>::dboxes::recv_ref_coarse, tl2s,rl-1, ml, &dh<D>::dboxes::send_ref_fine, time); } diff --git a/Carpet/CarpetLib/src/ggf.hh b/Carpet/CarpetLib/src/ggf.hh index 98c556800..b99c7aa10 100644 --- a/Carpet/CarpetLib/src/ggf.hh +++ b/Carpet/CarpetLib/src/ggf.hh @@ -1,4 +1,4 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.hh,v 1.16 2003/10/14 16:39:16 schnetter Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.hh,v 1.17 2003/11/05 16:18:39 schnetter Exp $ #ifndef GGF_HH #define GGF_HH @@ -80,7 +80,7 @@ public: // Modifiers // VGF - void recompose (const int initialise_upto = -1); + void recompose (const int initialise_from, const bool do_prolongate); // Cycle the time levels by rotating the data sets void cycle (int rl, int c, int ml); @@ -88,8 +88,11 @@ public: // Flip the time levels by exchanging the data sets void flip (int rl, int c, int ml); +#if 0 // Copy data from current time level to all previous time levels void copytoprevs (int rl, int c, int ml); +#endif + // Helpers @@ -105,39 +108,45 @@ protected: protected: // Copy a region - void copycat (int tl1, int rl1, int c1, int ml1, + void copycat (comm_state<D>& state, + int tl1, int rl1, int c1, int ml1, const ibbox dh<D>::dboxes::* recv_list, int tl2, int rl2, int ml2, const ibbox dh<D>::dboxes::* send_list); // Copy regions - void copycat (int tl1, int rl1, int c1, int ml1, + void copycat (comm_state<D>& state, + int tl1, int rl1, int c1, int ml1, const iblist dh<D>::dboxes::* recv_list, int tl2, int rl2, int ml2, const iblist dh<D>::dboxes::* send_list); // Copy regions - void copycat (int tl1, int rl1, int c1, int ml1, + void copycat (comm_state<D>& state, + int tl1, int rl1, int c1, int ml1, const iblistvect dh<D>::dboxes::* recv_listvect, int tl2, int rl2, int ml2, const iblistvect dh<D>::dboxes::* send_listvect); // Interpolate a region - void intercat (int tl1, int rl1, int c1, int ml1, + void intercat (comm_state<D>& state, + int tl1, int rl1, int c1, int ml1, const ibbox dh<D>::dboxes::* recv_list, const vector<int> tl2s, int rl2, int ml2, const ibbox dh<D>::dboxes::* send_list, CCTK_REAL time); // Interpolate regions - void intercat (int tl1, int rl1, int c1, int ml1, + void intercat (comm_state<D>& state, + int tl1, int rl1, int c1, int ml1, const iblist dh<D>::dboxes::* recv_list, const vector<int> tl2s, int rl2, int ml2, const iblist dh<D>::dboxes::* send_list, CCTK_REAL time); // Interpolate regions - void intercat (int tl1, int rl1, int c1, int ml1, + void intercat (comm_state<D>& state, + int tl1, int rl1, int c1, int ml1, const iblistvect dh<D>::dboxes::* recv_listvect, const vector<int> tl2s, int rl2, int ml2, const iblistvect dh<D>::dboxes::* send_listvect, @@ -154,25 +163,25 @@ public: // synchronised. They don't need to be prolongated. // Copy a component from the next time level - void copy (int tl, int rl, int c, int ml); + void copy (comm_state<D>& state, int tl, int rl, int c, int ml); // Synchronise the boundaries of a component - void sync (int tl, int rl, int c, int ml); + void sync (comm_state<D>& state, int tl, int rl, int c, int ml); // Prolongate the boundaries of a component - void ref_bnd_prolongate (int tl, int rl, int c, int ml, CCTK_REAL time); + void ref_bnd_prolongate (comm_state<D>& state, int tl, int rl, int c, int ml, CCTK_REAL time); // Restrict a multigrid level - void mg_restrict (int tl, int rl, int c, int ml, CCTK_REAL time); + void mg_restrict (comm_state<D>& state, int tl, int rl, int c, int ml, CCTK_REAL time); // Prolongate a multigrid level - void mg_prolongate (int tl, int rl, int c, int ml, CCTK_REAL time); + void mg_prolongate (comm_state<D>& state, int tl, int rl, int c, int ml, CCTK_REAL time); // Restrict a refinement level - void ref_restrict (int tl, int rl, int c, int ml, CCTK_REAL time); + void ref_restrict (comm_state<D>& state, int tl, int rl, int c, int ml, CCTK_REAL time); // Prolongate a refinement level - void ref_prolongate (int tl, int rl, int c, int ml, CCTK_REAL time); + void ref_prolongate (comm_state<D>& state, int tl, int rl, int c, int ml, CCTK_REAL time); diff --git a/Carpet/CarpetLib/src/gh.cc b/Carpet/CarpetLib/src/gh.cc index f8eeea4c0..b37fb3f15 100644 --- a/Carpet/CarpetLib/src/gh.cc +++ b/Carpet/CarpetLib/src/gh.cc @@ -1,4 +1,4 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.cc,v 1.22 2003/09/19 16:06:41 schnetter Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.cc,v 1.23 2003/11/05 16:18:39 schnetter Exp $ #include <assert.h> #include <stdlib.h> @@ -37,7 +37,8 @@ template<int D> void gh<D>::recompose (const rexts& exts, const rbnds& outer_bounds, const rprocs& procs, - const int initialise_upto) + const int initialise_from, + const bool do_prolongate) { DECLARE_CCTK_PARAMETERS; @@ -81,6 +82,9 @@ void gh<D>::recompose (const rexts& exts, assert (all(extents[rl][c][ml].stride() == extents[rl][0][ml].stride())); assert (extents[rl][c][ml].is_aligned_with(extents[rl][0][ml])); + for (int cc=c+1; cc<components(rl); ++cc) { + assert ((extents[rl][c][ml] & extents[rl][cc][ml]).empty()); + } } } } @@ -161,7 +165,7 @@ void gh<D>::recompose (const rexts& exts, } for (typename list<dh<D>*>::iterator d=dhs.begin(); d!=dhs.end(); ++d) { - (*d)->recompose(initialise_upto); + (*d)->recompose (initialise_from, do_prolongate); } } diff --git a/Carpet/CarpetLib/src/gh.hh b/Carpet/CarpetLib/src/gh.hh index c748e93ee..db046185a 100644 --- a/Carpet/CarpetLib/src/gh.hh +++ b/Carpet/CarpetLib/src/gh.hh @@ -1,4 +1,4 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.hh,v 1.13 2003/05/02 14:23:12 schnetter Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.hh,v 1.14 2003/11/05 16:18:39 schnetter Exp $ #ifndef GH_HH #define GH_HH @@ -86,7 +86,8 @@ public: // Modifiers void recompose (const rexts& exts, const rbnds& outer_bounds, const rprocs& procs, - const int initialise_upto = -1); + const int initialise_from, + const bool do_prolongate); // Helpers cexts make_reflevel_multigrid_boxes (const vector<ibbox>& exts, diff --git a/Carpet/CarpetLib/src/make.code.defn b/Carpet/CarpetLib/src/make.code.defn index ab3874e22..0b718747a 100644 --- a/Carpet/CarpetLib/src/make.code.defn +++ b/Carpet/CarpetLib/src/make.code.defn @@ -1,5 +1,5 @@ # Main make.code.defn file for thorn CarpetLib -*-Makefile-*- -# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/make.code.defn,v 1.9 2003/10/14 14:53:46 hawke Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/make.code.defn,v 1.10 2003/11/05 16:18:39 schnetter Exp $ # Source files in this directory SRCS = bbox.cc \ @@ -14,20 +14,27 @@ SRCS = bbox.cc \ gh.cc \ th.cc \ vect.cc \ + checkindex.F77 \ + copy_3d_int4.F77 \ copy_3d_real8.F77 \ prolongate_3d_real8.F77 \ + prolongate_3d_real8_rf2.F77 \ prolongate_3d_real8_o3.F77 \ + prolongate_3d_real8_o3_rf2.F77 \ prolongate_3d_real8_o5.F77 \ prolongate_3d_real8_2tl.F77 \ prolongate_3d_real8_2tl_o3.F77 \ prolongate_3d_real8_2tl_o5.F77 \ prolongate_3d_real8_3tl.F77 \ + prolongate_3d_real8_3tl_rf2.F77 \ prolongate_3d_real8_3tl_o3.F77 \ + prolongate_3d_real8_3tl_o3_rf2.F77 \ prolongate_3d_real8_3tl_o5.F77 \ prolongate_3d_real8_minmod.F77 \ prolongate_3d_real8_2tl_minmod.F77 \ prolongate_3d_real8_3tl_minmod.F77 \ - restrict_3d_real8.F77 + restrict_3d_real8.F77 \ + restrict_3d_real8_rf2.F77 # Subdirectories containing source files SUBDIRS = diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3_rf2.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3_rf2.F77 index daa130251..6c3cec79f 100644 --- a/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3_rf2.F77 +++ b/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3_rf2.F77 @@ -1,5 +1,5 @@ c -*-Fortran-*- -c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3_rf2.F77,v 1.3 2004/03/11 12:03:09 schnetter Exp $ +c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3_rf2.F77,v 1.1 2003/11/05 16:18:39 schnetter Exp $ #include "cctk.h" #include "cctk_Parameters.h" @@ -87,11 +87,6 @@ c bbox(:,3) is stride c This could be handled, but is likely to point to an error elsewhere call CCTK_WARN (0, "Internal error: region extent is empty") end if - if (mod(srcbbox(d,2) - srcbbox(d,1), srcbbox(d,3)).ne.0 - $ .or. mod(dstbbox(d,2) - dstbbox(d,1), dstbbox(d,3)).ne.0 - $ .or. mod(regbbox(d,2) - regbbox(d,1), regbbox(d,3)).ne.0) then - call CCTK_WARN (0, "Internal error: array extents are not integer multiples of the strides") - end if regkext = (regbbox(d,2) - regbbox(d,1)) / regbbox(d,3) + 1 srckoff = (regbbox(d,1) - srcbbox(d,1)) / dstbbox(d,3) offsetlo = regbbox(d,3) @@ -146,7 +141,7 @@ c Quadratic (second order) time interpolation call CCTK_WARN (0, "Internal error: arrays have same time") end if if (t.lt.min(t1,t2,t3)-eps .or. t.gt.max(t1,t2,t3)+eps) then - call CCTK_WARN (0, "Internal error: extrapolation in time in time") + call CCTK_WARN (0, "Internal error: extrapolation") end if s1fac = (t - t2) * (t - t3) / ((t1 - t2) * (t1 - t3)) diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_rf2.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_rf2.F77 index a77498fef..2395b0821 100644 --- a/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_rf2.F77 +++ b/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_rf2.F77 @@ -1,5 +1,5 @@ c -*-Fortran-*- -c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_rf2.F77,v 1.3 2004/03/11 12:03:09 schnetter Exp $ +c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_rf2.F77,v 1.1 2003/11/05 16:18:39 schnetter Exp $ #include "cctk.h" #include "cctk_Parameters.h" @@ -77,11 +77,6 @@ c bbox(:,3) is stride c This could be handled, but is likely to point to an error elsewhere call CCTK_WARN (0, "Internal error: region extent is empty") end if - if (mod(srcbbox(d,2) - srcbbox(d,1), srcbbox(d,3)).ne.0 - $ .or. mod(dstbbox(d,2) - dstbbox(d,1), dstbbox(d,3)).ne.0 - $ .or. mod(regbbox(d,2) - regbbox(d,1), regbbox(d,3)).ne.0) then - call CCTK_WARN (0, "Internal error: array extents are not integer multiples of the strides") - end if if (regbbox(d,1).lt.srcbbox(d,1) $ .or. regbbox(d,1).lt.dstbbox(d,1) $ .or. regbbox(d,2).gt.srcbbox(d,2) @@ -120,7 +115,7 @@ c Quadratic (second order) time interpolation call CCTK_WARN (0, "Internal error: arrays have same time") end if if (t.lt.min(t1,t2,t3)-eps .or. t.gt.max(t1,t2,t3)+eps) then - call CCTK_WARN (0, "Internal error: extrapolation in time") + call CCTK_WARN (0, "Internal error: extrapolation") end if s1fac = (t - t2) * (t - t3) / ((t1 - t2) * (t1 - t3)) diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_o3_rf2.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_o3_rf2.F77 index 8a3b2629d..ad101b6c6 100644 --- a/Carpet/CarpetLib/src/prolongate_3d_real8_o3_rf2.F77 +++ b/Carpet/CarpetLib/src/prolongate_3d_real8_o3_rf2.F77 @@ -1,5 +1,5 @@ c -*-Fortran-*- -c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_o3_rf2.F77,v 1.2 2004/03/11 12:03:09 schnetter Exp $ +c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_o3_rf2.F77,v 1.1 2003/11/05 16:18:39 schnetter Exp $ #include "cctk.h" #include "cctk_Parameters.h" @@ -74,11 +74,6 @@ c bbox(:,3) is stride c This could be handled, but is likely to point to an error elsewhere call CCTK_WARN (0, "Internal error: region extent is empty") end if - if (mod(srcbbox(d,2) - srcbbox(d,1), srcbbox(d,3)).ne.0 - $ .or. mod(dstbbox(d,2) - dstbbox(d,1), dstbbox(d,3)).ne.0 - $ .or. mod(regbbox(d,2) - regbbox(d,1), regbbox(d,3)).ne.0) then - call CCTK_WARN (0, "Internal error: array extents are not integer multiples of the strides") - end if regkext = (regbbox(d,2) - regbbox(d,1)) / regbbox(d,3) + 1 srckoff = (regbbox(d,1) - srcbbox(d,1)) / dstbbox(d,3) offsetlo = regbbox(d,3) diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_rf2.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_rf2.F77 index 556f4d092..6c59cc533 100644 --- a/Carpet/CarpetLib/src/prolongate_3d_real8_rf2.F77 +++ b/Carpet/CarpetLib/src/prolongate_3d_real8_rf2.F77 @@ -1,5 +1,5 @@ c -*-Fortran-*- -c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_rf2.F77,v 1.2 2004/03/11 12:03:09 schnetter Exp $ +c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_rf2.F77,v 1.1 2003/11/05 16:18:39 schnetter Exp $ #include "cctk.h" #include "cctk_Parameters.h" @@ -66,11 +66,6 @@ c bbox(:,3) is stride c This could be handled, but is likely to point to an error elsewhere call CCTK_WARN (0, "Internal error: region extent is empty") end if - if (mod(srcbbox(d,2) - srcbbox(d,1), srcbbox(d,3)).ne.0 - $ .or. mod(dstbbox(d,2) - dstbbox(d,1), dstbbox(d,3)).ne.0 - $ .or. mod(regbbox(d,2) - regbbox(d,1), regbbox(d,3)).ne.0) then - call CCTK_WARN (0, "Internal error: array extents are not integer multiples of the strides") - end if if (regbbox(d,1).lt.srcbbox(d,1) $ .or. regbbox(d,1).lt.dstbbox(d,1) $ .or. regbbox(d,2).gt.srcbbox(d,2) diff --git a/Carpet/CarpetLib/src/restrict_3d_real8_rf2.F77 b/Carpet/CarpetLib/src/restrict_3d_real8_rf2.F77 index b7f6a0454..71a8edfba 100644 --- a/Carpet/CarpetLib/src/restrict_3d_real8_rf2.F77 +++ b/Carpet/CarpetLib/src/restrict_3d_real8_rf2.F77 @@ -1,5 +1,5 @@ c -*-Fortran-*- -c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/restrict_3d_real8_rf2.F77,v 1.2 2004/03/11 12:03:09 schnetter Exp $ +c $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/restrict_3d_real8_rf2.F77,v 1.1 2003/11/05 16:18:39 schnetter Exp $ #include "cctk.h" #include "cctk_Parameters.h" @@ -43,7 +43,7 @@ c bbox(:,3) is stride call CCTK_WARN (0, "Internal error: strides disagree") end if if (dstbbox(d,3).ne.srcbbox(d,3)*2) then - call CCTK_WARN (0, "Internal error: destination strides are not twice the source strides") + call CCTK_WARN (0, "Internal error: destination strides are not twice the sourcd strides") end if if (mod(srcbbox(d,1), srcbbox(d,3)).ne.0 $ .or. mod(dstbbox(d,1), dstbbox(d,3)).ne.0 @@ -54,11 +54,6 @@ c bbox(:,3) is stride c This could be handled, but is likely to point to an error elsewhere call CCTK_WARN (0, "Internal error: region extent is empty") end if - if (mod(srcbbox(d,2) - srcbbox(d,1), srcbbox(d,3)).ne.0 - $ .or. mod(dstbbox(d,2) - dstbbox(d,1), dstbbox(d,3)).ne.0 - $ .or. mod(regbbox(d,2) - regbbox(d,1), regbbox(d,3)).ne.0) then - call CCTK_WARN (0, "Internal error: array extents are not integer multiples of the strides") - end if if (regbbox(d,1).lt.srcbbox(d,1) $ .or. regbbox(d,1).lt.dstbbox(d,1) $ .or. regbbox(d,2).gt.srcbbox(d,2) diff --git a/Carpet/CarpetLib/src/vect.hh b/Carpet/CarpetLib/src/vect.hh index 22d8ef054..dfb724694 100644 --- a/Carpet/CarpetLib/src/vect.hh +++ b/Carpet/CarpetLib/src/vect.hh @@ -1,4 +1,4 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/vect.hh,v 1.19 2003/08/28 21:07:27 schnetter Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/vect.hh,v 1.20 2003/11/05 16:18:39 schnetter Exp $ #ifndef VECT_HH #define VECT_HH @@ -58,7 +58,7 @@ public: elt[0]=x; elt[1]=y; } - /** Constructor for 3-element vectors from 4 elements. */ + /** Constructor for 3-element vectors from 3 elements. */ vect (const T x, const T y, const T z) { assert (D==3); elt[0]=x; elt[1]=y; elt[2]=z; @@ -149,7 +149,7 @@ public: /** Return a non-writable element of a vector. */ // Don't return a reference; *this might be a temporary - const T operator[] (const int d) const { + T operator[] (const int d) const { assert(d>=0 && d<D); return elt[d]; } diff --git a/Carpet/CarpetReduce/schedule.ccl b/Carpet/CarpetReduce/schedule.ccl index 82626dee2..f64908c6e 100644 --- a/Carpet/CarpetReduce/schedule.ccl +++ b/Carpet/CarpetReduce/schedule.ccl @@ -1,5 +1,5 @@ # Schedule definitions for thorn CarpetReduce -# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/schedule.ccl,v 1.1 2001/12/11 13:08:56 schnetter Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetReduce/schedule.ccl,v 1.2 2003/11/05 16:18:39 schnetter Exp $ schedule CarpetReduceStartup at STARTUP { diff --git a/Carpet/CarpetSlab/interface.ccl b/Carpet/CarpetSlab/interface.ccl index 96d182492..29a7298be 100644 --- a/Carpet/CarpetSlab/interface.ccl +++ b/Carpet/CarpetSlab/interface.ccl @@ -1,5 +1,5 @@ # Interface definition for thorn CarpetSlab -# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/interface.ccl,v 1.6 2003/09/04 16:23:22 tradke Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/interface.ccl,v 1.7 2003/11/05 16:18:39 schnetter Exp $ implements: Hyperslab @@ -16,3 +16,82 @@ uses include header: gdata.hh uses include header: dh.hh uses include header: ggf.hh uses include header: gh.hh + + + +CCTK_INT FUNCTION \ + Hyperslab_Get (CCTK_POINTER_TO_CONST IN cctkGH, \ + CCTK_INT IN mapping_handle, \ + CCTK_INT IN proc, \ + CCTK_INT IN vindex, \ + CCTK_INT IN timelevel, \ + CCTK_INT IN hdatatype, \ + CCTK_POINTER IN hdata) + +CCTK_INT FUNCTION \ + Hyperslab_GetList (CCTK_POINTER_TO_CONST IN cctkGH, \ + CCTK_INT IN mapping_handle, \ + CCTK_INT IN num_arrays, \ + CCTK_INT ARRAY IN procs, \ + CCTK_INT ARRAY IN vindices, \ + CCTK_INT ARRAY IN timelevels, \ + CCTK_INT ARRAY IN hdatatypes, \ + CCTK_POINTER ARRAY IN hdata, \ + CCTK_INT ARRAY OUT retvals) + +#CCTK_INT FUNCTION \ +# Hyperslab_DefineLocalMappingByIndex (CCTK_POINTER_TO_CONST IN cctkGH, \ +# CCTK_INT IN vindex, \ +# CCTK_INT IN hdim, \ +# CCTK_INT ARRAY IN direction, \ +# CCTK_INT ARRAY IN origin, \ +# CCTK_INT ARRAY IN extent, \ +# CCTK_INT ARRAY IN downsample, \ +# CCTK_INT IN table_handle, \ +# CCTK_INT CCTK_FPOINTER IN \ +# conversion_fn (CCTK_INT IN nelems, \ +# CCTK_INT IN src_stride, \ +# CCTK_INT IN dst_stride, \ +# CCTK_INT IN src_type, \ +# CCTK_INT IN dst_type, \ +# CCTK_POINTER_TO_CONST IN from, \ +# CCTK_POINTER IN to), \ +# CCTK_INT ARRAY OUT hsize_local, \ +# CCTK_INT ARRAY OUT hsize_global, \ +# CCTK_INT ARRAY OUT hoffset_global) + +CCTK_INT FUNCTION \ + Hyperslab_DefineGlobalMappingByIndex (CCTK_POINTER_TO_CONST IN cctkGH, \ + CCTK_INT IN vindex, \ + CCTK_INT IN hdim, \ + CCTK_INT ARRAY IN direction, \ + CCTK_INT ARRAY IN origin, \ + CCTK_INT ARRAY IN extent, \ + CCTK_INT ARRAY IN downsample, \ + CCTK_INT IN table_handle, \ + CCTK_INT CCTK_FPOINTER IN \ + conversion_fn (CCTK_INT IN nelems, \ + CCTK_INT IN src_stride, \ + CCTK_INT IN dst_stride, \ + CCTK_INT IN src_type, \ + CCTK_INT IN dst_type, \ + CCTK_POINTER_TO_CONST IN from, \ + CCTK_POINTER IN to), \ + CCTK_INT ARRAY OUT hsize) + +CCTK_INT FUNCTION Hyperslab_FreeMapping (CCTK_INT IN mapping_handle) + + + +PROVIDES FUNCTION Hyperslab_Get \ + WITH CarpetSlab_Get LANGUAGE C +PROVIDES FUNCTION Hyperslab_GetList \ + WITH CarpetSlab_GetList LANGUAGE C +PROVIDES FUNCTION Hyperslab_Get \ + WITH CarpetSlab_Fill LANGUAGE C +PROVIDES FUNCTION Hyperslab_DefineGlobalMappingByIndex \ + WITH CarpetSlab_DefineGlobalMappingByIndex LANGUAGE C +#PROVIDES FUNCTION Hyperslab_DefineLocalMappingByIndex \ +# WITH CarpetSlab_DefineLocalMappingByIndex LANGUAGE C +PROVIDES FUNCTION Hyperslab_FreeMapping \ + WITH CarpetSlab_FreeMapping LANGUAGE C diff --git a/Carpet/CarpetSlab/src/slab.cc b/Carpet/CarpetSlab/src/slab.cc index f6d613441..d72499b79 100644 --- a/Carpet/CarpetSlab/src/slab.cc +++ b/Carpet/CarpetSlab/src/slab.cc @@ -1,4 +1,4 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/slab.cc,v 1.12 2003/10/14 16:39:17 schnetter Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/slab.cc,v 1.13 2003/11/05 16:18:39 schnetter Exp $ #include <assert.h> #include <stdlib.h> @@ -8,6 +8,8 @@ #include "cctk.h" +#include "util_Table.h" + #include "bbox.hh" #include "bboxset.hh" #include "dh.hh" @@ -21,7 +23,7 @@ #include "slab.hh" extern "C" { - static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/slab.cc,v 1.12 2003/10/14 16:39:17 schnetter Exp $"; + static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/slab.cc,v 1.13 2003/11/05 16:18:39 schnetter Exp $"; CCTK_FILEVERSION(Carpet_CarpetSlab_slab_cc); } @@ -33,7 +35,249 @@ namespace CarpetSlab { - void* GetSlab (cGH* const cgh, + // Mapping object + // (just store the mapping) + struct mapping { + int vindex; + int hdim; + int * origin; // [vdim] + int * dirs; // [hdim] + int * stride; // [hdim] + int * length; // [hdim] + }; + + + + int StoreMapping (mapping * const mp) + { + int const table = Util_TableCreate (UTIL_TABLE_FLAGS_DEFAULT); + assert (table>=0); + int const ierr = Util_TableSetPointer (table, mp, "mapping"); + assert (ierr>=0); + return table; + } + + mapping * RetrieveMapping (int const table) + { + CCTK_POINTER mp; + int const ierr = Util_TableGetPointer (table, &mp, "mapping"); + assert (ierr>=0); + return (mapping *)mp; + } + + void DeleteMapping (int const table) + { + int const ierr = Util_TableDestroy (table); + assert (ierr>=0); + } + + + + void FillSlab (const cGH* const cgh, + const int dest_proc, + const int n, + const int ti, + const int hdim, + const int origin[/*vdim*/], + const int dirs[/*hdim*/], + const int stride[/*hdim*/], + const int length[/*hdim*/], + void* const hdata) + { + // Check Cactus grid hierarchy + assert (cgh); + + // Check destination processor + assert (dest_proc>=-1 && dest_proc<CCTK_nProcs(cgh)); + + // Check variable index + assert (n>=0 && n<CCTK_NumVars()); + + // Get info about variable + const int group = CCTK_GroupIndexFromVarI(n); + assert (group>=0); + const int n0 = CCTK_FirstVarIndexI(group); + assert (n0>=0); + const int var = n - n0; + assert (var>=0); + + // Get info about group + cGroup gp; + CCTK_GroupData (group, &gp); + assert (gp.dim<=dim); + assert (CCTK_QueryGroupStorageI(cgh, group)); + const int typesize = CCTK_VarTypeSize(gp.vartype); + assert (typesize>0); + + if (gp.grouptype==CCTK_GF && reflevel==-1) { + CCTK_WARN (0, "It is not possible to use hyperslabbing for a grid function in global mode"); + } + const int rl = gp.grouptype==CCTK_GF ? reflevel : 0; + + // Check dimension + assert (hdim>=0 && hdim<=gp.dim); + + // Check timelevel + const int num_tl = gp.numtimelevels; + assert (ti>=0 && ti<num_tl); + const int tl = -ti; + + // Check origin +// for (int d=0; d<dim; ++d) { +// assert (origin[d]>=0 && origin[d]<=sizes[d]); +// } + + // Check directions + for (int dd=0; dd<hdim; ++dd) { + assert (dirs[dd]>=1 && dirs[dd]<=dim); + } + + // Check stride + for (int dd=0; dd<hdim; ++dd) { + assert (stride[dd]>0); + } + + // Check length + for (int dd=0; dd<hdim; ++dd) { + assert (length[dd]>=0); + } + + // Check extent +// for (int dd=0; dd<hdim; ++dd) { +// assert (origin[dirs[dd]-1] + length[dd] <= sizes[dirs[dd]]); +// } + + // Get insider information about variable + const gh<dim>* myhh; + const dh<dim>* mydd; + const ggf<dim>* myff; + assert (group < (int)arrdata.size()); + myhh = arrdata[group].hh; + assert (myhh); + mydd = arrdata[group].dd; + assert (mydd); + assert (var < (int)arrdata[group].data.size()); + myff = arrdata[group].data[var]; + assert (myff); + + // Detemine collecting processor + const int collect_proc = dest_proc<0 ? 0 : dest_proc; + + // Determine own rank + const int rank = CCTK_MyProc(cgh); + + // Calculate global size + int totalsize = 1; + for (int dd=0; dd<hdim; ++dd) { + totalsize *= length[dd]; + } + + // Allocate memory + assert (hdata); + if (dest_proc==-1 || rank==dest_proc) { + memset (hdata, 0, totalsize * typesize); + } + + // Get sample data + const gdata<dim>* mydata; + mydata = (*myff)(tl, rl, 0, 0); + + // Stride of data in memory + const vect<int,dim> str = mydata->extent().stride(); + + // Stride of collected data + vect<int,dim> hstr = str; + for (int dd=0; dd<hdim; ++dd) { + hstr[dirs[dd]-1] *= stride[dd]; + } + + // Lower bound of collected data + vect<int,dim> hlb(0); + for (int d=0; d<gp.dim; ++d) { + hlb[d] = origin[d] * str[d]; + } + + // Upper bound of collected data + vect<int,dim> hub = hlb; + for (int dd=0; dd<hdim; ++dd) { + hub[dirs[dd]-1] += (length[dd]-1) * hstr[dirs[dd]-1]; + } + + // Calculate extent to collect + const bbox<int,dim> hextent (hlb, hub, hstr); + assert (hextent.size() == totalsize); + + // Create collector data object + void* myhdata = rank==collect_proc ? hdata : 0; + gdata<dim>* const alldata = mydata->make_typed(-1); + alldata->allocate (hextent, collect_proc, myhdata); + + // Done with the temporary stuff + mydata = 0; + + for (comm_state<dim> state; !state.done(); state.step()) { + + // Loop over all components, copying data from them + BEGIN_LOCAL_COMPONENT_LOOP (cgh, gp.grouptype) { + + // Get data object + mydata = (*myff)(tl, rl, component, mglevel); + + // Calculate overlapping extents + const bboxset<int,dim> myextents + = ((mydd->boxes[rl][component][mglevel].sync_not + | mydd->boxes[rl][component][mglevel].interior) + & hextent); + + // Loop over overlapping extents + for (bboxset<int,dim>::const_iterator ext_iter = myextents.begin(); + ext_iter != myextents.end(); + ++ext_iter) { + + // Copy data + alldata->copy_from (state, mydata, *ext_iter); + + } + + } END_LOCAL_COMPONENT_LOOP; + + } // for step + + // Copy result to all processors + if (dest_proc == -1) { + vector<gdata<dim>*> tmpdata(CCTK_nProcs(cgh)); + vector<comm_state<dim> > state; + + for (int proc=0; proc<CCTK_nProcs(cgh); ++proc) { + if (proc != collect_proc) { + void* myhdata = rank==proc ? hdata : 0; + tmpdata[proc] = mydata->make_typed(-1); + tmpdata[proc]->allocate (alldata->extent(), proc, myhdata); + tmpdata[proc]->copy_from (state[proc], alldata, alldata->extent()); + } + } + + for (int proc=0; proc<CCTK_nProcs(cgh); ++proc) { + if (proc != collect_proc) { + tmpdata[proc]->copy_from (state[proc], alldata, alldata->extent()); + } + } + + for (int proc=0; proc<CCTK_nProcs(cgh); ++proc) { + if (proc != collect_proc) { + tmpdata[proc]->copy_from (state[proc], alldata, alldata->extent()); + delete tmpdata[proc]; + } + } + + } // Copy result + + delete alldata; + } + + + + void* GetSlab (const cGH* const cgh, const int dest_proc, const int n, const int ti, @@ -134,6 +378,7 @@ namespace CarpetSlab { // Allocate memory void* hdata = 0; if (dest_proc==-1 || rank==dest_proc) { + assert (0); hdata = malloc(totalsize * typesize); assert (hdata); memset (hdata, 0, totalsize * typesize); @@ -176,45 +421,61 @@ namespace CarpetSlab { // Done with the temporary stuff mydata = 0; - // Loop over all components, copying data from them - BEGIN_LOCAL_COMPONENT_LOOP (cgh, gp.grouptype) { - - const int c = gp.grouptype==CCTK_GF ? component : 0; + for (comm_state<dim> state; !state.done(); state.step()) { - // Get data object - mydata = (*myff)(tl, rl, c, mglevel); - - // Calculate overlapping extents - const bboxset<int,dim> myextents - = ((mydd->boxes[rl][c][mglevel].sync_not - | mydd->boxes[rl][c][mglevel].interior) - & hextent); - - // Loop over overlapping extents - for (bboxset<int,dim>::const_iterator ext_iter = myextents.begin(); - ext_iter != myextents.end(); - ++ext_iter) { + // Loop over all components, copying data from them + BEGIN_LOCAL_COMPONENT_LOOP (cgh, gp.grouptype) { - // Copy data - alldata->copy_from (mydata, *ext_iter); + // Get data object + mydata = (*myff)(tl, rl, component, mglevel); - } + // Calculate overlapping extents + const bboxset<int,dim> myextents + = ((mydd->boxes[rl][component][mglevel].sync_not + | mydd->boxes[rl][component][mglevel].interior) + & hextent); + + // Loop over overlapping extents + for (bboxset<int,dim>::const_iterator ext_iter = myextents.begin(); + ext_iter != myextents.end(); + ++ext_iter) { + + // Copy data + alldata->copy_from (state, mydata, *ext_iter); + + } + + } END_LOCAL_COMPONENT_LOOP; - } END_LOCAL_COMPONENT_LOOP; + } // for step // Copy result to all processors if (dest_proc == -1) { + vector<gdata<dim>*> tmpdata(CCTK_nProcs(cgh)); + vector<comm_state<dim> > state; + for (int proc=0; proc<CCTK_nProcs(cgh); ++proc) { if (proc != collect_proc) { - void* myhdata = rank==proc ? hdata : 0; - gdata<dim>* const tmpdata = mydata->make_typed(-1); - tmpdata->allocate (alldata->extent(), proc, myhdata); - tmpdata->copy_from (alldata, alldata->extent()); - delete tmpdata; - + tmpdata[proc] = mydata->make_typed(-1); + tmpdata[proc]->allocate (alldata->extent(), proc, myhdata); + tmpdata[proc]->copy_from (state[proc], alldata, alldata->extent()); } } + + for (int proc=0; proc<CCTK_nProcs(cgh); ++proc) { + if (proc != collect_proc) { + tmpdata[proc]->copy_from (state[proc], alldata, alldata->extent()); + } + } + + for (int proc=0; proc<CCTK_nProcs(cgh); ++proc) { + if (proc != collect_proc) { + tmpdata[proc]->copy_from (state[proc], alldata, alldata->extent()); + delete tmpdata[proc]; + } + } + } // Copy result delete alldata; @@ -225,6 +486,271 @@ namespace CarpetSlab { + CCTK_INT CarpetSlab_Get (cGH const * const cctkGH, + CCTK_INT const mapping_handle, + CCTK_INT const proc, + CCTK_INT const vindex, + CCTK_INT const timelevel, + CCTK_INT const hdatatype, + void * const hdata) + { + // Check arguments + assert (cctkGH); + assert (mapping_handle>=0); + assert (proc==-1 || proc>=0 && proc<CCTK_nProcs(cctkGH)); + assert (vindex>=0 && vindex<CCTK_NumVars()); + assert (timelevel>=0); + assert (hdatatype>=0); + assert (hdata); + + // Get mapping + const mapping * const mp = RetrieveMapping (mapping_handle); + assert (mp); + + // Calculate total size + size_t size = 1; + for (size_t d=0; d<(size_t)mp->hdim; ++d) { + size *= mp->length[d]; + } + + // Get type size + size_t const sz = CCTK_VarTypeSize (hdatatype); + assert (sz>0); + + // Forward call + FillSlab (cctkGH, proc, vindex, timelevel, + mp->hdim, mp->origin, mp->dirs, mp->stride, mp->length, hdata); + + return 0; + } + + + + CCTK_INT CarpetSlab_GetList (cGH const * const cctkGH, + CCTK_INT const mapping_handle, + CCTK_INT const num_arrays, + CCTK_INT const * const procs, + CCTK_INT const * const vindices, + CCTK_INT const * const timelevels, + CCTK_INT const * const hdatatypes, + void * const * const hdata, + CCTK_INT * const retvals) + { + // Check arguments + assert (cctkGH); + assert (mapping_handle>=0); + assert (num_arrays>=0); + assert (procs); + assert (vindices); + assert (timelevels); + assert (hdatatypes); + assert (hdata); + assert (retvals); + + // Remember whether there were errors + bool everyting_okay = true; + + // Loop over all slabs + for (int n=0; n<num_arrays; ++n) { + // Forward call + retvals[n] = CarpetSlab_Get (cctkGH, mapping_handle, procs[n], + vindices[n], timelevels[n], hdatatypes[n], + hdata[n]); + everyting_okay = everyting_okay && retvals[n]; + } + + return everyting_okay ? 0 : -1; + } + + + + CCTK_INT CarpetSlab_DefineLocalMappingByIndex (cGH const * const cctkGH, + CCTK_INT const vindex, + CCTK_INT const hdim, + CCTK_INT const * const direction, + CCTK_INT const * const origin, + CCTK_INT const * const extent, + CCTK_INT const * const downsample, + CCTK_INT const table_handle, + CCTK_INT (* const conversion_fn) (CCTK_INT const nelems, + CCTK_INT const src_stride, + CCTK_INT const dst_stride, + CCTK_INT const src_type, + CCTK_INT const dst_type, + void const * const from, + void * const to), + CCTK_INT * const hsize_local, + CCTK_INT * const hsize_global, + CCTK_INT * const hoffset_global) + { + CCTK_WARN (0, "not implemented"); + return 0; + } + + + + CCTK_INT CarpetSlab_DefineGlobalMappingByIndex (cGH const * const cctkGH, + CCTK_INT const vindex, + CCTK_INT const hdim, + CCTK_INT const * const direction, + CCTK_INT const * const origin, + CCTK_INT const * const extent, + CCTK_INT const * const downsample, + CCTK_INT const table_handle, + CCTK_INT (* const conversion_fn) (CCTK_INT const nelems, + CCTK_INT const src_stride, + CCTK_INT const dst_stride, + CCTK_INT const src_type, + CCTK_INT const dst_type, + void const * const from, + void * const to), + CCTK_INT * const hsize) + { + // Check arguments + assert (cctkGH); + assert (vindex>=0 && vindex<CCTK_NumVars()); + assert (hdim>=0 && hdim<=dim); + assert (direction); + assert (origin); + assert (extent); + assert (downsample); + assert (table_handle>=0); + assert (hsize); + + // Get more information + int const vdim = CCTK_GroupDimFromVarI (vindex); + assert (vdim>=0 && vdim<dim); + assert (hdim<=vdim); + + // Not implemented + assert (! conversion_fn); + + // Allocate memory + mapping * mp = new mapping; + mp->origin = new int[vdim]; + mp->dirs = new int[hdim]; + mp->stride = new int[hdim]; + mp->length = new int[hdim]; + + // Calculate more convenient representation of the direction + int dirs[dim]; // should really be dirs[hdim] + // The following if statement is written according to the + // definition of "dir". + if (hdim==1) { + // 1-dimensional hyperslab + int mydir = 0; + for (int d=0; d<vdim; ++d) { + if (direction[d]!=0) { + mydir = d+1; + break; + } + } + assert (mydir>0); + for (int d=0; d<vdim; ++d) { + if (d == mydir-1) { + assert (direction[d]!=0); + } else { + assert (direction[d]==0); + } + } + dirs[0] = mydir; + } else if (hdim==vdim) { + // vdim-dimensional hyperslab + for (int d=0; d<vdim; ++d) { + dirs[d] = d+1; + } + } else if (hdim==2) { + // 2-dimensional hyperslab with vdim==3 + assert (vdim==3); + int mydir = 0; + for (int d=0; d<vdim; ++d) { + if (direction[d]==0) { + mydir = d+1; + break; + } + } + assert (mydir>0); + for (int d=0; d<vdim; ++d) { + if (d == mydir-1) { + assert (direction[d]==0); + } else { + assert (direction[d]!=0); + } + } + int dd=0; + for (int d=0; d<vdim; ++d) { + if (d != mydir-1) { + dirs[dd] = d+1; + ++dd; + } + } + assert (dd==hdim); + } else { + assert (0); + } + // Fill remaining length + for (int d=vdim; d<dim; ++d) { + dirs[d] = d+1; + } + + // Calculate lengths + for (int dd=0; dd<hdim; ++dd) { + if (extent[dd]<0) { + int gsh[dim]; + int ierr = CCTK_GroupgshVI(cctkGH, dim, gsh, vindex); + assert (!ierr); + const int totlen = gsh[dirs[dd]-1]; + assert (totlen>=0); + // Partial argument check + assert (origin[dirs[dd]-1]>=0); + assert (origin[dirs[dd]-1]<=totlen); + assert (downsample[dd]>0); + hsize[dd] = (totlen - origin[dirs[dd]-1]) / downsample[dd]; + } else { + hsize[dd] = extent[dd]; + } + assert (hsize[dd]>=0); + } + + // Store information + mp->vindex = vindex; + mp->hdim = hdim; + for (size_t d=0; d<(size_t)hdim; ++d) { + mp->origin[d] = origin[d]; + mp->dirs[d] = dirs[d]; + mp->stride[d] = downsample[d]; + mp->length[d] = hsize[d]; + } + + return 0; + } + + + + CCTK_INT CarpetSlab_FreeMapping (CCTK_INT const mapping_handle) + { + // Check arguments + assert (mapping_handle>=0); + + // Get mapping + mapping * mp = RetrieveMapping (mapping_handle); + assert (mp); + + // Delete storage + DeleteMapping (mapping_handle); + + // Free memory + delete [] mp->origin; + delete [] mp->dirs; + delete [] mp->stride; + delete [] mp->length; + delete mp; + + return 0; + } + + + int Hyperslab_GetHyperslab (cGH* const GH, const int target_proc, const int vindex, diff --git a/Carpet/CarpetSlab/src/slab.h b/Carpet/CarpetSlab/src/slab.h index c72222700..e57205aa0 100644 --- a/Carpet/CarpetSlab/src/slab.h +++ b/Carpet/CarpetSlab/src/slab.h @@ -1,4 +1,4 @@ -/* $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/slab.h,v 1.1 2002/10/24 10:53:48 schnetter Exp $ */ +/* $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/slab.h,v 1.2 2003/11/05 16:18:39 schnetter Exp $ */ #ifndef CARPETSLAB_H #define CARPETSLAB_H @@ -10,6 +10,65 @@ namespace CarpetSlab { extern "C" { #endif + CCTK_INT CarpetSlab_Get (cGH const * const cctkGH, + CCTK_INT const mapping_handle, + CCTK_INT const proc, + CCTK_INT const vindex, + CCTK_INT const timelevel, + CCTK_INT const hdatatype, + void * const hdata); + + CCTK_INT CarpetSlab_GetList (cGH const * const cctkGH, + CCTK_INT const mapping_handle, + CCTK_INT const num_arrays, + CCTK_INT const * const procs, + CCTK_INT const * const vindices, + CCTK_INT const * const timelevels, + CCTK_INT const * const hdatatypes, + void * const * const hdata, + CCTK_INT * const retvals); + + CCTK_INT CarpetSlab_DefineLocalMappingByIndex (cGH const * const cctkGH, + CCTK_INT const vindex, + CCTK_INT const hdim, + CCTK_INT const * const direction, + CCTK_INT const * const origin, + CCTK_INT const * const extent, + CCTK_INT const * const downsample, + CCTK_INT const table_handle, + CCTK_INT (* const conversion_fn) (CCTK_INT const nelems, + CCTK_INT const src_stride, + CCTK_INT const dst_stride, + CCTK_INT const src_type, + CCTK_INT const dst_type, + void const * const from, + void * const to), + CCTK_INT * const hsize_local, + CCTK_INT * const hsize_global, + CCTK_INT * const hoffset_global); + + CCTK_INT CarpetSlab_DefineGlobalMappingByIndex (cGH const * const cctkGH, + CCTK_INT const vindex, + CCTK_INT const hdim, + CCTK_INT const * const direction, + CCTK_INT const * const origin, + CCTK_INT const * const extent, + CCTK_INT const * const downsample, + CCTK_INT const table_handle, + CCTK_INT (* const conversion_fn) (CCTK_INT const nelems, + CCTK_INT const src_stride, + CCTK_INT const dst_stride, + CCTK_INT const src_type, + CCTK_INT const dst_type, + void const * const from, + void * const to), + CCTK_INT * const hsize); + + CCTK_INT CarpetSlab_FreeMapping (CCTK_INT const mapping_handle); + + + + /* Old interface -- don't use */ int Hyperslab_GetHyperslab (cGH* const GH, const int target_proc, const int vindex, diff --git a/CarpetAttic/CarpetIOFlexIO/src/ioflexio.cc b/CarpetAttic/CarpetIOFlexIO/src/ioflexio.cc index 7456091de..2a00ebe5c 100644 --- a/CarpetAttic/CarpetIOFlexIO/src/ioflexio.cc +++ b/CarpetAttic/CarpetIOFlexIO/src/ioflexio.cc @@ -305,7 +305,9 @@ namespace CarpetIOFlexIO { gdata<dim>* const tmp = data->make_typed (n); tmp->allocate (ext, 0); - tmp->copy_from (data, ext); + for (comm_state<dim> state; !state.done(); state.step()) { + tmp->copy_from (state, data, ext); + } // Write data if (CCTK_MyProc(cgh)==0) { @@ -600,7 +602,9 @@ namespace CarpetIOFlexIO { } // Copy into grid function - data->copy_from (tmp, ext); + for (comm_state<dim> state; !state.done(); state.step()) { + data->copy_from (state, tmp, ext); + } // Delete temporary copy delete tmp; diff --git a/CarpetDev/CarpetJacobi/par/charge_te_cjac.par b/CarpetDev/CarpetJacobi/par/charge_te_cjac.par index 3e77742dd..55e601b39 100644 --- a/CarpetDev/CarpetJacobi/par/charge_te_cjac.par +++ b/CarpetDev/CarpetJacobi/par/charge_te_cjac.par @@ -1,4 +1,4 @@ -# $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/par/charge_te_cjac.par,v 1.1 2003/09/02 14:35:58 schnetter Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/par/charge_te_cjac.par,v 1.2 2003/11/05 16:18:38 schnetter Exp $ # Erik Schnetter <schnetter@uni-tuebingen.de> ActiveThorns = "Boundary CartGrid3D CarpetIOASCII IOBasic IOUtil Time Carpet CarpetLib CarpetRegrid CarpetReduce NaNCatcher IDScalarWave WaveToyC IDSWTEsimple TATelliptic CarpetJacobi" @@ -31,7 +31,7 @@ CarpetRegrid::outerbounds = "[[ [[1,0],[1,0],[1,0]] ]]" grid::type = byspacing grid::domain = octant -grid::dxyz = 0.6 +grid::dxyz = 1.2 grid::avoid_origin = no IDScalarWave::initial_data = charge-TATelliptic-simple diff --git a/CarpetDev/CarpetJacobi/par/charge_te_cjac_2.par b/CarpetDev/CarpetJacobi/par/charge_te_cjac_2.par index 6f0131e74..8489ceab8 100644 --- a/CarpetDev/CarpetJacobi/par/charge_te_cjac_2.par +++ b/CarpetDev/CarpetJacobi/par/charge_te_cjac_2.par @@ -1,4 +1,4 @@ -# $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/par/charge_te_cjac_2.par,v 1.1 2003/09/02 14:35:58 schnetter Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/par/charge_te_cjac_2.par,v 1.2 2003/11/05 16:18:38 schnetter Exp $ # Erik Schnetter <schnetter@uni-tuebingen.de> ActiveThorns = "Boundary CartGrid3D CarpetIOASCII IOBasic IOUtil Time Carpet CarpetLib CarpetRegrid CarpetReduce NaNCatcher IDScalarWave WaveToyC IDSWTEsimple TATelliptic CarpetJacobi" @@ -31,7 +31,7 @@ CarpetRegrid::outerbounds = "[[ [[1,0],[1,0],[1,0]] ]]" grid::type = byspacing grid::domain = octant -grid::dxyz = 0.6 +grid::dxyz = 1.2 grid::avoid_origin = no IDScalarWave::initial_data = charge-TATelliptic-simple diff --git a/CarpetExtra/FOWaveToyF77/interface.ccl b/CarpetExtra/FOWaveToyF77/interface.ccl index 4c8131958..22f12f5db 100644 --- a/CarpetExtra/FOWaveToyF77/interface.ccl +++ b/CarpetExtra/FOWaveToyF77/interface.ccl @@ -1,5 +1,5 @@ # Interface definition for thorn WaveToyF77 -# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/FOWaveToyF77/interface.ccl,v 1.7 2003/06/27 16:44:03 schnetter Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/FOWaveToyF77/interface.ccl,v 1.8 2003/11/05 16:18:39 schnetter Exp $ implements: fowavetoy inherits: boundary grid idfoscalarwave @@ -21,7 +21,7 @@ CCTK_REAL scalarevolve_derivs type=GF timelevels=3 phiz, } "Time and space derivatives of phi" -CCTK_INT FUNCTION Boundary_SelectGroupForBC(CCTK_POINTER IN GH, \ +CCTK_INT FUNCTION Boundary_SelectGroupForBC(CCTK_POINTER_TO_CONST IN GH, \ CCTK_INT IN faces, CCTK_INT IN boundary_width, CCTK_INT IN table_handle, \ CCTK_STRING IN group_name, CCTK_STRING IN bc_name) USES FUNCTION Boundary_SelectGroupForBC diff --git a/CarpetExtra/IDScalarWaveFO/schedule.ccl b/CarpetExtra/IDScalarWaveFO/schedule.ccl index 4b68757d8..0a4a3407f 100644 --- a/CarpetExtra/IDScalarWaveFO/schedule.ccl +++ b/CarpetExtra/IDScalarWaveFO/schedule.ccl @@ -1,5 +1,5 @@ # Schedule definitions for thorn IDScalarWaveFO -# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/IDScalarWaveFO/schedule.ccl,v 1.1 2003/06/18 18:24:29 schnetter Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/IDScalarWaveFO/schedule.ccl,v 1.2 2003/11/05 16:18:40 schnetter Exp $ SCHEDULE IDScalarWaveFO_InitialData AT initial { diff --git a/CarpetExtra/IDScalarWaveFO/src/initialdata.F77 b/CarpetExtra/IDScalarWaveFO/src/initialdata.F77 index 59774e2b7..f65fc7084 100644 --- a/CarpetExtra/IDScalarWaveFO/src/initialdata.F77 +++ b/CarpetExtra/IDScalarWaveFO/src/initialdata.F77 @@ -1,7 +1,8 @@ -c $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/IDScalarWaveFO/src/initialdata.F77,v 1.2 2003/06/18 20:08:29 schnetter Exp $ +c $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/IDScalarWaveFO/src/initialdata.F77,v 1.3 2003/11/05 16:18:40 schnetter Exp $ #include "cctk.h" #include "cctk_Arguments.h" +#include "cctk_Functions.h" #include "cctk_Parameters.h" subroutine IDScalarWaveFO_InitialData (CCTK_ARGUMENTS) diff --git a/CarpetExtra/SpaceTimeToy/interface.ccl b/CarpetExtra/SpaceTimeToy/interface.ccl index 3d91b639b..66ac49c7a 100644 --- a/CarpetExtra/SpaceTimeToy/interface.ccl +++ b/CarpetExtra/SpaceTimeToy/interface.ccl @@ -1,12 +1,12 @@ # Interface definition for thorn SpaceTimeToy -# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/SpaceTimeToy/interface.ccl,v 1.6 2003/06/18 18:24:29 schnetter Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/SpaceTimeToy/interface.ccl,v 1.7 2003/11/05 16:18:40 schnetter Exp $ implements: spacetimetoy inherits: hydrotoy boundary driver grid -CCTK_INT FUNCTION Boundary_SelectGroupForBC(CCTK_POINTER IN GH, \ +CCTK_INT FUNCTION Boundary_SelectGroupForBC(CCTK_POINTER_TO_CONST IN GH, \ CCTK_INT IN faces, CCTK_INT IN boundary_width, CCTK_INT IN table_handle, \ CCTK_STRING IN group_name, CCTK_STRING IN bc_name) USES FUNCTION Boundary_SelectGroupForBC diff --git a/CarpetExtra/WaveToyExpl/interface.ccl b/CarpetExtra/WaveToyExpl/interface.ccl index f565d69c8..f05166930 100644 --- a/CarpetExtra/WaveToyExpl/interface.ccl +++ b/CarpetExtra/WaveToyExpl/interface.ccl @@ -1,5 +1,5 @@ # Interface definition for thorn WaveToyExpl -# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/WaveToyExpl/interface.ccl,v 1.1 2003/06/18 18:24:30 schnetter Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/WaveToyExpl/interface.ccl,v 1.2 2003/11/05 16:18:40 schnetter Exp $ IMPLEMENTS: WaveToyExpl @@ -7,7 +7,7 @@ INHERITS: boundary grid -CCTK_INT FUNCTION Boundary_SelectGroupForBC(CCTK_POINTER IN GH, \ +CCTK_INT FUNCTION Boundary_SelectGroupForBC(CCTK_POINTER_TO_CONST IN GH, \ CCTK_INT IN faces, CCTK_INT IN boundary_width, CCTK_INT IN table_handle, \ CCTK_STRING IN group_name, CCTK_STRING IN bc_name) diff --git a/CarpetExtra/WaveToyFO/schedule.ccl b/CarpetExtra/WaveToyFO/schedule.ccl index 1455dade7..7d48aba67 100644 --- a/CarpetExtra/WaveToyFO/schedule.ccl +++ b/CarpetExtra/WaveToyFO/schedule.ccl @@ -1,5 +1,5 @@ # Schedule definitions for thorn WaveToyFO -# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/WaveToyFO/schedule.ccl,v 1.3 2003/07/20 21:03:43 schnetter Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/WaveToyFO/schedule.ccl,v 1.4 2003/11/05 16:18:41 schnetter Exp $ STORAGE: scalarevolve[3] STORAGE: scalarevolvedot @@ -42,12 +42,12 @@ SCHEDULE GROUP ApplyBCs IN MoL_PostStep AFTER WaveToyFO_Boundaries -SCHEDULE WaveToyFO_Boundaries IN PostRestrict +SCHEDULE WaveToyFO_Boundaries AT POSTRESTRICT { LANG: Fortran SYNC: scalarevolve } "Select boundary conditions after restricting" -SCHEDULE GROUP ApplyBCs IN PostRestrict AFTER WaveToyFO_Boundaries +SCHEDULE GROUP ApplyBCs AT POSTRESTRICT AFTER WaveToyFO_Boundaries { } "Apply boundary conditions after restricting" diff --git a/CarpetExtra/WaveToyMoL/interface.ccl b/CarpetExtra/WaveToyMoL/interface.ccl index 4cf1fea82..9b9c4a769 100644 --- a/CarpetExtra/WaveToyMoL/interface.ccl +++ b/CarpetExtra/WaveToyMoL/interface.ccl @@ -1,5 +1,5 @@ # Interface definition for thorn WaveToyMoL -# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/WaveToyMoL/interface.ccl,v 1.1 2003/06/18 18:24:30 schnetter Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/WaveToyMoL/interface.ccl,v 1.2 2003/11/05 16:18:41 schnetter Exp $ IMPLEMENTS: WaveToyMoL @@ -7,7 +7,7 @@ INHERITS: boundary grid MethodOfLines -CCTK_INT FUNCTION Boundary_SelectGroupForBC(CCTK_POINTER IN GH, \ +CCTK_INT FUNCTION Boundary_SelectGroupForBC(CCTK_POINTER_TO_CONST IN GH, \ CCTK_INT IN faces, CCTK_INT IN boundary_width, CCTK_INT IN table_handle, \ CCTK_STRING IN group_name, CCTK_STRING IN bc_name) diff --git a/CarpetExtra/WaveToyMoL/schedule.ccl b/CarpetExtra/WaveToyMoL/schedule.ccl index 5404f562c..2b405ae15 100644 --- a/CarpetExtra/WaveToyMoL/schedule.ccl +++ b/CarpetExtra/WaveToyMoL/schedule.ccl @@ -1,5 +1,5 @@ # Schedule definitions for thorn WaveToyMoL -# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/WaveToyMoL/schedule.ccl,v 1.2 2003/07/20 21:03:43 schnetter Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/WaveToyMoL/schedule.ccl,v 1.3 2003/11/05 16:18:41 schnetter Exp $ STORAGE: scalarevolve[3] STORAGE: scalarevolvedot @@ -42,12 +42,12 @@ SCHEDULE GROUP ApplyBCs IN MoL_PostStep AFTER WaveToyMoL_Boundaries -SCHEDULE WaveToyMoL_Boundaries IN PostRestrict +SCHEDULE WaveToyMoL_Boundaries AT POSTRESTRICT { LANG: Fortran SYNC: scalarevolve } "Select boundary conditions after restricting" -SCHEDULE GROUP ApplyBCs IN PostRestrict AFTER WaveToyMoL_Boundaries +SCHEDULE GROUP ApplyBCs AT POSTRESTRICT AFTER WaveToyMoL_Boundaries { } "Apply boundary conditions after restricting" |