aboutsummaryrefslogtreecommitdiff
path: root/Carpet/CarpetSlab/src/slab.cc
diff options
context:
space:
mode:
Diffstat (limited to 'Carpet/CarpetSlab/src/slab.cc')
-rw-r--r--Carpet/CarpetSlab/src/slab.cc586
1 files changed, 556 insertions, 30 deletions
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,
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-06 10:30:47 +0000