We’ve taken the world apart but we have no idea what to do with the pieces.

― Chuck Palahniuk

There are a lot of ways to turn differential equations into a discrete system and numerical solve them. The choices usually come down to three different, but slightly ambiguous choices, finite differences, finite volumes, or finite elements. One of the important, but fuzzy pieces of knowledge is the actual meaning of the variables you’re solving in the first place. Some clarity regarding the variables detailed identity can come in handy, if not be essential to high fidelity solution. With any luck we can shed some light on this.

The place to start is writing down the governing (usually differential) equations. If your problem is geometrically complex, finite element methods have a distinct appeal. The finite element method has a certain “turn the crank” approach that takes away much of the explicit decision-making in discretization, but the decision and choices are so much deeper than this if you really care about the answer.

Tiny details imperceptible to us decide everything!

― W.G. Sebald

I’ve never been a fan of not thinking deeply about what you are doing at a very specific way. Finite element aficionados seem to be very keen on avoiding as much thought as possible in discretization with a philosophy of choosing your element and associated shape functions and letting the chips fall where they may. For simple problems with a lot of regularity (smoothness) this can work well and even be advantageous, but for difficult problems (i.e., hyperbolic PDEs) this can be disastrous. In the end people would be far better served by putting more thought into the transition from the continuous to the discrete.

My recent posts have been examples of the sorts of details that really matter a lot in determining the quality of results in computing. Ideas like convergence, limiting solutions, dissipation, and accuracy all matter a tremendous degree and make the difference between stability and instability, high and low fidelity, run of the mill and state of the art, and ultimately high quality. Failure to pay acute attention to the details of the discretization will result in mediocrity.

It should come as no surprise that I don’t particularly care for the finite element method, so in keeping with my tastes I’m focusing on the finite volume and finite difference methods today. There are significant differences between the two that should be taken into account when deriving discrete approximations. Perhaps more interestingly, there is a fairly well defined way to translate between the two points of view. This translation makes for a useful addition to anyone’s discretization “toolbox”.

Once upon a time there was no such thing as the “finite volume method” it was simply a special finite difference method. Some finite differences employed a discrete conservation principle, and could be thought of as directly updating conserved quantities. These distinctions are not terribly important until methods become higher than second-order in order of accuracy. Gradually the methodology became distinct enough that it was worth making the distinction. The term “finite volume” came into the vernacular in about 1973 and stuck (an earlier paper in 1971 coined “finite area” method in 2-D).

The start of the distinctions between the two approaches involves how the equations are updated. For a finite difference method the equations should be updated using the differential equation form of the equations at a point in space. For a finite volume method the equations should be updated in a manner consistent with an integral conservation principle for a well-defined region. So the variables in a finite difference method are defined at a point, and the values in a finite volume method are defined as the integral of that quantity over a region. Transfers between adjoining regions via fluxes are conserved, what leaves one volume should enter the next one. Nothing about the finite difference method precludes conservation, but nothing dictates it either. For a finite volume method conservation is far more intrinsic to its basic formulation.

Consider a conservation law, . One might be interested in approximating with either finite differences or finite volumes. Once a decision is made, the approximation approach falls out naturally. In the finite difference approach, one takes the solution at points in space (or in the case of this PDE, the fluxes, ) and interpolates these values in some sort of reasonable manner. Then the derivative of the flux is evaluated. The update equation is then . This then can be used to update the solution by treating the time like an ODE integration. This is often called the “method of lines”.

For a finite volume method the approach is demonstrably different. The update of the variable makes explicit contact with the notion that it is a quantity that is conserved in space. For the simple PDE like that above the update equation requires the updating of the variables with fluxes computed at the edge of the cell. Notice that the fluxes have to be evaluated at the edges, which can assure that the variable, is conserved discretely. The trick to doing the finite volume method properly is the conversion of the set of cell (element) conserved values of to the point values at the edges as fluxes. It isn’t entirely simple. Again an interpolation can be utilized to achieve this end, but the interpolation form must adhere to the character of the conservation of the variable. Thus when the interpolant is integrated over the cell it must return the conserved quantity in that cell precisely. This can be accomplished through the application of the classical Legendre polynomial basis for example.

Perhaps you’re now asking about how to translate between these two points of view. It turns out to be a piece of cake. It is another useful technique to place into the proverbial toolbox. It turns out that the translation between the two can be derived from the interpolations discussed about and rightfully mirrors each other. If one takes a set of control volume, integrally averaged, values and recovers the corresponding point values the formula is remarkably simple. Here I will refer to control volume values by and the point values by . Then we can transform to point values via . The inverse operation is and is derived by integrating the point wise interpolation over a cell. For higher order approximations these calculation are a bit more delicate than these formula imply!

For finite element methods we usually have the standard flavor of the continuous finite elements. Thinking about the variables in that case they are defined by nodal values, the “shape function” describing their variation in space, and its appropriately weighted integral. A more modern and exciting approach is discontinuous Galerkin, which does not require continuity of solution across element boundaries. The lowest order version of this method is equivalent to a low-order finite volume scheme. The variable is the zeroth moment of the solution over a cell. One way of looking at high order discontinuous Galerkin methods is taking successive moments of the solution over the cells (elements). This method holds great promise because of its high fidelity and great locality.

This is just the start of this exploration, but the key is know what your variables really mean.

Little details have special talents in creating big problems!

― Mehmet Murat ildan

Here, I will provide some concepts that may spur further development through attacking each of these three shortcomings in newer methods. Part of the route to analysis is utilization of evolution equations and the analysis of numerical methods using modified equation analysis. These concepts were essential in ideas that were key in the formulation of the original generation of (linear) monotone methods [HHL76] and revisiting these concepts should provide impetus for new advances. The approach is two-fold: explain why current methods are favored and provide a catalog of characteristics for new methods to follow.

In the past, I have used modified equation analysis in the past to unveil the form of the subgrid modeling provided by high-resolution methods for large eddy simulation (the so-called MILES or Implicit LES approach). This approach has provided a systematic approach to providing explanations for the effectiveness of high-resolution methods for modeling turbulence by revealing the implied subgrid model in implicit in the methods.

In solving hyperbolic conservation laws the work of Lax [Lax73] is preeminent. Lax identified the centrality of conservation form to computing weak solutions (along with Wendroff) [LW60]. Once weak solutions are achieved, it is essential to produce physically admissible solutions, which the concept of vanishing viscosity provides. Harten, Hyman and Lax [HHL76] provided the critical connection to numerical methods by showing that upwind differencing provided physically admissible solutions utilizing the method of modified equations to the analysis of the numerical method (the same can be shown for Lax-Friedrichs-based methods). These methods form the foundation upon which high-resolution methods are based.

Of course while these theoretical developments were made parallel discoveries were happening in the first high-resolution methods. Independently three different researchers discovered the basic principle of high-resolution methods, nonlinear selection of computational stencils so that the method becomes nonlinear even for linear problems [Boris71,VanLeer72,Kolgan72]. Each method allowed the barrier theorem of Godunov [God59] to be overcome and paved the way for further developments. Ultimately, Harten provided the work to tie these streams of work together by introducing the total variation concepts to formalize the mathematics [Harten83]. These methods form a foundation for a renaissance in numerical hyperbolic conservation laws.

Subsequently, there have been several attempts to move beyond the first generation of high resolution methods, but each has met with limited success [HEOC87,Shu87]. While this work has yielded some practical methods, none has shared the extreme success of the initial set of methods and their basis in mathematics. These methods have not yielded the same quantum leap in performance as the original methods, and not met the same success in being ubiquitous in applications. The reasons for this practical failure are multiple chiefly being cost, fragility, and failure of formal high-order accuracy to yield practical accuracy. I studied some of these basic tradeoffs with Jeff Greenough [GR04] between archetypical methods PLM [CG85] and fifth-order WENO [JS95].

The trick to forgetting the big picture is to look at everything close up.

I am proposing revisiting the foundational aspects of high-resolution methods by returning to modified equation analysis of numerical methods and providing connections to providing proof of physically admissible solutions, and their control of variation in the solution (i.e., connecting to variation-diminishing properties). To conduct such analysis I am suggesting the study of the energy evolution of equations, and/or the continuous variation evolution. Each equation will admit a vanishing viscosity solution, and the modified equation analysis can provide a connection of the numerical method to admissibility. The properties of the vanishing viscosity solution define the nature of desirable results. For a simple hyperbolic PDE, , the basic physical admissibility condition can be given . For higher order methods, the viscosity must be nonlinear and a function of the solution itself. Similarly the variation evolution may be derived and its vanishing viscosity solution can be similarly stated. In the rest of the discussion here, we will focus on the linear version of the equation for simplicity, .

When modified equation analysis is applied to upwind differencing, the classical result can be easily derived for the numerical viscosity, . For the variation equation the result is similar, but informative, . We can see that the term proportional to is positive negative and will reduce the variation where ever the sign of variation changes.

If we analyze classical TVD method using the minmod limiter, the connection to physical admissibility becomes clear for the kinetic energy evolution, . Similarly the variation evolution provides a similar connection with more texture on the nature of the solution with the leading truncation error term being a flux . When expanding the flux in the modified differential equation, two terms emerge as being non-conservative of variation, . The second term is negative definite while the first term could be either positive or negative. We would categorize the scheme as being quite likely to be variation diminishing, but not absolutely. Discretely, we know that the method is TVD. Note that the viscosity coefficient, , is now a nonlinear function of the solution itself and could be called a “hyper-viscosity”. It is important that the actual discrete scheme limits the magnitude of this ratio to one, so the ratio does not blow up in practice.

