Trust is a big word. 

The whole concept of trust or being trusted is a cornerstone of civilization, and presently our society has a frightfully short supply of it.  The crisis in trust has been forty years in the making and its side effects are profound and wide reaching.  I’m going to talk about one of these impacts as it pertains to the conduct of science.  To say that the public doesn’t trust science or scientists is an understatement even though polls would count scientists as one of the more trustworthy professions.  The cost to our society is immense with greater expense, less achievement, and ultimately horrible wastes of money, lives and opportunity.  All of this tends to create a viscous cycle that simply undermines trust further, and reinforces the attitudes that caused the problems in the first place.

Let us be clear, the cost of the lack of trust we have is massive.  It eats away at the very fabric of society like a corrosive substance.  Instead of applying trust to the conduct of activities, we demand checks and rechecks and accounting and progress reports all to provide the proof that business is being conducted in a trustworthy manner.  The lack of trust itself acts to corrode the environment and lead to an even greater degree of double-checks.  It is a classic death-spiral.  If the flow of good faith isn’t restored, when will this end?

As we look backwards in time, the early to mid 70’s seem to be the origin of this wholesale change in society’s attitudes.  Perhaps earlier during the 60’s the movement against the “establishment” laid the groundwork, with a general sense of widespread dishonesty and distortion of facts by the government in the conduct of the Vietnam War (a concomitant conflict as part of the Cold War).  The bellwether event of the 70’s was the Watergate scandal on the systematic, criminal activity in the White House ending with President Nixon’s resignation.  The unwinding of trust continued with the energy shock, Love Canal, Three Mile Island all etched on the National Psyche.  Science looms large in each of these areas and became part of the American experience of disappointment and mistrust.  People wanted accountability from their government and the science it funded.  These failures were viewed as a waste and embarrassment, and the system had to be changed.  So, driven by this the public and by proxy Congress demanded change.  Trust in our institutions and science as one of them had been revoked.

Where does this leave use today?  We have put numerous systems into place that enforce accountability.  The problem is that the accountability itself is not accountable and our laws demand that we literally spend $10 to save $1.  Because of all the issues with trust, I probably cost twice as much as I should to the taxpayers, and this isn’t the worst of it.  The work I do isn’t daring enough, and I don’t have nearly as much freedom as I should.  Ironically not enough responsibility is required of me either.   This undermines the value of what I do even more.  All this accountability and bean counting and progress report writing, and conduct of operations and training I do is a terrible return on investment for taxpayers.  In the end, all these things are their fault for their systematic lack of trust in the system and the cost of making sure I’m not wasting money by wasting a lot of money.

Scientists are terrified of failure and do not take risks.  Without the risk, much of the reward is gone.  In a recent address Bill Press noted the massive return on the investment in research associated with the Human Genome Project (on the order of 200:1).  This is true enough, but the government, and business can’t pick the winners that well.  We have to fund a lot of research that never pans out to get one of these massive payouts.  The reality is that we might have had to fund 30-100 times as much work to get the payout of the Human Genome Project.  This undermines the return on investment, but if you don’t fund all the research you never get the winner.  It is much like the rule of venture capital; you have to fund a lot of dogs to get that one winner.  The same applies to science.

Americans seem to be convinced we can manage ourselves into a sure thing.  We can’t, but no one, especially a Congressman likes to be told this simple truth.  Sometimes the best things you can do is put your trust in others and simply support good work done by good people.  Sometimes it works out, sometimes it doesn’t.  When it doesn’t, it doesn’t mean the people are bad, or the work is bad, it means they weren’t lucky enough to choose the thing that wins.  It isn’t fraud or incompetence either; its just bad luck.   

Simply put, success, learning and achievement depends completely on failures.  We must learn to accept outright failure as the glorious foundation of upon which success is built.  Ironic yes, but it’s the trust.  We must allow ourselves to fail in order to have the chance of success.

The people who win Nobel prizes are not always the best; they are just the luckiest.  Most of them are really good scientists, but lots of equally good or better scientists never get the fame, mostly because they weren’t lucky enough to be at the right place at the right time.  The problem is that won’t do that Nobel Prize-Winning work without risk and without funding a lot of dogs.  Realize that the Nobel Prize is a deeply lagging indicator, and the United States should start to see a deep falloff in our accumulation of this prize in the decades to come.  This future failure as a society has its origin in trust, or more properly the lack of it.

Lack of trust may be one way we can easily capture the malaise we find ourselves in.  To get out of this we need to go out and fail, and fail some more.  We might call it “glorious failure” because in its wake the seed of success are planted, and given the environment that allows you to fail once more success in the truest sense can be achieved.  To get to the point where we can allow this requires courage, faith in humanity and trust.

The core idea here is that adaptive mesh technology in codes seems to encourage a certain level of outright sloppiness because of the adaptivity engendering a false sense of confidence.

An interesting (and generally good) development in computation is the increasingly pervasive availability of mesh adaptivity particularly in computational fluid mechanics.  My comments are not only relevant to CFD although they might be the worst abusers.  Of course, things like verification are most often done in CFD, so it might all work out in the wash.  Hand-in-hand with the use of mesh adaptivity is a seeming laxness in attitudes toward solution quality.  In other words, adaptive meshes seem to fill people with the belief that they are inoculated from worrying about the quality of their calculations.  It is a bit ironic considering that the users of adaptive meshes are exhibiting an implicit commitment to mesh converged solutions, but often avoid doing any of the actions that might confirm this.  There seems to be an innate belief that the solution is sufficiently high quality if the mesh is refined enough.  They seem to be ignoring some of the deeper dangers associated with adaptivity.  In fact, an adaptive mesh is provides a substantially greater risk to solution quality and even greater care is necessary.

The first and most common incarnation of a misbegotten belief in intrinsic quality with adaptivity involves mesh sensitivity, or poor man’s verification.  The common statement is that we refined the mesh and the answer didn’t change, thus the mesh is fine enough.  This might be true, or it might simply indicate a stalled stagnant solution.  Without assessing the actual convergence of the solution towards a mesh independent answer these two situations cannot be distinguished easily from one another.  The adaptivity also add the danger that should the mesh refinement be chosen poorly, the mesh adaptive process may itself be causing the solution to stagnate.  I’ve actually seen this in the case where the problem is focused on a single region in the flow, and we bias the mesh refinement toward that region while ignoring the rest.   Instead of getting a better answer where we think we care about the result we actually hurt the quality by ignoring important parts of the solution that influence the result we care about.  For some people this may be counter-intuitive, but the honest truth is that the quality of the solution where we care about it is not only a function of the solution in that particular location.  The way to avoid this problem is simple enough, check the results for convergence and make certain that the solution is simply not stagnating.  It also may not be a good idea to de-refine the solution where something important is happening even if we don’t care about the details in that location.

