diff options
Diffstat (limited to 'src/apply.c')
| -rw-r--r-- | src/apply.c | 627 |
1 files changed, 627 insertions, 0 deletions
diff --git a/src/apply.c b/src/apply.c new file mode 100644 index 0000000..114c1f7 --- /dev/null +++ b/src/apply.c @@ -0,0 +1,627 @@ +#include <assert.h> +#include <stdlib.h> +#include <string.h> + +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" + +#include "util_ErrorCodes.h" +#include "util_Table.h" + +#include "reflection.h" + + + +static const char rcsid[] = "$Header$"; +CCTK_FILEVERSION(AEIDevelopment_ReflectionSymmetry_apply_c); + + + +#define COPY(VARTYPE) \ + static void \ + copy_##VARTYPE (VARTYPE const * restrict const srcvar, \ + VARTYPE * restrict const dstvar, \ + int const ni, int const nj, int const nk, \ + int const imin, int const jmin, int const kmin, \ + int const imax, int const jmax, int const kmax, \ + int const ioff, int const joff, int const koff, \ + int const idir, int const jdir, int const kdir, \ + int const parity) \ + { \ + int i, j, k; \ + \ + for (k=kmin; k<kmax; ++k) \ + { \ + for (j=jmin; j<jmax; ++j) \ + { \ + for (i=imin; j<imax; ++i) \ + { \ + int const dstind = i + ni * (j + nj * k); \ + int const ii = ioff + idir * i; \ + int const jj = joff + jdir * j; \ + int const kk = koff + 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;) \ + } \ + } \ + } \ + } + +#define REAL(x) x +#define IMAG(x) /* nothing */ +#define RE /* nothing */ +#define IM /* nothing */ + +#ifdef CCTK_INT1 +COPY(CCTK_INT1) +#endif + +#ifdef CCTK_INT2 + COPY(CCTK_INT2) +#endif + +#ifdef CCTK_INT4 + COPY(CCTK_INT4) +#endif + +#ifdef CCTK_INT8 + COPY(CCTK_INT8) +#endif + +#ifdef CCTK_REAL4 + COPY(CCTK_REAL4) +#endif + +#ifdef CCTK_REAL8 + COPY(CCTK_REAL8) +#endif + +#ifdef CCTK_REAL16 + COPY(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 CCTK_REAL4 + COPY(CCTK_COMPLEX8) +#endif + +#ifdef CCTK_REAL8 + COPY(CCTK_COMPLEX16) +#endif + +#ifdef CCTK_REAL16 + COPY(CCTK_COMPLEX32) +#endif + +#undef REAL +#undef IMAG +#undef RE +#undef IM + +#undef COPY + + + +static int +BndReflectVI (cGH const * restrict const cctkGH, + int const * restrict const stencil, + int const vi) +{ + DECLARE_CCTK_PARAMETERS; + + int gi; + cGroup group; + cGroupDynamicData data; + int firstvar, numvars; + char * restrict fullname; + char * restrict groupname; + + void * restrict varptr; + + int table; + char tensortypealias[1000]; + enum tensortype { UNKNOWN, SCALAR, VECTOR, TENSOR }; + enum tensortype ttype; + int tcomponent; + + int do_reflection[6]; + int do_stagger[6]; + + int dir, face; + + int lsh[3], imin[3], imax[3], ioff[3], idir[3]; + + int parity; + + int d; + + int ierr; + + + + /* Check arguments */ + if (! cctkGH) + { + CCTK_WARN (0, "Argument cctkGH is NULL"); + } + if (! stencil) + { + CCTK_WARN (0, "Argument stencil is NULL"); + } + if (vi < 0 || vi >= CCTK_NumVars()) + { + CCTK_WARN (0, "Illegal variable index"); + } + + if (verbose) + { + fullname = CCTK_FullName (vi); + if (! fullname) + { + CCTK_WARN (0, "Internal error in CCTK_FullName"); + } + CCTK_VInfo (CCTK_THORNSTRING, + "Applying reflection boundary conditions to \"%s\"", + fullname); + free (fullname); + } + + + + /* Get and check group information */ + gi = CCTK_GroupIndexFromVarI (vi); + if (gi < 0 || gi > CCTK_NumGroups()) + { + CCTK_WARN (0, "Internal error in CCTK_GroupIndexFromVarI"); + } + + ierr = CCTK_GroupData (gi, &group); + assert (!ierr); + assert (group.grouptype == CCTK_GF); + assert (group.disttype == CCTK_DISTRIB_DEFAULT); + assert (group.stagtype == 0); + + firstvar = CCTK_FirstVarIndexI (gi); + assert (firstvar>=0 && firstvar<CCTK_NumVars()); + numvars = CCTK_NumVarsInGroupI (gi); + assert (numvars>=0); + + ierr = CCTK_GroupDynamicData (cctkGH, gi, &data); + assert (!ierr); + + table = CCTK_GroupTagsTableI(gi); + assert (table>=0); + + varptr = CCTK_VarDataPtrI (cctkGH, 0, vi); + assert (varptr); + + + + /* Get and check tensor type information */ + ierr = Util_TableGetString + (table, sizeof tensortypealias, tensortypealias, "tensortypealias"); + if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY) { + groupname = CCTK_GroupName(gi); + assert (groupname); + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Tensor type alias not declared for group \"%s\" -- assuming a scalar", + groupname); + free (groupname); + strcpy (tensortypealias, "scalar"); + } else if (ierr<0) { + groupname = CCTK_GroupName(gi); + assert (groupname); + CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, + "Error in tensor type alias declaration for group \"%s\"", + groupname); + free (groupname); + } + + ttype = UNKNOWN; + tcomponent = 0; + if (CCTK_EQUALS (tensortypealias, "scalar")) + { + /* scalar */ + if (numvars != 1) { + groupname = CCTK_GroupName(gi); + assert (groupname); + CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING, + "Group \"%s\" has the tensor type alias \"scalar\", but contains more than 1 element", + groupname); + free (groupname); + } + ttype = SCALAR; + tcomponent = 0; + } + else if (CCTK_EQUALS (tensortypealias, "u")) + { + /* vector */ + assert (numvars == 3); + ttype = VECTOR; + tcomponent = vi - firstvar; + } else if (CCTK_EQUALS (tensortypealias, "dd_sym")) { + /* symmetric tensor */ + assert (numvars == 6); + ttype = TENSOR; + tcomponent = vi - firstvar; + } else { + char * groupname = CCTK_GroupName(gi); + assert (groupname); + CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING, + "Illegal tensor type alias for group \"%s\"", + groupname); + free (groupname); + } + + switch (ttype) + { + case SCALAR: + assert (tcomponent>=0 && tcomponent<1); + break; + case VECTOR: + assert (tcomponent>=0 && tcomponent<3); + break; + case TENSOR: + assert (tcomponent>=0 && tcomponent<6); + break; + default: + assert (0); + } + + + + /* Reflection symmetry information */ + do_reflection[0] = reflection_x; + do_reflection[1] = reflection_upper_x; + do_reflection[2] = reflection_y; + do_reflection[3] = reflection_upper_y; + 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; + + + + /* Loop over all directions and faces */ + for (dir=0; dir<3; ++dir) + { + for (face=0; face<2; ++face) + { + /* If there is a reflection symmetry on that face */ + if (do_reflection[2*dir+face]) + { + /* If we have the outer boundary of that face */ + if (cctkGH->cctk_bbox[2*dir+face]) + { + + /* Find parity */ + switch (ttype) + { + case SCALAR: + parity = +1; + break; + case VECTOR: + parity = dir == tcomponent ? -1 : +1; + break; + case TENSOR: + switch (tcomponent) + { + case 0: parity = +1; break; + case 1: parity = (dir == 0 || dir == 1) ? -1 : +1; break; + case 2: parity = (dir == 0 || dir == 2) ? -1 : +1; break; + case 3: parity = +1; break; + case 4: parity = (dir == 1 || dir == 2) ? -1 : +1; break; + case 5: parity = +1; break; + default: assert (0); + } + break; + default: + assert (0); + } + + /* Find region extent */ + for (d=0; d<3; ++d) + { + lsh[d] = cctkGH->cctk_lsh[d]; + imin[d] = 0; + imax[d] = cctkGH->cctk_lsh[d]; + ioff[d] = 0; + } + if (face == 0) + { + imax[d] = cctkGH->cctk_nghostzones[d]; + ioff[d] = (+ cctkGH->cctk_nghostzones[d] + + (do_stagger[2*d+face] ? 0 : 1)); + } + else + { + imin[d] = cctkGH->cctk_lsh[d] - cctkGH->cctk_nghostzones[d]; + ioff[d] = (- cctkGH->cctk_nghostzones[d] + - (do_stagger[2*d+face] ? 0 : 1)); + } + + /* Copy region */ + switch (group.vartype) + { + +#ifdef CCTK_INT1 + case CCTK_VARIABLE_INT1: +#ifdef CCTK_INTER_PRECISION_1 + case CCTK_VARIABLE_INT: +#endif + copy_CCTK_INT1 (varptr, varptr, + lsh[0], lsh[1], lsh[2], + imin[0], imin[1], imin[2], + imax[0], imax[1], imax[2], + ioff[0], ioff[1], ioff[2], + idir[0], idir[1], idir[2], + parity); + break; +#endif + +#ifdef CCTK_INT2 + case CCTK_VARIABLE_INT2: +#ifdef CCTK_INTER_PRECISION_2 + case CCTK_VARIABLE_INT: +#endif + copy_CCTK_INT2 (varptr, varptr, + lsh[0], lsh[1], lsh[2], + imin[0], imin[1], imin[2], + imax[0], imax[1], imax[2], + ioff[0], ioff[1], ioff[2], + idir[0], idir[1], idir[2], + parity); + break; +#endif + +#ifdef CCTK_INT4 + case CCTK_VARIABLE_INT4: +#ifdef CCTK_INTER_PRECISION_4 + case CCTK_VARIABLE_INT: +#endif + copy_CCTK_INT4 (varptr, varptr, + lsh[0], lsh[1], lsh[2], + imin[0], imin[1], imin[2], + imax[0], imax[1], imax[2], + ioff[0], ioff[1], ioff[2], + idir[0], idir[1], idir[2], + parity); + break; +#endif + +#ifdef CCTK_INT8 + case CCTK_VARIABLE_INT8: +#ifdef CCTK_INTER_PRECISION_8 + case CCTK_VARIABLE_INT: +#endif + copy_CCTK_INT8 (varptr, varptr, + lsh[0], lsh[1], lsh[2], + imin[0], imin[1], imin[2], + imax[0], imax[1], imax[2], + ioff[0], ioff[1], ioff[2], + idir[0], idir[1], idir[2], + parity); + break; +#endif + +#ifdef CCTK_REAL4 + case CCTK_VARIABLE_REAL4: +#ifdef CCTK_REAL_PRECISION_4 + case CCTK_VARIABLE_REAL: +#endif + copy_CCTK_REAL4 (varptr, varptr, + lsh[0], lsh[1], lsh[2], + imin[0], imin[1], imin[2], + imax[0], imax[1], imax[2], + ioff[0], ioff[1], ioff[2], + idir[0], idir[1], idir[2], + parity); + break; +#endif + +#ifdef CCTK_REAL8 + case CCTK_VARIABLE_REAL8: +#ifdef CCTK_REAL_PRECISION_8 + case CCTK_VARIABLE_REAL: +#endif + copy_CCTK_REAL8 (varptr, varptr, + lsh[0], lsh[1], lsh[2], + imin[0], imin[1], imin[2], + imax[0], imax[1], imax[2], + ioff[0], ioff[1], ioff[2], + idir[0], idir[1], idir[2], + parity); + break; +#endif + +#ifdef CCTK_REAL16 + case CCTK_VARIABLE_REAL16: +#ifdef CCTK_REAL_PRECISION_16 + case CCTK_VARIABLE_REAL: +#endif + copy_CCTK_REAL16 (varptr, varptr, + lsh[0], lsh[1], lsh[2], + imin[0], imin[1], imin[2], + imax[0], imax[1], imax[2], + ioff[0], ioff[1], ioff[2], + idir[0], idir[1], idir[2], + parity); + break; +#endif + +#ifdef CCTK_REAL4 + case CCTK_VARIABLE_COMPLEX8: +#ifdef CCTK_COMPLEX_PRECISION_8 + case CCTK_VARIABLE_COMPLEX: +#endif + copy_CCTK_COMPLEX8 (varptr, varptr, + lsh[0], lsh[1], lsh[2], + imin[0], imin[1], imin[2], + imax[0], imax[1], imax[2], + ioff[0], ioff[1], ioff[2], + idir[0], idir[1], idir[2], + parity); + break; +#endif + +#ifdef CCTK_REAL8 + case CCTK_VARIABLE_COMPLEX16: +#ifdef CCTK_COMPLEX_PRECISION_16 + case CCTK_VARIABLE_COMPLEX: +#endif + copy_CCTK_COMPLEX16 (varptr, varptr, + lsh[0], lsh[1], lsh[2], + imin[0], imin[1], imin[2], + imax[0], imax[1], imax[2], + ioff[0], ioff[1], ioff[2], + idir[0], idir[1], idir[2], + parity); + break; +#endif + +#ifdef CCTK_REAL16 + case CCTK_VARIABLE_COMPLEX32: +#ifdef CCTK_COMPLEX_PRECISION_32 + case CCTK_VARIABLE_COMPLEX: +#endif + copy_CCTK_COMPLEX32 (varptr, varptr, + lsh[0], lsh[1], lsh[2], + imin[0], imin[1], imin[2], + imax[0], imax[1], imax[2], + ioff[0], ioff[1], ioff[2], + idir[0], idir[1], idir[2], + parity); + break; +#endif + + default: + CCTK_WARN (0, "Unsupported variable type"); + } + + } /* if cctk_bbox */ + } /* if do_reflection */ + } /* for face */ + } /* for dir */ + + /* Success */ + return 0; +} + + + +void +ReflectionSymmetry_Apply (CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + + int nvars; + CCTK_INT * restrict indices; + CCTK_INT * restrict faces; + CCTK_INT * restrict widths; + CCTK_INT * restrict tables; + int vi; + int dim; + int * restrict stencil; + int i; + int istat; + int ierr; + + if (! cctkGH) + { + CCTK_WARN (0, "Argument cctkGH is NULL"); + } + + nvars = Boundary_SelectedGVs (cctkGH, 0, NULL, NULL, NULL, NULL, NULL); + if (nvars < 0) + { + CCTK_WARN (0, "Internal error in Boundary_SelectedGVs"); + } + + indices = malloc (nvars * sizeof *indices); + if (! indices) + { + CCTK_WARN (0, "Out of memory"); + } + faces = malloc (nvars * sizeof *faces); + if (! faces) + { + CCTK_WARN (0, "Out of memory"); + } + widths = malloc (nvars * sizeof *widths); + if (! widths) + { + CCTK_WARN (0, "Out of memory"); + } + tables = malloc (nvars * sizeof *tables); + if (! tables) + { + CCTK_WARN (0, "Out of memory"); + } + + istat = Boundary_SelectedGVs + (cctkGH, nvars, indices, faces, widths, tables, 0); + if (istat != nvars) + { + CCTK_WARN (0, "Internal error in Boundary_SelectedGVs"); + } + + for (i=0; i<nvars; ++i) + { + vi = indices[i]; + if (vi < 0 || vi >= CCTK_NumVars()) + { + CCTK_WARN (0, "Illegal variable index"); + } + + if (widths[i] < 0) + { + CCTK_WARN (0, "Illegal boundary width"); + } + + dim = CCTK_GroupDimFromVarI (vi); + if (dim < 0) + { + CCTK_WARN (0, "Illegal dimension"); + } + + stencil = malloc (dim * sizeof *stencil); + if (! stencil) + { + CCTK_WARN (0, "Out of memory"); + } + ierr = CCTK_GroupnghostzonesVI (cctkGH, dim, stencil, vi); + if (ierr) + { + CCTK_WARN (0, "Internal error in CCTK_GroupnghostzonesVI"); + } + + ierr = BndReflectVI (cctkGH, stencil, vi); + if (ierr) + { + CCTK_WARN (0, "Internal error in BndReflectVI"); + } + + free (stencil); + } + + free (indices); + free (faces); + free (widths); + free (tables); +} |