commit changes from 2 weeks of standalone programming on laptop

back into CVS:

main change is to Jacobian structure query API (code+docs),
also implement about 2/3 of this in the code


git-svn-id: http://svn.cactuscode.org/arrangements/CactusBase/LocalInterp/trunk@22 df1f8a13-aa1d-4dd4-9681-27ded5b42416

1 files changed, 101 insertions, 56 deletions

diff --git a/doc/documentation.tex b/doc/documentation.tex
index 4ae2b3d..fc18afe 100644
--- a/doc/documentation.tex
+++ b/doc/documentation.tex
@@ -142,7 +142,7 @@ We use \verb|pt|, \verb|in|, and \verb|out| as generic 0-origin integer
 subscripts into the arrays of interpolation points, input arrays, and
 output arrays respectively.  We use \verb|(i,j,k)| as a generic
 \verb|N_dims|-vector of integer subscripts into the input array
-\verb|input_arrays[in]|.  (Thus \defn{{\tt |(i,j,k)|} space} refers to
+\verb|input_arrays[in]|.  (Thus \defn{{\tt (i,j,k)} space} refers to
 the grid of data points.)  We usually only write array subscripting
 expressions for the 3-D case; the restrictions/generalizations to
 other dimensions should be obvious.
@@ -253,7 +253,8 @@ The semantics of this are as follows:  For each \verb|axis|,
 	the grid is of infinite extent), would require data
 	from outside the grid.
 \end{description}
-Other values of \verb|out_of_range_tolerance[axis]| are illegal.
+Other values of \verb|out_of_range_tolerance[axis]| are reserved for
+future use.
 
 To provide the ``fuzz'' noted above, \verb|out_of_range_tolerance[]|
 should default to having all elements set to a small positive value,
@@ -268,7 +269,11 @@ say $10^4 \epsilon$, where $\epsilon$ is the ``machine epsilon''.%%%
 	 they are about $1.19{\times}10^{-7}$ and
 	 $2.22{\times}10^{-16}$ respectively.
 	 }%%%
-{} but none of the interpolators do this at present.)
+{}  However, at present all interpolators actually set the default
+value to $10^4 \epsilon_d$, where $\epsilon_d$ is the machine epsilon
+for C \verb|double| values, even though \verb|CCTK_REAL| may actually
+be a different data type.  (This is a bug, and may get fixed some
+day\dots)
 
 If any interpolation points are out of range (as determined by the
 \verb|out_of_range_tolerance[]| critera described above), the
@@ -366,7 +371,7 @@ where
 offset = input_array_offsets[in]
 (stride_i,stride_j,stride_k) = input_array_strides[]
 \end{verbatim}
-\verb|(i,j,k)| run from \verb|input_array_min_subscripts[]|
+and where \verb|(i,j,k)| run from \verb|input_array_min_subscripts[]|
 to \verb|input_array_max_subscripts[]| inclusive
 (n.b.~this is an {\em inclusive\/} range, \ie{}
 $\verb|min| \le \verb|(i,j,k)| \le \verb|max|$).
@@ -387,6 +392,7 @@ contiguous 1-D arrays.
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 \subsection{Derivatives}
+\label{sect-generic-options/derivatives}
 
 Some interpolators can optionally take derivatives as part of the
 interpolation, \ie{} if we view the input data as being samples of
@@ -417,9 +423,9 @@ where some of the input arrays have many different operations applied
 to them.
 
 Negative \verb|operation_codes[out]| values  are reserved for future
-expansion.  An \verb|operation_codes[out]| value which is $\ge 0$ is
-taken as a decimal integer encoding a coordinate partial derivative:
-each decimal digit means to take the coordinate partial derivative along
+use.  An \verb|operation_codes[out]| value which is $\ge 0$ is taken
+as a decimal integer encoding a coordinate partial derivative: each
+decimal digit means to take the coordinate partial derivative along
 that (1-origin) axis; the order of the digits in a number is ignored.
 For example:
 \begin{flushleft}
@@ -475,10 +481,8 @@ questions:
 	 the Jacobian-query API would probably specify a cutoff
 	 to return an approximate Jacobian with reasonable sparsity.
 	 }%%%