We can apply the same technique to ENO or WENO methods and see the nature of its regularization at least for smooth solutions. For the famous fifth-order WENO we can derive the modified equation for the kinetic energy and find the leading order viscosity coefficient , a nonlinear and linear hyper-viscous regularization. For the variation evolution the variation changing terms are . The middle term is negative definite, but the first and third terms will be ambiguous for variation diminishment. This sort of behavior is observed for this scheme with the variation both shrinking and growing as time advances. Perhaps more importantly, the fully discrete method does not limit the size of the hyper-viscosity and it could blow up. The nonlinear hyperviscosity is itself ambiguous in terms of its variation diminishing character, which is distinctly different than the TVD scheme.

For the classical ENO scheme the result is quite similar. In the classical ENO scheme, the smoothest stencil is chosen hierarchically starting with first-order upwind. In other words the second order scheme is chosen to be the closest to the first-order scheme and the third-order stencil is chosen that is closest to the second-order stencil. We can analyze the third-order version using modified equation analysis. The effective numerical viscosity is a combination on linear and nonlinear hyperviscosity, . This can be analyzed for the variation equation and produces ambiguously signed non-conservative terms, . We find that ENO is quite similar to WENO, and the hyperviscosity magnitude is not controlled by the nonlinearity in the method. Practically speaking this is a deep concern when considering the robustness of the method.

The last method we analyze is one that I have suggested as a path forward. In a nutshell it is the high order method that is closest to a certain TVD method in some sense (I use the value of edge fluxes). Importantly, this method can still degenerate to the first-order upwind scheme while being formally third-order accurate. If the TVD method is chosen to be the first-order upwind method, we can analyze the result. The nonlinear dissipative flux is a third-order method, . For the variation equation we get a remarkably simple non-conservative term that is negative definite, . This is an extremely hopeful form indicating excellent variation properties while retaining high order.

The modified equation analysis is enabled by a combination of symbolic algebra (Mathematica) and certain transformations formerly utilized to analyze implicit large eddy simulation [MR02,RM05,GMR07].

Finally, we can utilize these methods to produce a new class of methods that offer the promise of greater resolution and more provable connections to physical admissibility. The basic concept is to use TVD methods (and their close relatives) to ground high-order discretizations. The higher order approximations will be chosen to be closest to the TVD method is a well-defined sense. Through these means we can assure that the high-order method simultaneously provides accuracy and variation control by analyzing the resulting schemes via modified equation analysis. These methods will typically have similar regularizations to the WENO methods and provide high resolution results that are both formally accurate when appropriate and higher resolution than the TVD methods they are grounded by. For a third-order method the analysis of this method provides the following results for the kinetic energy equation, and the variation evolution. I believe this approach will provide a beneficial path for method development and analysis.

Out of clutter, find simplicity.

If you’re not confused, you’re not paying attention.

[Harten83] Harten, Ami. “High resolution schemes for hyperbolic conservation laws.”Journal of computational physics 49, no. 3 (1983): 357-393.

[HEOC87] Harten, Ami, Bjorn Engquist, Stanley Osher, and Sukumar R. Chakravarthy. “Uniformly high order accurate essentially non-oscillatory schemes, III.” Journal of computational physics 71, no. 2 (1987): 231-303.

[HHL76] Harten, Amiram, James M. Hyman, Peter D. Lax, and Barbara Keyfitz. “On finite‐difference approximations and entropy conditions for shocks.”Communications on pure and applied mathematics 29, no. 3 (1976): 297-322.
[LW60] Lax, Peter, and Burton Wendroff. “Systems of conservation laws.”Communications on Pure and Applied mathematics 13, no. 2 (1960): 217-237.

[Lax73] Lax, Peter D. Hyperbolic systems of conservation laws and the mathematical theory of shock waves. Vol. 11. SIAM, 1973.

[Boris71] Boris, Jay P., and David L. Book. “Flux-corrected transport. I. SHASTA, A fluid transport algorithm that works.” Journal of computational physics 11, no. 1 (1973): 38-69.

(Boris, Jay P. A Fluid Transport Algorithm that Works. No. NRL-MR-2357. NAVAL RESEARCH LAB WASHINGTON DC, 1971.)

[VanLeer73] van Leer, Bram. “Towards the ultimate conservative difference scheme I. The quest of monotonicity.” In Proceedings of the Third International Conference on Numerical Methods in Fluid Mechanics, pp. 163-168. Springer Berlin Heidelberg, 1973.

[Kolgan72] Kolgan, V. P. “Application of the minimum-derivative principle in the construction of finite-difference schemes for numerical analysis of discontinuous solutions in gas dynamics.” Uchenye Zapiski TsaGI [Sci. Notes Central Inst. Aerodyn] 3, no. 6 (1972): 68-77.

(van Leer, Bram. “A historical oversight: Vladimir P. Kolgan and his high-resolution scheme.” Journal of Computational Physics 230, no. 7 (2011): 2378-2383.)

[God59] Godunov, Sergei Konstantinovich. “A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics.”Matematicheskii Sbornik 89, no. 3 (1959): 271-306.

[Shu87] Shu, Chi-Wang. “TVB uniformly high-order schemes for conservation laws.”Mathematics of Computation 49, no. 179 (1987): 105-121.

[GR04] Greenough, J. A., and W. J. Rider. “A quantitative comparison of numerical methods for the compressible Euler equations: fifth-order WENO and piecewise-linear Godunov.” Journal of Computational Physics 196.1 (2004): 259-281.

[CG85] Colella, Phillip, and Harland M. Glaz. “Efficient solution algorithms for the Riemann problem for real gases.” Journal of Computational Physics 59.2 (1985): 264-289. & Colella, Phillip. “A direct Eulerian MUSCL scheme for gas dynamics.” SIAM Journal on Scientific and Statistical Computing 6.1 (1985): 104-117.

[JS95] Jiang, Guang-Shan, and Chi-Wang Shu. “Efficient Implementation of Weighted ENO Schemes.” Journal of Computational Physics 1.126 (1996): 202-228.

[GLR07] Grinstein, Fernando F., Len G. Margolin, and William J. Rider, eds. Implicit large eddy simulation: computing turbulent fluid dynamics. Cambridge university press, 2007.

[RM05] Margolin, L. G., and W. J. Rider. “The design and construction of implicit LES models.” International journal for numerical methods in fluids 47, no. 10‐11 (2005): 1173-1179.

[MR02] Margolin, Len G., and William J. Rider. “A rationale for implicit turbulence modelling.” International Journal for Numerical Methods in Fluids 39, no. 9 (2002): 821-841.

Linearity Breeds Contempt

— Peter Lax

The development of high-resolution numerical methods for hyperbolic conservation laws has been an outstanding achievement for computational physics. These methods have provided an essential balance of accuracy (fidelity or resolution) with physical admissibility with computational tractability. While these methods were a tremendous achievement, their progress has stalled in several respects. After a flurry of development, the pace has slowed and adoption of new methods into “production” codes has come to a halt. There are good reasons for this worth exploring if we wish to see if progress can be restarted. This post builds upon the two posts from last week, which describes tools that may be used to develop methods.

Things have stalled for the very reasons that the adoption was so rapid and complete. The sort of radical success that the high-resolution (monotonicity-preserving) methods experienced may be difficult to impossible to replicate with new methods. The biggest thing that the monotonicity-preserving methods did was change the fundamental nature of the numerical viscosity from a form that inherently computed laminar, syrupy-looking flows to a numerical viscosity that allowed flows to take on an energetic character that gives fluid dynamics is power and beauty.

Why did the resolution simply stop at the rise from first- to second-order methods? One reason is probably deeply connected to the fundamental nature of nonlinear fluid systems and their inherent dissipation when viscosity is small. As it turns out the canonical small viscosity structures in flows, shocks and turbulence exhibit the same general form for the scaling of dissipation in these cases. Both show a dissipation that is proportional to the change in velocity in the direction of the flow cubed (the normal or longitudinal velocity). The fundamental numerical error for a second-order method in conservation (or control volume) form produces the correct asymptotic nonlinear dissipation of flow energy into heat. The monotonicity-preserving numerical methods allowed this term to achieve prominence in solutions helping to produce more physical appealing solutions.

The issue is that newer methods need to cure smaller ills in the current solutions. Sometimes monotonicity is not the right thing to enforce, and extreme values in the solution are damped severely, and shapes are greatly distorted. These problems are important, but not nearly as profound as the set of issues that monotonicity-preserving methods cured. We cannot expect the same level of transformative solution improvement. The case for adopting any newer methods is bound to be more subtle and less clear cut. Existing approaches that improve on monotonicity preservation are generally quite expensive without providing a commensurate improvement in solution quality.

Moreover the improvements in the solution from the new methods are generally highly subjective and lack a quantitative edge. At best, the case for using them is highly subjective and far from clear-cut. Adding to the problems are two key aspects in numerical analysis that further cloud and undermine the case for advancing. First, solutions are inherently lacking sufficient smoothness to give high-order accurate methods the full benefit of their properties. Secondly, the stability of the newer methods is substantially less than the adopted generation of methods, which provide almost as much robustness as the first-order methods they easily displaced. The newer methods are far too fragile for production codes. Part of the reason for this is the tendency to eschew every using the first-order approximation even locally. Thirdly the nonlinear stability principle associated with the newer methods is far weaker than monotonicity preservation and sometimes even precipitates new modes of instability in the solution.

Three strikes and you’re out

— A rule in Baseball, the “American Pastime”

Now we can revisit the three main issues for the new methods and provide suggestions that may allow progress to begin anew. The three shortcomings I identified above are the following:

1. Lack of smoothness, and full utility for high-order accuracy
2. Removal of the first-order method as a robustness mechanism
3. A much weaker nonlinear stability principle, and new potential instabilities

It might be useful to revisit some old advise from Bell, Colella and Trangenstein [BCT89] before unveiling my suggestions on principles to use in developing high resolution methods using monotoncity preservation,

