\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{graphicx} \AtBeginDocument{{\noindent\small Ninth MSU-UAB Conference on Differential Equations and Computational Simulations. \emph{Electronic Journal of Differential Equations}, Conference 20 (2013), pp. 27--37.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu} \thanks{\copyright 2013 Texas State University - San Marcos.} \vspace{9mm}} \begin{document} \setcounter{page}{27} \title[\hfilneg EJDE-2013/Conf/20/ \hfil Fast evaluation] {Fast evaluation of complex equations of state} \author[E. M. Collins, E. A. Luke \hfil EJDE-2013/conf/20 \hfilneg] {Eric M. Collins, Edward A. Luke} % in alphabetical order \address{Eric M. Collins \newline Center for Advanced Vehicular Systems, Mississippi State University, Box 9627, Mississippi State, MS 39762, USA} \email{emc@cavs.msstate.edu} \address{Edward A. Luke \newline Computer Science and Engineering, Mississippi State University, Box 9637, Mississippi State, MS 39762, USA} \email{luke@cse.msstate.edu} \thanks{Published October 31, 2013.} \subjclass[2000]{65D15} \keywords{Tabular approximation; adaptive sampling; Bezier surface; \hfill\break\indent spline; interpolation} \begin{abstract} One of the most common operations encountered in computational fluid dynamics (CFD) solvers is the evaluation of the caloric and thermal equations of state (EoS) which are required to compute thermodynamic state variables from the conserved values that are typically being advanced by the simulation. The complexity of these calculations can vary widely depending on the nature of the fluid under consideration. We present a method for generating an interpolated representation of the EoS that is fairly inexpensive to evaluate regardless of the complexity of the actual underlying state equations. This approach has the advantage of being agnostic towards the original representation; whether it be a complex analytic expression, expensive iterative method, or interpolated from empirical data. \end{abstract} \maketitle \numberwithin{equation}{section} \allowdisplaybreaks \section{Introduction} Computational efficiency has always been one of the most important issues facing developers of high-performance computing applications. However, as high-performance computing hardware continues to improve -- delivering ever-increasing FLOPs -- engineers continue to crave increasingly accurate simulations based on increasingly complex physical models. Obviously, this comes with additional computational cost, but what often goes unacknowledged is that different types of model complexity contribute to computational cost in significantly different ways. In this paper, we briefly look at some of the competing demands of efficiency versus complexity, particularly when that complexity is embedded deep within solver iteration hierarchy. We then propose a solution which could shift much of the expense required to evaluate certain aspects of the model to a pre-processing step. This step generates an approximation of the original model to some prescribed level of accuracy. A fairly simple algorithm is then employed evaluate the interpolated function within the context of the solver hierarchy. The present technique is largely derived from several previous efforts which have utilized tabular look-up techniques to avoid costly computations~\cite{1996:Swesty}\cite{2007:Xia.etal}. Our proposed solution extends these techniques by offering adaptive table sizes and resolutions to satisfy user specified accuracy constraints as well as support for enforcement of the continuity of the tabular approximation function and its derivatives. \subsection{Motivation} As high-performance computing hardware continues to improve, engineers continue to push the fidelity of computational simulations to accurately model increasingly complex physical phenomena. The increased complexity of the underlying physical models can have a drastic impact on the computational solvers that implement them. At a wide range of pressures and temperatures most gases behave according to the linear ideal gas law. The ideal gas law provides a simple closed form that relates pressure, temperature, and mixture to fluid enthalpy, density, and entropy. When subject to sufficiently high pressures and/or low temperatures, however, many gases will behave very differently. This change in behavior is often due to close range effects, such as Van der Waals forces. Under these conditions, the behavior of the gas can depart substantially from the linear model of the ideal gas, and closed form solutions may no longer exist that accurately describe the relationship between various thermodynamic variables. For fluid conditions in which simple closed form equations are not available, non-linear iterative solvers are often required to obtain values for enthalpy, density, or entropy at a given pressure and temperature. Luke and Cinnella~\cite{2007:Luke.etal} describe solving fluid equations using the Hirschfelder, Buehler, McGee and Sutton (HBMS) EoS. These relations involve three regions with nested non-linear functions. Another method which is often employed is the NIST REFPROP EoS that uses a modified Benedict-Webb-Rubin (MBWR) EoS~\cite{1994:Younglove.etal}. The MBWR EoS is described by 32 coefficients with an additional 20 coefficients used to describe saturation curve data and critical point properties. In the context of computational fluid dynamics solvers, in which each cell or mesh point requires a solution to these complex systems of equations, the use of non-linear solvers to obtain the required thermodynamic variables can easily become the dominant computational cost. In practice, many solvers will resort to utilizing simpler, less-accurate equations simply to make the simulation computationally tractable. In our approach, we have utilized an interpolated representation that is fairly inexpensive to evaluate, regardless of the complexity of the actual underlying state equations. This approach has the advantage of being agnostic towards the original representation, whether it be based on complex analytical expressions, expensive iterative methods, or interpolations of empirical data. The expensive evaluation operations are performed as a pre-processing step during which the underlying tabular approximation is generated. Once converted to our representation, all equations have nearly the same look-up and evaluation cost. While our initial inspiration was drawn from the adaptive Cartesian reconstructions of Xia, Li, and Merkle~\cite{2007:Xia.etal}, the underlying topological configuration of our table representation, is probably more similar to the polynomial splines over hierarchical T-meshes (PHT-Splines) of Deng et~al.~\cite{2008:Deng.etal}. However, rather than using B-spline based representations, we have chosen to utilize independent Bezier patches which have been carefully reconstructed to ensure $C^1$-continuity at the patch edges. Our algorithm allows for anisotropic refinement of the mesh depending on the features of the underlying surface data, while allowing, at most, only one hanging node per mesh element edge. In the following sections we describe a generic function approximation algorithm for real valued functions of two independent variables, and a method for efficiently evaluating the approximated surface. For the sake of generality, the independent variables will be referred to as $x$ and $y$, while the function to be approximated will be denoted as $F(x,y)$. The actual choice of independent variables will depend on the details of a particular application. Our immediate application for this approximation method is to determine an unknown thermodynamic state variable from two known variables which uniquely define the state. For example, density or energy can be computed from pressure and temperature. We may make occasional reference to this application in the following sections as it has been the primary motivation for some of our design choices. \section{Methodology} At the lowest level, the underlying representation is a piecewise $C^1$-continuous cubic Bezier surface. The cubic Bezier surface is a commonly used modeling element which is found in computer graphics (CG) and computer-aided geometry design (CAGD) applications. The mathematical description of the surface is both flexible and easy to evaluate~\cite{1997:Piegl.Tiller}. $$\label{BezierSurfDef} F(u,v) = \sum_{i=0}^3 \sum_{j=0}^3 B_i^3(u) B_j^3(v) b_{ij}$$ Here, $B_i^3(u)$ are the third-degree (fourth-order) Bernstein polynomials, $$\label{BernsteinPolynomial} B_i^3(u) = \frac{3!}{i!(3-i)!}(1-u)^{3-i} u^i$$ and $b_{ij}$ are the control points. The shape of the surface is completely determined by a linear combination of these control points. The blending coefficients can be determined by evaluating the Bernstein polynomials (Eq.~\ref{BernsteinPolynomial}) at the parametric location $(u,v) \in [0,1]\times[0,1]$, corresponding to the desired point on the surface, or by the evaluation of an equivalent algorithm such as the deCasteljau algorithm~\cite{1997:Piegl.Tiller}. Each cubic Bezier surface in the representation is referred to as a patch. Each patch has boundaries aligned with constant coordinate directions (e.g. constant $x$ and $y$). The surface is generated by performing a recursive subdivision of the desired domain until the resulting Bezier patches are capable of interpolating the original function to a specified error tolerance or a prescribed maximum recursion depth. Depending on the nature of the function to be approximated, sampling and subdivision can be accomplished in the space of the original independent variables ($x$ and $y$), or in the logarithmic space ($\log(x)$ and $\log(y)$). Once the patches are sufficiently refined to resolve the function, we establish $C^1$-continuity conditions at the interface between the patches. This process involves several steps to ensure that the resulting surface is both well defined, and well behaved. Details are provided in the following sections. \subsection{Lagrange approximation} During the first stage of approximation, all Bezier patches are generated by sampling the underlying function at 16 uniformly spaced locations ($4 \times 4$) within the respective sub-domain of the patch. The control points for the patch are then computed such that the cubic Bezier surface will pass through the sampled points. We refer to this as the Lagrange approximation. \begin{figure}[ht] \begin{center} \includegraphics[width=0.7\textwidth]{fig1} % sample_locations \end{center} \caption{Sampling sites for Lagrange interpolation} \label{LagrangeInterpolation} \end{figure} The control points may be found by solving a $16\times16$ linear system formed from the evaluation of the Bezier surface at the parametric locations of the sampled points. Let $B_{ij}(u,v)=B_i^3(u) B_j^3(v)$, and our system is formed as: $\mathbf{f} = \mathbf{B} \mathbf{b}$ where $\mathbf{f}$ is a vector of the 16 sampled values. The matrix $\mathbf{B}$ stores computed values of $B_{ij}(u,v)$. The parametric values for $(u,v)$ are fixed in each row according to the location at which the corresponding sample point $f\in\mathbf{f}$ was evaluated. And finally, $\mathbf{b}$ is a vector of 16 control points which uniquely define the surface. If we utilize a 2D indexing scheme to reflect the 2D nature of both the sample mesh and the control net, the matrix and vectors can be expanded as follows: $\begin{bmatrix} f_{00} \\ \vdots \\ f_{mn} \\ \vdots \\ f_{33} \\ \end{bmatrix} = \begin{bmatrix} B_{00}(u_0,v_0) & \cdots & B_{ij}(u_0,v_0) & \cdots & B_{33}(u_0,v_0) \\ \vdots & \vdots & \ddots & \vdots \\ B_{00}(u_m,v_n) & \cdots & B_{ij}(u_m,v_n) & \cdots & B_{33}(u_m,v_n) \\ \vdots & \vdots & \ddots & \vdots \\ B_{00}(u_3,v_3) & \cdots & B_{ij}(u_3,v_3) & \cdots & B_{33}(u_3,v_3) \\ \end{bmatrix} \begin{bmatrix} b_{00} \\ \vdots \\ b_{ij} \\ \vdots \\ b_{33} \\ \end{bmatrix}$ Here, the $f_{mn}$ values correspond to the 16 sampled values which are located at the parametric coordinates $(u_m, v_n)$: $u_m = \frac{m}{3}, \quad v_n = \frac{n}{3}$ and the mapping from parametric to state space is given by: $x = (1-u)*x_{lo} + u*x_{hi}, \; y = (1-v)*y_{lo} + v*y_{hi}$ with bounds of each patch from $(x_{lo},y_{lo})$ to $(x_{hi},y_{hi})$. To obtain the control points, we invert the system matrix $\mathbf{B}$ and apply it to both sides. $\mathbf{B}^{-1} \mathbf{f} = \mathbf{b}$ Since the Bernstein polynomials are independent of both the sampled values and the control points, the entries in the system matrix $\mathbf{B}$ are constant. Therefore, $\mathbf{B}$ and $\mathbf{B}^{-1}$ should only need to be evaluated once. \subsection{Recursive refinement} The recursive subdivision step is carried out by choosing a split-plane in the $x$ or $y$ direction. The splits are computed so as to divide the patch's subdomain in half. In the $x$-direction, the split is computed using: $\bar{x} = \frac{1}{2}(x_{hi} + x_{lo}).$ The split in the $y$ direction is computed in a similar manner. For each candidate split, two pairs of Bezier patches are generated - one on either side of the split plane. The pair which reduces the error in the approximation the most is the one selected for the refinement. This approach allows for anisotropic refinement of the domain in the direction where it is needed most first. \begin{figure}[ht] \begin{center} \includegraphics[width=0.7\textwidth]{fig2} % split_locations \end{center} \caption{Candidate patch splits with recomputed sample sites} \label{SplitLocations} \end{figure} Notice that the recomputed sample points along shared edges are evaluated from the original function definition. When these sampled values are used as input into the Lagrange approximation algorithm, the resulting surface definitions produce a common shared edge approximation. This patch splitting process continues until each patch satisfies a prescribed local error tolerance. Alternative stopping criteria can also be provided. For example, a maximum number of recursive refinements can be specified to prevent the algorithm from excessive refinement near features such as singularities or discontinuities. The recursion halts once a patch either satisfies the convergence criteria, or it exceeds the stopping criteria. The error in the patch approximations are evaluated by performing a numerical integration using tenth-order Gauss-Legendre quadrature. Technically, only a fifth or sixth order quadrature should be required to properly account for the leading error terms. However, we felt that the extra sampling resolution was justified to ensure that no important features were left unresolved during this refinement process. More importantly, care must be taken to ensure that the locations chosen to sample the error for the numerical integration do not coincide with the locations where the function was sampled to generate the approximation. It is in these sampled locations that the Bezier surface approximation will exactly interpolate to the sampled function, which would result in the integrated error norms being inappropriately biased. \subsection{Balance criteria} After the initial refinement process has completed, the resulting set of patches form a piecewise discontinuous surface as a collection of fully independent cubic Bezier surfaces (i.e. each patch has their own unique set of control points). As previously mentioned, any patches which share a common edge with exactly one other patch will have interpolated to the same four sampled values on that edge, and will therefore be $C^0$-continuous along that edge. However, any patch that shares an edge with more than one adjacent patch is not guaranteed to possess this property. It is likely that the approximated surface will be discontinuous along those edges. Before we can establish $C^1$ continuity throughout the surface, we must first attempt to minimize the effects of sampling bias by reducing the amount by which adjacent patches are allowed to differ in size. A balancing criteria has been established whereby all patches are required to be within one level of refinement of each of its neighbors. One consequence of this criteria is that all hanging nodes will only appear on the mid-points of patch edges. In other words, each patch will have at most two adjacent neighboring patches along each of its edges. Thus the sampling rate for adjacent patches will vary by a factor of two at most. \begin{figure}[ht] \begin{center} \includegraphics[width=0.4\textwidth]{fig3} % forced_split \end{center} \caption{Patch splits are required to satisfy the balance criteria.} \label{ForcedSplit} \end{figure} The set of patches which resulted from the initial recursive refinement steps are now further refined until all patches satisfy the balance criteria. Each of these subdivided patches are once again fit with cubic Bezier approximations to the original function using the Lagrange interpolation method described above. Since the balance step consists only of patches being subdivided, the accuracy of the approximation will most likely be improved. Thus, the ability of the patches to satisfy the error tolerance criteria is not effected. In addition, the most refined patches will, by definition, not have any adjacent neighbors which are more refined. Therefore, this process will not create any patches below the lowest refinement level, preserving the recursion limit if one has been specified. In other words, neither of our halting criteria is violated by this step. \subsection{Hermite approximation} In addition to the Lagrange approximation technique, we also make use of a second technique for generating cubic Bezier surface patches. We refer to this second method as the Hermite approximation. Rather than attempting to fit the patch to 16 uniformly sampled values, the function and its derivatives at the corner nodes are utilized. More specifically, this method uses the value of the function ($F$), its first derivatives with respect to each of the independent variables ($F_x$, $F_y$), and the cross-derivative ($F_{xy}$), evaluated at each corner node. \begin{figure}[ht] \begin{center} \includegraphics[width=0.7\textwidth]{fig4} % hermite_approx \end{center} \caption{Hermite approximation from the function value and its derivatives provided at the corner nodes} \label{HermiteApproximation} \end{figure} In our implementation, the values at the corner nodes are computed for each patch from the Bezier patches that were generated from the Lagrange approximations during the initial refinement steps. That is, rather than use derivatives derived from the original function (which may or may not be readily available), we utilize the derivative information encoded within the existing Bezier interpolations. \subsection{Reconciliation} Since the derivative data is computed independently for each patch, there may be a multiplicity of values at shared corner nodes. This node information is then reconciled amongst the adjacent patches in a two-step process. Step 1: For any node which is solely a corner node (i.e. is not a hanging node for any adjacent patch), we account for the following observations: (a) The function values should all already agree since the corner nodes are among the sixteen original sampled points in each patch. Therefore, they should all already be evaluating to the exact value of the original function. (b) To minimize the potential for the approximation to oscillate near rapidly changing features (e.g. discontinuities or singularities), we prefer to use the derivative with the smallest magnitude. This typically has the net effect of smoothing the approximated surface. (See figure \ref{SlopeLimiting}) \begin{figure}[htb] \begin{center} \includegraphics[width=0.7\textwidth]{fig5} % limiting_step \end{center} \caption{Slope limiting to reduce oscillations} \label{SlopeLimiting} \end{figure} Step 2: For the remaining nodes, we note that each hanging node is a hanging node for exactly one patch, and a corner node for the two neighboring patches. Since we desire the entire edge be $C^1$-continuous, we must insist that the two neighboring patches are generated in such a way that the edge shared with their less refined neighbor matches in the $C^1$ sense. To accomplish this we specify that the value of the function, and its three derivative values at the hanging node must be taken from the less refined patch. In other words, we evaluate the value and the derivatives of the less refined patch on the edge midpoint corresponding to the location of the hanging node. This information is then passed on to the two adjacent patches as corner node data for their respective interpolations. \begin{figure}[htb] \begin{center} \includegraphics[width=0.4\textwidth]{fig6} % hanging_node_dependencies \end{center} \caption{Hanging node dependency graph} \label{HangingNodeDependencies} \end{figure} Since it is likely that the shape of many of the patches will be necessarily altered by this approach, there is a specific order in which the updates to the patch shapes must be made. For this purpose, we form a data structure which records the data dependencies between all patches and hanging nodes. The data is then topologically sorted to obtain the proper ordering for the update operations. After all of the patches have been regenerated using the technique described above, the collective surface approximation should be $C^1$-continuous at all points within the specified domain. \subsection{Evaluation} Once we have established the piecewise $C^1$-continuous Bezier patch approximation to a given function, all that remains is to develop a fast method for evaluating an arbitrary point on this surface. There are two steps to the evaluation process: to search for the patch whose domain includes the point of interest, and then to evaluate the patch at that point. To these ends, we employ a $kd$-tree search algorithm, and the well-known deCasteljau evaluation algorithm. For the $kd$-tree search, we take advantage of the fact that the patches are organized hierarchically. Every time we split the domain during the refinement phases above, we also created a natural partitioning of the patches. This implies that for any node at any level of the tree, the sub-domain spanned by patches beneath that node is a contiguous rectangular region completely covered by the patches found in the child nodes beneath the current node. The $kd$-tree is generated by sorting the patches into a nested array based on the current split-direction and value. There are four possibilities regarding the distribution of child patches. If there is only one child cell, then it is considered a leaf node. If there are two or more child nodes, then there exists at least one split through the center of the sub-domain for which all child patches exist solely on one side or the other. The patches can therefore be partitioned by an x-coordinate split, a y-coordinate split, or both. If both directions are possible (i.e. no single patch spans the width or height of the remaining sub-domain), then the patches are sorted by splitting in the direction perpendicular to the orientation of the parent node in the tree. In other words, wherever possible, the orientation of the split will alternate from one level of the tree to the next. The look-up algorithm thus proceeds by simply comparing the coordinates of the desired point with the current split value and direction until a leaf node is reached. At that point, the local coordinates within the patch are computed, and then the Bezier patch is evaluated. The mapping from global to local coordinates is given by: \begin{align} u = \frac{x-x_L}{x_R-x_L} \\ v = \frac{y-y_L}{y_R-y_L} \end{align} where $(x_L,y_L)$ are the lower bounds of the patch, and $(x_R,y_R)$ are the upper bounds. The deCasteljau algorithm proceeds by recursively evaluating linear combinations of adjacent control points within the control net. In 2D, this evaluation has the form: \noindent{\tt for $k$ from 2 to 0}\\ \indent{\tt for all $(i,j) \in [0,k]\times[0,k]$} $b^k_{i,j} = (1-u)(1-v)b^{k+1}_{i,j} + u (1-v)b^{k+1}_{i+1,j} + (1-u) v b^{k+1}_{i,j+1} + u v b^{k+1}_{i+1,j+1}$ where $b^3_{i,j}$ are the sixteen original control points, and $b^0_{0,0}$ is the evaluation of the desired point on the surface. \section{Results} In high pressure environments -- such as those occurring in rocket engines, deep underground, or in certain industrial processes -- complex non-linear equations of state (EoS) are often required for the solver to make realistic predictions. Evaluation of these relations often requires an expensive iterative solver. The formulation of some EoS may also utilize expensive operations such as exponential and power functions. As a result, most of the computation time required to obtain solutions to these problems is spent evaluating the EoS. Using the method described above, we have been able to replace these expensive EoS evaluations with evaluations of our adaptive tabular surface approximations. We are able to obtain independent approximations the density and internal energy functions down to a specified error tolerance. Derived quantities, such as specific heats ($c_p$, $c_v$) and sound speed, can also be computed from these values and their derivatives obtained from the tabular fits. Since the evaluation of the cubic Bezier patch is generally much cheaper than the evaluation of complex equations of state, these surface approximations can dramatically improve the speed of computations with only modest loss in accuracy. To demonstrate the effectiveness of the proposed approach we compare the performance of our tabular equation of state against the NIST standard properties EoS known as REFPROP~\cite{1994:Younglove.etal} described earlier. Surface approximation tables are constructed for the super-critical region (with temperatures greater than the critical temperature) over the valid evaluation region for the REFPROP EoS of three different materials: $CO_2$, $O_2$, $H_2O$. The tables were refined to within about $0.1\%$ relative accuracy or better. Our method was able to achieve this goal with around $2,000$ patches for each of the species. An illustration of the tabular fit for $CO_2$ density is shown in figure~\ref{fig:CO2Tab}. The highly enriched region near the bottom right corner of this table is the region near the critical point of the fluid where density changes rapidly. \begin{figure}[htbp] \begin{center} \includegraphics[width=0.7\textwidth]{fig7} % CO2DensityTabCrop \end{center} \caption{Adapted table structure for density of super-critical $CO_2$.} \label{fig:CO2Tab} \end{figure} To compare the performance of our tabular EoS implementation to the REFPROP EoS, timing data was obtained for the evaluation of $200,000$ points sampled along the diagonals of the fit space. The results of this performance comparison are shown in table~\ref{tab:time}. For these tests, our adapted tabular EoS implementation achieves a performance speedup over the REFPROP EoS by approximately two orders of magnitude. \begin{table}[htbp] \caption{Comparison between NIST REFPROP and Tabular EoS} \label{tab:time} \begin{center} \begin{tabular}{|l||c|c|c|c|c|} \hline Model & Table & Error & REFPROP & Tabular & Speedup \\ & Size & & Time & Time & \\ \hline\hline $CO_2$ & 1968 & 0.011\% & 16.23s & 0.0786s & 206\\ $O_2$ & 2082 & 0.046\% & \phantom{1}7.10s & 0.0793s & \phantom{1}89 \\ $H_2O$ & 1949 & 0.110\% & 14.32s & 0.0783s & 183\\ \hline \end{tabular} \end{center} \end{table} It is important to note that the errors of the NIST REFPROP EoS are generally on the order of $0.1\%$ when compared to experimentally measured data. This highly efficient tabular implementation can achieve a significant performance increase with hardly any significant losses of accuracy. We also note that the isotropic tabular formulation of Xia et~al.~\cite{2007:Xia.etal} required $225,121$ patches to fit the super-critical region of $CO_2$ from the REFPROP database to an accuracy of only $0.1\%$. Our table was able to achieve a lower error bound with just $1,968$ patches. The improved efficiency of the tabular representation can be attributed to the anisotropic refinement afforded by our tabular fitting technique. Xia et~al. reported performance improvements over REFPROP for their method that were similar to our observed speedups. \subsection*{Conclusion} In this report, we have outlined an approach to generate interpolated look-up tables for functions of two independent variables. This technique may be deployed whenever the evaluation of one or more sufficiently complex functions begin to consume a significant fraction of the computational resources for a given simulation. By selecting a simple representation for the underlying surface approximation, the evaluation of complex functions and their first derivatives can be accelerated. In practice, we have observed that our technique can deliver a factor of two to three speed up in the context of multi-species chemically reacting flow simulations involving detonations and blast wave physics. As an additional benefit, this technique now provides a common interface for the flow solver to evaluate complex equations of state, regardless of their original source. We are now able to run simulations with EoS that were previously prohibitively expensive to evaluate or were in a form that would be difficult to incorporate into our solver in an effective manner. The extension of this technique to functions involving more independent variables is, at first glance, fairly straight-forward. However, the complexity of the refinement, the search algorithms for finding the desired patch, and the evaluation of the resulting $n$-dimensional hyper-surface may be rather more computationally challenging to implement efficiently. We are investigating application areas where such an approach may be warranted. \begin{thebibliography}{0} \bibitem{2008:Deng.etal} Jiansong Deng, Falai Chen, Xin Li, Changqi Hu, Weihua Tong, Zhouwang Yang, and Yuyu Feng. \newblock Polynomial splines over hierarchical {T}-meshes. \newblock {\em Graphical Models}, 70:76--86, 2008. \bibitem{2007:Luke.etal} Ed~Luke and Pasqualle Cinnella. \newblock Numerical simulations of mixtures of fluids using upwind algorithms. \newblock {\em Computers and Fluids}, 36(10):1547--1566, 2007. \bibitem{1997:Piegl.Tiller} Les Piegl and Wayne Tiller. \newblock {\em The {NURBS} book}. \newblock Monographs in visual communications. Springer, 1997. \bibitem{1996:Swesty} F.~Douglas Swesty. \newblock Thermodynamically consistent interpolation for equation of state tables. \newblock {\em Journal of Computational Physics}, 127(1):118--127, August 1996. \bibitem{2007:Xia.etal} Guoping Xia, Ding Li, and Charles~L. Merkle. \newblock Consistent properties reconstruction on adaptive {Cartesian} meshes for complex fluids computations. \newblock {\em Journal of Computational Physics}, 225:1175--1197, 2007. \bibitem{1994:Younglove.etal} Ben Younglove and Mark McLinden. \newblock An international standard equation of state for the thermodynamics properties of refrigerant 123 (2,2-dichloro-1,1,1-trifluoroethane. \newblock {\em Journal of Physical Chemistry Reference Data}, 23(5):731--779, 1994. \end{thebibliography} \end{document}