-\item	Does the interpolator always use the same molecule size
-	and shape (possibly depending on the interpolation order,
-	molecule family, or other such parameters), independent
-	of where the interpolation points are in the grid?%%%
+\item	Does the interpolation molecule size and/or shape depend
+	on where the interpolation points are in the grid?%%%
 \footnote{%%%
 	 We can always make the interpolation molecules be
 	 the ``same size and shape'' by padding them with
@@ -487,10 +491,18 @@ questions:
 	 but this would be very inefficient if many molecules
 	 were padded in this way.
 	 }%%%
-\item	Does the interpolator use the same interpolation molecule
-	size and shape for each output array?  (The answer to this
-	question may depend on whether and/or what derivatives are
-	being computed.)
+\item	If this interpolator supports computing derivatives
+	as described in section~\ref{sect-generic-options/derivatives},
+	does the interpolation molecule size and/or shape depend
+	on \verb|operation_codes[]|?
+\item	Does the interpolation molecule size and/or shape depend
+	on the actual floating-point values being interpolated?
+	(Examples of this might include ENO (essentially nonoscillatory)
+	and/or TVD (total-variation diminishing) interpolators for
+	hydrodynamics calculations.)
+\item	Do the actual floating-point values of the Jacobian depend
+	on the actual floating-point values being interpolated?
+	Equivalently, is the interpolation nonlinear?
 \end{itemize}
 Because the different cases differ so much in their complexity,
 we define several distinct APIs for obtaining the interpolator's
@@ -498,61 +510,87 @@ Jacobian and/or domain of dependence.
 
 %%%%%%%%%%%%%%%%%%%%
 