1. Use a good foundational high-order method, upstream centered, or high-order centered
2. Use an entropy satisfying Riemann solver
3. Add additional dissipation at strongly nonlinear or degenerate discontinutities.

These maxims are closely related to the reasons why this class of methods has been so successful:

1. Solutions are monotonicity preserving without inducing linear laminar viscosity and replacing it with a hyper viscosity that provides more energetic physically meaningful solutions,
2. It unveils the prominence of the control volume term, which provides an asymptotically appropriate stabilizing nonlinear dissipation term,
3. When the solution is significantly under-resolved or in-danger of unphysical solutions use a first-order methods as a “safety net”.

Putting all of these considerations into account and then thinking about how to move past the current state-of-the-art. We need to account for the weaknesses of current methods and advance beyond monotonicity preservation into account, and other advice, I can write down my three suggestions (keeping BCT’s advise as a base):

1. Base the nonlinear stability principle on monotonicity-preservation with detection of extrema and careful relaxation of those conditions
2. Continue to use the high-order base method unless it produces a distinct monotonicity violating solution.
3. Careful analysis of strongly nonlinear discontinuities to assure that the approximation is locally nonlinearly stable. If the solution is highly under-resolved give up high-order accuracy and degenerate to first order.

Next week I’ll unveil some new analysis that sheds light on how the methods work based on modified equation analysis. It also points to some details of what might work in the future.

[BCT89] Bell, John B., Phillip Colella, and John A. Trangenstein. “Higher order Godunov methods for general systems of hyperbolic conservation laws.” Journal of Computational Physics 82.2 (1989): 362-397.

[Hersh] Hersh, Reuben. Peter Lax, Mathematician: An Illustrated Memoir. Vol. 88. American Mathematical Soc., 2014.

[Lax1978] Lax, Peter D. “Accuracy and resolution in the computation of solutions of linear and nonlinear equations.” Selected Papers Volume I (2005): 184-194.

Boundary violations are deeply experienced.

― David W. Earle

The post yesterday is closely connected to this one.

Yesterday I talked about a general purpose form for “limiters” applied to polynomial reconstructions, . For methods of this type the reconstruction is important because it defines the value of the function where there is not data, at cell-edges or faces where fluxes are computed, . For a method in conservation form the edge-face values are essential to updating the solution. As I will discuss the value of these values and mindful modifications of the values are severely under-utilized in implementing methods. Various methods can be improved in terms of accuracy, resolution, dissipation and stability through the techniques discussed in this post.

All of these techniques will utilize an important observation for methods of this type; one can swap edge-face values without impacting the accuracy of the method. Consider the reconstructed edge value that can (usually) be two values with a reconstruction being applied to $p\left(\theta_j \right)$ and $ p\left(\theta_{j+1} \right)$. If both reconstructions have equal accuracy we could trade these values, or replace one value by the other without impacting accuracy. Other properties like stability, dissipation or phase error are effected by this change, but formal accuracy is preserved.

One of the key aspects of the use of these edge values is a model for what an upwind method looks like applied to the data. Consider a nonlinear flux where upwinding will be applied to the edge values from two neighboring cells to determine the physically admissible flux, . An upwind approximation will use edge values from cells and . The simplest and most useful approximation of the flux is where is an eigen-decomposition of the flux Jacobian. The numerical dissipation is proportional to the eigenvalues and the jump at the cell edge . Generally speaking the physical direction of dissipation is defined by the jump from the cell-center (average) values, $latex– u_j $. The extrapolated edge values can (often) change the sign of the jump meaning that the approximation is locally anti-diffusive. A picture of faces is worth showing so one can understand what each swap looks like and envision the impact.

Perhaps the first application of this idea was the contact steepener defined by Colella and Woodward in their PPM method. Material interfaces also known as contact discontinuities are infamously diffused by Eulerian finite volume schemes even ones as high-resolution as PPM. The method introduced with PPM removes most of the diffusion usually induced by the method. The method works by locating valid interfaces and in the cells with the interface swapping the interface values at the edges. Often a material interface will be dissipative and swapping the edge values will be anti-diffusive. It is essential to apply a nonlinear stability mechanism to the swapped interface values and the polynomial reconstruction so that the anti-dissipative effects are not destabilizing. A second word of caution is the ill-advised nature of applying this technique at a nonlinear discontinuity such as a shock or rarefaction. For those flow features this technique would invite disaster.

A second, far safer application of edge swapping would be increasing the dissipation. In cases where the high-order local dissipation is anti-dissipative because the extrapolated values at the edges have a differently signed jump across the interface, swapping the edge values would make this dissipative. This condition has recently been dubbed as “entropy stable”. The condition can be computed by checking whether the jump in edge values has the same sign as the jump in cell-centered values . This should be particularly important in stably and physically moving nonlinear structures on the grid.  A concern is that too much dissipation is created as the figure shows. In this case the jump is entropy stable, but larger than the first-order method. This could actually be unstable especially if a Lax-Friedrichs flux is used because it is more dissipative than upwinding.

The penultimate approach would be to use these techniques to remove dissipation entirely from the edge. This is enabled by simply setting the two edge values to be equal in value. Any convex average of the edge values would achieve this as well as simply over-writing one of the edge values by the other. Like the anti-dissipative methods this would need to be applied with care for the stability of the overall integration procedure. In several cases the resulting approximation is actually higher order accurate. For example by applying this to first-order upwinding, the approximation is promoted to second-order central differencing. This happens for odd-ordered upstream centered methods, i.e., the third-order upwind methods is promoted to fourth-order centered edge values, the fifth-order upwind to sixth-order centered, and so forth.

The final approach to discuss would simply take one of the edge values and completely replace it by the other. The most obvious way to achieve this would be to declare that the information at the edge is unambigiously one-directional, and use that direction to choose the edge values. In this case one would short-circuit the Riemann solver and would only be stable for super-sonic flow. Of course, the opposite would produce a demonstrably anti-diffusive effect and might be useful (and dangerous too).

Don’t live a normal life by default, push the boundaries of your potential.

― Steven Redhead

Colella, Phillip, and Paul R. Woodward. “The piecewise parabolic method (PPM) for gas-dynamical simulations.” Journal of computational physics 54.1 (1984): 174-201.

Yang, Huanan. “An artificial compression method for ENO schemes: the slope modification method.” Journal of Computational Physics 89.1 (1990): 125-160.

Balsara, Dinshaw S., and Chi-Wang Shu. “Monotonicity preserving weighted essentially non-oscillatory schemes with increasingly high order of accuracy.”Journal of Computational Physics 160.2 (2000): 405-452.

Fjordholm, Ulrik S., Siddhartha Mishra, and Eitan Tadmor. “Arbitrarily high-order accurate entropy stable essentially nonoscillatory schemes for systems of conservation laws.” SIAM Journal on Numerical Analysis 50.2 (2012): 544-573.

Banks, Jeffrey William, and Jeffrey Alan Furst Hittinger. “A new class of nonlinear finite-volume methods for Vlasov simulation.” Plasma Science, IEEE Transactions on 38.9 (2010): 2198-2207.

Dumbser, Michael, and Olindo Zanotti. “Very high order P N P M schemes on unstructured meshes for the resistive relativistic MHD equations.” Journal of Computational Physics 228.18 (2009): 6991-7006.

Munz, Claus-Dieter, et al. “Enhanced Accuracy for Finite-Volume and Discontinuous Galerkin Schemes via Non-intrusive Corrections.” Recent Developments in the Numerics of Nonlinear Hyperbolic Conservation Laws. Springer Berlin Heidelberg, 2013. 267-282.

You can only exceed your limits if you’ve discovered them.

― Roel van Sleeuwen

About 15 years ago I wrote a simple paper with Len Margolin about options for modifying limiters used in high resolution schemes. Generally the paper didn’t get the success it deserved like so much else. It was studied by folks at Brown University, and was repurposed as a general positivity limiter. I will grant that the researchers at Brown did some really good work with the limiter especially proving how it could preserve accuracy (the positivity preserving limiter also preserves the design accuracy of the original reconstruction). I returned to this form in one of my recent papers, “Reconsidering Remap Methods” published earlier this year and the topic of a number of blog posts. Last week I realized that even more can be done with it and the form is really convenient for expressing a wide variety of nonlinear stability mechanisms in a more uniform manner.

The moment that crystalized my thinking was a re-reading and analysis of a limiter devised my Phil Colella and Michael Sekora for the PPM algorithm. This limiter was a “competitor” to an entirely different “limiting” strategy that I worked with Jeff Greenough and Jim Kamm. I now realized that the two limiters could be written in very nearly identical ways. Moreover, both of these ideas can be viewed as rather direct decedents of work by HT Huynh and Suresh from the late 1990’s. Each method attempts to extend nonlinear stability ideas beyond simply preserving monotonicity, to preserving valid extrema in solutions. Monontonicity preservation has the problem of seriously damping extrema with the limiter. All of these limiters can be expressed as different choices for a common form allowing direct comparison or mix and match ideas to be utilized effectively.

So, what is this common form?

Given a polynomial (or more general functions perhaps) representation of a solution in a computational cell, where , the limiter is implemented as , is the general limiter and . The limiter has a general form, . The “Standard” is the test that defines the nonlinear control of the solution while “Reconstruction Choice” is the feature being controlled. I will give some examples below. It is important to note that is the reconstruction that produces an upwind scheme that is the canonical linear monotonicity preserving (positive) method, but only first-order accurate. Higher order polynomials can be defined for higher accuracy, but without the limiter they are invariably oscillatory.

The most common choice is the definition of a positivity preserving limiter. In this case where is the lowest value of that will be tolerated for physical reasons. Next our reconstruction is interrogated .

