diff options
author | schnetter <schnetter@082bdb00-0f4f-0410-b49e-b1835e5f2039> | 2004-06-01 16:09:35 +0000 |
---|---|---|
committer | schnetter <schnetter@082bdb00-0f4f-0410-b49e-b1835e5f2039> | 2004-06-01 16:09:35 +0000 |
commit | d02895af18ae9c3673202d6231da078cb3d74bac (patch) | |
tree | 8c8d99d5d3120685e345bf220dd416b76b86aabc | |
parent | 27ed7b1ede1863200d140d1f148a0199ac41df4f (diff) |
A thorn that applies reflection symmetry boundary conditions,
including symmetry interpolation.
git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/ReflectionSymmetry/trunk@2 082bdb00-0f4f-0410-b49e-b1835e5f2039
-rw-r--r-- | README | 10 | ||||
-rw-r--r-- | doc/documentation.tex | 144 | ||||
-rw-r--r-- | interface.ccl | 68 | ||||
-rw-r--r-- | param.ccl | 62 | ||||
-rw-r--r-- | schedule.ccl | 12 | ||||
-rw-r--r-- | src/apply.c | 627 | ||||
-rw-r--r-- | src/interpolate.c | 359 | ||||
-rw-r--r-- | src/make.code.defn | 8 | ||||
-rw-r--r-- | src/reflection.h | 29 | ||||
-rw-r--r-- | src/register.c | 66 |
10 files changed, 1385 insertions, 0 deletions
@@ -0,0 +1,10 @@ +CVS info : $Header$ + +Cactus Code Thorn ReflectionSymmetry +Thorn Author(s) : Erik Schnetter <schnetter@aei.mpg.de> +Thorn Maintainer(s) : Erik Schnetter <schnetter@aei.mpg.de> +-------------------------------------------------------------------------- + +Purpose of the thorn: + +Provide reflection symmetries, i.e. bitant, quadrant, and octant mode. diff --git a/doc/documentation.tex b/doc/documentation.tex new file mode 100644 index 0000000..fbc1f83 --- /dev/null +++ b/doc/documentation.tex @@ -0,0 +1,144 @@ +% *======================================================================* +% Cactus Thorn template for ThornGuide documentation +% Author: Ian Kelley +% Date: Sun Jun 02, 2002 +% $Header$ +% +% Thorn documentation in the latex file doc/documentation.tex +% will be included in ThornGuides built with the Cactus make system. +% The scripts employed by the make system automatically include +% pages about variables, parameters and scheduling parsed from the +% relevant thorn CCL files. +% +% This template contains guidelines which help to assure that your +% documentation will be correctly added to ThornGuides. More +% information is available in the Cactus UsersGuide. +% +% Guidelines: +% - Do not change anything before the line +% % START CACTUS THORNGUIDE", +% except for filling in the title, author, date, etc. fields. +% - Each of these fields should only be on ONE line. +% - Author names should be separated with a \\ or a comma. +% - You can define your own macros, but they must appear after +% the START CACTUS THORNGUIDE line, and must not redefine standard +% latex commands. +% - To avoid name clashes with other thorns, 'labels', 'citations', +% 'references', and 'image' names should conform to the following +% convention: +% ARRANGEMENT_THORN_LABEL +% For example, an image wave.eps in the arrangement CactusWave and +% thorn WaveToyC should be renamed to CactusWave_WaveToyC_wave.eps +% - Graphics should only be included using the graphicx package. +% More specifically, with the "\includegraphics" command. Do +% not specify any graphic file extensions in your .tex file. This +% will allow us to create a PDF version of the ThornGuide +% via pdflatex. +% - References should be included with the latex "\bibitem" command. +% - Use \begin{abstract}...\end{abstract} instead of \abstract{...} +% - Do not use \appendix, instead include any appendices you need as +% standard sections. +% - For the benefit of our Perl scripts, and for future extensions, +% please use simple latex. +% +% *======================================================================* +% +% Example of including a graphic image: +% \begin{figure}[ht] +% \begin{center} +% \includegraphics[width=6cm]{MyArrangement_MyThorn_MyFigure} +% \end{center} +% \caption{Illustration of this and that} +% \label{MyArrangement_MyThorn_MyLabel} +% \end{figure} +% +% Example of using a label: +% \label{MyArrangement_MyThorn_MyLabel} +% +% Example of a citation: +% \cite{MyArrangement_MyThorn_Author99} +% +% Example of including a reference +% \bibitem{MyArrangement_MyThorn_Author99} +% {J. Author, {\em The Title of the Book, Journal, or periodical}, 1 (1999), +% 1--16. {\tt http://www.nowhere.com/}} +% +% *======================================================================* + +% If you are using CVS use this line to give version information +% $Header$ + +\documentclass{article} + +% Use the Cactus ThornGuide style file +% (Automatically used from Cactus distribution, if you have a +% thorn without the Cactus Flesh download this from the Cactus +% homepage at www.cactuscode.org) +\usepackage{../../../../doc/latex/cactus} + +\begin{document} + +% The author of the documentation +\author{Erik Schnetter \textless schnetter@aei.mpg.de\textgreater} + +% The title of the document (not necessarily the name of the Thorn) +\title{} + +% the date your document was last changed, if your document is in CVS, +% please use: +% \date{$ $Date$ $} +\date{June 01 2004} + +\maketitle + +% Do not delete next line +% START CACTUS THORNGUIDE + +% Add all definitions used in this documentation here +% \def\mydef etc + +% Add an abstract for this thorn's documentation +\begin{abstract} + +\end{abstract} + +% The following sections are suggestive only. +% Remove them or add your own. + +\section{Introduction} + +\section{Physical System} + +\section{Numerical Implementation} + +\section{Using This Thorn} + +\subsection{Obtaining This Thorn} + +\subsection{Basic Usage} + +\subsection{Special Behaviour} + +\subsection{Interaction With Other Thorns} + +\subsection{Examples} + +\subsection{Support and Feedback} + +\section{History} + +\subsection{Thorn Source Code} + +\subsection{Thorn Documentation} + +\subsection{Acknowledgements} + + +\begin{thebibliography}{9} + +\end{thebibliography} + +% Do not delete next line +% END CACTUS THORNGUIDE + +\end{document} diff --git a/interface.ccl b/interface.ccl new file mode 100644 index 0000000..1a335bc --- /dev/null +++ b/interface.ccl @@ -0,0 +1,68 @@ +# Interface definition for thorn ReflectionSymmetry +# $Header$ + +IMPLEMENTS: ReflectionSymmetry + + + +CCTK_INT FUNCTION SymmetryRegister (CCTK_STRING IN sym_name) +REQUIRES FUNCTION SymmetryRegister + +CCTK_INT FUNCTION \ + SymmetryRegisterGrid \ + (CCTK_POINTER IN cctkGH, \ + CCTK_INT IN sym_handle, \ + CCTK_INT IN ARRAY which_faces, \ # array [N_FACES] + CCTK_INT IN ARRAY symmetry_zone_width) # array [N_FACES] +REQUIRES FUNCTION SymmetryRegisterGrid + +CCTK_INT FUNCTION \ + SymmetryRegisterGridInterpolator \ + (CCTK_POINTER IN cctkGH, \ + CCTK_INT IN sym_handle, \ + CCTK_INT CCTK_FPOINTER IN symmetry_interpolate \ + (CCTK_POINTER_TO_CONST IN cctkGH, \ + CCTK_INT IN N_dims, \ + CCTK_INT IN local_interp_handle, \ + CCTK_INT IN param_table_handle, \ + CCTK_INT IN coord_system_handle, \ + CCTK_INT IN N_interp_points, \ + CCTK_INT IN interp_coords_type, \ + CCTK_POINTER_TO_CONST ARRAY IN interp_coords, \ + CCTK_INT IN N_input_arrays, \ + CCTK_INT ARRAY IN input_array_indices, \ + CCTK_INT IN N_output_arrays, \ + CCTK_INT ARRAY IN output_array_types, \ + CCTK_POINTER ARRAY IN output_arrays, \ + CCTK_INT IN faces)) +REQUIRES FUNCTION SymmetryRegisterGridInterpolator + +CCTK_INT FUNCTION \ + SymmetryInterpolateFaces \ + (CCTK_POINTER_TO_CONST IN cctkGH, \ + CCTK_INT IN N_dims, \ + CCTK_INT IN local_interp_handle, \ + CCTK_INT IN param_table_handle, \ + CCTK_INT IN coord_system_handle, \ + CCTK_INT IN N_interp_points, \ + CCTK_INT IN interp_coords_type, \ + CCTK_POINTER_TO_CONST ARRAY IN interp_coords, \ + CCTK_INT IN N_input_arrays, \ + CCTK_INT ARRAY IN input_array_indices, \ + CCTK_INT IN N_output_arrays, \ + CCTK_INT ARRAY IN output_array_types, \ + CCTK_POINTER ARRAY IN output_arrays, \ + CCTK_INT IN faces) +REQUIRES FUNCTION SymmetryInterpolateFaces + + + +CCTK_INT FUNCTION Boundary_SelectedGVs \ + (CCTK_POINTER_TO_CONST IN GH, \ + CCTK_INT IN array_size, \ + CCTK_INT ARRAY OUT var_indicies, \ + CCTK_INT ARRAY OUT faces, \ + CCTK_INT ARRAY OUT boundary_widths, \ + CCTK_INT ARRAY OUT table_handles, \ + CCTK_STRING IN bc_name) +REQUIRES FUNCTION Boundary_SelectedGVs diff --git a/param.ccl b/param.ccl new file mode 100644 index 0000000..f6406ec --- /dev/null +++ b/param.ccl @@ -0,0 +1,62 @@ +# Parameter definitions for thorn ReflectionSymmetry +# $Header$ + +BOOLEAN verbose "Produce screen output while applying boundary conditions" +{ +} "no" + + + +BOOLEAN reflection_x "Reflection symmetry at the lower x boundary" +{ +} "no" + +BOOLEAN reflection_y "Reflection symmetry at the lower y boundary" +{ +} "no" + +BOOLEAN reflection_z "Reflection symmetry at the lower z boundary" +{ +} "no" + + + +BOOLEAN reflection_upper_x "Reflection symmetry at the upper x boundary" +{ +} "no" + +BOOLEAN reflection_upper_y "Reflection symmetry at the upper y boundary" +{ +} "no" + +BOOLEAN reflection_upper_z "Reflection symmetry at the upper z boundary" +{ +} "no" + + + +BOOLEAN avoid_origin_x "Stagger about the origin on the lower x boundary?" +{ +} "yes" + +BOOLEAN avoid_origin_y "Stagger about the origin on the lower y boundary?" +{ +} "yes" + +BOOLEAN avoid_origin_z "Stagger about the origin on the lower z boundary?" +{ +} "yes" + + + +BOOLEAN avoid_origin_upper_x "Stagger about the origin on the upper x boundary?" +{ +} "yes" + +BOOLEAN avoid_origin_upper_y "Stagger about the origin on the upper y boundary?" +{ +} "yes" + +BOOLEAN avoid_origin_upper_z "Stagger about the origin on the upper z boundary?" +{ +} "yes" diff --git a/schedule.ccl b/schedule.ccl new file mode 100644 index 0000000..8c5db0e --- /dev/null +++ b/schedule.ccl @@ -0,0 +1,12 @@ +# Schedule definitions for thorn ReflectionSymmetry +# $Header$ + +SCHEDULE ReflectionSymmetry_Register IN SymmetryRegister +{ + LANG: C +} "Register reflection symmetry boundaries" + +SCHEDULE ReflectionSymmetry_Apply IN BoundaryConditions +{ + LANG: C +} "Apply reflection symmetries" 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); +} diff --git a/src/interpolate.c b/src/interpolate.c new file mode 100644 index 0000000..8591fd7 --- /dev/null +++ b/src/interpolate.c @@ -0,0 +1,359 @@ +#include <assert.h> +#include <math.h> +#include <stdlib.h> +#include <string.h> + +#include "cctk.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_interpolate_c); + + + +CCTK_INT +ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH_, + CCTK_INT const N_dims, + CCTK_INT const local_interp_handle, + CCTK_INT const param_table_handle, + CCTK_INT const coord_system_handle, + CCTK_INT const N_interp_points, + CCTK_INT const interp_coords_type, + CCTK_POINTER_TO_CONST restrict const interp_coords[], + CCTK_INT const N_input_arrays, + CCTK_INT const input_array_indices[], + CCTK_INT const N_output_arrays, + CCTK_INT const output_array_types[], + CCTK_POINTER restrict const output_arrays[], + CCTK_INT const faces) +{ + cGH const * restrict const cctkGH = cctkGH_; + DECLARE_CCTK_PARAMETERS; + + int do_reflection[6]; + + CCTK_POINTER_TO_CONST restrict new_interp_coords[3]; + CCTK_INT newfaces; + + CCTK_INT * restrict operand_indices; + CCTK_INT * restrict operation_codes; + CCTK_INT * restrict output_array_indices; + + int dir; + + int iret; + + int m; + int n; + int d; + + int ierr; + + + + /* Get 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; + + newfaces = faces; + for (d=0; d<6; ++d) + { + if (do_reflection[d]) + { + assert (newfaces & (1 << d)); + newfaces &= ~ (1 << d); + } + } + + + + /* Fold coordinates */ + assert (interp_coords_type == CCTK_VARIABLE_REAL); + for (dir=0; dir<3; ++dir) + { + assert (! do_reflection[2*dir+1]); + + if (do_reflection[2*dir]) + { + new_interp_coords[dir] + = malloc (N_interp_points * sizeof *new_interp_coords[dir]); + assert (new_interp_coords[d]); + + for (n=0; n<N_interp_points; ++n) + { + CCTK_REAL const pos = ((CCTK_REAL const *)interp_coords[dir])[n]; + CCTK_REAL const newpos = fabs(pos); + ((CCTK_REAL *)new_interp_coords[dir])[n] = newpos; + } + } + else + { + new_interp_coords[dir] = interp_coords[dir]; + } + } + + + + /* Recursive call */ + iret = SymmetryInterpolateFaces + (cctkGH_, + N_dims, local_interp_handle, param_table_handle, coord_system_handle, + N_interp_points, interp_coords_type, new_interp_coords, + N_input_arrays, input_array_indices, + N_output_arrays, output_array_types, output_arrays, + newfaces); + + + + /* Free coordinates */ + for (dir=0; dir<3; ++dir) + { + if (do_reflection[2*dir]) + { + free (new_interp_coords[dir]); + } + } + + + /* Find output variable indices */ + operand_indices = malloc (N_output_arrays * sizeof *operand_indices); + assert (operand_indices); + ierr = Util_TableGetIntArray + (param_table_handle, N_output_arrays, operand_indices, "operand_indices"); + if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY) { + assert (N_output_arrays == N_input_arrays); + for (m=0; m<N_output_arrays; ++m) { + operand_indices[m] = m; /* set output index to input index */ + } + } else { + assert (ierr == N_output_arrays); + } + + operation_codes = malloc (N_output_arrays * sizeof *operation_codes); + assert (operation_codes); + ierr = Util_TableGetIntArray + (param_table_handle, N_output_arrays, operation_codes, "operation_codes"); + if (ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY) { + assert (N_output_arrays == N_input_arrays); + for (m=0; m<N_output_arrays; ++m) { + operation_codes[m] = 0; /* do not take derivatives */ + } + } else { + assert (ierr == N_output_arrays); + } + + output_array_indices + = malloc (N_output_arrays * sizeof *output_array_indices); + assert (output_array_indices); + for (m=0; m<N_output_arrays; ++m) { + assert (operand_indices[m]>=0 && operand_indices[m]<N_input_arrays); + output_array_indices[m] = input_array_indices[operand_indices[m]]; + assert (output_array_indices[m]>=0 + && output_array_indices[m]<CCTK_NumVars()); + } + + + + /* Unfold tensor types */ + for (m=0; m<N_output_arrays; ++m) + { + + int vi, gi; + int numvars, firstvar; + char * groupname; + + int table; + char tensortypealias[1000]; + enum tensortype { UNKNOWN, SCALAR, VECTOR, TENSOR }; + enum tensortype ttype; + int tcomponent; + + int parities[3]; + int check_dir[3]; + int needs_checking; + + + + vi = output_array_indices[m]; + assert (vi>=0 && vi<CCTK_NumVars()); + gi = CCTK_GroupIndexFromVarI (vi); + assert (gi>=0 && gi<CCTK_NumGroups()); + numvars = CCTK_NumVarsInGroupI(gi); + assert (numvars>0); + firstvar = CCTK_FirstVarIndexI(gi); + assert (firstvar>=0); + + + + /* 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); + } + + + + /* Calculate parities */ + parities[0] = parities[1] = parities[2] = +1; + switch (ttype) + { + case SCALAR: + /* do nothing */ + break; + case VECTOR: + parities[tcomponent] = -1; + break; + case TENSOR: + switch (tcomponent) + { + case 0: break; + case 1: parities[0] = parities[1] = -1; break; + case 2: parities[0] = parities[2] = -1; break; + case 3: break; + case 4: parities[1] = parities[2] = -1; break; + case 5: break; + default: assert (0); + } + break; + default: + assert (0); + } + + + + /* Take derivatives into account */ + { + int code = operation_codes[m]; + while (code) { + const int d = code % 10 - 1; + code /= 10; + assert (d>=0 && d<3); + parities[d] *= -1; + } + } + + + + /* Are there negative parities? */ + needs_checking = 0; + for (dir=0; dir<3; ++dir) + { + check_dir[d] = do_reflection[2*dir] && parities[dir] < 0; + needs_checking |= check_dir[d]; + } + + + + /* Loop over all points and unfold */ + if (needs_checking) + { + for (n=0; n<N_interp_points; ++n) + { + int parity = +1; + /* Is the point outside the domain? */ + for (dir=0; dir<3; ++dir) + { + if (check_dir[dir]) + { + CCTK_REAL const pos = ((CCTK_REAL const *)interp_coords[dir])[n]; + if (pos < 0) + { + /* Reflect the tensor component */ + parity *= -1; + } + } + } + ((CCTK_REAL *)output_arrays[m])[n] *= parity; + } + } + + } /* for m */ + + + + /* Free output variable indices */ + free (operand_indices); + free (operation_codes); + free (output_array_indices); + + + + /* Return */ + return iret; +} diff --git a/src/make.code.defn b/src/make.code.defn new file mode 100644 index 0000000..24a04f5 --- /dev/null +++ b/src/make.code.defn @@ -0,0 +1,8 @@ +# Main make.code.defn file for thorn ReflectionSymmetry +# $Header$ + +# Source files in this directory +SRCS = apply.c interpolate.c register.c + +# Subdirectories containing source files +SUBDIRS = diff --git a/src/reflection.h b/src/reflection.h new file mode 100644 index 0000000..fea4729 --- /dev/null +++ b/src/reflection.h @@ -0,0 +1,29 @@ +#include "cctk.h" +#include "cctk_Arguments.h" + +#ifndef REFLECTION_H +#define REFLECTION_H + +CCTK_INT +ReflectionSymmetry_Interpolate (CCTK_POINTER_TO_CONST restrict const cctkGH, + CCTK_INT const N_dims, + CCTK_INT const local_interp_handle, + CCTK_INT const param_table_handle, + CCTK_INT const coord_system_handle, + CCTK_INT const N_interp_points, + CCTK_INT const interp_coords_type, + CCTK_POINTER_TO_CONST restrict const interp_coords[], + CCTK_INT const N_input_arrays, + CCTK_INT const input_array_indices[], + CCTK_INT const N_output_arrays, + CCTK_INT const output_array_types[], + CCTK_POINTER restrict const output_arrays[], + CCTK_INT const faces); + +void +ReflectionSymmetry_Register (CCTK_ARGUMENTS); + +void +ReflectionSymmetry_Apply (CCTK_ARGUMENTS); + +#endif /* REFLECTION_H */ diff --git a/src/register.c b/src/register.c new file mode 100644 index 0000000..2ae4b58 --- /dev/null +++ b/src/register.c @@ -0,0 +1,66 @@ +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" + +#include "reflection.h" + + + +static const char rcsid[] = "$Header$"; +CCTK_FILEVERSION(AEIDevelopment_ReflectionSymmetry_register_c); + + + +void +ReflectionSymmetry_Register (CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + int do_reflection[6]; + CCTK_INT handle; + CCTK_INT faces[6]; + CCTK_INT width[6]; + int f; + CCTK_INT ierr; + + 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; + + for (f=0; f<6; ++f) + { + if (do_reflection[f]) + { + faces[f] = 1; + width[f] = cctk_nghostzones[f/2]; + } + else + { + faces[f] = 0; + width[f] = 0; + } + } + + handle = SymmetryRegister ("reflection_symmetry"); + if (handle < 0) + { + CCTK_WARN (0, "Could not register symmetry boundary condition"); + } + + ierr = SymmetryRegisterGrid (cctkGH, handle, faces, width); + if (ierr < 0) + { + CCTK_WARN (0, "Could not register the symmetry boundaries -- probably some other thorn has already registered the same boundary faces for a different symmetry"); + } + + ierr = SymmetryRegisterGridInterpolator + (cctkGH, handle, ReflectionSymmetry_Interpolate); + if (ierr < 0) + { + CCTK_WARN (0, "Could not register the symmetry interpolator"); + } +} |