The secondary issue with adaptivity relates to the overall attitude towards solution quality.   It has been said that adaptivity actually makes all verification unnecessary.  This is simply wrong-headed; adaptivity actually makes verification more necessary because the mesh selection procedure can be thought of as an additional approximation.  One is using the mesh selection to approximate the effect of having a fine mesh everywhere.  Adaptive code users often seem to implicitly expect that they are getting the fine grid everywhere solution at a discount by using adaptivity.  This might be approximately true under ideal circumstances, but when part of the mesh is selected to not be refined, one expects the solution to suffer.  One must also remember that the criteria for selecting the refined mesh (and where it is at lower resolution) is not entirely perfect, often basis on heuristics.  We need to check whether the degradation is significant.  The check for degradation is the verification process.  The faith placed in adaptivity is to do this where the solution is not important.  This faith can be seriously misplaced, and the entire procedure can go off the rails.  A large part of the cure for this is systematic code verification where the adaptivity is tested.

The second issue is the change in the discretization due to the variation in the mesh density.  One needs to make sure that the differencing of the equations in the presence of mesh adaption does not upset the basic properties of the original solution procedure.  Generally, one expects that the accuracy of the solution will be equal to the unadapted case AT BEST.  Usually the differencing in the face of adaption is degraded in terms of accuracy.  Again, the quality of the code and solution needs to be extensively tested.  In fact, the adaption places a greater burden for the testing of quality rather than lessening it as many assume.

Finally, I can speak of the upside of adaptivity.  Computational practitioners should consider adaptivity as providing two chief benefits: you can get refined meshes with less cost (but greater danger), and you can get solutions at multiple resolutions.  If you have done a sufficient job testing the quality of your code, then the refined mesh solutions are a massive benefit.  With solutions on multiple meshes you can test and estimate numerical errors more easily than standard approaches. 

In summary, adaptive meshes are wonderful techniques, but rather than lessen the burden on the computational code and analysts to prove quality, it increases the burden.  This is largely due to the greater complexity of adaptive mesh codes, and the potential to cut corners in the implementation.    As with most things, adaptivity is not a something for nothing technique, with great power comes, great responsibility! 

Use adaptivity, use it wisely and trust, but verify!

There are a lot of numerical methods out there and it might be hard to decide which one is right for your purpose.  A large part of the decision should be determined by the basic properties of the numerical method.  These properties can be determined by some basic analysis techniques.  I highly recommend doing these analyses personally rather than relying on what the authors of papers provide you.  A good place to start is reproducing the analysis found in textbooks, then research papers.  It provides the template to really learn how to do this.  These analyses have a number of benefits that I’ll outline below.

In the writing I will assume you’re pretty well versed in numerical methods, and you’ve had advanced mathematics classes at least an undergraduate level (and this will miss cause most of this is graduate level stuff).  On the other hand, why dumb it down? There is too much of that these days!

The classical starting point for numerical analysis is round-off error, which had importance back in the day when computers were crude.  Round-off error could be a real practical problem and even simple problems could be corrupted by the manner in which numbers were represented.  The analysis of Gaussian elimination is the archetype.  These issues have faded into the woodwork, but are staged to make a return to prominence as computers become less reliable.  This lack of reliability is related to the end of Moore’s law, and the physical limits of our present computing engineering.  As circuits become so small that the size of an atom becomes substantial, the difference between encoding a zero and one becomes fuzzy.  Cosmic rays can cause issues, and the roundoff problems will be back from the grave. 

Another topic you might cover in undergraduate classes is interpolation error.  This is another closely related topic.  In fact many numerical methods are designed via interpolation.  The interpolation problem is multifaceted

The most fundamental aspect of numerical analysis is the stability of an algorithm, and the simplest thing to analyze is the solution of ordinary differential equations.  The analysis is as simple as it gets, and is a good starting point.  Furthermore many algorithms can be split into two pieces, time and space differencing.  The ODE analysis can give results to the temporal part, and introduces the basic concepts. Like many of the basic analysis techniques, the analysis only formally applies to linear methods and models.  Below, a method that can cope with nonlineariy will be given (modified equation analysis).   The basic techniques can be found in many book, but Ascher and Petzold’s SIAM book is my favorite.

For partial differential equations Von Neumann stability analysis is the gold standard.  It was first introduced in the late 40’s at a lecture in Los Alamos, and the report on the topic was classified for decades.  Nevertheless Von Neumann shared his method with others on the outside and it rapidly spread to become the standard.  It was part of the paper published with Von Neumann that introduced artificial viscosity.  It is also called Fourier analysis because the technique involves replacing the spatial variables with a complex Fourier series.  Aside for this the conduct of Von Nuemann’s method is almost identical to ODE analysis, but produces a richer set of outcomes because the space and time aspects are characterized.  Trunctation error estimates, dissipation and dispersion (phase or propagation) aspects come out of the analysis.    Durran’s book on numerical geophysics contains a good introduction or Strikwerda’s numerical analysis book.

Pure dispersion analysis is closely related.  This most clearly can be applied to the pure spatial operator and used similarly to ODE analysis.  Unlike ODE analysis it also involves a Fourier series expansion, but only in space.  It is a simple, but limited alternative to Von Neumann’s method.

Modified equation analysis has less power than any of the methods above, but formally applies to nonlinear equations and importantly to nonlinear methods.  While it is limited in scope, it carries powerful intuitive capability.  One expands the method in terms of Taylor series, and looks at the remainder.  The remainder always involves more complex, higher order differential equations than the original equation.  For example, modified equations provided a very clear picture of implicit numerical diffusion arising from upwind differencing.  One derives the differential equation that you more closely approximate than the original equation, and in the case of upwind differncing you approximate the original advection law plus a diffusion term, and the solution looks like this.  This was important in the theoretical development because it showed how to connect numerical solutions to entropy, which determine what solutions are physical.

Energy methods are a form of analysis that I have little experience with.  I know there are some tests and conditions that exist with the finite element community, but I can’t personally vouch for the utility of them.

