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.