diff options
author | eschnett <> | 2001-03-10 19:54:00 +0000 |
---|---|---|
committer | eschnett <> | 2001-03-10 19:54:00 +0000 |
commit | 14c03fca5707f9b5954f84724e29d702955894a0 (patch) | |
tree | fe9a01ca3c2dea72092af7731244e25453acc607 | |
parent | 6b553d1f15dc4d2636ddb840313e0ff784ba29c3 (diff) |
Carpet now passes the Cart3d test suite. This required a rewrite of
Carpet now passes the Cart3d test suite. This required a rewrite of
CarpetSlabe, which now uses the data<> class instead of doing things
its own way. This in turn required some updates to data<>.
CarpetSlab now ignores (i. e. does not use) boundaries, except outer
boundaries. For that to work, the dh<> class has been given the
notion of an outer boundary.
In order to increase performance, the cctk_bbox[] information in
Carpet is not set correctly. As an approximation, it just reflects
whether the current grid component lies at the outer boundary, as
calculated using cctk_lbnd, cctk_ubnd, and cctk_gsh.
darcs-hash:20010310195459-f6438-8ee874d722de4df717e5eb3002722ceb8d4a9bb7.gz
-rw-r--r-- | Carpet/Carpet/param.ccl | 12 | ||||
-rw-r--r-- | Carpet/Carpet/src/carpet.cc | 60 | ||||
-rw-r--r-- | Carpet/Carpet/src/carpet.hh | 6 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/bbox.hh | 4 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/bboxset.cc | 4 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/bboxset.hh | 4 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/data.cc | 54 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/data.hh | 11 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/dh.cc | 40 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/dh.hh | 8 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/gdata.cc | 8 | ||||
-rw-r--r-- | Carpet/CarpetLib/src/gdata.hh | 17 | ||||
-rw-r--r-- | Carpet/CarpetSlab/src/carpetslab.cc | 678 | ||||
-rw-r--r-- | Carpet/CarpetSlab/src/carpetslab.h | 50 | ||||
-rw-r--r-- | Carpet/CarpetSlab/src/carpetslab.hh | 46 | ||||
-rw-r--r-- | CarpetAttic/Cart3dTest/src/Initial.c | 53 |
16 files changed, 504 insertions, 551 deletions
diff --git a/Carpet/Carpet/param.ccl b/Carpet/Carpet/param.ccl index 95ce595f1..13eb6b396 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.2 2001/03/07 12:59:40 eschnett Exp $ +# $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/param.ccl,v 1.3 2001/03/10 20:54:59 eschnett Exp $ shares: Cactus @@ -27,6 +27,11 @@ CCTK_INT global_nz "Grid size in z direction" 0:* :: "must be nonnegative" } 10 +CCTK_INT global_nsize "Grid size in each spatial direction" +{ + -1:* :: "must be nonnegative, or -1 to use the individual parameters" +} -1 + CCTK_INT ghost_size_x "Ghost zones in x direction" @@ -44,6 +49,11 @@ CCTK_INT ghost_size_z "Ghost zones in z direction" 0:* :: "must be nonnegative" } 1 +CCTK_INT ghost_size "Ghost zones in each spatial direction" +{ + -1:* :: "must be nonnegative, or -1 to use the individual parameters" +} -1 + CCTK_INT max_refinement_levels "Maximum number of refinement levels (including the base level)" diff --git a/Carpet/Carpet/src/carpet.cc b/Carpet/Carpet/src/carpet.cc index 7edb60644..9e1ae2cfb 100644 --- a/Carpet/Carpet/src/carpet.cc +++ b/Carpet/Carpet/src/carpet.cc @@ -1,4 +1,4 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Attic/carpet.cc,v 1.4 2001/03/07 12:59:53 eschnett Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Attic/carpet.cc,v 1.5 2001/03/10 20:55:03 eschnett 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 @@ -32,6 +32,8 @@ #include "carpet.hh" +static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Attic/carpet.cc,v 1.5 2001/03/10 20:55:03 eschnett Exp $"; + namespace Carpet { @@ -67,6 +69,9 @@ namespace Carpet { // data for grid functions vector<gfdesc> gfdata; // [group] + // active time level + int activetimelevel; // 0 for current, 1 for next + // current position on the grid hierarchy int mglevel; int reflevel; @@ -114,13 +119,24 @@ namespace Carpet { Checkpoint ("starting SetupGH..."); // ghost zones - const vect<int,dim> lghosts(ghost_size_x, ghost_size_y, ghost_size_z); - const vect<int,dim> ughosts(ghost_size_x, ghost_size_y, ghost_size_z); + vect<int,dim> lghosts, ughosts; + if (ghost_size == -1) { + lghosts = vect<int,dim>(ghost_size_x, ghost_size_y, ghost_size_z); + ughosts = vect<int,dim>(ghost_size_x, ghost_size_y, ghost_size_z); + } else { + lghosts = vect<int,dim>(ghost_size, ghost_size, ghost_size); + ughosts = vect<int,dim>(ghost_size, ghost_size, ghost_size); + } // grid size const int stride = (int)floor(pow(refinement_factor, max_refinement_levels) + 0.5); - const vect<int,dim> npoints(global_nx, global_ny, global_nz); + vect<int,dim> npoints; + if (global_nsize == -1) { + npoints = vect<int,dim>(global_nx, global_ny, global_nz); + } else { + npoints = vect<int,dim>(global_nsize, global_nsize, global_nsize); + } const vect<int,dim> str(stride); const vect<int,dim> lb(lghosts * str); @@ -207,6 +223,9 @@ namespace Carpet { } } + // active time level + activetimelevel = 0; + // current position reflevel = 0; mglevel = 0; @@ -342,6 +361,10 @@ namespace Carpet { Checkpoint ("%*sstarting RecursiveEvolve on level %d...", 2*reflevel, "", reflevel); + // Move activity to next time level + assert (activetimelevel == 0); + activetimelevel = 1; + // Prestep CCTK_ScheduleTraverse ("CCTK_PRESTEP", cgh, CallFunction); @@ -359,11 +382,15 @@ namespace Carpet { cgh->cctk_time += cgh->cctk_delta_time; } + // Restrict + Restrict (cgh); + // Cycle time levels CycleTimeLevels (cgh); - // Restrict - Restrict (cgh); + // Move activity back to current time level + assert (activetimelevel == 1); + activetimelevel = 0; // Poststep CCTK_ScheduleTraverse ("CCTK_POSTSTEP", cgh, CallFunction); @@ -391,8 +418,8 @@ namespace Carpet { RecursiveShutdown (cgh); - fprintf (stdout, "--------------------------------------------------------------------------------\n" - "Done.\n"); + CCTK_PRINTSEPARATOR; + printf ("Done.\n"); // earlier checkpoint before calling finalising MPI Checkpoint ("done with Shutdown."); @@ -575,11 +602,10 @@ namespace Carpet { cgh->cctk_gsh[d] = ((base.shape() / base.stride() + dd->lghosts + dd->ughosts)[d] * cgh->cctk_levfac[d]); - // TODO: set boundary information correctly - cgh->cctk_bbox[2*d ] = 1; - cgh->cctk_bbox[2*d+1] = 1; cgh->cctk_lbnd[d] = (ext.lower() / ext.stride())[d]; cgh->cctk_ubnd[d] = (ext.upper() / ext.stride())[d]; + cgh->cctk_bbox[2*d ] = cgh->cctk_lbnd[d] == 0; + cgh->cctk_bbox[2*d+1] = cgh->cctk_ubnd[d] == cgh->cctk_gsh[d]-1; for (int stg=0; stg<CCTK_NSTAGGER; ++stg) { // TODO: support staggering cgh->cctk_lssh[CCTK_LSSH_IDX(stg,d)] = cgh->cctk_lsh[d]; @@ -750,7 +776,8 @@ namespace Carpet { } const int n0 = CCTK_FirstVarIndexI(group); - const int tl = max(0, CCTK_NumTimeLevelsFromVarI(n0) - 1); + const int tl = max(0, (CCTK_NumTimeLevelsFromVarI(n0) + + activetimelevel - 2)); switch (CCTK_GroupTypeI(group)) { @@ -1001,6 +1028,7 @@ namespace Carpet { const int reflevels = max_refinement_levels; // arbitrary value const int mglevels = 1; // arbitrary value vector<vector<bbox<int,dim> > > bbss(reflevels); + // note: what this routine calls "ub" is "ub+str" elsewhere vect<int,dim> rstr = hh->baseextent.stride(); vect<int,dim> rlb = hh->baseextent.lower(); vect<int,dim> rub = hh->baseextent.upper() + rstr; @@ -1018,8 +1046,9 @@ namespace Carpet { vect<int,dim> cstr = rstr; vect<int,dim> clb = rlb; vect<int,dim> cub = rub; - const int zstep = ((rub[dim-1] - rlb[dim-1] + rstr[dim-1]) - / cstr[dim-1] + nprocs) / nprocs * cstr[dim-1]; + const int glonpz = (rub[dim-1] - rlb[dim-1]) / cstr[dim-1]; + const int locnpz = (glonpz + nprocs - 1) / nprocs; + const int zstep = locnpz * cstr[dim-1]; clb[dim-1] = rlb[dim-1] + zstep * c; cub[dim-1] = rlb[dim-1] + zstep * (c+1); if (c == nprocs-1) cub[dim-1] = rub[dim-1]; @@ -1103,7 +1132,8 @@ namespace Carpet { for (int group=0; group<CCTK_NumGroups(); ++group) { const int n0 = CCTK_FirstVarIndexI(group); - const int tl = max(0, CCTK_NumTimeLevelsFromVarI(n0) - 1); + const int tl = max(0, (CCTK_NumTimeLevelsFromVarI(n0) + + activetimelevel - 2)); switch (CCTK_GroupTypeI(group)) { diff --git a/Carpet/Carpet/src/carpet.hh b/Carpet/Carpet/src/carpet.hh index d94acdea4..a3cadb4c5 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.2 2001/03/05 14:30:03 eschnett Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/carpet.hh,v 1.3 2001/03/10 20:55:03 eschnett Exp $ #include <vector> @@ -47,6 +47,9 @@ namespace Carpet { }; extern vector<gfdesc> gfdata; // [group] + // active time level + extern int activetimelevel; // 0 for current, 1 for next + // current position on the grid hierarchy extern int mglevel; extern int reflevel; @@ -76,7 +79,6 @@ namespace Carpet { int EnableGroupComm (cGH *cgh, const char *groupname); int DisableGroupComm (cGH *cgh, const char *groupname); int Barrier (cGH *cgh); -// int ParallelInit (cGH *cgh); int Exit (cGH *cgh, int retval); int Abort (cGH *cgh, int retval); int myProc (cGH *cgh); diff --git a/Carpet/CarpetLib/src/bbox.hh b/Carpet/CarpetLib/src/bbox.hh index 6787b9240..95afb992e 100644 --- a/Carpet/CarpetLib/src/bbox.hh +++ b/Carpet/CarpetLib/src/bbox.hh @@ -5,7 +5,7 @@ copyright : (C) 2000 by Erik Schnetter email : schnetter@astro.psu.edu - $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bbox.hh,v 1.3 2001/03/07 13:00:57 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bbox.hh,v 1.4 2001/03/10 20:55:06 eschnett Exp $ ***************************************************************************/ @@ -96,6 +96,8 @@ public: bbox contracted_for (const bbox& b) const; // Set operations + // TODO: rename these; they clash with the bboxset operations + // (and name them & and | instead) // Smallest bbox containing both boxes bbox operator* (const bbox& b) const; bbox& operator*= (const bbox& b); diff --git a/Carpet/CarpetLib/src/bboxset.cc b/Carpet/CarpetLib/src/bboxset.cc index 1dd685178..b76dc67b6 100644 --- a/Carpet/CarpetLib/src/bboxset.cc +++ b/Carpet/CarpetLib/src/bboxset.cc @@ -5,7 +5,7 @@ copyright : (C) 2000 by Erik Schnetter email : schnetter@astro.psu.edu - $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bboxset.cc,v 1.1 2001/03/01 13:40:10 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bboxset.cc,v 1.2 2001/03/10 20:55:06 eschnett Exp $ ***************************************************************************/ @@ -38,7 +38,7 @@ bboxset<T,D>::bboxset () { template<class T, int D> bboxset<T,D>::bboxset (const box& b) { - bs.insert(b); + if (!b.empty()) bs.insert(b); assert (invariant()); } diff --git a/Carpet/CarpetLib/src/bboxset.hh b/Carpet/CarpetLib/src/bboxset.hh index 009d9afdd..98aa3fb65 100644 --- a/Carpet/CarpetLib/src/bboxset.hh +++ b/Carpet/CarpetLib/src/bboxset.hh @@ -5,7 +5,7 @@ copyright : (C) 2000 by Erik Schnetter email : schnetter@astro.psu.edu - $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bboxset.hh,v 1.1 2001/03/01 13:40:10 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bboxset.hh,v 1.2 2001/03/10 20:55:06 eschnett Exp $ ***************************************************************************/ @@ -34,6 +34,8 @@ // Forward definition template<class T, int D> class bboxset; +// TODO: add more functions for the other operators +// (but cf. bbox<T,D> first!) template<class T,int D> bboxset<T,D> operator- (const bbox<T,D>& b1, const bbox<T,D>& b2); template<class T,int D> diff --git a/Carpet/CarpetLib/src/data.cc b/Carpet/CarpetLib/src/data.cc index 3f8222f6c..b2c190865 100644 --- a/Carpet/CarpetLib/src/data.cc +++ b/Carpet/CarpetLib/src/data.cc @@ -5,7 +5,7 @@ copyright : (C) 2000 by Erik Schnetter email : schnetter@astro.psu.edu - $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.cc,v 1.3 2001/03/05 21:48:38 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.cc,v 1.4 2001/03/10 20:55:06 eschnett Exp $ ***************************************************************************/ @@ -56,13 +56,14 @@ data<T,D>::~data () { // Pseudo constructors template<class T, int D> -data<T,D>* data<T,D>::make_typed (const ibbox& extent, const int proc) const { - return new data(extent, proc); +data<T,D>* data<T,D>::make_typed () const { + return new data(); } // Storage management template<class T, int D> -void data<T,D>::allocate (const ibbox& extent, const int proc) { +void data<T,D>::allocate (const ibbox& extent, const int proc, + void* const mem=0) { assert (!_has_storage); _has_storage = true; // data @@ -77,13 +78,20 @@ void data<T,D>::allocate (const ibbox& extent, const int proc) { int rank; MPI_Comm_rank (dist::comm, &rank); if (rank==_proc) { - _storage = new T[_size]; + _owns_storage = !mem; + if (!mem) { + _storage = new T[_size]; + } else { + _storage = (T*)mem; + } + } else { + assert (!mem); } } template<class T, int D> void data<T,D>::free () { - if (_storage) delete [] _storage; + if (_storage && _owns_storage) delete [] _storage; _storage = 0; _has_storage = false; } @@ -99,9 +107,12 @@ void data<T,D>::transfer_from (generic_data<D>* gsrc) { // Processor management template<class T, int D> -void data<T,D>::change_processor (const int newproc) { - if (newproc == _proc) return; - +void data<T,D>::change_processor (const int newproc, void* const mem=0) { + if (newproc == _proc) { + assert (!mem); + return; + } + if (_has_storage) { int rank; MPI_Comm_rank (dist::comm, &rank); @@ -109,7 +120,12 @@ void data<T,D>::change_processor (const int newproc) { // copy from other processor assert (!_storage); - _storage = new T[_size]; + _owns_storage = !mem; + if (!mem) { + _storage = new T[_size]; + } else { + _storage = (T*)mem; + } T dummy; MPI_Status status; @@ -118,16 +134,20 @@ void data<T,D>::change_processor (const int newproc) { } else if (rank == _proc) { // copy to other processor - + + assert (!mem); assert (_storage); T dummy; MPI_Send (_storage, _size, dist::datatype(dummy), newproc, dist::tag, dist::comm); - delete [] _storage; + if (_owns_storage) { + delete [] _storage; + } _storage = 0; - + } else { + assert (!mem); assert (!_storage); } } @@ -166,7 +186,7 @@ void data<T,D>::copy_from (const generic_data<D>* gsrc, const ibbox& box) { } else { // copy to different processor - data* tmp = new data(box, src->_proc); + data* const tmp = new data(box, src->_proc); tmp->copy_from (src, box); tmp->change_processor (_proc); copy_from (tmp, box); @@ -228,7 +248,7 @@ void data<T,D>::interpolate_from (const generic_data<D>* gsrc, } else { // interpolate from other processor - data* tmp = new data(box, src->_proc); + data* const tmp = new data(box, src->_proc); tmp->interpolate_from (src, box); tmp->change_processor (_proc); copy_from (tmp, box); @@ -300,9 +320,9 @@ void data<T,D>::interpolate_from (const generic_data<D>* gsrc, } else { // interpolate from other processors - data* smp = new data(box, src->_proc); + data* const smp = new data(box, src->_proc); smp->copy_from (src, box); - data* tmp = new data(box, trc->_proc); + data* const tmp = new data(box, trc->_proc); tmp->copy_from (trc, box); interpolate_from (smp, sfact, tmp, tfact, box); delete smp; diff --git a/Carpet/CarpetLib/src/data.hh b/Carpet/CarpetLib/src/data.hh index 7b4336b97..84e4c92b7 100644 --- a/Carpet/CarpetLib/src/data.hh +++ b/Carpet/CarpetLib/src/data.hh @@ -5,7 +5,7 @@ copyright : (C) 2000 by Erik Schnetter email : schnetter@astro.psu.edu - $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.hh,v 1.2 2001/03/05 14:31:03 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.hh,v 1.3 2001/03/10 20:55:06 eschnett Exp $ ***************************************************************************/ @@ -27,7 +27,6 @@ #include "defs.hh" #include "dist.hh" #include "bbox.hh" -#include "bboxset.hh" #include "gdata.hh" #include "vect.hh" @@ -49,7 +48,6 @@ class data: public generic_data<D> { // Types typedef vect<int,D> ivect; typedef bbox<int,D> ibbox; - typedef bboxset<int,D> ibset; // Fields T* restrict _storage; // the data (if located on this processor) @@ -64,15 +62,16 @@ public: virtual ~data (); // Pseudo constructors - virtual data* make_typed (const ibbox& extent, const int proc) const; + virtual data* make_typed () const; // Storage management - virtual void allocate (const ibbox& extent, const int proc); + virtual void allocate (const ibbox& extent, const int proc, + void* const mem=0); virtual void free (); virtual void transfer_from (generic_data<D>* gsrc); // Processor management - virtual void change_processor (const int newproc); + virtual void change_processor (const int newproc, void* const mem=0); // Accessors virtual const T* storage () const { diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc index 5dc7ba0e7..756f37cad 100644 --- a/Carpet/CarpetLib/src/dh.cc +++ b/Carpet/CarpetLib/src/dh.cc @@ -6,7 +6,7 @@ copyright : (C) 2000 by Erik Schnetter email : schnetter@astro.psu.edu - $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.cc,v 1.3 2001/03/07 13:00:57 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.cc,v 1.4 2001/03/10 20:55:06 eschnett Exp $ ***************************************************************************/ @@ -66,18 +66,18 @@ void dh<D>::recompose () { boxes[rl][c].resize(h.mglevels(rl,c)); for (int ml=0; ml<h.mglevels(rl,c); ++ml) { const ibbox intr = h.extents[rl][c][ml]; - + // Interior // (the interior of the grid has the extent as specified by // the user) boxes[rl][c][ml].interior = intr; - + // Exterior (add ghost zones) // (the content of the exterior is completely determined by // the interior of this or other components; the content of // the exterior is redundant) boxes[rl][c][ml].exterior = intr.expand(lghosts, ughosts); - + // Boundaries (ghost zones only) // (interior + boundaries = exterior) boxes[rl][c][ml].boundaries = boxes[rl][c][ml].exterior - intr; @@ -118,8 +118,9 @@ void dh<D>::recompose () { for (int ml=0; ml<h.mglevels(rl,c); ++ml) { const ibbox intr = boxes[rl][c][ml].interior; const ibbox extr = boxes[rl][c][ml].exterior; + const ibset bnds = boxes[rl][c][ml].boundaries; - + // Sync boxes for (int cc=0; cc<h.components(rl); ++cc) { assert (ml<h.mglevels(rl,cc)); @@ -208,6 +209,31 @@ void dh<D>::recompose () { } // for c } // for rl + for (int rl=0; rl<h.reflevels(); ++rl) { + for (int c=0; c<h.components(rl); ++c) { + for (int ml=0; ml<h.mglevels(rl,c); ++ml) { + + // Boundaries that are neither synced nor prolonged from + // coarser grids (outer boundaries) + ibset& sync_not = boxes[rl][c][ml].sync_not; + + // The whole boundary + sync_not = boxes[rl][c][ml].boundaries; + + // Subtract boxes received during synchronisation + const iblistvect& recv_sync = boxes[rl][c][ml].recv_sync; + for (iblistvect::const_iterator lvi=recv_sync.begin(); + lvi!=recv_sync.end(); ++lvi) { + for (iblist::const_iterator li=lvi->begin(); + li!=lvi->end(); ++li) { + sync_not -= *li; + } + } + + } // for ml + } // for c + } // for rl + bases.resize(h.reflevels()); for (int rl=0; rl<h.reflevels(); ++rl) { if (h.components(rl)==0) { @@ -221,6 +247,8 @@ void dh<D>::recompose () { bases[rl][ml].exterior *= boxes[rl][c][ml].exterior; bases[rl][ml].interior *= boxes[rl][c][ml].interior; } + bases[rl][ml].boundaries + = bases[rl][ml].exterior - bases[rl][ml].interior; } } } @@ -247,6 +275,7 @@ void dh<D>::recompose () { cout << "boundaries=" << boxes[rl][c][ml].boundaries << endl; cout << "recv_sync=" << boxes[rl][c][ml].recv_sync << endl; cout << "recv_ref_bnd_coarse=" << boxes[rl][c][ml].recv_ref_bnd_coarse << endl; + cout << "sync_not=" << boxes[rl][c][ml].sync_not << endl; } } } @@ -258,6 +287,7 @@ void dh<D>::recompose () { cout << "rl=" << rl << " ml=" << ml << endl; cout << "exterior=" << bases[rl][ml].exterior << endl; cout << "interior=" << bases[rl][ml].interior << endl; + cout << "boundaries=" << bases[rl][ml].boundaries << endl; } } } diff --git a/Carpet/CarpetLib/src/dh.hh b/Carpet/CarpetLib/src/dh.hh index 8960e60d1..91634d1ca 100644 --- a/Carpet/CarpetLib/src/dh.hh +++ b/Carpet/CarpetLib/src/dh.hh @@ -6,7 +6,7 @@ copyright : (C) 2000 by Erik Schnetter email : schnetter@astro.psu.edu - $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.hh,v 1.2 2001/03/07 13:00:57 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.hh,v 1.3 2001/03/10 20:55:06 eschnett Exp $ ***************************************************************************/ @@ -73,13 +73,15 @@ class dh { iblistvect send_ref_bnd_fine; ibset boundaries; // boundaries - iblistvect recv_sync; // receive while syncing - iblistvect recv_ref_bnd_coarse; + iblistvect recv_sync; // received while syncing + iblistvect recv_ref_bnd_coarse; // received from coarser grids + ibset sync_not; // not received while syncing (outer boundary) }; struct dbases { ibbox exterior; // whole region (including boundaries) ibbox interior; // interior (without boundaries) + ibset boundaries; // boundaries }; typedef vector<dboxes> mboxes; // ... for each multigrid level diff --git a/Carpet/CarpetLib/src/gdata.cc b/Carpet/CarpetLib/src/gdata.cc index 36889de97..288eb9aa6 100644 --- a/Carpet/CarpetLib/src/gdata.cc +++ b/Carpet/CarpetLib/src/gdata.cc @@ -5,7 +5,7 @@ copyright : (C) 2000 by Erik Schnetter email : schnetter@astro.psu.edu - $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.cc,v 1.4 2001/03/07 13:00:57 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.cc,v 1.5 2001/03/10 20:55:06 eschnett Exp $ ***************************************************************************/ @@ -20,6 +20,7 @@ #include <cassert> #include <fstream> +#include <iomanip> #include "bbox.hh" #include "defs.hh" @@ -67,6 +68,8 @@ void generic_data<D>::write_ascii (const string name, const int time, ofstream file(name.c_str(), ios::app); assert (file.good()); + file << setprecision(15); + file << "#" << endl << "# iteration " << time << endl << "# time level " << tl << " refinement level " << rl @@ -111,7 +114,8 @@ void generic_data<D>::write_ascii (const string name, const int time, } else { // copy to processor 0 and output there - generic_data* tmp = make_typed(_extent, 0); + generic_data* const tmp = make_typed(); + tmp->allocate(_extent, 0); tmp->copy_from (this, _extent); tmp->write_ascii (name, time, org, dirs, tl, rl, c, ml); delete tmp; diff --git a/Carpet/CarpetLib/src/gdata.hh b/Carpet/CarpetLib/src/gdata.hh index 1d245be7c..4e77071d1 100644 --- a/Carpet/CarpetLib/src/gdata.hh +++ b/Carpet/CarpetLib/src/gdata.hh @@ -5,7 +5,7 @@ copyright : (C) 2000 by Erik Schnetter email : schnetter@astro.psu.edu - $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.hh,v 1.2 2001/03/05 14:31:03 eschnett Exp $ + $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.hh,v 1.3 2001/03/10 20:55:06 eschnett Exp $ ***************************************************************************/ @@ -30,7 +30,6 @@ #include "defs.hh" #include "dist.hh" #include "bbox.hh" -#include "bboxset.hh" #include "vect.hh" @@ -51,12 +50,12 @@ class generic_data { // Types typedef vect<int,D> ivect; typedef bbox<int,D> ibbox; - typedef bboxset<int,D> ibset; protected: // should be readonly // Fields bool _has_storage; // has storage associated (on some processor) + bool _owns_storage; // owns the storage ivect _shape, _stride; // shape and index order int _size; // size @@ -73,21 +72,25 @@ public: virtual ~generic_data (); // Pseudo constructors - virtual generic_data* make_typed (const ibbox& extent, const int proc) const - = 0; + virtual generic_data* make_typed () const = 0; // Storage management - virtual void allocate (const ibbox& extent, const int proc) = 0; + virtual void allocate (const ibbox& extent, const int proc, + void* const mem=0) = 0; virtual void free () = 0; virtual void transfer_from (generic_data* src) = 0; // Processor management - virtual void change_processor (const int newproc) = 0; + virtual void change_processor (const int newproc, void* const mem=0) = 0; // Accessors bool has_storage () const { return _has_storage; } + bool owns_storage () const { + assert (_has_storage); + return _owns_storage; + } virtual const void* storage () const = 0; diff --git a/Carpet/CarpetSlab/src/carpetslab.cc b/Carpet/CarpetSlab/src/carpetslab.cc index 15cb8a874..599b31bf3 100644 --- a/Carpet/CarpetSlab/src/carpetslab.cc +++ b/Carpet/CarpetSlab/src/carpetslab.cc @@ -1,25 +1,24 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/Attic/carpetslab.cc,v 1.3 2001/03/07 13:01:11 eschnett Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/Attic/carpetslab.cc,v 1.4 2001/03/10 20:55:09 eschnett Exp $ #include <cassert> #include <cstdlib> -#include <cstring> - -#include <mpi.h> #include "cctk.h" #include "Carpet/CarpetLib/src/bbox.hh" -#include "Carpet/CarpetLib/src/gdata.hh" +#include "Carpet/CarpetLib/src/bboxset.hh" #include "Carpet/CarpetLib/src/dh.hh" -#include "Carpet/CarpetLib/src/dist.hh" -#include "Carpet/CarpetLib/src/ggf.hh" +#include "Carpet/CarpetLib/src/gdata.hh" #include "Carpet/CarpetLib/src/gh.hh" +#include "Carpet/CarpetLib/src/ggf.hh" #include "Carpet/CarpetLib/src/vect.hh" #include "Carpet/Carpet/src/carpet.hh" #include "carpetslab.hh" +static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/Attic/carpetslab.cc,v 1.4 2001/03/10 20:55:09 eschnett Exp $"; + namespace CarpetSlab { @@ -28,81 +27,85 @@ namespace CarpetSlab { - int Hyperslab_GetLocalHyperslab (cGH* cgh, - const int n, - const int tl, - const int hdim, - const int origin [/*vdim*/], - const int dir [/*vdim*/], - const int len [/*hdim*/], - const int downsample [/*hdim*/], - void** const hdata, - int hsize [/*hdim*/], - int ghsize [/*hdim*/], - int hoffset [/*hdim*/]) + void* GetSlab (cGH* const cgh, + const int dest_proc, + const int n, + const int tl, + const int hdim, + const int origin[/*vdim*/], + const int dirs[/*hdim*/], + const int stride[/*hdim*/], + const int length[/*hdim*/]) { - // check current status - assert (mglevel>=0); + // Check global state assert (reflevel>=0); - assert (component>=0); // local := this component + assert (mglevel>=0); + + // Save global state + int saved_component = component; + component = -1; - // check arguments + // Check Cactus grid hierarchy + assert (cgh); + + // Check destination processor + assert (dest_proc>=-1 && dest_proc<CCTK_nProcs(cgh)); + + // Check dimension + assert (hdim>=0 && hdim<=dim); + + // Check variable index assert (n>=0 && n<CCTK_NumVars()); - assert (tl>=0); - assert (hdim>0 && hdim<=dim); - // the following assertion is too strict (better allow arbitrary values) - { - const int group = CCTK_GroupIndexFromVarI(n); - assert (group>=0); - switch (CCTK_GroupTypeFromVarI(n)) { - case CCTK_SCALAR: - abort(); - case CCTK_ARRAY: { - assert (group<(int)arrdata.size()); -// assert (arrdata[group].hh->is_local(reflevel, component)); -// const bbox<int,dim> ext = -// arrdata[group].dd->boxes[reflevel][component][mglevel].exterior; -// for (int d=0; d<dim; ++d) { -// assert (origin[d] >= ext.lower()[d] / ext.stride()[d] -// && origin[d] <= ext.upper()[d] / ext.stride()[d]); -// } - break; - } - case CCTK_GF: { - assert (group<(int)gfdata.size()); -// assert (hh->is_local(reflevel, component)); -// const bbox<int,dim> ext = -// dd->boxes[reflevel][component][mglevel].exterior; -// for (int d=0; d<dim; ++d) { -// assert (origin[d] >= ext.lower()[d] / ext.stride()[d] -// && origin[d] <= ext.upper()[d] / ext.stride()[d]); -// } - break; - } - default: - abort(); - } - } - // the following assertion is too strict (better allow arbitrary values) - for (int d=0; d<dim; ++d) assert (dir[d]==0 || dir[d]==1); - // the following assertion is too strict (better allow arbitrary values) - for (int d=0; d<hdim; ++d) assert (downsample[d]>0); - assert (hdata); - assert (hsize); - assert (ghsize); - assert (hoffset); - // get variable info + // Get info about variable const int group = CCTK_GroupIndexFromVarI(n); assert (group>=0); - const int var = n - CCTK_FirstVarIndexI(group); + const int n0 = CCTK_FirstVarIndexI(group); + assert (n0>=0); + const int var = n - n0; assert (var>=0); - // get grid hierarchy data hierarchy and grid function - gh<dim>* myhh; - dh<dim>* mydd; - generic_gf<dim>* myff; - switch (CCTK_GroupTypeFromVarI(n)) { + // 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); + + // Check timelevel + assert (tl>=0 && tl<gp.numtimelevels); + + // 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 generic_gf<dim>* myff; + switch (gp.grouptype) { case CCTK_SCALAR: abort(); case CCTK_ARRAY: @@ -114,7 +117,7 @@ namespace CarpetSlab { break; case CCTK_GF: myhh = hh; - mydd = dd; + mydd = Carpet::dd; // necessary for pre-ANSI C++ compilers assert (group < (int)gfdata.size()); assert (var < (int)gfdata[group].data.size()); myff = gfdata[group].data[var]; @@ -122,388 +125,243 @@ namespace CarpetSlab { default: abort(); } + assert (myhh); + assert (mydd); + assert (myff); - // get data - const generic_data<dim>* mydata - = (*myff)(tl, reflevel, component, mglevel); - - // get local and global bounding boxes - assert (reflevel < (int)mydd->boxes.size()); - assert (component < (int)mydd->boxes[reflevel].size()); - assert (mglevel < (int)mydd->boxes[reflevel][component].size()); - const bbox<int,dim> locbox - = mydd->boxes[reflevel][component][mglevel].exterior; - const bbox<int,dim> globox - = mydd->bases[reflevel][mglevel].exterior; - - // calculate more convenient representation of the direction - vect<int,dim> stride[hdim]; - // the following if statement is written according to the - // definition of "dir". - if (hdim==1) { - for (int d=0; d<dim; ++d) stride[0][d] = dir[d] * downsample[0]; - } else if (hdim==dim) { - for (int dd=0; dd<hdim; ++dd) { - for (int d=0; d<dim; ++d) stride[dd][d] = d==dd ? downsample[dd] : 0; - } - } else if (hdim==2) { - assert (dim==3); - if (dir[0]==0) { - assert (dir[1]!=0 && dir[2]!=0); - stride[0] = vect<int,dim>::dir(1); - stride[1] = vect<int,dim>::dir(2); - } else if (dir[1]==0) { - assert (dir[0]!=0 && dir[2]!=0); - stride[0] = vect<int,dim>::dir(0); - stride[1] = vect<int,dim>::dir(2); - } else if (dir[2]==0) { - assert (dir[0]!=0 && dir[1]!=0); - stride[0] = vect<int,dim>::dir(0); - stride[1] = vect<int,dim>::dir(1); - } else { - abort(); - } - for (int dd=0; dd<hdim; ++dd) stride[dd] *= downsample[dd]; - } else { - abort(); - } - for (int dd=0; dd<hdim; ++dd) stride[dd] *= locbox.stride(); - - // local lower bound - vect<int,dim> lbound; - for (int d=0; d<dim; ++d) lbound[d] = origin[d]; - lbound *= locbox.stride(); - - // local upper bound - vect<int,dim> ubound = lbound; - for (int dd=0; dd<hdim; ++dd) { - if (len[dd]<0) { - assert (any(stride[dd]>0)); - while (all(ubound + stride[dd] <= locbox.upper())) { - ubound += stride[dd]; - } - } else { - ubound += stride[dd] * (len[dd]-1); - } - } - - lbound = max(lbound, locbox.lower()); - ubound = min(ubound, locbox.upper()); - - // local size - int total_hsize = 1; - for (int dd=0; dd<hdim; ++dd) { - hsize[dd] = 0; - assert (any(stride[dd]>0)); - while (all(lbound + stride[dd] * hsize[dd] <= ubound)) ++hsize[dd]; - total_hsize *= hsize[dd]; - } - assert (total_hsize>=0); - - // sanity check - if (total_hsize>0) { - vect<int,dim> tmp = lbound; - for (int dd=0; dd<hdim; ++dd) { - tmp += stride[dd] * (hsize[dd]-1); - } - assert (all(tmp == ubound)); - } + // Detemine collecting processor + const int collect_proc = dest_proc<0 ? 0 : dest_proc; - // global lower bound - vect<int,dim> glbound; - for (int d=0; d<dim; ++d) glbound[d] = origin[d]; - glbound *= globox.stride(); + // Determine own rank + const int rank = CCTK_MyProc(cgh); - // global upper bound - vect<int,dim> gubound = glbound; + // Calculate global size + int totalsize = 1; for (int dd=0; dd<hdim; ++dd) { - if (len[dd]<0) { - assert (any(stride[dd]>0)); - while (all(gubound + stride[dd] <= globox.upper())) { - gubound += stride[dd]; - } - } else { - gubound += stride[dd] * (len[dd]-1); - } + totalsize *= length[dd]; } - glbound = max(glbound, globox.lower()); - gubound = min(gubound, globox.upper()); - - // global size - int total_ghsize = 1; - for (int dd=0; dd<hdim; ++dd) { - ghsize[dd] = 0; - assert (any(stride[dd]>0)); - while (all(glbound + stride[dd] * ghsize[dd] <= gubound)) { - ++ghsize[dd]; - } - total_ghsize *= ghsize[dd]; + // Allocate memory + void* hdata = 0; + if (dest_proc==-1 || rank==dest_proc) { + hdata = malloc(totalsize * typesize); + assert (hdata); } - assert (total_ghsize>=0); - // sanity check - if (total_ghsize>0) { - vect<int,dim> tmp = glbound; + if (hh->components(reflevel) > 0) { + + // Only temporarily + component = 0; + + // Get sample data + const generic_data<dim>* mydata; + mydata = (*myff)(tl, reflevel, component, mglevel); + + // 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) { - tmp += stride[dd] * (ghsize[dd]-1); - } - assert (all(tmp == gubound)); - } - - // local to global offset - for (int dd=0; dd<hdim; ++dd) { - hoffset[dd] = 0; - assert (any(stride[dd]>0)); - while (all(glbound + stride[dd] * (hoffset[dd]+1) <= lbound)) { - ++hoffset[dd]; + hstr[dirs[dd]-1] *= stride[dd]; } - } - - // sanity check - { - vect<int,dim> tmp = glbound; - for (int dd=0; dd<hdim; ++dd) { - tmp += stride[dd] * hoffset[dd]; + + // Lower bound of collected data + vect<int,dim> hlb; + for (int d=0; d<dim; ++d) { + hlb[d] = origin[d] * str[d]; } - assert (all(tmp == lbound)); - } - - // bail out if this component is on another processor - if (! myhh->is_local(reflevel, component)) { - *hdata = 0; + + // Upper bound of collected data + vect<int,dim> hub = hlb; for (int dd=0; dd<hdim; ++dd) { - hsize[dd] = 0; - hoffset[dd] = 0; + hub[dirs[dd]-1] += (length[dd]-1) * hstr[dirs[dd]-1]; } - return -1; - } - - // allocate the memory - *hdata = malloc(total_hsize * CCTK_VarTypeSize(CCTK_VarTypeI(n))); - assert (*hdata); - - if (total_hsize>0) { - // copy the data to user memory - char* const dest = (char*)*hdata; - const char* const src = (const char*)mydata->storage(); - const int sz = CCTK_VarTypeSize(CCTK_VarTypeI(n)); + // Calculate extent to collect + const bbox<int,dim> hextent (hlb, hub, hstr); + assert (hextent.num_points() == totalsize); - int dest_index[hdim]; - for (int dd=0; dd<hdim; ++dd) dest_index[dd] = 0; - for (;;) { - - vect<int,dim> src_index = lbound; - for (int dd=0; dd<hdim; ++dd) src_index += stride[dd] * dest_index[dd]; + // Create collector data object + void* myhdata = rank==collect_proc ? hdata : 0; + generic_data<dim>* const alldata = mydata->make_typed(); + alldata->allocate (hextent, collect_proc, myhdata); + + // Done with the temporary stuff + mydata = 0; + component = -1; + + // Loop over all components, copying data from them + assert (component == -1); + for (component=0; component<hh->components(reflevel); ++component) { - int di = 0; - for (int dd=hdim-1; dd>=0; --dd) di = di * hsize[dd] + dest_index[dd]; + // Get data object + mydata = (*myff)(tl, reflevel, component, mglevel); - const int si = mydata->offset(src_index); + // Calculate overlapping extents + const bboxset<int,dim> myextents + = ((mydd->boxes[reflevel][component][mglevel].sync_not + | mydd->boxes[reflevel][component][mglevel].interior) + & hextent); - memcpy(dest + sz*di, src + sz*si, sz); + // 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 (mydata, *ext_iter); + + } - for (int dd=0; dd<hdim; ++dd) { - ++dest_index[dd]; - if (dest_index[dd]<hsize[dd]) break; - dest_index[dd]=0; - if (dd==hdim-1) goto done; + } // Loop over components + component = -1; + + // Copy result to all processors + if (dest_proc == -1) { + for (int proc=0; proc<CCTK_nProcs(cgh); ++proc) { + if (proc != collect_proc) { + + void* myhdata = rank==proc ? hdata : 0; + generic_data<dim>* const tmpdata = mydata->make_typed(); + tmpdata->allocate (alldata->extent(), proc, myhdata); + tmpdata->copy_from (alldata, alldata->extent()); + delete tmpdata; + + } } - } - done: ; + } // Copy result + + delete alldata; - } // if total_hsize>0 + } // if components>0 - return 0; + // Restore global state + component = saved_component; + + // Success + return hdata; } - int Hyperslab_GetHyperslab (cGH* cgh, + int Hyperslab_GetHyperslab (cGH* const GH, const int target_proc, - const int n, - const int tl, + const int vindex, + const int vtimelvl, const int hdim, - const int origin [/*vdim*/], - const int dir [/*vdim*/], - const int len [/*hdim*/], + const int global_startpoint [/*vdim*/], + const int directions [/*vdim*/], + const int lengths [/*hdim*/], const int downsample [/*hdim*/], void** const hdata, int hsize [/*hdim*/]) { - // check current status - assert (mglevel>=0); - assert (reflevel>=0); - - // in order to work, this requires that all processors have the - // same number of components - const int saved_component = component; - component = -1; + // Check some arguments + assert (hdim>=0 && hdim<=dim); - // check arguments - assert (n>=0 && n<CCTK_NumVars()); - assert (tl>=0); - assert (hdim>0 && hdim<=dim); - // the following assertion is too strict (better allow arbitrary values) - { - // TODO: make sure that origin is within the extent of this - // refinement / multigrid level - // (but no such extent is stored in dh) - // (it is now; use it!) - const int group = CCTK_GroupIndexFromVarI(n); - assert (group>=0); - switch (CCTK_GroupTypeFromVarI(n)) { - case CCTK_SCALAR: - abort(); - case CCTK_ARRAY: - assert (group<(int)arrdata.size()); - break; - case CCTK_GF: - assert (group<(int)gfdata.size()); - break; - default: - abort(); - } - } - // the following assertion is too strict (better allow arbitrary values) - for (int d=0; d<dim; ++d) assert (dir[d]==0 || dir[d]==1); - // the following assertion is too strict (better allow arbitrary values) - for (int d=0; d<hdim; ++d) assert (downsample[d]>0); + // Check output arguments assert (hdata); assert (hsize); - int collect_proc = target_proc; - if (collect_proc<0) collect_proc = 0; - - assert (hh->components(reflevel)>0); - *hdata = 0; - for (int dd=0; dd<hdim; ++dd) hsize[dd] = 0; - int totalhsize = 0; - - const int sz = CCTK_VarTypeSize(CCTK_VarTypeI(n)); - - MPI_Datatype type; - switch (CCTK_VarTypeI(n)) { -#define TYPECASE(N,T) \ - case N: { \ - assert (sz == sizeof(T)); \ - T dummy; \ - type = dist::datatype(dummy); \ - break; \ - } -#include "Carpet/Carpet/src/typecase" -#undef TYPECASE - default: - abort(); - } - - int rank; - MPI_Comm_rank (dist::comm, &rank); - - // loop over all components - for (component=0; component<hh->components(reflevel); ++component) { - - void* myhdata; - int myhsize[hdim], ghsize[hdim], hoffset[hdim]; - - const int retval = Hyperslab_GetLocalHyperslab - (cgh, n, tl, hdim, origin, dir, len, downsample, - &myhdata, myhsize, ghsize, hoffset); - - if (hh->is_local(reflevel,component)) { - assert (retval == 0); - } else { - assert (retval == -1); + // Calculate more convenient representation of the direction + int 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<dim; ++d) { + if (directions[d]!=0) { + mydir = d+1; + break; + } } - - int mytotalsize = 1; - for (int dd=0; dd<hdim; ++dd) mytotalsize *= myhsize[dd]; - - if (component==0) { - if (target_proc<0 || rank == target_proc) { - - for (int dd=0; dd<hdim; ++dd) hsize[dd] = ghsize[dd]; - - totalhsize = 1; - for (int dd=0; dd<hdim; ++dd) totalhsize *= hsize[dd]; - - if (rank == collect_proc) { - *hdata = malloc(totalhsize * sz); - assert (*hdata); - } else { - *hdata = 0; - } - + assert (mydir>0); + for (int d=0; d<dim; ++d) { + if (d == mydir-1) { + assert (directions[d]!=0); + } else { + assert (directions[d]==0); } } - - if (!myhdata && rank == collect_proc) { - MPI_Status status; - assert (mytotalsize==0); - const int src = hh->proc(reflevel, component); - MPI_Recv (&mytotalsize, 1, MPI_INT, src, 2001, - dist::comm, &status); - myhdata = malloc(mytotalsize * sz); - assert (myhdata); - MPI_Recv (myhdata, mytotalsize, type, src, 2001, - dist::comm, &status); - } else if (myhdata && rank != collect_proc) { - MPI_Send (&mytotalsize, 1, MPI_INT, collect_proc, 2001, dist::comm); - MPI_Send (myhdata, mytotalsize, type, collect_proc, 2001, dist::comm); - free (myhdata); - myhdata = 0; - mytotalsize = 0; + dirs[0] = mydir; + } else if (hdim==dim) { + // dim-dimensional hyperslab + for (int dd=0; dd<hdim; ++dd) { + dirs[dd] = dd+1; } - - if (myhdata>0) { - assert (rank == collect_proc); - - if (mytotalsize>0) { - int dest_index[hdim], src_index[hdim]; - for (int dd=0; dd<hdim; ++dd) dest_index[dd] = hoffset[dd]; - for (int dd=0; dd<hdim; ++dd) src_index[dd] = 0; - for (;;) { - int di=0; - for (int dd=hdim-1; dd>=0; --dd) { - di = di * hsize[dd] + dest_index[dd]; - } - int si=0; - for (int dd=hdim-1; dd>=0; --dd) { - si = si * myhsize[dd] + src_index[dd]; - } - - memcpy ((char*)*hdata + sz*di, (char*)myhdata + sz*si, sz); - - for (int dd=0; dd<hdim; ++dd) { - ++dest_index[dd]; - ++src_index[dd]; - if (src_index[dd] < myhsize[dd]) break; - dest_index[dd] = hoffset[dd]; - src_index[dd] = 0; - if (dd==hdim-1) goto done; - } - } - done: ; - - } // if mytotalsize>0 - - free (myhdata); - myhdata = 0; - } else { - assert (rank != collect_proc); + } else if (hdim==2) { + // 2-dimensional hyperslab with dim==3 + assert (dim==3); + int mydir = 0; + for (int d=0; d<dim; ++d) { + if (directions[d]==0) { + mydir = d+1; + break; + } } - + assert (mydir>0); + for (int d=0; d<dim; ++d) { + if (d == mydir-1) { + assert (directions[d]==0); + } else { + assert (directions[d]!=0); + } + } + int dd=0; + for (int d=0; d<dim; ++d) { + if (d != mydir-1) { + dirs[dd] = d+1; + ++dd; + } + } + assert (dd==hdim); + } else { + abort(); } - if (target_proc<0) { - MPI_Bcast (hdata, totalhsize*sz, MPI_BYTE, - collect_proc, CarpetMPICommunicator()); +#if 0 + // Invert directions + int invdir[dim]; + for (int d=0; d<dim; ++d) { + invdir[d] = 0; } + for (int dd=0; dd<hdim; ++dd) { + invdir[dirs[dd]-1] = dd+1; + } +#endif - component = saved_component; + // Calculate lengths + for (int dd=0; dd<hdim; ++dd) { + if (lengths[dd]<0) { + const int totlen = *CCTK_ArrayGroupSizeI(GH, dirs[dd]-1, vindex); + assert (totlen>=0); + // Partial argument check + assert (global_startpoint[dirs[dd]-1]>=0); + assert (global_startpoint[dirs[dd]-1]<=totlen); + assert (downsample[dd]>0); + hsize[dd] = (totlen - global_startpoint[dirs[dd]-1]) / downsample[dd]; + } else { + hsize[dd] = lengths[dd]; + } + assert (hsize[dd]>=0); + } - return 0; + // Get the slab + *hdata = GetSlab (GH, + target_proc, + vindex, + vtimelvl, + hdim, + global_startpoint, + dirs, + downsample, + hsize); + + // Return with success + return 1; } + + } // namespace CarpetSlab diff --git a/Carpet/CarpetSlab/src/carpetslab.h b/Carpet/CarpetSlab/src/carpetslab.h index dd3f86136..6e3fd7f00 100644 --- a/Carpet/CarpetSlab/src/carpetslab.h +++ b/Carpet/CarpetSlab/src/carpetslab.h @@ -1,28 +1,28 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/Attic/carpetslab.h,v 1.1 2001/03/07 13:01:11 eschnett Exp $ +/* $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/Attic/carpetslab.h,v 1.2 2001/03/10 20:55:09 eschnett Exp $ */ -#include "cctk.h" +#ifndef CARPETSLAB_H +#define CARPETSLAB_H -int Hyperslab_GetLocalHyperslab (cGH* GH, - int vindex, - int vtimelvl, - int hdim, - const int global_startpoint [/*vdim*/], - const int directions [/*vdim*/], - const int lengths [/*hdim*/], - const int downsample [/*hdim*/], - void** hdata, - int hsize [/*hdim*/], - int ghsize [/*hdim*/], - int hoffset [/*hdim*/]); +#ifdef __cplusplus +namespace CarpetSlab { + extern "C" { +#endif + + int Hyperslab_GetHyperslab (cGH* const GH, + const int target_proc, + const int vindex, + const int vtimelvl, + const int hdim, + const int global_startpoint [/*vdim*/], + const int directions [/*vdim*/], + const int lengths [/*hdim*/], + const int downsample [/*hdim*/], + void** const hdata, + int hsize [/*hdim*/]); + +#ifdef __cplusplus + } /* extern "C" */ +} /* namespace CarpetSlab */ +#endif -int Hyperslab_GetHyperslab (cGH* GH, - int target_proc, - int vindex, - int vtimelvl, - int hdim, - const int global_startpoint [/*vdim*/], - const int directions [/*vdim*/], - const int lengths [/*hdim*/], - const int downsample [/*hdim*/], - void** hdata, - int hsize [/*hdim*/]); +#endif /* !defined(CARPETSLAB_H) */ diff --git a/Carpet/CarpetSlab/src/carpetslab.hh b/Carpet/CarpetSlab/src/carpetslab.hh index 444d14cf1..48fa52c31 100644 --- a/Carpet/CarpetSlab/src/carpetslab.hh +++ b/Carpet/CarpetSlab/src/carpetslab.hh @@ -1,36 +1,22 @@ -// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/Attic/carpetslab.hh,v 1.2 2001/03/07 13:01:11 eschnett Exp $ +// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/Attic/carpetslab.hh,v 1.3 2001/03/10 20:55:09 eschnett Exp $ -#include "cctk.h" +#ifndef CARPETSLAB_HH +#define CARPETSLAB_HH + +#include "carpetslab.h" namespace CarpetSlab { - extern "C" { - - int Hyperslab_GetLocalHyperslab (cGH* GH, - int vindex, - int vtimelvl, - int hdim, - const int global_startpoint [/*vdim*/], - const int directions [/*vdim*/], - const int lengths [/*hdim*/], - const int downsample [/*hdim*/], - void** hdata, - int hsize [/*hdim*/], - int ghsize [/*hdim*/], - int hoffset [/*hdim*/]); - - int Hyperslab_GetHyperslab (cGH* GH, - int target_proc, - int vindex, - int vtimelvl, - int hdim, - const int global_startpoint [/*vdim*/], - const int directions [/*vdim*/], - const int lengths [/*hdim*/], - const int downsample [/*hdim*/], - void** hdata, - int hsize [/*hdim*/]); - - } // extern "C" + void* GetSlab (cGH* const cgh, + const int dest_proc, + const int n, + const int tl, + const int hdim, + const int origin[/*vdim*/], + const int dirs[/*hdim*/], + const int stride[/*hdim*/], + const int length[/*hdim*/]); } // namespace CarpetSlab + +#endif // !defined(CARPETSLAB_HH) diff --git a/CarpetAttic/Cart3dTest/src/Initial.c b/CarpetAttic/Cart3dTest/src/Initial.c index f88c25b28..384dd3453 100644 --- a/CarpetAttic/Cart3dTest/src/Initial.c +++ b/CarpetAttic/Cart3dTest/src/Initial.c @@ -17,7 +17,6 @@ int Cart3dTest_Initial (CCTK_ARGUMENTS) int r; int i,j,k; - int ind; int vi[9]; @@ -80,36 +79,42 @@ int Cart3dTest_Initial (CCTK_ARGUMENTS) for (j=0; j<cctk_lsh[1]; ++j) { for (i=0; i<cctk_lsh[0]; ++i) { - ind = CCTK_GFINDEX3D(cctkGH, i,j,k); + const int ii = cctk_lbnd[0] + i; + const int jj = cctk_lbnd[1] + j; + const int kk = cctk_lbnd[2] + k; + + const int ind = CCTK_GFINDEX3D(cctkGH, i,j,k); + + /* store the position in the components */ - s[ind] = (( 1 * 100 + i) * 100 + j) * 100 + k; + s[ind] = (( 1 * 100 + ii) * 100 + jj) * 100 + kk; - vx[ind] = ((11 * 100 + i) * 100 + j) * 100 + k; - vy[ind] = ((12 * 100 + i) * 100 + j) * 100 + k; - vz[ind] = ((13 * 100 + i) * 100 + j) * 100 + k; + vx[ind] = ((11 * 100 + ii) * 100 + jj) * 100 + kk; + vy[ind] = ((12 * 100 + ii) * 100 + jj) * 100 + kk; + vz[ind] = ((13 * 100 + ii) * 100 + jj) * 100 + kk; - txx[ind] = ((21 * 100 + i) * 100 + j) * 100 + k; - txy[ind] = ((22 * 100 + i) * 100 + j) * 100 + k; - txz[ind] = ((23 * 100 + i) * 100 + j) * 100 + k; - tyy[ind] = ((24 * 100 + i) * 100 + j) * 100 + k; - tyz[ind] = ((25 * 100 + i) * 100 + j) * 100 + k; - tzz[ind] = ((26 * 100 + i) * 100 + j) * 100 + k; + txx[ind] = ((21 * 100 + ii) * 100 + jj) * 100 + kk; + txy[ind] = ((22 * 100 + ii) * 100 + jj) * 100 + kk; + txz[ind] = ((23 * 100 + ii) * 100 + jj) * 100 + kk; + tyy[ind] = ((24 * 100 + ii) * 100 + jj) * 100 + kk; + tyz[ind] = ((25 * 100 + ii) * 100 + jj) * 100 + kk; + tzz[ind] = ((26 * 100 + ii) * 100 + jj) * 100 + kk; - ax[ind] = ((31 * 100 + i) * 100 + j) * 100 + k; - ay[ind] = ((32 * 100 + i) * 100 + j) * 100 + k; - az[ind] = ((33 * 100 + i) * 100 + j) * 100 + k; + ax[ind] = ((31 * 100 + ii) * 100 + jj) * 100 + kk; + ay[ind] = ((32 * 100 + ii) * 100 + jj) * 100 + kk; + az[ind] = ((33 * 100 + ii) * 100 + jj) * 100 + kk; - fxx[ind] = ((41 * 100 + i) * 100 + j) * 100 + k; - fxy[ind] = ((42 * 100 + i) * 100 + j) * 100 + k; - fxz[ind] = ((43 * 100 + i) * 100 + j) * 100 + k; - fyx[ind] = ((44 * 100 + i) * 100 + j) * 100 + k; - fyy[ind] = ((45 * 100 + i) * 100 + j) * 100 + k; - fyz[ind] = ((46 * 100 + i) * 100 + j) * 100 + k; - fzx[ind] = ((47 * 100 + i) * 100 + j) * 100 + k; - fzy[ind] = ((48 * 100 + i) * 100 + j) * 100 + k; - fzz[ind] = ((49 * 100 + i) * 100 + j) * 100 + k; + fxx[ind] = ((41 * 100 + ii) * 100 + jj) * 100 + kk; + fxy[ind] = ((42 * 100 + ii) * 100 + jj) * 100 + kk; + fxz[ind] = ((43 * 100 + ii) * 100 + jj) * 100 + kk; + fyx[ind] = ((44 * 100 + ii) * 100 + jj) * 100 + kk; + fyy[ind] = ((45 * 100 + ii) * 100 + jj) * 100 + kk; + fyz[ind] = ((46 * 100 + ii) * 100 + jj) * 100 + kk; + fzx[ind] = ((47 * 100 + ii) * 100 + jj) * 100 + kk; + fzy[ind] = ((48 * 100 + ii) * 100 + jj) * 100 + kk; + fzz[ind] = ((49 * 100 + ii) * 100 + jj) * 100 + kk; } } |