-\subsubsection{Querying the Interpolator about the Jacobian's structure}
+\subsubsection{Determining the Jacobian's structure}
+\label{sect-generic-options/Jacobian/structure}
 
 To allow generic code to determine which of the different Jacobian-structure
 cases applies, (and thus which APIs to use), an interpolator which
-supports Jacobian operations should support using the following
-parameter-table entries to query the interpolator about the Jacobian
-structure:
+supports Jacobian operations should report its Jacobian structure
+using the following parameter-table entries:
 
 The parameter \verb|molecule_family| may may be used to query what
 molecule family is being used.  This is described in detail in
 section~\ref{sect-generic-options/molecule-family}.
 
-This parameter may be used to query whether the interpolation molecule's
-size and/or shape varies with the interpolation coordinates:
+If the interpolation molecule size and/or shape vary with the
+interpolation coordinates, the interpolator should set the parameter
 \begin{verbatim}
 CCTK_INT MSS_is_fn_of_interp_coords;
 \end{verbatim}
-The semantics of this are that if this key is present%%%
-\footnote{%%%
-	 ``MSS'' abbreviates ``{\bf m}olecule {\bf s}ize
-	 and/or {\bf s}hape''.
-	 }%%%
-{} (the value doesn't matter), then the interpolator will (re)set
-the value to 0~if the molecule size and shape do {\em not\/} vary
-with the interpolation coordinates, or 1~if the molecule size and/or
-shape {\em do\/} vary with the interpolation coordinates.
-
-This parameter may be used to query whether the interpolation molecule's
-size and/or shape varies from one output array to another:
+to 1.  Otherwise (\ie{} if the interpolation molecule size and shape
+are independent of the interpolation coordinates) it should set this
+parameter to 0.
+
+If the interpolator supports computing derivatives as described
+in section~\ref{sect-generic-options/derivatives}, and if the
+interpolation molecule's size and/or shape varies with
+\verb|operation_codes[]|, the interpolator should set the
+parameter
+\begin{verbatim}
+CCTK_INT MSS_is_fn_of_which_operation;
+\end{verbatim}
+to 1.  Otherwise (\ie{} if the interpolator doesn't support computing
+derivatives, or if the interpolator does support computing derivatives
+but the interpolation molecule size and shape are independent of the
+\verb|operation_code[]| values), it should set this parameter to 0.
+Note that this query tests whether the molecule size and/or shape
+depend on \verb|operation_codes[]| in general, independent of whether
+there are in fact any distinct values (or even any values at all) passed
+in \verb|operation_codes[]| in this particular interpolator call.  In
+other words, this is a query about the basic design of the interpolator,
+not about this particular call.
+
+If the interpolation molecule's size and/or shape varies with the
+actual floating-point values of the input arrays, the interpolator
+should set the parameter
+\begin{verbatim}
+CCTK_INT MSS_is_fn_of_input_array_values;
+\end{verbatim}
+to 1.  Otherwise (\ie{} if the interpolation molecule size and shape
+are independent of the input array values; this is a necessary, but not
+sufficient, condition for the interpolation to be linear), it should
+set this parameter to 0.
+
+If the actual floating-point values of the Jacobian~\eqref{eqn-Jacobian}
+(for a given \verb|out|, \verb|in|, and \verb|pt|) depend on the actual
+floating-point values of the input arrays (\ie{} if the interpolation
+is nonlinear), the interpolator should set the parameter
 \begin{verbatim}
-CCTK_INT MSS_is_fn_of_which_output;
+CCTK_INT Jacobian_is_fn_of_input_array_values;
 \end{verbatim}
-The semantics of this are that if this key is present (the value
-doesn't matter), then the interpolator will (re)set the value to
-0~if the molecule size and shape do {\em not\/} vary from one output
-array to another (this includes the case where there is only 1~output
-array!), or 1~if the molecule size and/or shape {\em do\/} vary from
-one output array to another.  Note that since the answer to this may
-depend on whether and/or what derivatives are being computed, setting
-this key may force the interpolator to scan completely through the
-derivative specifications.  This probably isn't {\em very\/} costly,
-but you may want to avoid unnecessarily paying that cost by deleting
-the query key (\verb|MSS_is_fn_of_which_output|) from the table
-on future interpolator calls.
+to 1.  Otherwise (\ie{} if the interpolation is linear) it should set
+this parameter to 0.
 
 %%%%%%%%%%%%%%%%%%%%
 
 \subsubsection{Fixed-Size Hypercube-Shaped Molecules}
 \label{sect-generic-options/Jacobian/fixed-sized-hypercube}
 
-The simplest case (and the only one for which we have defined an API
-at present) is when the molecules are hypercube-shaped and of (small)
-fixed size, independent of the interpolation coordinates (though likely
-depending on the interpolation order or other such parameters), and
-not varying from one output array to another.
+The simplest case (and the only one for which we have defined an
+API at present) is when the molecules are hypercube-shaped and of
+(typically small) fixed size, independent of the interpolation
+coordinates and the actual floating-point values in the input arrays
+(though presumably depending on the interpolation order and on
+\verb|operation_code|).  In other words, this case applies if
+(and only if) the Jacobian structure information described in
+section~\ref{sect-generic-options/Jacobian/structure}
+returns
+\begin{verbatim}
+MSS_is_fn_of_interp_coords = 0
+MSS_is_fn_of_which_operation = 0
+MSS_is_fn_of_input_array_values = 0
+Jacobian_is_fn_of_input_array_values = 0
+\end{verbatim}
 
 The following parameters may be used to query the molecule size:
 \begin{verbatim}
@@ -578,22 +616,29 @@ The following parameters may be used to query the
 Jacobian~\eqref{eqn-Jacobian} itself:
 \begin{verbatim}
 CCTK_REAL *const Jacobian_pointer[N_output_arrays];
-const CCTK_INT   Jacobian_offset [N_output_arrays];
+const CCTK_INT   Jacobian_offset [N_output_arrays]; /* optional */
 
-/* the next 2 table entries are shared by all Jacobians */
+/* the next 3 table entries are shared by all Jacobians */
 const CCTK_INT Jacobian_interp_point_stride;
 const CCTK_INT Jacobian_m_strides[N_dims];
+const CCTK_INT Jacobian_part_stride;                /* optional */
 \end{verbatim}
-Then the interpolator would store the Jacobian~\eqref{eqn-Jacobian} in
+All the optional entries default to 0 if omitted.
+Then for each \verb|out| where \verb|Jacobian_pointer[out] != NULL|,
+the interpolator would store the Jacobian~\eqref{eqn-Jacobian} in
 \begin{verbatim}
-Jacobian_pointer[out][offset + pt*Jacobian_interp_point_stride
-                             + mi*stride_i + mj*stride_j + mk*stride_k]
+Jacobian_pointer[out][offset
+                      + pt*Jacobian_interp_point_stride
+                      + mi*stride_i + mj*stride_j + mk*stride_k
+                      + part*Jacobian_part_stride]
 \end{verbatim}
 where
 \begin{verbatim}
 offset = Jacobian_offset[out]
 (stride_i,stride_j,stride_k) = Jacobian_m_strides[]
 \end{verbatim}
+and where \verb|part| is 0 for real values and the real parts of complex
+values, and 1 for the imaginary parts of complex values.
 
 By appropriately setting the various stride parameters, this allows
 a fairly wide variety of possible storage layouts for the Jacobian.