From 3e40b5b2ccc932dff681a9e9b2c1666740d20d08 Mon Sep 17 00:00:00 2001 From: eschnett Date: Fri, 8 Mar 2013 20:31:14 +0000 Subject: Replace Cactus complex number type with C/C++ complex numbers Map CCTK_COMPLEX to "double complex" in C, and "complex" in C++. (It is already mapped to "double complex" in Fortran.) Update type definitions. Re-implement Cactus complex number math functions by calling the respective C functions. Update thorn that access real and imaginary parts of complex numbers to use standard-conforming methods instead. git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/ReflectionSymmetry/trunk@55 082bdb00-0f4f-0410-b49e-b1835e5f2039 --- src/apply.c | 153 +++++++++++++++++++++++++++++++++++++++++------------------- 1 file changed, 106 insertions(+), 47 deletions(-) diff --git a/src/apply.c b/src/apply.c index f179c0e..0ef2674 100644 --- a/src/apply.c +++ b/src/apply.c @@ -56,20 +56,14 @@ int const jj = jjoff + jdir * j; \ int const kk = kkoff + kdir * k; \ int const srcind = ii + ni * (jj + nj * kk); \ - REAL(dstvar[dstind] RE = parity * srcvar[srcind] RE;) \ - IMAG(dstvar[dstind] IM = parity * srcvar[srcind] IM;) \ + dstvar[dstind] = parity * srcvar[srcind]; \ } \ } \ } \ -#define COPY_POST(VARTYPE) \ +#define COPY_POST(VARTYPE) \ } -#define REAL(x) x -#define IMAG(x) /* nothing */ -#define RE /* nothing */ -#define IM ERROR ERROR ERROR - #ifdef HAVE_CCTK_INT1 COPY_PRE(CCTK_INT1) #pragma omp parallel for @@ -119,16 +113,6 @@ COPY_LOOP(CCTK_REAL16) COPY_POST(CCTK_REAL16) #endif -#undef REAL -#undef IMAG -#undef RE -#undef IM - -#define REAL(x) x -#define IMAG(x) x -#define RE .Re -#define IM .Im - #ifdef HAVE_CCTK_COMPLEX8 COPY_PRE(CCTK_COMPLEX8) #pragma omp parallel for @@ -150,11 +134,6 @@ COPY_LOOP(CCTK_COMPLEX32) COPY_POST(CCTK_COMPLEX32) #endif -#undef REAL -#undef IMAG -#undef RE -#undef IM - #undef COPY_PRE #undef COPY_LOOP #undef COPY_POST @@ -177,21 +156,35 @@ BndReflectVI (cGH const * restrict const cctkGH, int table; char tensortypealias[1000]; - enum tensortype { UNKNOWN, + enum tensortype { TT_UNKNOWN, SCALAR, VECTOR, SYMTENSOR, SYMTENSOR3, TENSOR, WEYLSCALARS_REAL, MANUALCARTESIAN }; enum tensortype ttype; CCTK_INT tensorparity; int tcomponent; + /* The stagger type is defined for grid variable groups, and + depending on this type, each variable in the group may be + staggered differently. This implies that groups with FACE or EDGE + staggering have well-defined tensor types. */ + /* The stagger type is defined assuming a cell-centered grid, since + this seems to be the common case where different stagger types + are used. If the grid is vertex centered, then we currently abort + when staggered grid variables are encountered. (We determine the + centering of the grid via the avoid_origin_* parameters.) */ + char staggertype[1000]; + enum staggertype { ST_UNKNOWN, CELL, FACE, EDGE, VERTEX }; + enum staggertype stype; + int do_reflection[6]; - int do_stagger[6]; + int do_stagger_grid[6]; int dir, face; - int lsh[3], imin[3], imax[3], ioff[3], idir[3]; + int ash[3], imin[3], imax[3], ioff[3], idir[3]; int parity; + int do_stagger; int manual_parities[3]; int d; @@ -287,7 +280,7 @@ BndReflectVI (cGH const * restrict const cctkGH, free (groupname); } - ttype = UNKNOWN; + ttype = TT_UNKNOWN; tcomponent = 0; if (CCTK_EQUALS (tensortypealias, "scalar")) { /* scalar */ @@ -414,6 +407,45 @@ BndReflectVI (cGH const * restrict const cctkGH, + /* Get and check stagger type information */ + ierr = Util_TableGetString + (table, sizeof staggertype, staggertype, "staggertype"); + if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY) { + /* assume cell centering */ + strcpy (staggertype, "cell"); + } else if (ierr<0) { + char * groupname = CCTK_GroupName(gi); + assert (groupname); + CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, + "Error in stagger type declaration for group \"%s\": %d", + groupname, ierr); + free (groupname); + } + + stype = ST_UNKNOWN; + if (CCTK_EQUALS (staggertype, "cell")) { + /* cell */ + stype = CELL; + } else if (CCTK_EQUALS (staggertype, "face")) { + /* face */ + stype = FACE; + } else if (CCTK_EQUALS (staggertype, "edge")) { + /* edge */ + stype = EDGE; + } else if (CCTK_EQUALS (staggertype, "vertex")) { + /* vertex */ + stype = VERTEX; + } else { + char * groupname = CCTK_GroupName(gi); + assert (groupname); + CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, + "Illegal stagger type \"%s\" for group \"%s\"", + staggertype, groupname); + free (groupname); + } + + + /* Reflection symmetry information */ do_reflection[0] = reflection_x; do_reflection[1] = reflection_upper_x; @@ -422,12 +454,12 @@ BndReflectVI (cGH const * restrict const cctkGH, do_reflection[4] = reflection_z; do_reflection[5] = reflection_upper_z; - do_stagger[0] = avoid_origin_x; - do_stagger[1] = avoid_origin_upper_x; - do_stagger[2] = avoid_origin_y; - do_stagger[3] = avoid_origin_upper_y; - do_stagger[4] = avoid_origin_z; - do_stagger[5] = avoid_origin_upper_z; + do_stagger_grid[0] = avoid_origin_x; + do_stagger_grid[1] = avoid_origin_upper_x; + do_stagger_grid[2] = avoid_origin_y; + do_stagger_grid[3] = avoid_origin_upper_y; + do_stagger_grid[4] = avoid_origin_z; + do_stagger_grid[5] = avoid_origin_upper_z; @@ -512,9 +544,36 @@ BndReflectVI (cGH const * restrict const cctkGH, assert (0); } + /* Find staggering */ + do_stagger = do_stagger_grid[2*dir+face]; + switch (stype) { + case CELL: + /* do nothing */ + break; + case FACE: + assert (face == 0); /* assume lower face */ + assert (do_stagger); /* assume cell-centered grid */ + assert (ttype == VECTOR); /* TODO: support other tensor types */ + if (dir == tcomponent) do_stagger = !do_stagger; + break; + case EDGE: + assert (face == 0); /* assume lower face */ + assert (do_stagger); /* assume cell-centered grid */ + assert (ttype == VECTOR); /* TODO: support other tensor types */ + if (dir != tcomponent) do_stagger = !do_stagger; + break; + case VERTEX: + assert (face == 0); /* assume lower face */ + assert (do_stagger); /* assume cell-centered grid */ + do_stagger = !do_stagger; + break; + default: + assert (0); + } + /* Find region extent */ for (d=0; d<3; ++d) { - lsh[d] = cctkGH->cctk_lsh[d]; + ash[d] = cctkGH->cctk_ash[d]; imin[d] = 0; imax[d] = cctkGH->cctk_lsh[d]; ioff[d] = 0; @@ -539,11 +598,11 @@ BndReflectVI (cGH const * restrict const cctkGH, if (face == 0) { imax[dir] = cctkGH->cctk_nghostzones[dir]; ioff[dir] = (+ 2*cctkGH->cctk_nghostzones[dir] - 1 - + (do_stagger[2*dir+face] ? 0 : 1)); + + (do_stagger ? 0 : 1)); idir[dir] = -1; } else { imin[dir] = cctkGH->cctk_lsh[dir] - cctkGH->cctk_nghostzones[dir]; - ioff[dir] = - 1 - (do_stagger[2*dir+face] ? 0 : 1); + ioff[dir] = - 1 - (do_stagger ? 0 : 1); idir[dir] = -1; } @@ -554,7 +613,7 @@ BndReflectVI (cGH const * restrict const cctkGH, int const have_points = cctkGH->cctk_gsh[dir]; int const need_points = 3 * cctkGH->cctk_nghostzones[dir] - + !do_stagger[2*dir] + !do_stagger[2*dir+1]; + + !do_stagger_grid[2*dir] + !do_stagger_grid[2*dir+1]; if (need_points > have_points) { CCTK_VWarn (CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING, "Cannot apply symmetry boundary zones in the %s %c direction, since there seem to be more symmetry zones than interior zones", @@ -567,7 +626,7 @@ BndReflectVI (cGH const * restrict const cctkGH, switch (group.vartype) { #define ARGS varptr, varptr, \ - lsh[0], lsh[1], lsh[2], \ + ash[0], ash[1], ash[2], \ imin[0], imin[1], imin[2], \ imax[0], imax[1], imax[2], \ ioff[0], ioff[1], ioff[2], \ @@ -706,7 +765,7 @@ CheckBoundaryParameters (cGH const * restrict const cctkGH, CCTK_INT * restrict shiftout; int do_reflection[6]; /* This thorn's parameters */ - int do_stagger[6]; + int do_stagger_grid[6]; int d; int ierr; @@ -751,12 +810,12 @@ CheckBoundaryParameters (cGH const * restrict const cctkGH, do_reflection[4] = reflection_z; do_reflection[5] = reflection_upper_z; - do_stagger[0] = avoid_origin_x; - do_stagger[1] = avoid_origin_upper_x; - do_stagger[2] = avoid_origin_y; - do_stagger[3] = avoid_origin_upper_y; - do_stagger[4] = avoid_origin_z; - do_stagger[5] = avoid_origin_upper_z; + do_stagger_grid[0] = avoid_origin_x; + do_stagger_grid[1] = avoid_origin_upper_x; + do_stagger_grid[2] = avoid_origin_y; + do_stagger_grid[3] = avoid_origin_upper_y; + do_stagger_grid[4] = avoid_origin_z; + do_stagger_grid[5] = avoid_origin_upper_z; /* Check the boundary sizes */ for (d=0; d<6; ++d) { @@ -777,13 +836,13 @@ CheckBoundaryParameters (cGH const * restrict const cctkGH, "The %s %c face is a symmetry boundary. The corresponding CoordBase boundary must not be internal.", face, dir); } - if (do_stagger[d] != is_staggered[d]) { + if (do_stagger_grid[d] != is_staggered[d]) { CCTK_VWarn (CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING, "The %s %c face is a symmetry boundary. The symmetry condition and the corresponding CoordBase boundary must either be both staggered or both not staggered.", face, dir); } - if ((do_stagger[d] ? 0 : 1) != shiftout[d]) { + if ((do_stagger_grid[d] ? 0 : 1) != shiftout[d]) { CCTK_VWarn (CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING, "The %s %c face is a symmetry boundary. If the symmetry condition is staggered, then the corresponding CoordBase shiftout must be 0; otherwise it must be 1.", -- cgit v1.2.3