diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/CollectData1D.c | 240 | ||||
| -rw-r--r-- | src/GetHyperslab.c | 198 | ||||
| -rw-r--r-- | src/GetLocalLine.c | 334 | ||||
| -rw-r--r-- | src/Hyperslab.c | 99 | ||||
| -rw-r--r-- | src/Mapping.c | 51 | ||||
| -rw-r--r-- | src/make.code.defn | 2 |
6 files changed, 255 insertions, 669 deletions
diff --git a/src/CollectData1D.c b/src/CollectData1D.c deleted file mode 100644 index 38b5f60..0000000 --- a/src/CollectData1D.c +++ /dev/null @@ -1,240 +0,0 @@ - -#include <stdio.h> -#include <stdlib.h> - -#include "cctk.h" -#include "CactusPUGH/PUGH/src/include/pugh.h" -#include "PUGHSlab.h" - -static const char *rcsid = "$Header$"; - -CCTK_FILEVERSION(CactusPUGH_PUGHSlab_CollectData1D_c) - -/* #define DEBUGME 1 */ - - -#define PICKUP_LINE_DATA(GH, cctk_type, vsize, start, directions, vdata, hdata)\ - { \ - int i, dim, point [3]; \ - cctk_type *typed_data = (cctk_type *) hdata; \ - \ - \ - for (i = 0; i < vsize; i++) \ - { \ - for (dim = 0; dim < 3; dim++) \ - point [dim] = start [dim] + i*directions [dim]; \ - *typed_data++ = ((cctk_type *) vdata) \ - [CCTK_GFINDEX3D (GH, point [0], point [1], point [2])]; \ - } \ - } - -/* local function prototypes */ -int Hyperslab_CollectData1D(const cGH *GH, int vindex, int vtimelvl, - const int *origin, const int *directions, - int downsample, int length, - void **hdata, int *hsize, int proc); -int GetLocalLine(const cGH *GH, pGExtras *extras, - int vardim, const int *S_l, const int *u_l, int ds, int len_l, - int *locstart, int *locend, int *nlocpoints); - - -static int CollectLocalData1D(const cGH *GH, int vindex, int vtimelvl, - const int *origin, const int *directions, - int downsample, int length, - void **hdata, int *hsize) -{ - int idim, ierr; - cGroup groupinfo; - int locstart [3], locend [3]; - void *vdata = CCTK_VarDataPtrI (GH, vtimelvl, vindex); - pGExtras *extras = ((pGA *) PUGH_pGH (GH)->variables [vindex][vtimelvl])->extras; - - - - CCTK_GroupData (CCTK_GroupIndexFromVarI (vindex), &groupinfo); - - if (groupinfo.dim != 3) - { - CCTK_WARN (1, "CollectLocalData1D() processes 3D variables only"); - return (-1); - } - - ierr = GetLocalLine(GH, extras, 3, origin, directions, downsample, length, - locstart, locend, hsize); - if (ierr < 0) - { - CCTK_WARN (1, "GetLocalLine() failed"); - return (ierr); - } - - /* Transform from global to local index coordinate */ - for(idim=0;idim<3;idim++) - { - locstart[idim]-=GH->cctk_lbnd[idim]; - locend[idim] -=GH->cctk_lbnd[idim]; - } - -#ifdef DEBUGME - printf("CollectLocal1D: %d npoints lb: %d %d %d\n", - *hsize,GH->cctk_lbnd[0],GH->cctk_lbnd[1],GH->cctk_lbnd[2]); -#endif - - *hdata= malloc (*hsize * CCTK_VarTypeSize (groupinfo.vartype)); - - /* CCTK_VARIABLE_INT */ - switch (groupinfo.vartype) - { - case CCTK_VARIABLE_CHAR: - PICKUP_LINE_DATA(GH, CCTK_BYTE, *hsize, locstart, directions, vdata, *hdata); - break; - - case CCTK_VARIABLE_INT: - PICKUP_LINE_DATA(GH, CCTK_INT, *hsize, locstart, directions, vdata, *hdata); - break; - - case CCTK_VARIABLE_REAL: - PICKUP_LINE_DATA(GH, CCTK_REAL, *hsize, locstart, directions, vdata, *hdata); - break; - - case CCTK_VARIABLE_COMPLEX: - PICKUP_LINE_DATA(GH, CCTK_COMPLEX, *hsize, locstart, directions, vdata, *hdata); - break; - - default: - CCTK_WARN (1, "Unsupported variable type"); - ierr = -1; - break; - } - - return (ierr); -} - - -int Hyperslab_CollectData1D(const cGH *GH, int vindex, int vtimelvl, - const int *origin, const int *directions, - int downsample, int length, - void **hdata, int *hsize, int proc) -{ - - void *loc_dptr; - int loc_dcount; - void *all_dptr; - int all_dcount=0; - int vtype, vtypesize; - int myproc,nprocs; - int ierr; - -#ifdef CCTK_MPI - int iproc; - const pGH *pughGH; - CCTK_INT tmp, *rem_dcount; - MPI_Datatype vmpitype; - int *recvcnts, *displs; -#endif - - - vtype = CCTK_VarTypeI(vindex); - vtypesize= CCTK_VarTypeSize(vtype); - - myproc = CCTK_MyProc(GH); - nprocs = CCTK_nProcs(GH); - - - ierr = CollectLocalData1D(GH, vindex, vtimelvl, - origin, directions, downsample, length, - &loc_dptr, &loc_dcount); - - if (ierr) - { - CCTK_WARN (1, "CollectLocalData1D() failed"); - return (ierr); - } - -#ifdef CCTK_MPI - - pughGH = PUGH_pGH(GH); - - rem_dcount = NULL; - recvcnts = NULL; - displs = NULL; - vmpitype = PUGH_MPIDataType (pughGH, vtype); - - /* Send local number of data elements */ - tmp = (CCTK_INT) loc_dcount; - - if (proc < 0 || proc == myproc) - { - rem_dcount= (CCTK_INT *) malloc (nprocs * sizeof (CCTK_INT)); - recvcnts = (int *) malloc (nprocs * sizeof (int)); - displs = (int *) malloc (nprocs * sizeof (int)); - } - - if (proc < 0) - { - CACTUS_MPI_ERROR (MPI_Allgather (&tmp, 1, PUGH_MPI_INT, - rem_dcount, 1, PUGH_MPI_INT, - pughGH->PUGH_COMM_WORLD)); - } else { - CACTUS_MPI_ERROR (MPI_Gather (&tmp, 1, PUGH_MPI_INT, - rem_dcount, 1, PUGH_MPI_INT, - proc, pughGH->PUGH_COMM_WORLD)); - } - - if (proc < 0 || proc == myproc) - { - displs[0] = 0; - all_dcount = recvcnts[0] = (int) rem_dcount[0]; - for (iproc=1;iproc<nprocs;iproc++) - { - recvcnts[iproc] = (int) rem_dcount[iproc]; - displs[iproc] = displs[iproc-1] + recvcnts[iproc-1]; - all_dcount += recvcnts[iproc]; - } - -#ifdef DEBUGME - printf (" local hsize: %d\n", loc_dcount); - for (iproc=0;iproc<nprocs;iproc++) - printf(" DISPLS[%d] rem: %d displs: %d all %d \n", - iproc, recvcnts[iproc],displs[iproc],all_dcount); - fflush(stdout); -#endif - - all_dptr = malloc (all_dcount * vtypesize); - } - else - { - all_dptr = NULL; - } - - if (proc < 0) - { - CACTUS_MPI_ERROR (MPI_Allgatherv(loc_dptr, loc_dcount, vmpitype, - all_dptr, recvcnts, displs, vmpitype, - pughGH->PUGH_COMM_WORLD)); - } else { - CACTUS_MPI_ERROR (MPI_Gatherv(loc_dptr, loc_dcount, vmpitype, - all_dptr, recvcnts, displs, vmpitype, - proc, pughGH->PUGH_COMM_WORLD)); - } - - if (proc < 0 || proc == myproc) - { - free (rem_dcount); - free (recvcnts); - free (displs); - } - -#else - all_dptr = loc_dptr; - all_dcount = loc_dcount; -#endif - - /* assign return values */ - if (proc < 0 || proc == myproc) - { - *hdata = all_dptr; - *hsize = all_dcount; - } - - return (0); -} diff --git a/src/GetHyperslab.c b/src/GetHyperslab.c index b6e0af0..b1206e2 100644 --- a/src/GetHyperslab.c +++ b/src/GetHyperslab.c @@ -39,19 +39,12 @@ static int GetLocalHyperslab (const cGH *GH, int timelevel, int hdatatype, void *hdata); -/******************************************************************** - ******************** External Routines ************************ - ********************************************************************/ -/* Gerd's routine to get 1D lines from a 3D array */ -int Hyperslab_CollectData1D (const cGH *GH, int vindex, int vtimelvl, - const int *origin, - const int *directions, - int downsampling, - int length, - void **hdata, - int *hsize, - int proc); - +static int GetDiagonalFromFrom3D (const cGH *GH, + const hslab_mapping_t *mapping, + int vindex, + int timelevel, + int hdatatype, + void *hdata); CCTK_INT Hyperslab_Get (const cGH *GH, @@ -274,27 +267,11 @@ static int GetLocalHyperslab (const cGH *GH, return (0); } - /* FIXME: hack for getting diagonals from 3D variables - This is calling Gerd's CollectData1D() routine which can - extract non-axis-parallel lines too but is fixed to 3D data. */ + /* diagonals from 3D variables are treated special */ if (mapping->is_diagonal_in_3D) { - const int origin[] = {0, 0, 0}; - const int directions[] = {1, 1, 1}; - int dummy_hsize; - void *dummy_hdata; - - - retval = Hyperslab_CollectData1D (GH, vindex, timelevel, origin, directions, - 1, -1, &dummy_hdata, &dummy_hsize, - mapping->target_proc); - if (! retval) - { - memcpy (hdata, dummy_hdata, - mapping->totals * CCTK_VarTypeSize (vinfo.vartype)); - free (dummy_hdata); - } - + retval = GetDiagonalFromFrom3D (GH, mapping, vindex, timelevel, hdatatype, + hdata); return (retval); } @@ -427,7 +404,7 @@ static int GetLocalHyperslab (const cGH *GH, /* get the byte pointer into the source array */ typed_vdata = (const char *) vdata + point[0]; -#if 1 +#if 0 fprintf (stderr, "***** base vdata %p offset %d '%s'\n", vdata, point[0], CCTK_FullName (vindex)); #endif for (i = 1; i < vinfo.dim; i++) @@ -445,7 +422,7 @@ fprintf (stderr, "***** base vdata %p offset %d '%s'\n", vdata, point[0], CCTK_F } else { -#if 1 +#if 0 fprintf (stderr, "***** copying %d bytes from %p tp %p\n", dim0_hsize, typed_vdata, typed_hdata); #endif memcpy (typed_hdata, typed_vdata, dim0_hsize); @@ -489,3 +466,156 @@ fprintf (stderr, "***** copying %d bytes from %p tp %p\n", dim0_hsize, typed_vda return (0); } + + +static int GetDiagonalFromFrom3D (const cGH *GH, + const hslab_mapping_t *mapping, + int vindex, + int timelevel, + int hdatatype, + void *hdata) +{ + int i, j, k, myproc, nprocs, linear_idx; + CCTK_INT local_npoints, npoints; + int vdatatype, htypesize, vtypesize; + const char *vdata; + char *local_hdata; + const pGH *pughGH; + const pGExtras *extras; +#ifdef CCTK_MPI + CCTK_INT *global_npoints; + int *recvcnts, *displs; + MPI_Datatype mpidatatype; +#endif + + + + /* get the pGH pointer and the variable's GA structure */ + pughGH = PUGH_pGH (GH); + extras = ((const pGA *) pughGH->variables[vindex][timelevel])->extras; + + vdatatype = CCTK_VarTypeI (vindex); + if (hdatatype < 0) + { + hdatatype = vdatatype; + } + htypesize = CCTK_VarTypeSize (hdatatype); + vtypesize = CCTK_VarTypeSize (vdatatype); + vdata = (const char *) CCTK_VarDataPtrI (GH, timelevel, vindex); + + myproc = CCTK_MyProc (GH); + nprocs = CCTK_nProcs (GH); + if (nprocs == 1) + { + local_hdata = (char *) hdata; + } + else + { + local_hdata = (char *) malloc (mapping->global_hsize[0] * htypesize); + } + + i = mapping->global_startpoint[0] - extras->lb[myproc][0]; + j = mapping->global_startpoint[1] - extras->lb[myproc][1]; + k = mapping->global_startpoint[2] - extras->lb[myproc][2]; + local_npoints = 0; + for (npoints = 0; npoints < mapping->global_hsize[0]; npoints++) + { + if (i >= extras->ownership[0][0][0] && i < extras->ownership[0][1][0] && + j >= extras->ownership[0][0][1] && j < extras->ownership[0][1][1] && + k >= extras->ownership[0][0][2] && k < extras->ownership[0][1][2]) + { + linear_idx = i + j*extras->hyper_volume[1] + k*extras->hyper_volume[2]; + if (vdatatype != hdatatype) + { + mapping->conversion_fn (vdata + linear_idx*vtypesize, local_hdata, + 1, 1, 1); + } + else + { + memcpy (local_hdata, vdata + linear_idx*vtypesize, htypesize); + } + local_hdata += htypesize; + local_npoints++; + } + i += mapping->downsample[0]; + j += mapping->downsample[1]; + k += mapping->downsample[2]; + } + local_hdata -= local_npoints * htypesize; + +#ifdef CCTK_MPI + if (nprocs > 1) + { + /* allocate communication buffers */ + myproc = CCTK_MyProc (GH); + if (mapping->target_proc < 0 || mapping->target_proc == myproc) + { + global_npoints = (CCTK_INT *) malloc (nprocs * sizeof (CCTK_INT)); + recvcnts = (int *) malloc (2 * nprocs * sizeof (int)); + displs = recvcnts + nprocs; + } + else + { + global_npoints = NULL; recvcnts = displs = NULL; hdata = NULL; + } + + /* gather the number of local points on each processor */ + if (mapping->target_proc < 0) + { + CACTUS_MPI_ERROR (MPI_Allgather (&local_npoints, 1, PUGH_MPI_INT, + global_npoints, 1, PUGH_MPI_INT, + pughGH->PUGH_COMM_WORLD)); + } + else + { + CACTUS_MPI_ERROR (MPI_Gather (&local_npoints, 1, PUGH_MPI_INT, + global_npoints, 1, PUGH_MPI_INT, + mapping->target_proc, + pughGH->PUGH_COMM_WORLD)); + } + + /* compute the receive count and displacement vectors */ + if (mapping->target_proc < 0 || mapping->target_proc == myproc) + { + for (i = 0; i < nprocs; i++) + { + recvcnts[i] = (int) global_npoints[i]; + displs[i] = i == 0 ? 0 : displs[i-1] + recvcnts[i-1]; + } + } + + /* get the MPI datatype for the hyperslab data */ + mpidatatype = PUGH_MPIDataType (pughGH, hdatatype); + + /* gather the local hyperslab data from each processor */ + if (mapping->target_proc < 0) + { + CACTUS_MPI_ERROR (MPI_Allgatherv (local_hdata, local_npoints, + mpidatatype, hdata, recvcnts, displs, + mpidatatype, pughGH->PUGH_COMM_WORLD)); + } + else + { + CACTUS_MPI_ERROR (MPI_Gatherv (local_hdata, local_npoints, + mpidatatype, hdata, recvcnts, displs, + mpidatatype, mapping->target_proc, + pughGH->PUGH_COMM_WORLD)); + } + + /* free communication buffers */ + if (mapping->target_proc < 0 || mapping->target_proc == myproc) + { + free (global_npoints); + free (recvcnts); + } + } +#endif /* CCTK_MPI */ + + /* free allocated resources */ + if (nprocs > 1) + { + free (local_hdata); + } + + return (0); +} diff --git a/src/GetLocalLine.c b/src/GetLocalLine.c deleted file mode 100644 index 78a9b57..0000000 --- a/src/GetLocalLine.c +++ /dev/null @@ -1,334 +0,0 @@ -#include <stdio.h> -#include <stdlib.h> -#include <math.h> - -#include "cctk.h" -#include "CactusPUGH/PUGH/src/include/pugh.h" - -static const char *rcsid = "$Header$"; - -CCTK_FILEVERSION(CactusPUGH_PUGHSlab_GetLocalLine_c) - -#define MAX_DIM 3 -#define MAX_FACE 6 -#define BAD -42 -#define EPS 1e-10 - -#define ABS(a) ((a)<0 ? -(a) : (a)) -#define SIGN(a) ((a)<0.0 ? (-1.0) : (1.0)) -#define LIN_INDEX3D(extras,P) ((P[0])+ (extras)->nsize[0]*((P[1]) + (extras)->nsize[1]*(P[2]))) - -#define MIN(a,b) ((a)<(b) ? (a) : (b)) -#define MAX(a,b) ((a)>(b) ? (a) : (b)) - -/*$#define LINE_DEBUG$*/ - -int IsPointOnProc_INT (int dim, int *gownership, const int *pnt); -int IsPointOnProc_REAL(int dim, int *gownership, CCTK_REAL *pnt); - -int IsPointOnProc_INT(int dim, int *gownership, const int *pnt) -{ - int retval=1; /*assume YES */ - int idim; - for (idim=0;idim<dim;idim++) - { - /* out of range -> set retval to zero;*/ - if ((gownership[2*idim ]>pnt[idim]) || - (gownership[2*idim+1]<pnt[idim])) - retval&=0; /* incremental AND on the return value */ - } - return(retval); -} - -int IsPointOnProc_REAL(int dim, int *gownership, CCTK_REAL *pnt) -{ - int retval=1; /*assume YES */ - int idim; - for (idim=0;idim<dim;idim++) - { - /* out of range -> set retval to zero;*/ - if ((gownership[2*idim ]>pnt[idim]) || - (gownership[2*idim+1]<pnt[idim])) - retval&=0; /* incremental AND on the return value */ - } - return(retval); -} - -int GetLocalLine(cGH *GH, pGExtras *extras, - int vardim, const int *S_l, const int *u_l, int ds, int len_l, - int *locstart, int *locend, int *nlocpoints) -{ - - const int n_s[MAX_DIM][MAX_FACE] = - {{ -1, 1, 0, 0, 0, 0}, - { 0, 0,-1, 1, 0, 0}, - { 0, 0, 0, 0,-1, 1}}; - - - /* variables for proc. surfaces */ - CCTK_REAL lmd_s[MAX_FACE], lmd[2], t_lmd; - int P_s[MAX_DIM][MAX_FACE]; - int gownership[MAX_FACE]; - - int notcoll[MAX_FACE]; - CCTK_REAL t_pnt1[MAX_DIM]; - - int prochasline=0; - - int tmp1,tmp2; - CCTK_REAL rtmp1; - - int idim, iface, linsec; - int myproc; - - - -#ifdef LINE_DEBUG - printf("\n+++++++++ GET LOCAL LINE \n"); - printf("S_l [XYZ] : %d %d %d \n", - S_l[0], S_l[1], S_l[2]); - printf("u_l [XYZ] : %d %d %d ds %d \n", - u_l[0], u_l[1], u_l[2], ds); - - printf("n_s: face 0/1: [%d %d %d], [%d %d %d] \n", - n_s[0][0],n_s[1][0],n_s[2][0], - n_s[0][1],n_s[1][1],n_s[2][1]); - printf("n_s: face 2/3: [%d %d %d], [%d %d %d] \n", - n_s[0][2],n_s[1][2],n_s[2][2], - n_s[0][3],n_s[1][3],n_s[2][3]); - printf("n_s: face 4/5: [%d %d %d], [%d %d %d] \n", - n_s[0][4],n_s[1][4],n_s[2][4], - n_s[0][5],n_s[1][5],n_s[2][5]); -#endif - - myproc = CCTK_MyProc (GH); - - for (idim=0;idim<vardim;idim++) - { - /* GLOBAL OWNERSHIP: lower/upperbnd are INCLUSIVE */ - gownership[2*idim] = extras->lb[myproc][idim]; - if (GH->cctk_bbox[2*idim]==0) - gownership[2*idim] += extras->nghostzones[idim]; - - gownership[2*idim+1]= extras->ub[myproc][idim]; - if (GH->cctk_bbox[2*idim+1]==0) - gownership[2*idim+1]-= (extras->nghostzones[idim]); - } - - - for (idim=0;idim<vardim;idim++) - { - P_s[idim][0] = gownership[2*idim]; /* Aufpunkt, minus faces */ - P_s[idim][2] = gownership[2*idim]; /* Aufpunkt, minus faces */ - P_s[idim][4] = gownership[2*idim]; /* Aufpunkt, minus faces */ - P_s[idim][1] = gownership[2*idim+1]; /* Aufpunkt, plus faces */ - P_s[idim][3] = gownership[2*idim+1]; /* Aufpunkt, plus faces */ - P_s[idim][5] = gownership[2*idim+1]; /* Aufpunkt, plus faces */ - } - -#ifdef LINE_DEBUG - printf("P_s (face 0): %d %d %d \n", - P_s[0][0],P_s[1][0],P_s[2][0]); - printf("P_s (face 1): %d %d %d \n\n", - P_s[0][1],P_s[1][1],P_s[2][1]); -#endif - - /* Lamda for each surface */ - for (iface=0;iface<2*vardim;iface++) - { - tmp1=0.0; tmp2=0.0; - - for (idim=0;idim<vardim;idim++) - { - tmp1+= (S_l[idim]-P_s[idim][iface])*n_s[idim][iface]; - tmp2+= (ds*u_l[idim]*n_s[idim][iface]); - } - if (tmp2==0) - { - lmd_s[iface]=BAD; - notcoll[iface]=0; - } - else - { - notcoll[iface]=1; - rtmp1 = (CCTK_REAL)(-1.0)*(tmp1)/(tmp2); - if (rtmp1<0.0) lmd_s[iface]=BAD; - else lmd_s[iface]=rtmp1; -#ifdef LINE_DEBUG - printf("lmd_s (face%d) %f \n",iface,rtmp1); -#endif - } - } - - /* Find 2 valid intersection with the proc surfaces */ - /* tricky: duplicate points corners and the startpoint */ - /* We have to restrict the line to positive lambda, otherwise - it is getting extremly ugly in finding out the local - points */ - - lmd[0] = BAD; - lmd[1] = BAD; - linsec = BAD; - -#ifdef LINE_DEBUG - printf("testing for local Startpt: %d %d %d", - S_l[0],S_l[1],S_l[2]); -#endif - - if (IsPointOnProc_INT(vardim,gownership,S_l)==1) - { -#ifdef LINE_DEBUG - printf(" ... YES ... valid first\n"); -#endif - lmd[0] = 0; - linsec = LIN_INDEX3D(extras,S_l); - } - else - { -#ifdef LINE_DEBUG - printf(" ... NO \n"); -#endif - } - - - for (iface=0;iface<MAX_FACE; iface++) - { -#ifdef LINE_DEBUG - printf("checking face %d: coll: %d lmd_s %d ", - iface, notcoll[iface], lmd_s[iface]); -#endif - if (notcoll[iface]==1 && lmd_s[iface]!=BAD) - { - for (idim=0;idim<vardim;idim++) - t_pnt1[idim]=S_l[idim]+lmd_s[iface]*ds*u_l[idim]; - -#ifdef LINE_DEBUG - printf(" testing for local intersec. (face %d lamda %+4.2f): %4.2f %4.2f %4.2f\n", - iface, lmd_s[iface], t_pnt1[0],t_pnt1[1],t_pnt1[2]); -#endif - - /* check if point is local */ - if (IsPointOnProc_REAL(vardim,gownership,t_pnt1)==1) - { -#ifdef LINE_DEBUG - printf(" ... YES ... "); -#endif - - /* CASE1: first valid face: this point is good */ - if (linsec<0) - { -#ifdef LINE_DEBUG - printf(" valid first \n"); -#endif - lmd[0] = lmd_s[iface]; - linsec = LIN_INDEX3D(extras,t_pnt1); - } - /* CASE2: we already have intersec. - check for duplicates: eg. a diagonal through the n-space, with - one exception: make sure that we allow the endpoint to be located - on the startpoint if lmd == 0 (which is the case if we start at the last - point of a proc) */ - else if (linsec==LIN_INDEX3D(extras,t_pnt1) && lmd_s[iface]>EPS) - { -#ifdef LINE_DEBUG - printf(" but duplicate, skip \n"); -#endif - } - else - { -#ifdef LINE_DEBUG - printf(" valid second \n"); -#endif - lmd[1] =lmd_s[iface]; - } - } - else - { -#ifdef LINE_DEBUG - printf(" ... NO \n"); -#endif - } - } - else - { -#ifdef LINE_DEBUG - printf(" ... skipping face: %d %4.2f \n",iface,lmd_s[iface]); -#endif - } - } - - /* flag if a proc has intersections, order the intersections - if we have any */ - if (lmd[0]!=BAD && lmd[1]!=BAD) - { - prochasline=1; - - /* switch: lmd[0] has to be the closer one to S_l */ - if (lmd[0]>lmd[1]) - { - t_lmd = lmd[0]; - lmd[0] = lmd[1]; - lmd[1] = t_lmd; - } - lmd[0]=ceil (lmd[0]); - lmd[1]=floor(lmd[1]); -#ifdef LINE_DEBUG - printf("LOCAL INTERSECTIONS: (ordered in lmd) %4.2f %4.2f \n",lmd[0],lmd[1]); -#endif - } - else if (lmd[0]==BAD && lmd[1]==BAD) - { - prochasline = 0; -#ifdef LINE_DEBUG - printf("NO INTERSECTION ON THIS PROC\n"); -#endif - } - else if ((lmd[0]==BAD && lmd[1]!=BAD) || - (lmd[1]==BAD && lmd[0]!=BAD)) - { - printf("ERROR: only one local intersection: you are screwed. lmd 0/1 %4.2f %4.2f\n", - lmd[0],lmd[1]); - } - - /* patch if lambdas are greater specified line length */ - if (len_l>=0) - { - if (len_l>lmd[0] && len_l<lmd[1]) - lmd[1]=floor(len_l); - else if (len_l<lmd[0]) - { - lmd[0]=BAD; - lmd[1]=BAD; - prochasline=0; - } - } - - /* Calculate the intersections points in global index space*/ - if (prochasline) - { - nlocpoints[0]= ABS(lmd[0]-lmd[1])+1; - for (idim=0;idim<vardim;idim++) - { - locstart[idim] = (S_l[idim]+lmd[0]*ds*u_l[idim]); - locend[idim] = (S_l[idim]+lmd[1]*ds*u_l[idim]); - } - } - else - { - nlocpoints[0] = 0; - for (idim=0;idim<vardim;idim++) - { - locstart[idim] = BAD; - locend[idim] = BAD; - } - } - -#ifdef LINE_DEBUG - printf("\n>>>>> GLOBAL START/END: (%d %d %d) -> (%d %d %d) = %d\n", - locstart[0],locstart[1],locstart[2], - locend[0],locend[1],locend[2],nlocpoints[0]); -#endif - - return(0); - -} diff --git a/src/Hyperslab.c b/src/Hyperslab.c index 9c482e6..dd5f840 100644 --- a/src/Hyperslab.c +++ b/src/Hyperslab.c @@ -131,16 +131,6 @@ CCTK_FILEVERSION(CactusPUGH_PUGHSlab_Hyperslab_c) } -/*** local function prototypes ***/ -/* Gerd's routine to get 1D lines from a 3D array */ -int Hyperslab_CollectData1D (const cGH *GH, int vindex, int vtimelvl, - const int *origin, - const int *directions, - int downsampling, - int length, - void **hdata, - int *hsize, - int proc); /* routine to check parameters passed to the Hyperslab routines */ static const char *checkParameters (const cGH *GH, int vindex, int vtimelvl, int hdim, @@ -746,46 +736,6 @@ int Hyperslab_GetHyperslab (const cGH *GH, #endif - /* do some plausibility checks */ - errormsg = checkParameters (GH, vindex, vtimelvl, hdim, - global_origin, directions, extents, - downsample_, hdata, hsize); - - /* FIXME: hack for getting diagonals from 3D variables - This is calling Gerd's CollectData1D() routine which can - extract non-axis-parallel lines too but is fixed to 3D data. */ - if (errormsg && ! strcmp (errormsg, "Given normal vector isn't axis-parallel")) - { - int length; - - - /* get the info on the variable to extract a hyperslab from */ - CCTK_GroupData (CCTK_GroupIndexFromVarI (vindex), &vinfo); - - if (hdim != 1 || vinfo.dim != 3) - { - CCTK_WARN (1, "Non-axis-parallel hyperslabs are supported as 1D lines " - "from 3D variables only"); - return (-1); - } - - length = extents[0]; - if (length > 0) - { - length--; - } - return (Hyperslab_CollectData1D (GH, vindex, vtimelvl, global_origin, - directions, downsample_[0], length, - hdata, hsize, target_proc)); - } - - /* immediately return in case of errors */ - if (errormsg) - { - CCTK_WARN (1, errormsg); - return (-1); - } - /* get processor information */ myproc = CCTK_MyProc (GH); nprocs = CCTK_nProcs (GH); @@ -812,6 +762,55 @@ int Hyperslab_GetHyperslab (const cGH *GH, /* get the info on the variable to extract a hyperslab from */ CCTK_GroupData (CCTK_GroupIndexFromVarI (vindex), &vinfo); + /* do some plausibility checks */ + errormsg = checkParameters (GH, vindex, vtimelvl, hdim, + global_origin, directions, extents, + downsample_, hdata, hsize); + + /* FIXME: hack for getting diagonals from 3D variables */ + if (errormsg && ! strcmp (errormsg, "Given normal vector isn't axis-parallel")) + { + int mapping; + + + if (hdim != 1 || vinfo.dim != 3) + { + CCTK_WARN (1, "Non-axis-parallel hyperslabs are supported as diagonals " + "from non-staggered 3D variables only"); + return (-1); + } + + mapping = Hyperslab_DefineGlobalMappingByIndex (GH, vindex, 1, directions, + global_origin, extents, + downsample_, -1, + target_proc, NULL, hsize); + if (mapping < 0) + { + CCTK_WARN (1, "Failed to define hyperslab mapping for 3D diagonal"); + return (-1); + } + if (hsize[0] > 0) + { + *hdata = malloc (hsize[0] * CCTK_VarTypeSize (vinfo.vartype)); + retval = Hyperslab_Get (GH, mapping, vindex, 0, vinfo.vartype, *hdata); + } + else + { + *hdata = NULL; + retval = 0; + } + Hyperslab_FreeMapping (mapping); + + return (retval); + } + + /* immediately return in case of errors */ + if (errormsg) + { + CCTK_WARN (1, errormsg); + return (-1); + } + /* set the pointers to pass to Hyperslab_GetLocalHyperslab() */ if (nprocs == 1) { diff --git a/src/Mapping.c b/src/Mapping.c index c74be93..8e80893 100644 --- a/src/Mapping.c +++ b/src/Mapping.c @@ -218,6 +218,17 @@ CCTK_INT Hyperslab_DefineGlobalMappingByIndex ( return (-7); } } + + /* diagonals can be extracted from non-staggered 3D variables only */ + if (mapping->is_diagonal_in_3D && vinfo.stagtype != 0) + { + free (mapping->vectors); free (mapping); + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Hyperslab_DefineGlobalMappingByIndex: diagonals can be " + "extracted from non-staggered 3D variables only"); + return (-7); + } + for (vdim = 0; vdim < (unsigned int) vinfo.dim; vdim++) { mapping->do_dir[vdim] = 0; @@ -256,6 +267,15 @@ CCTK_INT Hyperslab_DefineGlobalMappingByIndex ( } hdim++; } + else if (mapping->is_diagonal_in_3D && + origin[vdim] + extent[0] > GA->extras->nsize[vdim]) + { + free (mapping->vectors); free (mapping); + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Hyperslab_DefineGlobalMappingByIndex: extent in " + "%d-direction exceeds grid size", vdim); + return (-8); + } } /* now fill out the hyperslab mapping structure */ @@ -281,6 +301,23 @@ CCTK_INT Hyperslab_DefineGlobalMappingByIndex ( } hdim++; } + else if (mapping->is_diagonal_in_3D) + { + mapping->totals = extent[0] / mapping->downsample[0]; + if (extent[0] % mapping->downsample[0]) + { + mapping->totals++; + } + /* subtract ghostzones for periodic BC */ + if (GA->connectivity->perme[vdim]) + { + mapping->totals -= 2 * GA->extras->nghostzones[vdim]; + } + if ((unsigned int) mapping->global_hsize[0] > mapping->totals) + { + mapping->global_hsize[0] = mapping->totals; + } + } } /* check whether the full local data patch was requested as hyperslab */ @@ -293,17 +330,12 @@ CCTK_INT Hyperslab_DefineGlobalMappingByIndex ( } else if (mapping->is_diagonal_in_3D) { - /* diagonals in 3D are special: global_hsize is minimum of grid extents */ - mapping->global_hsize[0] = GH->cctk_gsh[0]; - if (mapping->global_hsize[0] > GH->cctk_gsh[1]) - { - mapping->global_hsize[0] = GH->cctk_gsh[1]; - } - if (mapping->global_hsize[0] > GH->cctk_gsh[2]) + /* just initialize the downsample and global_startpoint vectors */ + for (vdim = 0; vdim < (unsigned int) vinfo.dim; vdim++) { - mapping->global_hsize[0] = GH->cctk_gsh[2]; + mapping->downsample[vdim] = mapping->downsample[0]; + mapping->global_startpoint[vdim] = origin[vdim]; } - mapping->totals = mapping->global_hsize[0]; } else { @@ -371,7 +403,6 @@ CCTK_INT Hyperslab_DefineGlobalMappingByIndex ( mapping->local_endpoint[vdim] = -1; } -#define DEBUG 1 #ifdef DEBUG printf ("direction %d: global ranges [%d, %d), local ranges[%d, %d)\n", vdim, diff --git a/src/make.code.defn b/src/make.code.defn index 7d322e8..d3c05c7 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -2,4 +2,4 @@ # $Header$ # Source files in this directory -SRCS = Mapping.c GetHyperslab.c Hyperslab.c NewHyperslab.c DatatypeConversion.c CollectData1D.c GetLocalLine.c +SRCS = Mapping.c GetHyperslab.c Hyperslab.c NewHyperslab.c DatatypeConversion.c |