Symbolic arithmetic makes it easy to do the analysis.  The course I took was to look at a very simple method (e.g., upwind differencing), and analyze it by hand.  I then took on a more complex method, but rapidly discovered that the analysis is terribly tedious and error prone.  For me, Mathematica came to the rescue!  I could do simple methods easily, and extend to more complicated methods, and still more complicated methods.  You can start to analyze systems of equations or solutions in multiple dimensions, or multi-level methods.  At some point you can push current symbolic capabilities to the breaking point, but by then you’re examining quite exotic methods.  The point is that you can do a great deal with relative ease.  It also really makes you appreciate the skill and patience of those that successfully analyzed complex schemes by hand!

Here are the reasons in no particular order:

1. You can do it easily, so there really isn’t an excuse for not.  If you have the basic analysis at hand, it can be reused to understand a particular case you’re interested in.
2. You learn by actually doing, not simply reading, not by simply coding and implementing (although that part of it is recommended too).  Holistically speaking one really knows their topic when you’ve read about, thought about, analyzed personally and implemented personally.  At this point you’ve successfully learned.  Learning is the heart of excellence.  Learning is the heart of expertise.
3. You’ll understand your calculations behave in ways you never antisipated.
4. You’ll start to see how to improve things and you’ll be equiped to test, understand and prove that your improvements are better.
5. Your view of numerical methods will become nuanced, and you’ll be a better user of them.
6. The difficulty of solving differential equations will become clear, and your appreciation for progress in the field will be more complete.
7. You will get a feel for numerical stability, dispersion, dissipation and accuracy.  This understanding will allow you to diagnose the nature of your simulations like a boss.
8. Your simulations will start to tell you new things that you never appreciated or even saw before.
9. You will learn something new almost every time you analyze a new method, or analyze an old method in a new way.

10. It will push you to learn new math, and expand your knowledge to areas you might not have thought you needed.

11. Your curiousity will be piqued.

12. You won’t rely on the published literature for knowledge.  It’s the give a man a fish and you feed him, teach him to fish and you feed him for life…

The greatest thing about the human mind is the ability to detect patterns.  This pattern detection and connection capability is in essence the holy grail of artificial intelligence.  I’ve noticed such a pattern across a broad spectrum of technical endeavors associated with the 1-norm.  The 1-norms latest foray into my consciousness comes from “compressed sensing,” but it resonates with a bunch of other things.  These patterns are always trying to tell us, “there is something here worth paying attention to”.

This is all the more amazing since the pattern detection capability that humans have is probably more related to our staying alive as nomadic hunters, 10’s of thousands of years ago.  In those days the pattern would say, “there’s something to eat” or “something is over there that might kill me” or “here is a good place to sleep” and so on.  Now the same wiring is providing increased acuity applied toward abstract constructs like functional spaces. Go figure.

I’m trying to write something up about compressed sensing and related fields, but don’t feel up to the task quite yet.  So instead of nursing that as a bonus this week I’ll just give my brain dump of what is intriguing me, and what tendrils are connected to other areas of science that are perhaps far from compressed sensing at first blush.

1. L1 is the sharpest and least convex norm.
2. “When a traveler reaches a fork in the road, the ℓ1 -norm tells him to take either one way or the other, but the ℓ2 -norm instructs him to head off into the bushes. ” – John F. Claerbout and Francis Muir, 1973. I If you think about the rest of the post, it might be clear what this means.
3. Error and convergence for discontinuous (hyperbolic) problems should be measured in the L1 norm.  The L1 norm is very closely related to compactness arguments and the total variation norm.  Total variation concepts have been essential in numerical approximations to hyperbolic conservation equations, and image processing.  At the close of this post it might be worth keeping this firmly in mind.
4. The least squares problem is ubiquitous in mathematics, physics and engineering largely because it is so easy to compute.  L2 is related to energy, which adds to the appeal.  People apply L2 almost without thought as it is built into numerous computing packages and it is an introductory topic.  This ubiquity is unfortunate because of all the bad things about L2 minimization framework.  Among these bad things are the numerous unstated assumptions that come along for the ride and aren’t even known (or acknowledged) by those that use it. In many applications L1 is more appropriate and useful, but computing solutions to the L1 minimization is inherently nonlinear and there aren’t any non-trivial algorithms.  Trivial introductory algorithms exist for L2.  Some minimization is better than none, but too much begins and ends with least squares.
5. Until recently I pondered what sort of math takes place in fractional “norms”.  I’ve received a partial answer with the compressed sensing work.  Little “L” norms aren’t norms because of lack of convexity, but can be computed.  L0 tells one how sparse a set is, that is it measures all the non-zero entries.  Because of the lack of convexity it isn’t very computable (NP-hard!).  The cool thing is that with “high probability” minimizing in L1 can give the answer to the same question.  This is the key to the utility of compressed sensing.
6. The L1 norm is related to robust statistics.  In other words, if you construct estimates using L1 minimization you are not susceptible to outliers in the data.  If instead you do the standard least squares estimator (L2), you are very susceptible to outliers in the data, that is if you have a corrupt data point and do a L2 fit, your result can be influenced greatly, L1 will cope by ignoring the outlier.  It turns out that most of the outlier detection schemes are predicated on assumptions from least squares, and assumptions associated with Gaussian statistics.
7. I found that a L1 fit to local data was a nice way of differencing hyperbolic PDEs.  It gave a much-improved solution compared to L2.  When in doubt my suggestion is solve the problem using L1 minimization first.  You’ll probably have to do it outside of MS Excel!
8. In analyzing the convergence of numerical approximations to PDE (i.e., solution verification), the L1 norm gives by far the best solutions.  If you were to use on norm to fit your data I believe L1 is the best.
9. Moving onto the ideas in the frame of compressed sensing, I came across the LASSO, which is a way of regularizing least squares minimization (which you might want or need to do if your system is under-determined or terribly conditioned).  The classical regularization is Tikhonov’s where a L2 penalty term is added to the matrix you solve (basically adding a term to the diagonal).  LASSO is instead a L1 based penalty.  Like most L1 things it works better once you solve the system.  There is more, and this is where things start to get really really cool.
10. LASSO is really a way to select models.  Usually the regularization as discussed in the previous point uses a small penalty, because you want the problem you solve to be very close to the original problem.  What if instead you make the penalty large?  With the L1 regularization you start to find that the minimized solution starts to set most of the solution to zero.  Ultimately as the penalty becomes large enough everything is zero except a single value.  This usually happens in a well-defined sequence and the ordering compared to penalty size tells you what variables are most important to the data.  You can then choose how complex you want to your model to be and which variables to keep and which to throw away.  I’m playing around with this concept right now.
11. Compressed sensing arose from signal processing work using wavelets.  A couple of techniques arose “matching pursuit” basically a greedy algorithm to find out what basis was the best representation of an image, and “basis pursuit” which is an all-at-once algorithm that finds the sparsest representation.  A certain version of basis pursuit that applies to noisy data, “basis pursuit denoising (BPDN)” is basically identical to LASSO.  The big difference is that LASSO is generally applied to over-determined problems and BPDN is applied to under-determined problems.  Both methods apply the magic of the L1 norm to do the selections.  These pursuits are just another model selection problem.
12. Finally compressed sensing can be discussed albeit briefly.    This concept arose from basis pursuit and has a phenomenal mathematical basis from the assurances associated with using L1 to get the L0 solution of the sparsest basis.  Much is made of the RIP, the restricted isometry property, which characterizes systems that compressed sensing can tackle although it has been noted that the theory is a bit pessimistic compared to what can be achieved in practice.  It looks like an awesome bit of math that has massive implications on signal processing, graphics, and data analysis.