The same can be done for monotonicity-preserving limiters. Here the most common application is “slope-limiting” where a common choice is here is the famous minimum modulus function that returns the smallest magnitude argument if the arguments have the same sign. where is the “high-order” slope that you desire to use and control the stability of.

For extrema preservation I will provide a couple of choices from my own work and Colella & Sekora. In both cases the extrema preservation will only be applied if a local extrema is detected by the scheme. Generally speaking a monotonicity preserving limiter will help do this by setting .

The idea that I worked on had a couple of key concepts: don’t give up on the original high-order scheme yet, at an extrema a ENO or WENO scheme might be useful, and the monotonicity preserving method has good nonlinear stability properties. My limiter worked on edge values that then helped define the reconstruction polynomial. It turns out that the general limiter can work here too, , and everything else follows. The parts of the limiter are with being the monotonicity preserving edge value, is the high-order edge values and is the WENO edge value. Next, the limiter is completed with . Part of the key idea is to use two complimentary nonlinear stability properties from monotonicity preservation and WENO to assure the stability of the method should the high-order method be chosen. In this case the high-order method is bounded by two nonlinearly stable approximation choices. I thought it actually worked pretty well.

Colella and Sekora had a different approach based on looking at the local curvature (second derivative) of the parabolic reconstruction used in PPM. If the second derivatives of the reconstruction and the neighborhood of the reconstruction all have the same sign and roughly the same magnitude, the extrema is viewed as valid and worth propagating undistorted by the sort of damping monotonicity preserving methods produce. A simple way to implement this approach can be implemented using the general limiter I discuss. First we need to define the second derivative of the polynomial reconstruction . In this case the specifics of the limiter are where is the classical second derivative on the mesh and . The choice of method to define the denominator in the limiter is . The resulting limiter can either be applied to the edge values or equivalently to the reconstructing polynomial.

In my next post I will talk about how edge values, can be used productively to improve the robustness, resolution and stability of these methods.

Set Goals, not Limits.

― Manoj Vaz

Rider, William J., and Len G. Margolin. “Simple modifications of monotonicity-preserving limiter.” Journal of Computational Physics 174.1 (2001): 473-488.

Qiu, Jianxian, and Chi-Wang Shu. “A Comparison of Troubled-Cell Indicators for Runge–Kutta Discontinuous Galerkin Methods Using Weighted Essentially Nonoscillatory Limiters.” SIAM Journal on Scientific Computing 27.3 (2005): 995-1013.

Zhang, Xiangxiong, and Chi-Wang Shu. “Maximum-principle-satisfying and positivity-preserving high-order schemes for conservation laws: survey and new developments.” Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences. The Royal Society, 2011.

Zhang, Xiangxiong, and Chi-Wang Shu. “On maximum-principle-satisfying high order schemes for scalar conservation laws.” Journal of Computational Physics229.9 (2010): 3091-3120.

Colella, Phillip, and Michael D. Sekora. “A limiter for PPM that preserves accuracy at smooth extrema.” Journal of Computational Physics 227.15 (2008): 7069-7076.

Suresh, A., and H. T. Huynh. “Accurate monotonicity-preserving schemes with Runge–Kutta time stepping.” Journal of Computational Physics 136.1 (1997): 83-99.

Huynh, Hung T. “Accurate upwind methods for the Euler equations.” SIAM Journal on Numerical Analysis 32.5 (1995): 1565-1619.

Rider, William J., Jeffrey A. Greenough, and James R. Kamm. “Accurate monotonicity-and extrema-preserving methods through adaptive nonlinear hybridizations.” Journal of Computational Physics 225.2 (2007): 1827-1848.

The beginning of wisdom is the definition of terms.

― Socrates

Accuracy, Fidelity, Resolution and Stability, Robustness, Reliability

These six words matter a lot in the discussion of numerical methods yet their meanings are quite poorly understood or even regularly conflict. Even worse they are poorly articulated in the literature despite their prominent importance in the character of numerical methods and their instantiation in code. Maybe I can make some progress toward rectifying this issue here. Consider this post to be a follow-on and compliment to the earlier post on improving CFD codes, https://williamjrider.wordpress.com/2015/07/10/cfd-codes-should-improve-but-wont-why/.

One of the biggest issues in providing clarity in these matters is the ultimate inability of classical mathematics to provide clear guidance as we move from the abstract and ideal world to the real and messy reality. In using methods to solve problems with a lot of reality things that work well or even optimally in an idealized version of the World fall apart. In the ideal World very high-order approximations are wonderful, powerful and efficient. Reality tears this apart and all of these characteristics change to problematic, fragile and expensive. The failure to bridge methods over to the real world undermines progress and strands innovation behind.

Accuracy is a simple thing to define, the deviation of an approximation from reality. For an idealized mathematical representation this measurement is “easy”. The problem is measuring this for any real world circumstance is quite difficult to impossible. One way to define accuracy is related to the asymptotic rate of convergence as the degree of approximation is changed this is “order-of-accuracy”. The connection of order of accuracy, and numerical precision is rather ill defined and fuzzy. Moreover reality makes this even harder.

A clear way that the World becomes real is shock waves (and other similar nonlinear structures). A shock wave or similarly discontinuous solution renders most numerical methods first-order accurate at best. The current knowledge about how formal ideal high-order accuracy connects to measurable accuracy in the presence of discontinuities is sketchy. This manifests itself in having research in methods focus on high-order accuracy without any understanding of how this would translate to accuracy for real World problems. Today the efficacy of high-order methods is merely an article of faith.

We have alternative definitions of accuracy focus on the terms fidelity and resolution. Both of these terms are even more fuzzy that accuracy. These both get applied to methods that provide their value to (more) real World circumstances where formal accuracy is diminished. Thus important classes of methods are defined as “high-fidelity” or “high-resolution”. Both of these definitions are used to imply capability to provide good solutions when the ugliness of reality intrudes into our idealized reality.

