From c895e5b15a407c5fc3f9109958e41a48844942f7 Mon Sep 17 00:00:00 2001 From: schnetter Date: Fri, 28 May 2004 15:44:34 +0000 Subject: Provide infrastructure to take symmetries into account for global interpolations. git-svn-id: http://svn.cactuscode.org/arrangements/CactusBase/SymBase/trunk@11 906471b6-c639-44d1-9ea0-3e3d6879f074 --- doc/documentation.tex | 183 +++++++++++++++++++++++- interface.ccl | 95 ++++++++++++- param.ccl | 4 + schedule.ccl | 11 +- src/Faces.c | 165 ++++++++++++++++++++++ src/Interpolation.c | 382 ++++++++++++++++++++++++++++++++++++++++++++++++++ src/Startup.c | 17 +++ src/Statistics.c | 134 ++++++++++++++++++ src/SymBase.h | 87 +++++++++++- src/make.code.defn | 2 +- 10 files changed, 1068 insertions(+), 12 deletions(-) create mode 100644 src/Interpolation.c create mode 100644 src/Statistics.c diff --git a/doc/documentation.tex b/doc/documentation.tex index c37694e..8e1808b 100644 --- a/doc/documentation.tex +++ b/doc/documentation.tex @@ -45,11 +45,11 @@ % % 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} +% \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: @@ -296,6 +296,179 @@ end do +\section{Symmetry Interpolation} + +When the computational domain has symmetries, then it is often very +convenient to be able to interpolate at points that are not present on +the actual computational grid, but can be mapped into the grid through +the symmetries. Thorn SymBase provides a mechanism by which symmetry +conditions can register routines that handle this mapping when a +global interpolator is called. Additionally, the driver has to be +aware that is calls thorn SymBase's mapping routine before it actually +interpolates. The whole mechanism is transparent for the user. + +The mechanism by which the grid points are mapped into the domain +works as follows: +\begin{enumerate} +\item The user calls \texttt{CCTK\_Interp} with a list of coordinates. +\item The flesh forwards this call to the driver. +\item The driver calls SymBase's aliased function. + \texttt{SymmetryInterpolate}, passing along all arguments. +\item SymBase sets a flag for all boundaries that have a symmetry + condition associated with it, and then calls + \texttt{SymmetryInterpolateFaces}, passing along all arguments. + This is the beginning of a chain of recursive calls. +\item \texttt{SymmetryInterpolateFaces} checks whether any faces are + flagged. +\item If no faces are flagged, SymBase calls the driver's aliased + function \texttt{DriverInterpolate}, which performs the actual + interpolation. This ends the chain of recursive calls. +\item If there are faces with symmetry conditions flagged, SymBase + chooses one such face, and then calls this face's symmetry + condition's ``symmetry interpolation'' routine, passing along all + arguments. +\item The ``symmetry interpolation'' routine maps the coordinates into + the domain by applying the symmetry condition for this face. It + then removes the flag for the corresponding face, and calls + \texttt{SymmetryInterpolateFaces}, passing along the arguments with + the changed interpolation locations. +\item After the actual interpolation has happened in the driver, the + recursive call will return. The ``symmetry interpolation'' routine + then examines the tensor types of the interpolated quantities and + un-maps the values back onto their original locations. That is, + e.g., after a reflection on the lower $x$-boundary, $x$-components + of vectors need their sign changed. +\item The chain of recursive calls unravels until the call to + \texttt{CCTK\_Interp} returns. +\end{enumerate} + +This mechanism has thus four players: +\begin{itemize} +\item The driver forwards the call to SymBase so that the list of + interpolation points can be mapped into the domain. +\item Thorn SymBase controls which symmetry conditions perform this + mapping on which faces. +\item Each symmetry boundary condition has to register a ``symmetry + interpolation'' routine that first maps the points into the domain, + then calls SymBase recursively, and finally corrects the tensor + types of the interpolated quantities. +\item Finally, the user calls \texttt{CCTK\_Interp} as before, and + everything happens transparently for him. +\end{itemize} + + + +\subsection{Driver Interaction} + +The driver has to call SymBase's aliased function +\texttt{SymmetryInterpolate}, and has to provide an aliased function +\texttt{DriverInterpolate}. Both functions have essentially the same +prototype as \texttt{CCTK\_Interp}, which are: + +\begin{verbatim} +CCTK_INT FUNCTION + SymmetryInterpolate + (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) +\end{verbatim} + +\begin{verbatim} +CCTK_INT FUNCTION + DriverInterpolate + (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) +\end{verbatim} + + + +\subsection{Interaction With Symmetry Conditions} + +The symmetry conditions have to register their ``symmetry +interpolation'' routines by calling SymBase's aliased function +\texttt{SymmetryRegisterGridInterpolator}. The ``symmetry +interpolation'' routine must use C linkage and must have the prototype + +\begin{verbatim} +CCTK_INT 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) +\end{verbatim} + +These arguments are essentially the same as those for +\texttt{CCTK\_Interp}, except that the bit field \texttt{faces} flags +those faces that still their symmetry boundary condition applied to +the interpolation points. The aliased function +\texttt{SymmetryRegisterGridInterpolator} has the prototype + +\begin{verbatim} +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)) +\end{verbatim} + +which takes a function pointer to the aforementioned ``symmetry +interpolation'' routine, while \texttt{sym\_handle} specifies which +symmetry condition this routine is for. This handle must have been +obtained from \texttt{SymmetryRegister}. + +\begin{quote} + The routine \texttt{SymmetryRegisterGridInterpolator} must be called + \emph{after} the symmetry faces have been selected by the call to + \texttt{SymmetryRegisterGrid}. +\end{quote} + + + % Do not delete next line % END CACTUS THORNGUIDE diff --git a/interface.ccl b/interface.ccl index 2c07bd8..c8ddacc 100644 --- a/interface.ccl +++ b/interface.ccl @@ -28,10 +28,10 @@ PROVIDES FUNCTION SymmetryNameOfHandle \ 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] + (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] PROVIDES FUNCTION SymmetryRegisterGrid \ WITH SymBase_SymmetryRegisterGrid \ LANGUAGE C @@ -60,6 +60,33 @@ PROVIDES FUNCTION SymmetryRegisterGN \ +# Register a symmetry interpolator: + +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)) +PROVIDES FUNCTION SymmetryRegisterGridInterpolator \ + WITH SymBase_SymmetryRegisterGridInterpolator \ + LANGUAGE C + + + # Get the symmetry table handle for a grid or grid array: CCTK_INT FUNCTION \ @@ -83,3 +110,63 @@ CCTK_INT FUNCTION \ PROVIDES FUNCTION SymmetryTableHandleForGN \ WITH SymBase_SymmetryTableHandleForGN \ LANGUAGE C + + + +# Interpolation + +CCTK_INT FUNCTION \ + SymmetryInterpolate \ + (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) +PROVIDES FUNCTION SymmetryInterpolate \ + WITH SymBase_SymmetryInterpolate \ + LANGUAGE C + +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) +PROVIDES FUNCTION SymmetryInterpolateFaces \ + WITH SymBase_SymmetryInterpolateFaces \ + LANGUAGE C + +CCTK_INT FUNCTION \ + DriverInterpolate \ + (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) +USES FUNCTION DriverInterpolate diff --git a/param.ccl b/param.ccl index e82ab1e..7803536 100644 --- a/param.ccl +++ b/param.ccl @@ -1,2 +1,6 @@ # Parameter definitions for thorn SymBase # $Header$ + +BOOLEAN verbose "Output symmetry boundary face descriptions after registration" +{ +} "yes" diff --git a/schedule.ccl b/schedule.ccl index f3d7f54..350d1f9 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -6,6 +6,15 @@ SCHEDULE SymBase_Startup AT CCTK_STARTUP LANG: C } "Register GH Extension for SymBase" -SCHEDULE GROUP SymmetryRegister AT CCTK_WRAGH +SCHEDULE GROUP SymBase_Wrapper AT CCTK_WRAGH +{ +} "Wrapper group for SymBase" + +SCHEDULE GROUP SymmetryRegister IN SymBase_Wrapper { } "Register your symmetries here" + +SCHEDULE SymBase_Statistics IN SymBase_Wrapper AFTER SymmetryRegister +{ + LANG: C +} "Print symmetry boundary face descriptions" diff --git a/src/Faces.c b/src/Faces.c index c1a668f..844dfc2 100644 --- a/src/Faces.c +++ b/src/Faces.c @@ -360,3 +360,168 @@ SymBase_SymmetryRegisterGN (CCTK_POINTER const cctkGH_, return SymBase_SymmetryRegisterGI (cctkGH_, sym_handle, which_faces, new_symmetry_zone_width, group_index); } + + + +/*@@ + @routine SymBase_SymmetryRegisterInterpolatorFaces + @author Erik Schnetter + @date 2004-05-25 + @desc + Register a symmetry interpolator for certain faces + @enddesc + @var sym_table + @vtype CCTK_INT + @vdesc Table which describes the grid or grid array + @vio in + @endvar + @var group_dim + @vtype CCTK_INT + @vdesc Dimension of the grid or grid array + @vio in + @endvar + @var sym_handle + @vtype CCTK_INT + @vdesc Symmetry handle + @vio in + @endvar + @var symmetry_interpolate + @vtype CCTK_FPOINTER + @vdesc Routine that applies symmetries to the interpolation points + and then calls back + (may be NULL) + @vio in + @endvar + @returntype CCTK_INT + @returndesc + status code + 0 for success + -1 if sym_table has an illegal value + -9 if group_dim has an illegal value + -2 if sym_handle has an illegal value + -10 symmetry_interpolate is NULL + @endreturndesc +@@*/ + +CCTK_INT +SymBase_SymmetryRegisterInterpolatorFaces (CCTK_INT const sym_table, + CCTK_INT const group_dim, + CCTK_INT const sym_handle, + SymmetryInterpolateFPointer const new_symmetry_interpolate) +{ + CCTK_INT symmetry_handle[100]; + CCTK_FPOINTER symmetry_interpolate[100]; + int face; + int ierr; + + /* Check arguments */ + if (sym_table < 0) + { + return -1; /* illegal argument */ + } + if (group_dim < 0) + { + return -9; /* illegal argument */ + } + if (sym_handle < 0 || sym_handle >= SymBase_num_symmetries) + { + return -2; /* illegal argument */ + } + if (! new_symmetry_interpolate) + { + return -10; + } + + /* Get table entries */ + if (2 * group_dim > 100) + { + CCTK_WARN (0, "internal error"); + } + ierr = Util_TableGetIntArray + (sym_table, 2 * group_dim, symmetry_handle, "symmetry_handle"); + if (ierr != 2 * group_dim) + { + CCTK_WARN (0, "internal error"); + } + ierr = Util_TableGetFPointerArray + (sym_table, 2 * group_dim, symmetry_interpolate, "symmetry_interpolate"); if (ierr != 2 * group_dim) { + CCTK_WARN (0, "internal error"); + } + + /* Update table entries */ + for (face = 0; face < 2 * group_dim; ++face) + { + if (symmetry_handle[face] == sym_handle) + { + symmetry_interpolate[face] = (CCTK_FPOINTER)new_symmetry_interpolate; + } + } + + /* Set table entries */ + ierr = Util_TableSetFPointerArray + (sym_table, 2 * group_dim, symmetry_interpolate, "symmetry_interpolate"); + if (ierr != 1) + { + CCTK_WARN (0, "internal error"); + } + + return 0; +} + + + +/*@@ + @routine SymBase_SymmetryRegisterGridInterpolator + @author Erik Schnetter + @date 2004-05-25 + @desc + Register a symmetry interpolator + for certain faces of the grid hierarchy + @enddesc + @var cctkGH + @vtype cGH * + @vdesc The grid hierarchy + @vio in + @endvar + @var sym_handle + @vtype CCTK_INT + @vdesc Symmetry handle + @vio in + @endvar + @var symmetry_interpolate + @vtype CCTK_FPOINTER + @vdesc Routine that applies symmetries to the interpolation points + and then calls back + @vio in + @endvar + @returntype CCTK_INT + @returndesc + status code + 0 for success + Error codes of SymBase_SymmetryRegisterInterpolatorFaces + @endreturndesc +@@*/ + +CCTK_INT +SymBase_SymmetryRegisterGridInterpolator (CCTK_POINTER const cctkGH_, + CCTK_INT const sym_handle, + SymmetryInterpolateFPointer const new_symmetry_interpolate) +{ + cGH const *const cctkGH = cctkGH_; + struct SymBase const *symdata; + + if (!cctkGH) + { + CCTK_WARN (0, "internal error"); + } + + symdata = CCTK_GHExtension (cctkGH, "SymBase"); + if (!symdata) + { + CCTK_WARN (0, "internal error"); + } + + return SymBase_SymmetryRegisterInterpolatorFaces + (symdata->sym_table, cctkGH->cctk_dim, sym_handle, + new_symmetry_interpolate); +} diff --git a/src/Interpolation.c b/src/Interpolation.c new file mode 100644 index 0000000..6873ebe --- /dev/null +++ b/src/Interpolation.c @@ -0,0 +1,382 @@ +/*@@ + @file Interpolation.c + @author Erik Schnetter + @date 2004-05-11 + @desc + Apply symmetry conditions during interpolation + @version $Header$ + @enddesc +@@*/ + +#include + +#include "cctk.h" +#include "cctk_Arguments.h" + +#include "util_ErrorCodes.h" +#include "util_Table.h" + +#include "SymBase.h" + + + +/* the rcs ID and its dummy function to use it */ +static const char *const rcsid = "$Header$"; +CCTK_FILEVERSION (CactusBase_SymBase_Interpolation_c); + + + +/*@@ + @routine SymBase_SymmetryInterpolate + @author Erik Schnetter + @date 2004-05-11 + @desc + Adjust the coordinates of the interpolation points, + call the driver's interpolation function, + and then adjust the interpolated tensor components. + All arguments are equivalent to the CCTK_InterpGridArrays + function call. + @enddesc + @var cctkGH + @vtype CCTK_POINTER_TO_CONST + @vdesc Grid hierarchy + @vio in + @endvar + @var N_dims + @vtype CCTK_INT + @vdesc same as for the CCTK_InterpGridArrays call + @vio in + @endvar + @var local_interp_handle + @vtype CCTK_INT + @vdesc same as for the CCTK_InterpGridArrays call + @vio in + @endvar + @var param_table_handle + @vtype CCTK_INT + @vdesc same as for the CCTK_InterpGridArrays call + @vio in + @endvar + @var coord_system_handle + @vtype CCTK_INT + @vdesc same as for the CCTK_InterpGridArrays call + @vio in + @endvar + @var N_interp_points + @vtype CCTK_INT + @vdesc same as for the CCTK_InterpGridArrays call + @vio in + @endvar + @var interp_coords_type + @vtype CCTK_INT + @vdesc same as for the CCTK_InterpGridArrays call + @vio in + @endvar + @var interp_coords + @vtype CCTK_POINTER_TO_CONST[] + @vdesc same as for the CCTK_InterpGridArrays call + @vio in + @endvar + @var N_input_arrays + @vtype CCTK_INT + @vdesc same as for the CCTK_InterpGridArrays call + @vio in + @endvar + @var input_array_indices + @vtype CCTK_INT[] + @vdesc same as for the CCTK_InterpGridArrays call + @vio in + @endvar + @var N_output_arrays + @vtype CCTK_INT + @vdesc same as for the CCTK_InterpGridArrays call + @vio in + @endvar + @var output_array_types + @vtype CCTK_INT[] + @vdesc same as for the CCTK_InterpGridArrays call + @vio in + @endvar + @var output_arrays + @vtype CCTK_POINTER[] + @vdesc same as for the CCTK_InterpGridArrays call + @vio in + @endvar + @returntype CCTK_INT + @returndesc + status code + 0 for success + -91 GH is NULL + -92 N_dims is not equal to cctkGH->cctk_dim + @endreturndesc +@@*/ + +CCTK_INT +SymBase_SymmetryInterpolate (CCTK_POINTER_TO_CONST 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 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 const output_arrays[]) +{ + cGH const * restrict const cctkGH = cctkGH_; + CCTK_INT sym_table; + CCTK_FPOINTER symmetry_interpolate[100]; + CCTK_INT faces; + int f; + int ierr; + + /* Check arguments */ + if (! cctkGH) + { + CCTK_WARN (1, "cctkGH is NULL"); + return -91; + } + + if (N_dims != cctkGH->cctk_dim) + { + CCTK_WARN (1, "The number of dimensions is not equal to the GH's number of dimensions"); + return -92; + } + + /* Get table */ + sym_table = SymmetryTableHandleForGrid (cctkGH); + if (sym_table < 0) CCTK_WARN (0, "internal error"); + + /* Get table entries */ + if (2 * cctkGH->cctk_dim > 100) + { + CCTK_WARN (0, "internal error"); + } + ierr = Util_TableGetFPointerArray + (sym_table, 2 * cctkGH->cctk_dim, + symmetry_interpolate, "symmetry_interpolate"); + if (ierr != 2 * cctkGH->cctk_dim) + { + CCTK_WARN (0, "internal error"); + } + + /* Find all faces that need symmetries applied */ + faces = 0; + for (f=0; f<2*cctkGH->cctk_dim; ++f) { + if (symmetry_interpolate[f]) { + faces |= (1 << f); + } + } + + /* Forward the call */ + return SymBase_SymmetryInterpolateFaces + (cctkGH, N_dims, + local_interp_handle, param_table_handle, coord_system_handle, + N_interp_points, interp_coords_type, interp_coords, + N_input_arrays, input_array_indices, + N_output_arrays, output_array_types, output_arrays, + faces); +} + + + +/*@@ + @routine SymBase_SymmetryInterpolateFaces + @author Erik Schnetter + @date 2004-05-11 + @desc + Adjust the coordinates of the interpolation points, + call the driver's interpolation function, + and then adjust the interpolated tensor components, + but only for the selected faces. + When no faces are selected, forward the call to the + driver's interpolator throug the aliased function + "DriverInterpolate". + All arguments except "faces" are equivalent to the + CCTK_InterpGridArrays function call. + @enddesc + @var cctkGH + @vtype CCTK_POINTER_TO_CONST + @vdesc Grid hierarchy + @vio in + @endvar + @var N_dims + @vtype CCTK_INT + @vdesc same as for the CCTK_InterpGridArrays call + @vio in + @endvar + @var local_interp_handle + @vtype CCTK_INT + @vdesc same as for the CCTK_InterpGridArrays call + @vio in + @endvar + @var param_table_handle + @vtype CCTK_INT + @vdesc same as for the CCTK_InterpGridArrays call + @vio in + @endvar + @var coord_system_handle + @vtype CCTK_INT + @vdesc same as for the CCTK_InterpGridArrays call + @vio in + @endvar + @var N_interp_points + @vtype CCTK_INT + @vdesc same as for the CCTK_InterpGridArrays call + @vio in + @endvar + @var interp_coords_type + @vtype CCTK_INT + @vdesc same as for the CCTK_InterpGridArrays call + @vio in + @endvar + @var interp_coords + @vtype CCTK_POINTER_TO_CONST[] + @vdesc same as for the CCTK_InterpGridArrays call + @vio in + @endvar + @var N_input_arrays + @vtype CCTK_INT + @vdesc same as for the CCTK_InterpGridArrays call + @vio in + @endvar + @var input_array_indices + @vtype CCTK_INT[] + @vdesc same as for the CCTK_InterpGridArrays call + @vio in + @endvar + @var N_output_arrays + @vtype CCTK_INT + @vdesc same as for the CCTK_InterpGridArrays call + @vio in + @endvar + @var output_array_types + @vtype CCTK_INT[] + @vdesc same as for the CCTK_InterpGridArrays call + @vio in + @endvar + @var output_arrays + @vtype CCTK_POINTER[] + @vdesc same as for the CCTK_InterpGridArrays call + @vio in + @endvar + @var faces + @vtype CCTK_INT + @vdesc Bit mask selecting which faces' symmetries should be treated. + @vio in + @endvar + @returntype CCTK_INT + @returndesc + status code + 0 for success + -91 GH is NULL + -92 N_dims is not equal to cctkGH->cctk_dim + -93 faces is NULL + or the return value from DriverInterpolate + @endreturndesc +@@*/ + +CCTK_INT +SymBase_SymmetryInterpolateFaces (CCTK_POINTER_TO_CONST 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 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 const output_arrays[], + CCTK_INT const faces) +{ + cGH const * restrict const cctkGH = cctkGH_; + CCTK_INT sym_table; + CCTK_FPOINTER symmetry_interpolate[100]; + int face; + int ierr; + + /* Check arguments */ + + if (! cctkGH) + { + CCTK_WARN (1, "cctkGH is NULL"); + return -91; + } + + if (N_dims != cctkGH->cctk_dim) + { + CCTK_WARN (1, "The number of dimensions is not equal to the GH's number of dimensions"); + return -92; + } + + /* Get table */ + sym_table = SymmetryTableHandleForGrid (cctkGH); + if (sym_table < 0) CCTK_WARN (0, "internal error"); + + /* Get table entries */ + if (2 * cctkGH->cctk_dim > 100) + { + CCTK_WARN (0, "internal error"); + } + ierr = Util_TableGetFPointerArray + (sym_table, 2 * cctkGH->cctk_dim, + symmetry_interpolate, "symmetry_interpolate"); + if (ierr != 2 * cctkGH->cctk_dim) + { + CCTK_WARN (0, "internal error"); + } + + /* Find a face that needs a symmetry applied */ + for (face=0; face<2*N_dims; ++face) + { + if (faces & (1 << face)) + { + break; + } + } + + if (face == 2*N_dims) + { + /* All faces are done */ + + /* Call the real interpolator */ + ierr = CCTK_IsFunctionAliased ("DriverInterpolate"); + if (! ierr) + { + CCTK_WARN (0, "The aliased function \"DriverInterpolate\" has not been provided"); + } + ierr = DriverInterpolate + (cctkGH, N_dims, + local_interp_handle, param_table_handle, coord_system_handle, + N_interp_points, interp_coords_type, interp_coords, + N_input_arrays, input_array_indices, + N_output_arrays, output_array_types, output_arrays); + + } + else + { + /* Treat the face */ + + /* Recursive call to a symmetry condition */ + if (! symmetry_interpolate[face]) + { + CCTK_WARN (0, "internal error"); + } + ierr = ((SymmetryInterpolateFPointer)symmetry_interpolate[face]) + (cctkGH, N_dims, + local_interp_handle, param_table_handle, coord_system_handle, + N_interp_points, interp_coords_type, interp_coords, + N_input_arrays, input_array_indices, + N_output_arrays, output_array_types, output_arrays, + faces); + + } + + return ierr; +} diff --git a/src/Startup.c b/src/Startup.c index 298d0ee..886b4fa 100644 --- a/src/Startup.c +++ b/src/Startup.c @@ -95,6 +95,7 @@ SymBase_Setup (tFleshConfig * const config, struct SymBase *symdata; CCTK_INT symmetry_handle[100]; CCTK_INT symmetry_zone_width[100]; + CCTK_FPOINTER symmetry_interpolate[100]; int group; int face; int ierr; @@ -111,6 +112,7 @@ SymBase_Setup (tFleshConfig * const config, { symmetry_handle[face] = -1; symmetry_zone_width[face] = 0; + symmetry_interpolate[face] = NULL; } /* Create grid symmetry table */ @@ -139,6 +141,13 @@ SymBase_Setup (tFleshConfig * const config, { CCTK_WARN (0, "Internal Error: Failed to create symmetry_zone_width array"); } + ierr = Util_TableSetFPointerArray + (symdata->sym_table, + 2 * cctkGH->cctk_dim, symmetry_interpolate, "symmetry_interpolate"); + if (ierr) + { + CCTK_WARN (0, "Internal Error: Failed to create symmetry_fold_points array"); + } /* Create grid array symmetry tables */ symdata->array_sym_tables @@ -190,6 +199,14 @@ SymBase_Setup (tFleshConfig * const config, { CCTK_WARN (0, "Internal Error: Failed to create symmetry_zone_width array"); } + ierr = Util_TableSetFPointerArray + (symdata->array_sym_tables[group], + 2 * CCTK_GroupDimI (group), symmetry_interpolate, + "symmetry_interpolate"); + if (ierr) + { + CCTK_WARN (0, "Internal Error: Failed to create symmetry_fold_points array"); + } break; default: diff --git a/src/Statistics.c b/src/Statistics.c new file mode 100644 index 0000000..e31cac6 --- /dev/null +++ b/src/Statistics.c @@ -0,0 +1,134 @@ +/*@@ + @file Statistics.c + @author Erik Schnetter + @date 2004-05-25 + @desc + Provide information about the registered symmetry conditions + @version $Header$ + @enddesc +@@*/ + +#include + +#include "cctk.h" +#include "cctk_Parameters.h" + +#include "util_Table.h" + +#include "SymBase.h" + + + +/* the rcs ID and its dummy function to use it */ +static const char *const rcsid = "$Header$"; +CCTK_FILEVERSION (CactusBase_SymBase_Statistics_c); + + + +/*@@ + @routine SymBase_Statistics + @author Erik Schnetter + @date 2004-05-25 + @desc + If verbose, print the symmetry boundary conditions for each + face. + Warn if symmetries do not have interpolation operators + associated with them. + @enddesc +@@*/ + +void +SymBase_Statistics (CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + CCTK_INT sym_table; + CCTK_INT symmetry_handle[100]; + CCTK_FPOINTER symmetry_interpolate[100]; + char const * symmetry_name; + int face; + int ierr; + + + + /* Get table */ + sym_table = SymmetryTableHandleForGrid (cctkGH); + if (sym_table < 0) CCTK_WARN (0, "internal error"); + + /* Get table entries */ + if (2 * cctkGH->cctk_dim > 100) + { + CCTK_WARN (0, "internal error"); + } + ierr = Util_TableGetIntArray + (sym_table, 2 * cctkGH->cctk_dim, symmetry_handle, "symmetry_handle"); + if (ierr != 2 * cctkGH->cctk_dim) + { + CCTK_WARN (0, "internal error"); + } + ierr = Util_TableGetFPointerArray + (sym_table, 2 * cctkGH->cctk_dim, + symmetry_interpolate, "symmetry_interpolate"); + if (ierr != 2 * cctkGH->cctk_dim) + { + CCTK_WARN (0, "internal error"); + } + + + + /* Print information about the registered symmetry faces, if desired */ + + if (verbose) + { + + /* Loop over all faces */ + for (face=0; face<2*cctkGH->cctk_dim; ++face) { + if (symmetry_handle[face] >= 0) { + symmetry_name = SymBase_SymmetryNameOfHandle (symmetry_handle[face]); + if (! symmetry_name) + { + CCTK_WARN (0, "internal error"); + } + if (face >= 6) + { + CCTK_WARN (0, "only 3 dimensions are supported so far"); + } + CCTK_VInfo (CCTK_THORNSTRING, + "Symmetry on %s %c-face: %s", + face % 2 == 0 ? "lower" : "upper", + "xyz"[face/2], + symmetry_name); + } + } + + } + + + + /* Warn about symmetries without interpolators */ + + /* Loop over all registered symmetries */ + for (face=0; face<2*cctkGH->cctk_dim; ++face) { + if (symmetry_handle[face] >= 0) { + if (symmetry_interpolate[face] == NULL) + { + symmetry_name = SymBase_SymmetryNameOfHandle (symmetry_handle[face]); + if (! symmetry_name) + { + CCTK_WARN (0, "internal error"); + } + if (face >= 6) + { + CCTK_WARN (0, "only 3 dimensions are supported so far"); + } + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "The symmetry %s on the %s %c-face has not registered a symmetry interpolator", + symmetry_name, + face % 2 == 0 ? "lower" : "upper", + "xyz"[face/2]); + } + } + } + +} diff --git a/src/SymBase.h b/src/SymBase.h index 397c0a1..b4e118e 100644 --- a/src/SymBase.h +++ b/src/SymBase.h @@ -14,21 +14,60 @@ #include "cctk.h" +#include "cctk_Arguments.h" /* SymBase's GH extension */ struct SymBase { + /* The table handles below identify tables that have the following + entries: + + CCTK_INT symmetry_handle[2*dim] + CCTK_INT symmetry_zone_width[2*dim] + CCTK_FPOINTER symmetry_interpolate[2*dim] + + dim is here the dimension of the group containing the variable, + i.e. there is one entry per face. The true type of the + symmetry_interpolate entries is SymmetryInterpolateFPointer. + + (Given that user code is not allowed to modify these tables, and + that this structure is defined here in this header file, there is + no real need to use tables for this. It would be possible to put + these data directly into this structure, thereby simplifying user + code and not make anything less flexible.) + */ + /* Grid symmetry table handle */ int sym_table; - /* Grid array symmetry table handles */ + /* Grid array symmetry table handles, one handle per grid + variable */ int *array_sym_tables; }; +typedef +CCTK_INT +(* SymmetryInterpolateFPointer) (CCTK_POINTER_TO_CONST 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 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 const output_arrays[], + CCTK_INT const faces); + + + /* Number of registered symmetries */ extern size_t SymBase_num_symmetries; @@ -73,6 +112,16 @@ SymBase_SymmetryRegisterGN (CCTK_POINTER const cctkGH_, CCTK_INT const *const new_symmetry_zone_width, CCTK_STRING const group_name); +CCTK_INT +SymBase_SymmetryRegisterInterpolatorFaces (CCTK_INT const sym_table, + CCTK_INT const group_dim, + CCTK_INT const sym_handle, + SymmetryInterpolateFPointer const new_symmetry_interpolate); +CCTK_INT +SymBase_SymmetryRegisterGridInterpolator (CCTK_POINTER const cctkGH_, + CCTK_INT const sym_handle, + SymmetryInterpolateFPointer const new_symmetry_interpolate); + /* Table.c */ CCTK_INT SymBase_SymmetryTableHandleForGrid (CCTK_POINTER_TO_CONST const cctkGH_); @@ -83,6 +132,42 @@ CCTK_INT SymBase_SymmetryTableHandleForGN (CCTK_POINTER_TO_CONST const cctkGH_, CCTK_STRING const group_name); +/* Interpolate.c */ +CCTK_INT +SymBase_SymmetryInterpolate (CCTK_POINTER_TO_CONST 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 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 const output_arrays[]); +CCTK_INT +SymBase_SymmetryInterpolateFaces (CCTK_POINTER_TO_CONST 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 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 const output_arrays[], + CCTK_INT const faces); + + +/* Statistics.c */ +void +SymBase_Statistics (CCTK_ARGUMENTS); + #endif /* ! defined SYMBASE_H */ diff --git a/src/make.code.defn b/src/make.code.defn index 84b3c24..53c838c 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -2,7 +2,7 @@ # $Header$ # Source files in this directory -SRCS = Faces.c Handles.c Startup.c Table.c +SRCS = Faces.c Handles.c Interpolation.c Startup.c Statistics.c Table.c # Subdirectories containing source files SUBDIRS = -- cgit v1.2.3