I came across the area from a blog post on the most influential papers in statistics in the last decade.  David Donoho’s paper introducing the concept has about 8000 citations in less than 10 years!  I met David last year at a workshop, and had no idea about why he was so well known, he just gave a really nice talk on reproducible software in an educational setting.

If you really want an introduction to compressed sensing read Terry Tao’s wonderful, wonderful post (http://terrytao.wordpress.com/2007/04/13/compressed-sensing-and-single-pixel-cameras/) on the single pixel camera.  His math papers on related topics are pretty great too!

I hope this was interesting and perhaps inspiring.  There is cool new math out there, and better yet it is powerful and useful.

There are two standard ways of thinking about differential equations, the strong and the weak form.  The strong form is usually the most familiar involving derivatives that all exist.  By the same token it is less useful because derivatives don’t always exist, or a least not in a form that is familiar (i.e., singularities form or exist and evolve).  The weak form is less familiar, but more general and more generally useful.  The weak form involves integrals of the strong form differential equations, and describes solutions that are more general because of this.  This is because we can integrate across the areas that undermine the strong solutions and access the underlying well-behaved solution.  Functions that have derivatives that are undefined have perfectly well defined integral.  One of the keys to the utility of the weak form is the existence of many more solutions for the weak form than the associated strong form.   With these many solutions we have additional burden to find solutions that are meaningful.  By meaningful I mean physical, or more plainly those that can be found in the natural universe.*

Just as the differential equations can be cast in the strong and weak form, initial and boundary conditions can do the same.  Just as with the differential equations these differences are consequential in how the solutions to differential equations perform.  I have found that both ideas are useful in the context of successfully running problems.  If you run problems then you have to deal with boundary conditions and boundary conditions are essential, and essentially ignored in most discussions.  Perhaps not conincidently the boundary conditions are the ugliest and kludgiest part of most codes.  If you want to look at bad code, look at how boundary conditions are implemented.

The classic approach for setting boundary conditions in a finite volume code is the use of “ghost” cells.  These ghost cells are added onto the grid in layers outside the domain where the solution is obtained.  The values in the ghost cells is set in order to achieve the boundary effect when the stencil is applied to the entire method (i.e., the finite difference method).  For a finite volume method, the values in the cells are updated through the application of fluxes in and out of each of the cells.  This step is where issues occur namely that the fluxes are not necessarily the same fluxes one would get by applying the boundary condition correctly.  Ghost cells use the strong form of the PDE in their mindset.  This is the imposition of the continuously differentiable equation outside the domain.  One proviso is that a control volume method is inherently a weak form of the PDE concept, so there is a sort of mixed boundary condition once you decide how to compute the flux at the boundary.

Fluxing is where the weak boundary condition comes in.  In the weak boundary condition, the flux itself or the process of the flux calculation is used to impose the boundary condition.   If one is using a solver based upon a Riemann solver then the state going into the Riemann solution on the boundary from outside the domain imposes the boundary condition.  The best thing about this approach is that the update to the values in the cells is updated in a manner that is consistent with the boundary condition. 

If you are using a first-order method for integrating the equations of motion, the weak boundary condition is the only one that matters.  For most modern shock capturing methods, the first-order method is essential for producing quality results.  In this case the weak boundary condition determines the wave interaction at the boundary because it defines the state of the fluid at the boundary on the outside of the domain.

Just to be more concrete regarding the discussion I will provide a bit of detail on one particular boundary condition, a reflection (or inviscid wall).  At a wall the normal velocity is zero, but the velocity extrapolated to the boundary from cell-centeres is rarely identically zero.  The other variables all have a zero gradient.  The normal velocity outside the boundary in the ghost cell is the mirror of the cells internal, thus taking on a value of negative the value on the physical grid cells.  By the same token, the other variables take identical values to those inside the domain.  More particularly for the strong boundary conditions, the first ghost cell takes the same value as the first interior cell except for normal velocity where the value is negative the interior value.  If you have a second ghost cell then the same is applied to values taken from the second interior cell, and so on.  For the weak boundary condition, the same procedure is taken, but applied to the values at the boundary extrapolated from the last interior cell.

In practice, I use both strong and weak boundary conditions to minimize problems in each step of the finite volume method.  In this manner the boundary conditions are reinforced by each step and the solution is more consistent with the desired conditions.  This is clearest when a scheme with a wide stencil is used where non-boundary cells access ghost cells.  In practice the discrepancy between the strongly and weakly imposed boundary conditions is small; however, when limiters and other nonlinear numerical mechanisms are used in a solver these differences can become substantial.

 If you want to read more I would suggest the chapter on boundary conditions from Culbert Laney’s excellent “Computational Gasdynamics” that touches upon some of these issues.

 * Mathematicians can be interested in more general unphysical solutions and do lots of beautiful mathematics.  The utility of this work is may be questionable, or at least scrutized more than it often is. The beauty of it is a matter of aesthetics, much like art.  An example would be the existence of solutions to the incompressible Euler equations where the lack of connection to the physical world is two fold**: incompressibility isn’t entirely physical (infinite sound speeds! No thermodynamics) and setting viscosity equal to zero isn’t physical either (physical solutions come from viscosity being positive definite, which is a different limiting process).  These distinctions seem to be lost in some math papers.

 **Remarkably the lack of connection to the physical World doesn’t limit things to being useless.  For many applications the unphysical approximation of incompressibility is useful.  Many engineering applications profitably use incompressible flow because it gets rid of sound waves that are not particularly important or useful to get right.  The same is true for inviscid flows such as potential flow, which can be used for aircraft design at a basic level.  The place where this lack of physical connection should be more worrisome is in physics.  Take turbulent flow, which is thought of as a great-unsolved physics problem.  The fact that turbulence is believed to be associated with the mathematics of incompressible flows, yet remains largely unsolved should come as no surprise.  I might posit that the ability to make progress on physical problems with unphysical approximation might be dubious.

For the compressible Euler equations one requires some sort of dissipation mechanism to get stable numerical solutions.  Without dissipation you are limited to solving useless and uninteresting problems.  For almost any choice of initial data, shocks will automatically form, and dissipation must arise. This is the “regularized” in a singularity of this blog’s title.  The original dissipation is the Richtmyer-Von Neumann artificial viscosity.

A few years later (1951-1952) Lax recognized the sensibility of using the conservation form of the equations to compute shocks.  This was the “Lax-Friedrich’s” method written as a finite difference scheme, but in fact was a finite volume method (with a sneaky Riemann solver hidden underneath, but no one really realized this for a long time).  Lax’s method was first-order accurate and very dissipative (in fact one can show that it is the most dissipative stable method!).  It is notable that Lax did his original work in Los Alamos, which is the same place where the Richtmyer-Von Neumann method arose.

Around the same time in Russia, Sergei Godunov was working on his PhD thesis.  He was asked to develop a method to demonstrate a real computer that his institute would soon install by solving the Euler equations.  Finding methods provided to him unworkable and unstable, he developed something new.  This became Godunov’s method, which is distinguished by using a Riemann solver to evolve the solution in time.  The Riemann solution was developed in Germany in 1859 as the solution of an idealized initial value problem, the evolution of two semi-infinite discontinuous states.  Godunov realized it could be used as for a finite period of time (a time step size) as a semi-analytical portion of a numerical solution.

Things percolated along for a decade or so.  Artificial viscosity became the workhorse method across the World at nuclear weapons’ Labs.  The basic method developed relatively little save a handful of improvements in the early 50’s in Los Alamos, and extensions to two dimensions (this is almost certainly an overly harsh assessment).  Lax worked with Wendroff to introduce a second-order method that improved on Lax-Friedrichs, but gave oscillations too.  Lax also sorted through a general theory on the mathematics of hyperbolic conservation laws.  This work formed much of the basis for his being awarded the Abel prize.  Godunov, on the other hand, received little support and eventually quit working on fluid dynamics, and moved into numerical linear algebra while moving to Siberia in 1969.

Then everything changed.  At virtually the same time Godunov gave up the importance of Godunov’s work was rediscovered by two young researchers in the West. In fact it was never lost as noted by the coverage of Godunov’s method in Richtmyer and Morton’s book, “Finite Difference Methods for Initial Value Problems”.  He really was ignored in the Soviet Union. First, Jay Boris overcame Godunov’s barrier theorem with flux corrected transport, then by Bram Van Leer who did the same while embracing Godunov’s method as a framework as well.   Actually, Boris’ papers don’t mention Godunov, but upon reflection it is clear that he was aware of Godunov’s work and its importance.  Finally Kolgan did little appreciated work that also extended Godunov’s work directly, but went largely unnoticed until Van Leer’s recent work to bring his work to the attention of the International community.  Unfortunately Kolgan died before his work was discovered or completed.  Boris’ work was popular in the (plasma) physics community, while Van Leer’s approach had more traction with astrophysicists and aerospace engineers.  With Van Leer’s work, Riemann solvers became an important part of numerical methods and spawned what James Quirk described as a “cottage industry” in the creation of them.

Artificial viscosity kept on advancing at the Labs. Moving to three dimensions and incorporating innovations from other fields such as the limiters that Boris and Van Leer used to overcome Godunov’s barrier theorem.   Even more recently a couple of brilliant French researchers (Bruno Depres and Pierre-Henri Maire) have brought Riemann solvers to multidimensional Lagrangian codes, basically allowing Godunov’s method to improve upon the methods evolved from Richtmyer and Von Neumann’s original efforts.

It is largely a matter of philosophy and the same results can be achieved via either approach.  Nonetheless, the philosophy for developing a method can greatly influence the method developed and the results achieved.  Think about the relative success for finite volume and finite element methods in inventing new approaches for compressible flows.  The finite volume philosophy has been far more successful even though finite element aficionados would claim that they can be completely equivalent.  The same applies here, philosophy matters.

So what methodology is better?  Neither, but I do favor on balance Riemann solvers.  Why?  Here are 12 reasons:

1. Riemann solvers force you analyze the equations deeply; this is a good thing, right?  The better your analysis, the better your Riemann solver and the better your method will be.
2. Riemann solvers have exact solutions, and a sequence of well-controlled approximations
3. Riemann solvers allow one to control the dissipation quite well using your deep knowledge of mathematics and the physics.
4. Riemann solvers have a deep level of nonlinearity
5. Some of the approximations are really cool.  The paper by Harten, Lax and Van Leer (HLL, SIAM Review 1983) is wonderful and the solvers derived from that work have a wonderful form.
6. Multi-dimensional approaches are around the corner courtesy of the French work mentioned at the end of the narrative
7. They can use results from artificial viscosity to make your Riemann solver better especially in the case of strong shocks (see my Computers & Fluids paper from 1999)
8. You can design your level of detail, and dissipation wave-by-wave with enormous flexibility meaning you can ignore details, but in a mindful planned manner.
9. The Riemann solver is abstracted from the details of the mesh (they are mesh independent).  This is a huge problem for artificial viscosity where most forms have an explicit choice regarding length scales.

10. All first-order monotone schemes can be placed into a single unified framework.  Godunov’s method is the least dissipative one, and Lax-Friedrich’s is the most dissipative one.

11. Riemann solvers work well in both physics and engineering settings.

12. Riemann solvers naturally work with the thermodynamics of real materials while artificial viscosity appears to have arbitrary coefficients that can be tuned by uses (not a good thing!).  Riemann solver work can provide a clear story about the true non-arbitrary nature of these coefficients (ironically recognized for the first time by another Soviet scientist, Kurapatenko).

Artificial viscosity has advantages too; here is my list of these:

1. The nonlinearity of the dissipation for shocks is clear from the get-go, with the quadratic viscosity arising from an analysis of the Riemann solution by Richtmyer.   He didn’t recognize this, but it is true.
2. The method extends to multiple dimensions sensibly
3. They perform well with multiple physical effects
4. They have a low computational overhead
5. They can use results from Riemann solvers to make themselves better (point 12 above).
6. They are robust
7. They have a direct analogy to physical dissipation and you have control over the magnitude of dissipation directly through the coefficients.
8. The basic concepts can be extended to hyperviscosity, which is a viscosity consistent with a higher-order solver.
9. Artificial viscosity was one of the first turbulence models!

In the final analysis you should probably use a combination.  My preferred mode of operation is to make the Riemann solver primary and use artificial viscosity to clean up loose ends.

Code and calculation verification are an overlooked step for assuring the quality of codes and calculations.  Perhaps people don’t have enough good reasons to do this work.  It can often be labor intensive, frustrating and disappointing.  That is a real awful sales pitch! I really believe in doing this sort of work, but the need to provide better reasoning is clear.  I’ll get at the core of the issue right away, and then expound on ancillary benefits.

In the final analysis, not doing code verification is simply asking for trouble.  Doing code verification well is another matter altogether, but doing it well is not as necessary as simply doing it.  Many developers, students, and researchers simply compare the solutions to benchmarks using visual means (i.e., comparing to the benchmark solution in the infamous eyeball norm, or the “viewgraph” norm if its presented).  This is better than nothing, but not by much at all.  It is certainly easier and more convenient to simply not verify calculations.

Very few of us actually create codes that are free of bugs (in fact I would posit none of us). To not verify is to commit an act of extreme hubris.  Nevertheless, admitting one’s intrinsic fallibility is difficult; dealing with the impact of one’s fallibility is inevitable.

So, without further philosophizing, here is my list:

1. Don’t assert that your code is correct; prove it’s correct.  This is the scientific method; respect it, and apply it appropriately.  Treat a code as you would an experiment, and apply many of the same procedures and measures to ensure quality results.
2. Mistakes found later in the code development process are harder and more expensive to fix.  There is vast evidence of this and I recommend reading the material on the Capability Maturity Model for Software, or better yet Steve McConnell’s book, “Code Complete,” which is a masterpiece.
3. Once completed (and you aren’t ever really done), you will be confident in how your code will preform.  You will be confident that you’ve done things correctly.  You can project this confidence to others and the results you and your handiwork produce.
4. You will find the problems with your code instead of someone else like a code customer or user. As much as we dislike finding problems, some one else finding our problems is more painful.
5. Verification will allow you to firmly establish the accuracy properties of your method if you look at error and convergence rates.  You might have to confront the theory behind you method and problem, and this might help you learn something.  All of this is really good for you.
6. Doing the same thing with your calculations will allow you to understand the error associated with solving the equation approximately.  Again, confront the available theory, or its lack of availability.  It will provide you much needed humility.
7. It is embarrassing when a numerical error is influencing or hiding a model deficiency.  Worse yet, it is badly conducted science and engineering.  Don’t be responsible for more embarrassments.
8. When you are calibrating your model (admit it, you do it!), you might just be calibrating the model to the numerical error.  You want to model physics, not truncation error, right?
9. Verification results will force you to confront really deep questions that will ultimately make you a better scientist or engineer.  Science is about asking good questions, and verification is about asking good numerical questions.
10. You are a professional, right?  Doing verification is part of due diligence, it is being the best you can be.   Adopting a personal quality mentality is important in one’s development.  If you are in the business of numerical solutions, verification is a key part of the quality arsenal.
11. You won’t really understand your code until you look really hard at the results, and verification helps you understand the details you are examining.  You will looks at your work more deeply than before and that is a good thing.
12. Conducting real error analysis can help you make sure and prove that your mesh is adequate for the problem you are solving.  Just because a mesh looks good, or looks like the thing you’re simulating isn’t evidence that it actually allows you to simulate that thing.
13. It is 2013, so come on!  Please do things in a way that reflects a modern view of how to do computational science.

The excuses for not doing verification of numerical solutions are myriad.  One of the best although it is unstated, is that verification just doesn’t work all the time.  The results are not “robust” even though the scientist “knows” the results are OK.  A couple things happen to produce this outcome:

1. The results are actually not OK,
2. The mesh isn’t in the “asymptotic” range of convergence
3. The analysis of verification is susceptible to numerical problems.

I’ll comment on each of these in turn and suggest some changes to mindset and analysis approach.

First, I’d like to rant about a few more pernicious ways people avoid verification.  These are especially bad because they often think they are doing it.  For example, let’s say you have an adaptive mesh code (this is done without AMR, but more common with AMR because changing resolution is often so easy).   You get solutions at several resolutions, and you “eyeball” the solutions.  The quantity you or your customer cares about isn’t changing much as you refine the mesh.  You declare victory.

What have you done? It is not verification, it is mesh sensitivity.  You could actually have a convergent result, but you actually don’t have proof.  The solution could be staying the same because it isn’t changing, and is in fact, mesh-insensitive.  What would verification bring to the table?  It would provide a convergence rate, and an error estimate.  In fact, error estimates are the true core of verification.  The error estimate is a much stronger statement about the influence of mesh on your solution, and you almost never see it.

Why? Quite often the act of computing the error estimate actually undermines your faith in the seemingly wonderful calculation and leaves you with questions you can’t answer.  It is much easier to simply exist in comfortable ignorance and believe that your calculations are awesome.  This state of ignorance is the most common way that people almost do verification of calculations, but fail at the final hurtle.  Even at the highest level of achievement in computational science, the last two paragraphs describe the state of affairs.  I think its rather pathetic.

OK, rant off.

Let’s get the zeroth point out of the way first.  Verification is first and foremost about estimating numerical errors.  This counters the oft-stated purpose associated with “order” verification where the rate of convergence for a numerical method is computed.  Order verification is used exclusively with code verification where a solution is known to be sufficiently differentiable to provide proof that a method achieves the right rate of convergence as a function of the mesh parameters.  It is an essential part of the verification repertoire, and it squashes a more important reason to verify, error quantification whether true or estimated.

The first point is the reason for doing verification in the first place.   You want to make sure that you understand how large the impact of numerical error is on your numerical solution.  If the numerical errors are large they can overwhelm the modeling you are interested in doing.  If the numerical errors grow as a function of the mesh parameters, something is wrong.  It could be the code, or it could be the model, or the mesh, but something is amiss and the solution isn’t trustworthy.  If it isn’t checked, you don’t know.

The second point is much more subtle.  So let’s get the elephant in the room identified, meshes used in practical numerical calculations are almost never asymptotic.  This is true even in the case of what is generously called “direct numerical simulation (DNS)” where it is claimed that the numerical effects are small.  Rarely is there an error estimate in sight.  I’ve actually looked into this and the errors are much larger than scientists would have you believe, and the rates of convergence are clearly not asymptotic.

All of this is bad enough, but there is more and it is not easy to understand.  Unless the mesh parameters are small the rate of convergence should systematically deviate from the theoretical rate in a non-trivial way.  Depending on the size of the parameter and the nature of the equation being solved, the correct convergence rate could be smaller or larger than expected.  All of this can be analyzed for ideal equations such as a linear ordinary differential equation.  Depending on the details of the ODE method, and the solution one can get radically different rates of convergence.

The third point is the analysis of numerical solutions.  Usually we just take our sequence of solution and apply standard regression to solve for the convergence rate and estimated converged solution.  This simple approach is the heart of many unstated assumptions that we shouldn’t be making without at least thinking about them.   Standard least squares relies a strong assumption about the data and its errors to begin with.  It assumes that the errors from the regression are normally distributed (i.e., Gaussian).  Very little about numerical error leads one to believe this is true.  Perhaps in the case where the errors are dominated by dissipative dynamics a Gaussian would be plausible, but again this analysis itself only holds in the limit where the mesh is asymptotic.  If one is granted the luxury of analyzing such data, the analysis methodology, frankly, matters little.

What would I suggest as an alternative?

One of the problems that plagues verification is bogus results associated with either bad data, changes in convergence behavior, or outright failure of the (nonlinear) regression.  Any of these should be treated as an outlier and disregarded.  Most common outlier analysis itself relies on the assumption of Gaussian statistics.   Again, making use of this assumption is unwarranted.  Standard statistics using the mean, and standard deviation is the same thing.  Instead one should use median statistics, which can withstand the presence of up to half the data being outliers without problems.  This is the definition of robust and this is what you should use.

Do not use a single regression to analyze data, but instead do many regressions using different formulations of the regression problem, and apply constraints to the solution using your knowledge.  If you have the luxury of many mesh solutions run the regression over various subsets of your data.  For example, you know that a certain quantity is positive, or better yet must take on values between certain limits.  The same applies to convergence rates, you generally have an idea what would be reasonable from analysis; i.e., a first-order method might converge at a rate between one-half and two.   Use these constraints to make your regression fits better and more guaranteed to produce results you can use.  Make sure you throw out results that show that your second-order method is producing 25^th order convergence.  This is simply numerical garbage and there is no excuse.

At the end of this you will have a set of numerical estimates for the error and convergence rate.  Use median statistics to choose the best result and the variation from the best so that outliers along the way are disregarded naturally.

Verification should (almost) always produce a useful result and injecting ideas from robust statistics can do this.

Going beyond this point leads us to some really beautiful mathematics that is hot property now (L1 norms leading to compressed sensing, robust statistics, …).   The key is not using the standard statistical toolbox without at least thinking about it and justifying its use.  Generally in verification work it is not justifiable.  For general use a constrained L1 regression would be the starting point, maybe it will be more generally available soon.

Verification is an important activity for numerical computation that is too seldom done. Almost anybody with a lick of common sense will acknowledge how important it is and vital to achieve.  On the other hand very few actually do it.

Why?

 It is where arrogant hubris intersects laziness.  It is where theory goes to die.  Verification produces an unequivocal measure of your work.  Sometimes the result you get is unpleasant and sends you back to the drawing board.  It basically checks whether you are full of shit or not.  Most people who are full of shit would rather not admit it. 

 What is actually worse is to do verification of someone else’s work.  This puts one in the position of telling someone else that they are full of shit.  This never makes you popular or welcome.  For this reason I’ve compared V&V work to the “dark side” of the force, and by virtue of this analogy, myself to a Sith lord.   At least it’s a pretty cool bad guy.

 Areas where I’ve worked a good bit are methods for shock physics.   There is a standard problem almost everyone working on shock physics uses to test their method, Sod’s shock tube.  Gary Sod introduced his problem is a 1978 paper in the Journal of Computational Physics to compare existing methods.  By modern standards the methods available in 1978 were all bad. This was only presented as comparisons between the numerical solution, and the exact solution graphically.  This is important to be sure, but so much more could be done. Sod’s problem is a Riemann problem, which can be solved and evaluated exactly (actually solved accurately using Newton’s method).  The cool thing with an exact solution is that you can compute the errors in your numerical solution.   Unfortunately Sod didn’t do this. 

 Again, why?

 Instead he reported run time for the methods, as if to say, “all of these suck equally, so what matters is getting the bad answer the fastest”.  A more technical response to this query is that shock captured solutions to problems with discontinuous solutions (Sod’s problem, has a rarefaction, contact discontinuity and a shock wave) are intrinsically limited to first-order accuracy.  At a detail level the solution is to Sod’s shock tube is typically less than first-order accurate because the contact discontinuity’s solution convergences at less than first-order and eventually this dominates the numerical error under mesh refinement (covered in a paper by Banks, Aslam, and Rider(me)).  The mindset is “well its going to be first order anyway so what’s the point of computing convergence rates?” 

 I’m going to tell you what the point is.  The point is that the magnitude of the errors in the solutions is actually a lot different even among a whole bunch of methods of similar construction and accuracy.  Not every method is created equal and there is a good bit of difference between the performance of Godunov’s original method, Van Leer’s method, Colella’s refined PLM method, WENO, PPM, etc.  Ultimately what someone really cares about is the overall accuracy of the solution for amount of computational resource, as well as scheme robustness and extensibility. 

 So going back to the reality, people continued to use Sod’s problem and presented results with a graphical comparison usually at a single grid, and a comparison at runtime.  I never saw anyone present the numerical error.  People would only provide numerical errors for problems where they expected their method to achieve its full order of accuracy, i.e., smooth problems like advecting a Gaussian function or Sine wave.   The only exception came from the flux-corrected transport world where they did present numerical accuracy for the propagation of square waves, but not systems of equations. 

 In 2004-2005 Jeff Greenough and I attacked this low hanging piece of fruit.  Jeff’s Lab was using a weighted ENO (WENO) code developed at Brown University (where the inventor of WENO is a professor).  Jeff worked on a code that used Colella’s PLM method (basically an improved version of Van Leer’s second-order Godunov method).  I had my own implementations of both methods in a common code base.  We undertook a comparison of the two methods head-to-head, with the following question in mind: “what method is the most efficient for different classes of problems?”  We published the answer in the Journal of Computational Physics. 

 Before giving the answer to the question a little more background is needed.  WENO methods are enormously popular in the research community with the fifth-order method being the archetype.  I would say that the implementation of a WENO method is much more elegant and beautiful than the second-order Godunov method.  It is compact and modular, and actually fun to implement.   It is also expensive.  

 We used a set of 1-D problems: Sod’s shock tube, a stronger shock tube, a popular blast wave problem, and a shock hitting smooth density perturbations.  For the first three problems, the 2^nd order method either out performed or was even with the 5^th order method at the same mesh.  On the final problem WENO was better, until we asked the efficiency question.  Accounting for runtime to achieve the same accuracy, the 2^nd order method outperformed WENO on every problem.  If the problem was shock dominated, the difference in efficiency wasn’t even close.  WENO was only competitive on the problem with lots of smooth structure.  On the advection of a Gaussian pulse WENO wins at almost any level of error.

 This asks a deeper question, why would a second-order method be more efficient?

 This answer gets to the heart of method design.  The 2^nd order method cheaply introduces some elements of high-order methods, i.e., it uses 4^th order approximations to the slope.  Its nonlinear stability mechanism is cheap, and it uses a single step time integrator rather than a multistep integrator.  The single step integrator provides a more generous time step size.  The combination of the larger time step, single step integrator and simpler nonlinear stability mechanism equates to a factor of six in runtime.  In addition to accommodate the fifth-order formal accuracy, the WENO method actually increases the dissipation used near shocks, which increases the error.

 The bottom line is that verification is about asking good questions, answering these questions, and then going back to the well… it is in essence the scientific method. 

The name of the original successful shock capturing method is “artificial viscosity,” and it is terrible.  John Von Neumann and Robert Richtmyer were both mathematical geniuses, perhaps, at least in this case, poor at marketing.   To be fair, they struggled with the name, and artificial viscosity was better than “mock” or “fictitious” viscosity, which Richtmyer considered in his earlier Los Alamos reports (LA-671 and LA-699).  Nonetheless, we are left with this less than ideal name. 

I’ve often wished I could replace this name with “shock viscosity” or better “shock dissipation” because the impact of artificial viscosity is utterly real.  It is physical and necessary to compute the solutions of shocks automatically without resolving the ridiculously small length and time scales associated with the physical dissipation.  In a very true sense, the magnitude of viscosity used in artificial viscosity is the correct amount of dissipation for the manner in which the shock is being represented.

The same issue comes about in the models of turbulence.  Again, the effects of nonlinearity enormously augment the impact of physical dissipation.  The nonlinearity is associated with complex small-scale structures whose details are unimportant to the large-scale evolution of the flow.  When simulating these circumstances numerically we are often only interested in the large-scale flow field and we can use an enhanced (nonlinear) numerical viscosity to stably compute the flow.  This approach in both shocks and turbulence has been ridiculously successful.  In turbulence this approach goes by the name implicit large eddy simulation.  It has been extremely controversial, but in the last decade this approach has grown in acceptance.

The key to this whole line of thinking is that dissipation is essential for physical behavior of numerical simulations.  While too much dissipation is undesirable, too little dissipation is a disaster.  The most dangerous situation is when a simulation is stable (that is, runs to completion), and produces seemingly plausible results, but has less dissipation than physically called for.  In this case the simulation will produce an “entropy-violating” solution.  In other words, the result will be unphysical, that is, not achievable in reality.  This is truly dangerous and far less desirable than the physical, but overly dissipated result (IMHO).  I’ve often applied a maxim to simulations, “when in doubt, diffuse it out”.   In other words, dissipation while not ideal (a play on words!) is better than too little dissipation, which allows physically unachievable solutions to persist.

Too often numerical practitioners seek to remove numerical dissipation without being mindful of the delicate balance between, the dissipation that is excessive, and the necessity of dissipation for guaranteeing physically relevant (or admissible) results.  It is a poorly appreciated aspect of nonlinear solutions for physical systems that the truly physical solutions are not those that have no dissipation, but rather produce a finite amount of dissipation.  This finite amount is determined by the large-scale variations in the flow, and is proportional to the third power in these variations. 

Very similar scaling laws exist for ideal incompressible and compressible flows.  In the incompressible case, Kolmogorov discovered the scaling law in the 1940’s in the form of his refined similarity hypothesis, or the “4/5ths” law.  Remarkably, Kolmogorov did this work in the middle of the Nazi invasion of the Soviet Union, and during the darkest days of the war for the Soviets.  At nearly the same time in the United States, Hans Bethe discovered a similar scaling law for shock waves.  Both can be written in stunningly similar forms where the time rate of change of kinetic energy due to dissipative processes (or change in entropy), is proportional to the third power of large-scale velocity differences, and completely independent of the precise value of viscosity. 

These scaling laws are responsible for the success of artificial viscosity, and large eddy simulation.  In fact, the origin of large eddy simulation is artificial viscosity.  The original suggestion that led to the development of the first large eddy simulation by Smagorinsky was made by Von Neumann’s collaborator, Jules Charney to remove numerical oscillations from early weather simulations in 1956.  Smagorinsky implemented this approach in three dimensions in what became the first large eddy simulation, and the first global circulation model.  This work was the origin of a major theme in turbulence modeling, and climate modeling.  The depth of this common origin has rarely been elaborated upon, but I believe the commonality has profound implications.  

At the very least, it is important to realize that these different fields have a common origin.  Dissipation isn’t something to be avoided at all costs, but rather something to manage carefully.  Better yet, dissipation is something to be modeled carefully, and controlled.  In numerical simulations, stability is the most important thing because without stability, the simulation is worthless.  The second priority is producing a physically meaningful (or realizable) simulation.  In other words, we want a simulation that produces a result that matches a situation achievable in the real world.  The last condition is accuracy.  We want a solution that is as accurate as possible, without sacrificing the previous two conditions. 

Too often, the researcher gets hung up on these conditions in the wrong order (like prizing accuracy above all else).   The practitioner applying simulations in an engineering context often does not prize accuracy enough.   The goal is to apply these constraints in balance and in the right order (stability then realizability then accuracy).  Getting the order and balance right is the key to high quality simulations.