Peter Lax provided a definition of resolution in an unfortunately obscure source. There and in a fascinating interview with Lax by Phil Colella (http://history.siam.org/oralhistories/lax.htm) the concept was discussed with an astounding proposition, perhaps higher than second-order was too accurate and produced solutions that did not have enough “room” to capture true solutions. It is not a final answer, but it does yield a direction to thinking about such things.

Here is the key passage on the topic of resolution:

“COLELLA: So you have shocks, why second-order? What were you thinking?

LAX: Well, you want greater accuracy, but even more you want greater resolution. I defined a concept of resolution. If you take a difference method and you consider a set of initial value problems of interest, which in practice could be some ball in L1-space, anything will do, and then you look at the states into which it develops after the unit-time, any given time, that’s another set. The first surprise is that this set is much smaller for nonlinear…for linear equations where time is reversible, the size of this set is roughly the same as the original set. For nonlinear equations, which are not reversible and where the wave information is actually destroyed, it’s a much smaller set. And the measure of the set that is relevant is what’s called entropy or capacity with respect to some given scale delta. So the first thing to look at is what is the capacity or entropy of this set of exact solutions. Then you take a numerical method, you start, you discretize the same set of initial data, then you look at what you get after time t goes to whatever the test time was. A method has a proper resolving power if the size of this set is comparable to the size of the exact solution; if it’s very much smaller it clearly cannot resolve. And first-order methods have resolution that is too low, and many details are just washed out. Second-order methods have better resolution. In fact, I was trying to – well, I want to bring up the question: could it be that methods that are even higher order (third, fourth) have perhaps too much resolution, more resolution than is needed? I just bring this up as a question.”

I might offer a bit of support for that concept in the case of genuinely nonlinear problems below. In a nutshell, the second-order methods with conservation form provide truncation error that matches important aspects of the true physics. Higher order methods will not capture this aspect of the physics. I’ll also note that Len Margolin and I have followed a similar, but different line of thinking looking at implicit large eddy simulation (ILES). ILES is an observation that high-resolution methods appear to provide effective turbulence modeling without the benefit of explicitly added subgrid modeling.

So let’s talk about the archetype of nonlinear, real World, messy computations, shock waves. In some ways shocks are really nice, they are inherently dissipative even in the case where the system is free of explicit molecular viscosity. Dissipation in the limit of zero viscosity is one of the most profound aspects of our mathematical description of reality. For physical systems with a quadratic nonlinearity including shocks and turbulence, this dissipation scales, with being the velocity and being a constant. At its core is the imposition of reality on idealized math to describe reality, and provide a useful, utilitarian description of mathematically singular structures. This character is present in turbulence as well. Both have basically the same scaling law and deep philosophical implications.

This form of nonlinear dissipation comes directly from the application of the conservation form to methods with second-order accuracy. For energy this term is precisely the form of the asymptotic law except for its connection to the discrete system. If the method achieves a formally higher than second-order accuracy this term disappears. For very simple second-order schemes there are truncation errors that compete with this fortuitous term, but if the linear accuracy of the method is higher order, this term is the leading and dominant truncation error. This may explain why schemes like PPM, and FCT methods produce high quality turbulence simulations without explicit modeling, but methods like minmod or WENO do not. The minmod scheme has a nonlinear truncation error that dominates the control volume term. For WENO method the higher order accuracy means the dissipation is dominated by a combination hyperviscous terms.

These deep philosophical implications are ignored by the literature for the most part, with shocks and turbulence defining a separation of focus. The connections between these topics are diffuse and unfocused. A direct connection would be a stunning breakthrough, but entrenched interests in both areas conspire against this. This remarkable similarity of the limiting dissipation in the absence of viscosity have been systematically ignored by scientists. I see it as utterly compelling or simply brutally obvious for a quadratic nonlinearity. Either way the similarity is meaningful. One of the key problems is that turbulence is almost completely grounded in the belief that it can be completely described by incompressible flow. No one seems to ever question this assumption.

Incompressibility is a physically limited approximation of reality, but not reality. It renders the equations to be intractable in some ways (see the Clay prize for proving the existence of solutions!). The unphysical nature of the equations is two-fold: sound speeds are infinite and thermodynamics are removed (especially harmful is the loss of the second law). Perhaps more problematically is the loss of the very nonlinearity known to be the source of dissipation without viscosity for shock waves, that is the steepening of arbitrary disturbances into discontinuous shock waves.

I’ve written before about stability and robustness with a focus on the commonality of their definition https://williamjrider.wordpress.com/2014/12/03/robustness-is-stability-stability-is-robustness-almost/. The default basic methodology for stability analysis was discussed too https://williamjrider.wordpress.com/2014/07/15/conducting-von-neumann-stability-analysis/. If we add the term “reliable” the situation is quite analogous to the issues with accuracy. We ultimately don’t have the right technical definitions for the useful character of practical reliability and robustness of numerical methods and their executable instantiation in code. Stability is necessary for robustness and reliable, but robustness and reliable imply even more. Typically the concept of robustness applies to practical computational methods used for real World (i.e., messy) problems.

The key issue for high order methods is the inherently non-smoothness and lack of clean structure in the real world. This messiness renders high-order methods of questionable utility. Showing that high-order methods improve real world, practical, pragmatic calculation is challenge for the research community working in this area. Generally high-order methods show a benefit, but at a cost that makes their viability in production software questionable. In addition the high-order methods tend to be more fragile than their lower order cousins. The two questions of robustness in use and efficiency are the keys to progress.

Given all of these considerations, what is a path forward to improving existing production codes with higher order methods?

I will close with a set of proposals on how we might see our way clear to improving methods in codes by balancing requirements for high-order accuracy, high-resolution, robustness and stability. The goal is to improve the solution of “real” “practical” problems, not idealized problems associated with publishing research papers.

1. For practical accuracy high-order only matters for linear modes in the problem. Therefore seek high-order only for the leading order terms in the expansion. Full nonlinear accuracy is a waste of effort. Full nonlinear accuracy only matters if the flow is fully resolved and the fields contain the level of smoothness equal to the scheme (they never do!). This would allow the quadratures usually invoked by formally high-order methods could be reduced along with the costs.
2. For nonlinear structures, you just need second-order accuracy, which gives you a truncation error that matches the asymptotic structure of dissipation analytically. Removing this term may actually harm the solution rather than improve. The reasoning follows Lax’s comments above.
3. Nonlinear stability is more important than linear stability in fact nonlinear stability will allow you to use methods that are locally linearly unstable. Extending useful nonlinear stability beyond monotonicity is one of the keys to improving codes.
4. Developing nonlinear stability principles beyond monotonicity preservation is one of the keys to progress. A test of a good principle is its ability to allow the use of linearly unstable methods without catastrophe. The principle should not create too much dissipation outside of shocked regions (this is why ENO and WENO are not good enough principles). In a key way monotonicity-preserving methods naturally extended the linear monotone and first order methods. The next step beyond monotonicity preservation has not built upon this foundation, but rather introduced an entirely different concept. A more hierarchical approach may be needed to achieve something of more practical utility.
5. The option of fully degrading the accuracy of the method to first-order accuracy must always be in play. This step is the key to robustness. Methods that do not allow this to happen will never be robust enough for production work. This is another reason why ENO and WENO don’t work for production codes.

Logically, all things are created by a combination of simpler less capable components

– Scott Adams (or the Dogbert Principle that applies to high resolution schemes in Laney, Culbert B. Computational gasdynamics. Cambridge University Press, 1998).

References

Lax, Peter D. “Accuracy and resolution in the computation of solutions of linear and nonlinear equations.” Selected Papers Volume I (2005): 184-194.

Margolin, Len G., and William J. Rider. “A rationale for implicit turbulence modelling.” International Journal for Numerical Methods in Fluids 39.9 (2002): 821-841.

Grinstein, Fernando F., Len G. Margolin, and William J. Rider, eds. Implicit large eddy simulation: computing turbulent fluid dynamics. Cambridge university press, 2007.

Rider, William J., Jeffrey A. Greenough, and James R. Kamm. “Accurate monotonicity-and extrema-preserving methods through adaptive nonlinear hybridizations.” Journal of Computational Physics 225.2 (2007): 1827-1848.

Broadwell, James E. “Shocks and energy dissipation in inviscid fluids: a question posed by Lord Rayleigh.” Journal of Fluid Mechanics 347 (1997): 375-380.

Bethe, Hans Albrecht. Report on” The Theory of Shock Waves for an Arbitrary Equation of State”… 1942.

It is dangerous to be right in matters on which the established authorities are wrong.

― Voltaire

Last week I received a question via email that prompted this post. It proposed that the title of the post is true. It talked about the benefits of pushing the envelope with high performance computing. The gist of the thought is that by pushing the envelope with computing we can effectively use the mainstream high end computing resources better. This is true without a doubt. It is a benefit of having a bleeding edge research program in any area. A better question is whether it is the most beneficial way to allocate our current efforts.

Everything in excess is opposed by nature.

― Hippocrates

The ability of such a program to impact the World positively is a much more difficult and nuanced question. For high performance computing we have seen decades long focus on the power and speed of the hardware that has fueled a growth in peak computing speed consistent with Moore’s law. Unfortunately a host of essential capabilities for realizing this computing power as a scientific capability have not been similarly supported. Without these other capabilities such as physical models, solution methods, algorithms, the computing hardware is nothing more than a very expensive way to use electricity. The very things that make computers really useful for the purposes of modeling and simulation are the things we have not invested in for these same decades.

The distance between insanity and genius is measured only by success

― Ian Fleming

The issue with it is that this benefit does not exist in a vacuum. There are limits to the financial and human resources that may be devoted to the objective of “predictive” modeling and simulation. My politically incorrect assertion is that the focus on high performance computing hardware is a suboptimal approach to achieving the end result of maximizing the capability for modeling and simulation. The devotion to progress in hardware is sapping the resources that might be applied to attacking this problem in a more balanced manner.

If we have no heretics we must invent them, for heresy is essential to health and growth.

― Yevgeny Zamyatin

This imbalance is primarily exemplified by the failure to invest in people, experiments and models. When I speak of investing in people, it goes far beyond simply paying people. Investing in people means creating systems where people can develop and grow in their capability while feeling safe and secure to take huge risks. Talented people who take risks are necessary for progress, and without such risk-taking progress stagnates. Without taking risks we cannot develop talent, the two are intertwined.

An expert is someone who knows some of the worst mistakes that can be made in his (her) subject, and how to avoid them.

– Werner Heisenberg

We have destroyed the vitality of our experimental sciences, which further amplifies the destruction of our scientific staff. Experimental science is absolutely necessary to advance science. This has the knock-on effect of undermining the creation of new, better models for science. Having the World’s fastest computer cannot replace any of these shortcomings.

It doesn’t matter how beautiful your theory is, it doesn’t matter how smart you are. If it doesn’t agree with experiment, it’s wrong.

― Richard Feynman

A big part of the problem is that aspects of modeling and simulation closer to the modeling are starved. The modeling associated with high performance computing is static and relatively free of progress. We are not progressing toward introducing new models into codes and ultimately practical use in a healthy way. The same can be said for solution methods that provide either more powerful or effective ways of solving models. These methods then provide the impetus for algorithms that systematically provide means of solution. Our current computing emphasis is geared toward efficiently delivering existing algorithmic solutions for existing methods for existing models.

History warns us … that it is the customary fate of new truths to begin as heresies and to end as superstitions.

― Thomas Henry Huxley

Inadequacies in models are simply allowed to persist while the fiction of their rescue by more powerful computing platforms a false promise. The real scientific answer to this issue is that a more powerful computer, better software, better algorithms or better methods cannot save a model that is incorrect. Despite this maxim we keep attacking modeling and simulation dominantly through computing hardware.

…a more powerful computer, better software, better algorithms or better methods cannot save a model that is incorrect.

The fastest computer has even taken a page from Cold War nationalism with the “missile gap” being replaced by the “supercomputer gap”. Chinese prowess in computing is demonstrated by their supercomputer and funding is supposed to follow the American effort to regain the crown of fastest. I, for one, am far more worried about the Chinese and Russian investments in human resources in modeling, methods and algorithms than computers. By all appearances these investments are significant. The American public should be far more worried by the encroachment of mediocrity into the research staff at our National Labs than our lack of the fastest computers.

In the republic of mediocrity, genius is dangerous.

― Robert G. Ingersoll

The management of our National research programs and devotion to intellectually questionable priorities is by far a greater threat to our National security than anything our adversaries are doing. A perfect example of the problem is the increasingly legacy nature of our computer codes. We looked at the local history of one type of code and noted that up until 1990 we had a new code every five years. Since 1990 we have simply kept the same old codes. This is 25 years with the same code with the same methods and same approach. We have basically lost an entire generation of staff. We have lost a generation of progress and research. The models and methods are frozen in time. This is a recipe for mediocrity at best, disaster at worst.

Progress is born of doubt and inquiry. The Church never doubts, never inquires. To doubt is heresy, to inquire is to admit that you do not know—the Church does neither.

― Robert G. Ingersoll

So, yes, having leading edge computing is a great, wonderful and important thing for the country. It’s true for any country desiring international leadership. Getting to properly defining what leading edge computing is actually comprised of becomes difficult. A completely naive and incorrect way to define leading edge is having the fastest or most powerful computer on Earth. Everyone knows what a fast computer is, but a powerful computer is a more subtle question. I would argue that in many respects my iPhone is more powerful and useful than virtually any supercomputer I’ve used. The problem in defining powerful is the limited utility of supercomputers. Supercomputers are important for solving scientific problems, which are necessarily limited in context.

The riskiest thing we can do is just maintain the status quo.

― Bob Iger
Moreover, this effort for leading edge computing lies in a resource constrained trade space and the focus on hardware leaves other efforts starved for funding or focus.  Even this discussion leaves most of the important nuance untouched, the dependence on people and their talent. The issues around the efficacy of the HPC efforts are subtle and far more nuanced than the mere power of the computer. A powerful supercomputer is useless without talented people to use it. The problem in the United States is that people are something we have chronically and systematically under-invested in. Our universities are in decline and part of a vastly corrupt system that underserves the public at a massive cost. The consequences of this decline in education are then amplified by the destruction of the social contracts associated with post-educational work.

Unhappiness lies in that gap between our talents and our expectations.

― Sebastian Horsley

Employees are viewed as commodities and infinitely replicable, even at National Labs. The lifetime employment necessary for deep sustainable expertise has been replaced by an attitude more appropriate for Wal-Mart. Over the past couple of decades the sort of strong scientific leadership once provided by the National Labs has been replaced by Lab employees who are little more than “sheeple” who bend to political will rather than speak up and offer their expertise instead of politically correct pablum. Today the Labs simply do what they are told. Their spirit has been beat out of them. I might even be so bold and to say that the attempt to lead in scientific computing says more about our lack of scientific leadership than our commitment to it.

The smart way to keep people passive and obedient is to strictly limit the spectrum of acceptable opinion, but allow very lively debate within that spectrum….

― Noam Chomsky
If the models, methods, and codes used on our fastest computer are lacking, the computer’s value is diminished to the point of being worthless. A computer will provide results that are as good as the codes running on it. If the people running the problems on the computer are similarly lacking, the computer’s value is diminished as well. The computer is only answering questions as good as the people asking the questions. I believe that we have systematically failed to make investments in the models; methods, codes and people sufficient to make the focus on computing power pay off. We have created a new generation of legacy codes (just because its written in C++ does not keep it from being a legacy code!), with legacy models and methods and a staff that cannot fully understand the codes or calculations they are running. The fiction that all we need to do is refine the mesh and run calculations on bigger computers to predict nature continues to hold sway.

Every truth in this world stretched beyond its limits will become a false doctrine.

― K.P. Yohannan
This is a situation that the mismanagement of the labs has created.  The DOE has done the same thing to its Labs that DOD did. DOD foolishly destroyed its research labs 30 years ago, and over the last 20 years, DOE has followed a similar path toward mediocrity. The national resource of these Labs is being allowed to fade and wither.  We have allowed the Labs to atrophy. Our approach to high performance computing is but one example of this. The situation is even worse when you look at what we have done to our experimental sciences. Under these conditions having the lead in computing hardware will do little actually support our national security because we have failed to keep competence in the fundamentals necessary for efficacy.

Change almost never fails because it’s too early. It almost always fails because it’s too late.

― Seth Godin
In other words for the hardware to really matter in the delivery of predictive simulation and modeling, everything upstream of the machine needs to be right.  This includes the people.  We have failed to invest in leading edge technology in the very things that make the supercomputer valuable. We have a skewed and unbalanced view of how computing works, which allows the justification of our current programmatic path, but fails to deliver true progress.

Heretics are the new leaders. The ones who challenge the status quo, who get out in front of their tribes, who create movements.

― Seth Godin

Happiness is the pleasantest of emotions; because of this, it is the most dangerous. Having once felt happiness, one will do anything to maintain it, and losing it, one will grieve.

― Kij Johnson

Rarely has any movie stuck with me like Inside-Out. It is hard to imagine how an animated movie about an 11-year-old girl has left such a lasting impression on me, but it has. It is advertised as a kid’s movie with hidden themes for their parents to enjoy, but that is wrong. It is a movie for adults that will entertain your kids. The movie is well done like you’d expect from Pixar and it packs a massive punch. It’s a really good movie, but it’s a better adult movie than kid’s movie.

What has given this movie such staying power?

It gave me a lot to consider that can apply to my own life. It provided an open playful theme to deeply think about extremely deep themes in how you personally manage your response to reality through your emotions. The main characters in Inside-Out are the emotions of the 11-year-old girl, Riley. The dominant character is Joy, with other key emotions Sadness, Anger, Fear and Disgust all vying to respond to her World. When the story begins, Joy dominates Riley’s life, but a major life event throws everything into turmoil. The main conflict is how much should you allow sadness into your life and what happens when you don’t.

Riley has spent 11 years growing up in Minnesota living a life of fun, family, friends, and especially hockey. Her family picks up, and moves to San Francisco and all hell breaks loose inside Riley. Everything she has known and loved in her life has been taken away. Sadness begins to creep into her reactions, and Joy tries to reign in this reaction. Push comes to shove and disaster ensues. Riley’s emotions fall apart along with her. Riley begins to make bad decisions based on inappropriate emotional responses, and her life spirals out of control. The foundations of her entire psyche begin to fall apart. Meanwhile Joy and Sadness try to get back into her emotionalwheelhouse (along with all the hijinks that the child movie viewers will enjoy).

Emotions are not irrational, but how we rationalize the world. The key is to use the right ones for the information we are being handed. We need to use our cognitive tools represented in emotions to process our own world and respond appropriately. Irrational responses are applying the wrong emotion to the situation, which can result in bad outcomes. The world throws a lot at us and the reactions need to be proper. We have a myriad of events that range from wonderful to horrible, and our emotions need to match the event. This is rational. Therefore sadness, fear, anger or disgust despite being negative are important and proper for many things. They are necessary for our survival.

In the end, Joy and Sadness return to save the day with a key realization. Sadness is appropriate and important as a response to traumatic events. Joy is not. Trying to make everything joyous and happy is simply wrong and inappropriate. Sadness adds a texture to formerly unambiguously happy memories that now have complexity the younger Riley lacked. Riley is growing up and rebuilding her emotional world. Joy is no longer so dominant; the other emotions now have much more sway in her emotional makeup. Letting Sadness in was the key to her response to crisis.

The lesson will always repeat itself, unless you see yourself as the problem–not others.

― Shannon L. Alder

A key lesson for us today is that unhappiness and sadness are the right response to many things in life. We have a culture that drives home the message that we should always be happy and satisfied, and feel joy. Many situations call for sadness, or anger, or fear, or disgust. If we don’t feel the right emotion our reaction to the situation is inappropriate and harmful to our long-term well-being. The movie has the deep and impactful message for all of us that “negative” emotions are more than simply okay, they are powerful and proper responses to many of life’s events.

Even more powerfully there are circumstances where happiness or joy is utterly wrong and harmful. The attempt to imprint joy onto these situations hurts us and provides an improper personal response to life’s travails. Life is about balance and the push and pull of events. For us to learn, grow and develop correctly, we must process and respond to events in a way that fits the events. When we do not respond properly we can hurt our future. Moreover the full repertoire of emotions is needed to live our lives with the tools to deal with what life throws at us. The so-called negative emotions should be embraced when they fit the circumstances. In the case of Riley sadness was proper and healthy, and the attempt to feel nothing but joy left her in a tailspin of almost disastrous magnitude.

Don’t cry because it’s over, smile because it happened.

― Dr. Seuss

If failure is not an option, then neither is success.

― Seth Godin

Alternate titles

Progress in CFD has stalled. Why?

Why are methods in CFD codes so static?

Why is the status quo in CFD so persistent?

Status quos are made to be broken.

― Ray Davis

There are a lot of reasons for lack of progress in CFD codes, and here I will examine one particular issue. The reality is that there is a myriad of issues plaguing modern codes. I’ve written about issues with our modeling and its lack of suitability for tackling modern simulation questions. One of the major issues is the declaration that success won’t be reached until computers are far more powerful. This is also testimony to the lack of faith in innovation and creativity in research (risk aversion and fear of failure being key). As a result funding and focus for improving the fundamentals of CFD codes has dried up. It like the community has collectively thrown up it hands and said, “its not worth it!”

The riskiest thing we can do is just maintain the status quo.

― Bob Iger

We have an overly focused research program toward utilizing the next generation of computing hardware. The major overarching issue is a general lack of risk taking in our research programs spanning from government funded pure research, through applied research programs and extending to industrially focused research. Without a tolerance for failure and hence a risk, the ability to make progress is utterly undermined. This more than anything explains why the codes are generally vehicles of status quo practice rather than dynamos of innovation.

Yesterday’s adaptations are today’s routines.

― Ronald A. Heifetz

If one travels back into the mid-1980’s there was a massive revolution in numericalmethods in CFD codes. Methods that were introduced at that time remain at the core of CFD codes today. The reason was the development of new methods that were so unambiguously better than the previous alternatives that the change was a fait accompli. Codes produced results with the new methods that were impossible to achieve with previous methods. At that time a broad and important class of physical problems in fluid dynamics were suddenly open to successful simulation. Simulation results were more realistic and physically appealing and the artificial and unphysical results of the past were no longer a limitation.  

These methods were high-resolution methods such as flux corrected transport (FCT), high-order Godunov, total variation diminishing (TVD), and other formulations for solving hyperbolic conservation laws. These terms are in other words the convective or inertial terms in the governing equations transporting quantities through waves most typically through the bulk motion of the fluid. These new (at that time) methods produced results that when compared with preceding options were simply superior by virtually any conceivable standard. In addition, the new methods were not either overly complex or expensive to use. The principles associated with their approach to solving the equations combined the best, most appealing aspects of previous methods in a novel fashion. They became the standard method almost overnight.

Novelty does not require intelligence, but ignorance, which is why the young excel in this branch.

― Anthony Marais

This was accomplished because the methods were nonlinear even for linear equations meaning that the domain of dependence for the approximation is a function of the solution itself. Earlier methods were linear meaning that the approximation was the same without regard for the solution. Before the high-resolution methods you had two choices either a low-order method that would wash out the solution, or a high-order solution that would have unphysical solutions. Theoretically the low-order solution is superior in a sense because the solution could be guaranteed to be physical. This happened because the solution was found using a great deal of numerical or artificial viscosity. The solutions were effectively laminar (meaning viscously dominated) thus not having energetic structures that make fluid dynamics so exciting, useful and beautiful.

When your ideas shatter established thought, expect blowback.

― Tim Fargo

The new methods would use higher accuracy approximations as much as possible (orsafe to do so), and only use the lower accuracy, dissipative method when absolutely necessary. Making these choices on the fly is the core of the magic of these methods. The new methods alleviated the bulk of this viscosity, but did not entirely remove it. This is good and important because some viscosity in the solution is essential to connect the results to the real world. Real world flows all have some amount of viscous dissipation. This fact is essential for success in computing shock waves where having dissipation allows the selection of the correct solution.

The status quo is never news, only challenges to it.

― Malorie Blackman

The dissipation is the essence of important phenomena such as turbulence as well. The viscous nature of things can be seen through a technique known as the method of modified equations. This method of numerical analysis derives the equations that the numerical method effectively solves. Because of numerical error when you solve an equation numerically, the solution more closely matches a more complex equation.

In the case of simple hyperbolic conservation laws that define the inertial part of fluid dynamics, the low order accuracy methods solve an equation with classical viscous terms that match those seen in reality although generally the magnitude of viscosity is much larger than the real world. Thus these methods produce laminar (syrupy) flows as a matter of course. This makes these methods unsuitable for simulating most conditions of interest to engineering and science. It also makes these methods very safe to use and virtually guarantee a physically reasonable (if inaccurate) solution.

The new methods get rid of these large viscous terms and replace it with a smaller viscosity that depends on the structure of the solution. The results with the new methods are stunningly different and produce the sort of rich nonlinear structures found in nature (or something closely related). Suddenly codes produced solutions that matched reality far more closely. It was a night and day difference in method performance, once you tried the new methods there was no going back.

Negative results are just what I want. They’re just as valuable to me as positive results. I can never find the thing that does the job best until I find the ones that don’t.

― Thomas A. Edison

This is the crux of the issue with moving on to even more advanced methods, the quantum leap in performance to be had then simply won’t be repeated. The newer methods will not yield a change like the initial movement to high-resolution methods. The newer methods will be better and more accurate, but not Earth-shatteringly so. In today’s risk adverse world making a change for the sake of continual improvement is almost impossible to sell. The result is stagnation and lack of progress.

The problems don’t end there by a long shot. Because of the massive improvement in solutions to be had with the first generation of high resolution methods to a very large extent cost wasn’t an issue. With the next generation of methods, the improvements are far more modest and the cost of using them is an issue. So far, these methods are simply too expensive to displace the older methods.

The issues don’t even stop there. The new methods also tend to have relatively large errors compared to their cost. In addition the newer methods tend to be fragile and may not handle difficult situations robustly. The demands of maintaining formally high-order accuracy are quite expensive (time, space integration demands are costly whereas the first generation high resolution methods are simple and cheap). The result is that the newer approaches are methods that “do not pay their way.”

The balance of accuracy and cost has not been negotiated well. This whole dynamic is worth a good bit of discussion.

The key to this issue is the lack of capacity for high-order accuracy to be achieved in practical problems. To get high-order accuracy the solution needs to be smooth and differentiable. Real problems conspire against this sort of character at virtually every turn with singular structures both in the solution itself, not to mention geometry or physical properties. Real objects are rough and imperfect, which tends to breed more structure in solutions. Shock waves are the archetype of the problem that undermines high-order accuracy, but the problem hardly stops there.

The measure of intelligence is the ability to change.

― Albert Einstein

All of these factors conspire to produce in real problems, results only improve their accuracy at first-order (or worse), which means that double the mesh produces half the error. In other words, the accuracy is linearly proportional to the mesh spacing. This is a big deal as the second-order means that halving the mesh yields a four times reduction in error. Third-order would yield an eight times reduction. The reality is everything gives first-order accuracy or worse. The key for high-order working at all is that the high-order methods give a lower starting point for the error, which it sometimes does. The problem is that high-order methods are too expensive to justify the improvements they provide. The question is whether the benefits of practical accuracy can be achieved without incurring the costs typical for such methods.

Sometimes a clearly defined error is the only way to discover the truth

― Benjamin Wiker

The higher costs of the high-order methods are associated with a multitude of the characteristics of these methods. The basic steps associated with creating the high-order approximations use more data and involve many more operations than existing methods. If this wasn’t bad enough these methods often require a multiple of evaluations to integrate their approximations using quadratures. In cases of time-dependent methods, these methods often require more steps and require smaller time steps than the standard methods. To make matters even worse these methods are often not applicable to complex geometries associated with real problems. If you add on relative fragility and small gains in practical accuracy, you get the state of affairs we see today.

Restlessness is discontent — and discontent is the first necessity of progress. Show me a thoroughly satisfied man — and I will show you a failure.

― Thomas A. Edison

Meanwhile the theoretical and mathematical communities will tie themselves to high formal order of accuracy even when the methods are inefficient. The very communities that we should depend on to break this log jam are not motivated to deal with the actual problem. We are left in a lurch where no progress is being made toward improving the work horse methods in our codes.

To improve is to change; to be perfect is to change often.

― Winston S. Churchill

The cost part is almost a uniformly disappointing part of these methods most of which is dedicated to achieving formally high-order results. The irony is that the formal order of accuracy is immaterial to their practical and pragmatic utility. Almost no effort has been devoted to understanding how this cost accuracy dynamic can be negotiated. Without progress and understanding of these issues, the older methods, which now are standard will simply not move forward. Thus we had a great leap forward 25-30 years ago followed by stasis and stagnation.

Change almost never fails because it’s too early. It almost always fails because it’s too late.

― Seth Godin

Here are some “fun” research papers to read on these topics.

[Harten83] Harten, Ami. “High resolution schemes for hyperbolic conservation laws.”Journal of computational physics 49, no. 3 (1983): 357-393.

[HEOC87] Harten, Ami, Bjorn Engquist, Stanley Osher, and Sukumar R. Chakravarthy. “Uniformly high order accurate essentially non-oscillatory schemes, III.” Journal of computational physics 71, no. 2 (1987): 231-303.

[Lax73] Lax, Peter D. Hyperbolic systems of conservation laws and the mathematical theory of shock waves. Vol. 11. SIAM, 1973.

[LW60] Lax, Peter, and Burton Wendroff. “Systems of conservation laws.”Communications on Pure and Applied mathematics 13, no. 2 (1960): 217-237.

[Boris71] Boris, Jay P., and David L. Book. “Flux-corrected transport. I. SHASTA, A fluid transport algorithm that works.” Journal of computational physics 11, no. 1 (1973): 38-69.

(Boris, Jay P. A Fluid Transport Algorithm that Works. No. NRL-MR-2357. NAVAL RESEARCH LAB WASHINGTON DC, 1971.)

[VanLeer73] van Leer, Bram. “Towards the ultimate conservative difference scheme I. The quest of monotonicity.” In Proceedings of the Third International Conference on Numerical Methods in Fluid Mechanics, pp. 163-168. Springer Berlin Heidelberg, 1973.

[Shu87] Shu, Chi-Wang. “TVB uniformly high-order schemes for conservation laws.”Mathematics of Computation 49, no. 179 (1987): 105-121.

[GLR07] Grinstein, Fernando F., Len G. Margolin, and William J. Rider, eds. Implicit large eddy simulation: computing turbulent fluid dynamics. Cambridge university press, 2007.

[RM05] Margolin, L. G., and W. J. Rider. “The design and construction of implicit LES models.” International journal for numerical methods in fluids 47, no. 10‐11 (2005): 1173-1179.

[MR02] Margolin, Len G., and William J. Rider. “A rationale for implicit turbulence modelling.” International Journal for Numerical Methods in Fluids 39, no. 9 (2002): 821-841.

I suppose it is tempting, if the only tool you have is a hammer, to treat everything as if it were a nail.

― Abraham Maslow

A collection of my thoughts on issues with modeling and simulation for the future with an emphasis on modeling. Without modeling improvements, the promise of simulation cannot be achieved. Our current approach to high performance computing focused on faster computers without balance is intellectually bankrupt. Without changes in our fundamental philosophy on modeling in computational simulation, the investments in hardware will yield little benefit for society.

A common point of view these days regards the existing code base as a massive investment to be preserved and ported to the new generation of computers. What is often not articulated is the antiquated nature of the modeling available in these codes. The approach used to model material has been in use for over 50 years and has quite a lot in common with a “textbook” approach to material description found in undergraduate courses in engineering. These undergraduate courses define the basic approach to analysis of systems, which are by their very nature macroscopic. Are these descriptions appropriate for highly detailed computer models? And at what point must the description change to reflect the physics of the scales being explored?

Computers allow a far more detailed description of these systems discretizing the macroscopic system into pieces that should resolve ever more of the microstructure and its response. The problem is that description of the materials is still almost isomorphic to the philosophy expressed in the undergraduate textbook. The material’s description is invariant all the way from the macroscopic view to scales that are clearly uncovering microscopic details. This highlights several clear scientific problems that extending the current code base to exascale computation will only exacerbate.

A key to moving forward is the recognition that the class of problems that can be simulated has grown. The older homogeneous average response modeling is still useful and valid, but only for a restricted class of problems. New capabilities and models will enrich simulation’s impact through providing avenues to new classes of problem solving. The new class of problems is defined by creating simulations that more faithfully provide the role of experiments and device testing. The simulations should be able to selectively probe the cases where off-normal response of devices arises. This will allow the analysis to assist in the determination of the limits of operation and safety for engineered systems.

1. At a macroscopic level, systems are not deterministic, yet the models we rely upon are. The models are exercised in an overly deterministic manner.
2. The material descriptions are invariant to the scales, which they are described at.
3. The questions answered by the codes do not match the asked about these systems any more.
4. A scientific vibrant field would not tolerate the level of inflexibility implied by current modeling practice. Vibrant science would demand that the models evolve to better match reality.

Addressing this set of issues is going to require deep scientific investigations, and perhaps more deep cultural evolution. We have a wealth of approaches investigating and solving problems associated with multiscale, multiphysics which bridge detailed microstructural simulation to macroscopic scales. The problem is that these approaches are not being utilized to solve the applied problems we currently tackle with our code base. None of these methods are being forced to displace the ancient techniques we rely upon today. As a result the state of practice is stuck in quicksand and remains static.

The importance of modeling as a driver for simulation capability should be obvious as well as its role as the essence of utility for the entire enterprise. This importance is not as obvious when looking at the balance of efforts in simulation science. For example no amount of accuracy, or computer power, or software quality can rescue a model that is inadequate or wrong. Only a focus on improving the model itself can rescue it. Today improving models is far down on the list of priorities for simulation despite its primal role in the quality of the enterprise. Nearby issues of solution methods and algorithms for model are also poorly funded. Most of the emphasis is tilted toward high performance computing and are implicitly predicated on the models themselves being correct.

Even if the models were judged to be correct, the advances in experimental science should be providing pressure to improve models. Improvements in detection, information and analysis will all yield ever better experimental measurements and access to uniquely innovative experimental investigations. These should provide a constant impetus to advance models beyond their current state. This tension is essential to the conduct of high quality science. If science is healthy there is a push and pull with theory and experiment where a theoretical advance will drive experiments, or new experimental observations will drive theory to explain. Without modeling being allowed to advance in response to experimental evidence, the fundamental engine of science is broken.

Furthermore the culture of analysis in engineering and science reinforces these approaches. First and foremost is the commitment to deterministic outcomes in simulation. Experimental science makes it very clear that our everyday macroscopic world has stochastic elements. There is a deterministic aspect to events, but the non-deterministic aspects are equally essential. By and large our analysis of experiments and simulations works steadfastly to remove the stochastic. Usually this is the adoption of averaging (or regression fits to data). These average properties or events then become the working model of our systems. In the past this approach allowed great progress, but more and more our engineered systems are defined more properly by the extremes of behavior that they can exhibit.

Our entire modeling approach especially that used in simulation are completely ill suited to address these extreme behaviors. A fundamental change in modeling and simulation philosophy is necessary to advance our understanding. Our models do not produce actual physically realizable simulations because not system actually acts like an average system everywhere. Instead the average behavior results from variations in behavior throughout the system. Sometimes these variations produce effects that are exactly associated with the newer questions being asked about extreme behavior.

The new methods do not displace the need for the old methods, indeed the new methods should appropriately limit to the solutions found by the old methods. The new methods allow the resolution of scale-dependent behavior, and off-average behavior of the system, but need to be self-consistent with traditional methods for simulation. Perhaps just as importantly, those conducting simulations should be deeply aware of when the old methods lose validity both in terms of scale-dependent behavior, and the questions being addressed through the simulation.

This brings the idea of questions to the forefront of the discussion. What questions are being addressed via simulation? There is a set of questions that older simulation methods are distinctly capable of answering. These questions are not the same questions driving the need for simulation capability today. In providing new models for simulation, the proper questions are primal in importance.

The current simulation capability is tied to answering old questions, which are valid today, but less important as new topics are crowding them out. Examples of the older questions are “what is the performance of this system under average conditions?” “What is the yield of this production process?” “How large is the average margin of performance beyond the requirements for the system?” With the key aspect of the questions being answered being the capacity of the modeling to attack the average properties and performance of engineered systems. By the same token the uncertainty we can assess via modeling today via simulation is the lack of knowledge about the average behavior of these systems, which is not the same as the uncertainty in the behavior of the actual system.

This mindset influences the experimental comparisons done for the purposes of validation as well. Experimental data is often processed into an average, and then compared to the simulation. No single experiment is appropriately simulated, but rather the simulation is modeling the average of the experiments. As such, the simulations are not truly modeling reality because for many physical systems, the average response of the system is never produced in a single experiment. As discussed below, this mindset infects the interpretation of experiments in a deeply pernicious manner.

The new questions being asked of simulations are subtly different, and require different models, methods, algorithms and codes to answer. Key among these questions is “how much variation in the behavior of the system can be expected?” “How often will the system manifest certain extreme behavior?” “How will the entire population of the system behave under certain conditions?” “What is the worst behavior to expect from the system and how likely is it to happen?”

Ideally, the calculation should be the same as observations from a physical experiment (validation), not the average of all experiments. In this way our simulations do not model any reality today because they are almost invariably too homogeneous and deterministic in character. Experiments, on the other hand, are heterogeneous and variable yielding some degree of stochastic response. Systems truly have both a variable stochastic character, which usually acts as a non-deterministic component around a major homogeneous and deterministic aspect of the system. Today our models are predominantly focused on the homogeneous, deterministic aspect of these systems. This aspect is the focus of traditional models and the older questions. The new questions are clearly focused on the secondary stochastic aspects that we average away today. The result is a strong tendency to treat single experiments inappropriately as instances of average response when they are simply a single instance from a population of possible experiments. When the deterministic calculation is forced to compare too closely to the non-deterministic aspects of an experiment, problems ensue.

Of course this decomposition is only approximate. For nonlinear systems the separation between stochastic and deterministic is dependent on the circumstances, and the nature of the system itself. Some instances of the system will yield a different decomposition because of the coupling of the system’s response to variability. Examples of the newer questions to be addressed by simulation abound in areas such as device engineering, stockpile stewardship and weather/climate modeling. For example, a key aspect of an engineered device is the portion of the population of devices that can be expected to fail under the extreme conditions associated with normal use. This may have significant reliability consequences and economic side effects. Similar questions are key in stockpile stewardship in part to address shortcomings in the degree of testing in the field or as populations of devices diminish and reduce statistical method’s effectiveness. Extreme weather events such a rain, wind or snowfall have extreme consequences on the mortality and economic impacts on society. The degree to which climate change causes an increase in such occurrences has significant policy consequences. Simulations are being relied up to an ever-greater degree to estimate this issue.

In many cases the modeling in our workhorse engineering analysis codes is quite recognizable from our undergraduate engineering textbooks. Rather than forming a distinct field of study as the modeling unveils more mesoscopic and ultimately microscopic details, the modeling is still couched in terms of the macroscopic methods used in classical desk modeling and calculations. The modeling does not account for the distinct aspects of being applied to a discretized system where smaller scales are available. Many of these models are clearly associated with average, steady state behavior of the full macroscopic system. The multiscale modeling is simply short-circuited by the traditional view of modeling defined in many codes. For example continuum codes for fluids, solids, heat transfer and mechanics all use uniform, homogenized properties for solving problems. The philosophy is virtually identical to the macroscopic material description that would be familiar to undergraduate engineering students.

This is madness! This was a reasonable fifty years ago as these methods first came into use and the number of computational elements was small and the elements were large. Today these methods are quite mature, and the number elements is huge and their size is clearly separated from the large scale. The scale separation dictates that a model that more properly describes the material at the scale of simulation overturns the homogenized models. A homogenized material can only describe the homogenized outcome, or the average solution for the material. Furthermore this homogeneous model will not match any actual circumstances from reality.

One of the key aspects of real experiments is the ever-present random component of results. The initial and boundary conditions all have a random uncontrolled variability that yields the variation in results. In homogenized simulations, this aspect of reality is washed out and for this reason the simulation is unrealizable in the real World. At times the random component is significant enough that the result of the experiment will radically depart from the average response. In these cases, however small in probability, the simulations fall completely short of serving to replace experiments and testing. This aspect of simulation is woefully lacking from current plans despite in centrality to the role of the simulation in providing a transformative scientific tool.

Another place where current simulation approaches fall demonstrably short of serving the modeling of reality are ideal models. These models are often mathematically beautiful evoking Hamiltonian structure and deep provable properties that breed devotion by the mathematically inclined. All of this simply detracts from the lack of physical reality bound up in this idealization. These models lack dissipative forces, which define the presence of the second law of thermodynamics, a necessary element for continua associated with reality. By too greatly focusing on the beauty and majesty of the ideal model, the primal focus of modeling reality is ultimately sacrificed. This is simply too great a price to pay for beauty. More perniciously the approach produces models with seemingly wonderful properties and rigor that seduce the unwary into modeling the World in utterly unphysical manners. In many cases the modeling is constructed as the solution to the ideal model plus an explicit model for the non-ideal effects. It should be a focus of modeling to assess whether the intrinsically unphysical aspects of the ideal model are polluting the objective of modeling reality.

In computing there is a chain of activities that provide value to the World. Modeling is the closest thing to reality. No amount of computing speed, algorithmic efficiency, and methodological accuracy can rescue a model that is inadequate. Once a model is defined in needs to be solved on the computer via a method. The method can be streamlined and made more efficient via algorithmic advances. Finally we must consider that all of these need to have software for implementation and as well as mapping to the computing hardware. At the end of the chain the computing hardware is dependent on everything above it for its capacity to impact our reality. Again, modeling is the absolute key to any value at all in simulation.

Your assumptions are your windows on the world. Scrub them off every once in a while, or the light won’t come in.

― Isaac Asimov