diff options
| author | schnetter <> | 2003-11-05 15:18:00 +0000 |
|---|---|---|
| committer | schnetter <> | 2003-11-05 15:18:00 +0000 |
| commit | f1bfa75fa4425d8cf55937fa26b33b2c6ba8f27a (patch) | |
| tree | db14971237fb532d89eff53dde64282084b73e6f | |
| parent | 1072312e58ca8b5d10bcb738a117a5fdf9cd6415 (diff) | |
Many changes that accumulated while Cactus and Carpet diverged.
Many changes that accumulated while Cactus and Carpet diverged.
Add processor splitting mechanism "along-dir" that splits along a
specified direction.
Rename group PostRestrict to bin POSTRESTRICT.
Prolongate initial data only when desired. This saves much time.
Sorry, Ian.
Fix bug in time level cycling of grid arrays. (Note: grid arrays
should not have time levels.)
Fix time_t bug on IRIX.
Make sure that there is no integer overflow when there are many
refinement levels. Always put parentheses around
(maxreflevelfact/reflevelfact).
Fix typo in Carpet verbose output.
Add debug output in processor splitting.
Communicate in three stages: Irecv, (work), Isend, Wait. This might
be more efficient. Much more, potentially.
Fix bug in processor layout of grid arrays. Sorry, Ian.
Make the interpolator interpolate between time levels. Untested.
Fix bug in processor communication in interpolator. Sorry, Ian.
Rewrite prolongation operators to make them twice as fast. There you
are, Ian.
Move prolongation operator kind handling from data to gdata.
Add official hyperslabbing interfaces to CarpetSlab.
Adapt to new cGH * handling.
darcs-hash:20031105151837-07bb3-758a87ff0355dba053269df4b7d7d79bea018669.gz
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" |