“The quantitative influence of these terms can be made as small as one wishes by choice of a sufficiently fine mesh for the numerical integrations.” – John von Neumann and Robert D. Richtmyer

The implementation of such a method is one of my most enduring achievements and contributions to national security. It was put into a major production code for Los Alamos. The methodology remains important and may be a path forward into a development of improved capability across a whole spectrum of codes. It can be explained with some simple principles. I hope to do that here for a broader audience.

“A new algorithm for the volume tracking of interfaces in two dimensions is presented. The algorithm is based upon a well-defined, second-order geometric solution of a volume evolution equation.” – William J. Rider and Douglas B. Kothe

Another significant professional accomplishment was the paper I wrote with Doug Kothe on volume tracking using volume-of-fluid methods. It’s got nearly 2,500 citations for the Journal of Computational Physics article we wrote. There was a related Annual Reviews Fluid Mechanics article that came out a couple of years after our article. I’m sure our article would probably have 4,000 citations were it not for that. It has a similar citation count to ours. Nonetheless, it is a great accomplishment, and I’ve been told that the code that I wrote (almost as a hobby and part of an unpaid project at Los Alamos) is still being used. It’s the same production code I contributed to for steepened methods.

I didn’t actually put the code that I wrote into the production code. This was the work of a colleague, Ed Dendy, who did the heavy lifting. I was the idea man, and we worked in concert with one another.

In the late ’90s, we discovered that the use of adaptive mesh refinement was insufficient to retain sharp interfaces. This in physical circumstances where those sharp interfaces should not diffuse. The situation where this becomes most acute is when the flow is expanding, but there’s a material interface embedded. In this situation, the material interface grows almost without bound and becomes extremely diffuse. For some simulations, this is unacceptable and must be mitigated.

This shows up in a test problem that is quite common. The Woodward-Colella interacting blast wave problem at late time has a contact discontinuity that becomes embedded in an expansion wave. For most methods, the contact becomes extremely distended. This ends up being a very critical test for the ability of the code to handle a sharp interface.

Thus, we needed to add the ability to retain sharp interfaces to this code. There were two means to do this:

1. The volume tracking that I discussed in the paper with Kothe. I wrote the toolbox code to demonstrate how to do this in an object-oriented way. I’ve discussed this before, where I learned that the standard way of implementing volume tracking with a series of nested if-then loops had too much cyclomatic complexity to be easily maintained or extended.

2. There is a second way to maintain a sharp interface. This is the notion of what’s often called steepeners. This ends up being a controlled way of combating numerical diffusion. This often involves anti-diffusion, and there are a host of methods for doing this. Some are physically motivated, effectively looking like surface tension; these use an analogy to physics to achieve sharp interfaces.

“Simplicity is prerequisite for reliability.” – Edsger W. Dijkstra

This idea originally came from Ami Harten. He introduced a method called the artificial compression method (ACM). This method worked to remove numerical dissipation. It would remove it in a safe way that removed the dissipation where it was unnecessary. ACM was amongst the set of first methods designed to overcome Godunov’s theorem. It would turn out to be useful with later generations of methods.

My mind, when developing numerical methods, works in a Godunov framework. A geometric reconstruction of the profile in a cell is used. That reconstruction informs how fluxes are computed. Below, I’m going to describe and show the performance of methods using this formulation. I’ll also show some plots to illustrate the effect. I will note that this was the project that started to show me what Claude could do that blew my fucking mind.

The key idea in all of this is that, in order to get anti-diffusion, you reverse the way you get diffusion. Numerical diffusion is often viewed as the bane of numerical simulation. This is doubly so in places like Los Alamos, where the canonical numerical method for solving problems is the Lagrangian method, where interfaces are naturally preserved as invariants. Thus, the numerical diffusion necessary in Eulerian computations is viewed as pathological. I’ve spoken before about how wrong-headed this is. The solution of hyperbolic PDEs, which include the linearly degenerate interfaces such as contacts and shears, requires a modicum of numerical diffusion in order to ensure numerically and physically stable and consistent interfaces. Without numerical diffusion, solutions are a mess. Nonetheless, it is around interfaces that much of the annoyance and distress of code users revolves.

I will just note that Lagrangian material interfaces are themselves filled with some pathologies. The most famous of these pathologies is known as wall heating, which arises most classically in terms of the Noh problem, where the heating that results from it is quite extreme. In planar geometry, this is an infinitely strong shock bouncing off a wall, and there’s a dip in the density and a rise in internal energy that is called “wall heating.”

The desire for this Lagrangian material identity-invariant point of view can be seen in our physics education. By and large, when we learn physics in college, the Lagrangian frame of reference is taken. In some cases, this is merely due to simplicity. In other cases, it has a genuine character. In terms of relationships like solid mechanics, the Lagrangian frame is king. Most of more complex dynamics has its foundation with mechanics. Radiation transport is the same. All of this matters greatly in Los Alamos.

The problem with the Lagrangian frame is that phenomena such as shear and turbulence (and not to mention diffusion of a real physical kind) will always induce the destruction of the Lagrangian frame over time. This becomes unavoidable in virtually any multidimensional complex simulation. In this case, some degree of Eulerian representation is necessary if we want to compute the physical circumstances we care about. This is where interface tracking comes in.

Volume tracking methods are generally thought of as ways to keep the interfaces completely sharp in a nearly Lagrangian approach. David Youngs introduced the idea of viewing this geometrically. Youngs introduced a linear planar interface for tracking. This built on a representation known as SLIC from Lawrence Livermore that represented the interfaces as a series of straight lines. In some ways, this was analogous to the Godunov method, where constant values represent the reconstruction in a cell. One of the issues is that this method can produce side effects that are quite profound, and you trade numerical surface tension for numerical diffusion.

With numerical surface tension, the interfaces bead and form substructures or islands that are equally physically incorrect. You trade this for excessive diffusion, so pathologies are unavoidable. Wall heating is another one you get. The key is to be cognizant of the issues and manage them.

Steepened or anti-diffusive transport methods are an alternative. It was this alternative that I divined. Ultimately, my friend Ed implemented it in the production code at Los Alamos. As fate would have it, the production code people gave the method a completely different name than you would find in the literature. It produced a sense that it was unique and separate from what was there. They called the method “interface preserving.” As with most labels, it’s not correct in any way, but it also incorrectly divorces the approach from its origin.

“Numerical analysis is the study of algorithms for the problems of continuous mathematics.” – Lloyd N. Trefethen

In the remainder of this article, I’m going to describe how these methods work and how they perform. In showing how they perform, I’m going to go on a slight tangent to present a sidebar discussion of how much current AI work using Claude specifically assisted me. It’s an example of the sort of incredible boost to productivity these tools can provide.

As they discovered with Claude, I could reproduce almost everything I could do with Mathematica in numerical analysis of numerical methods. I used to do all this work with Mathematica, which carries a heavy price tag that a retired person such as myself isn’t very eager to bear.

I started to experiment with Claude and began to go through a host of numerical analysis tests and techniques that I had used over the years to see what the AI could produce. Ultimately, I want to use Claude to reproduce the results we already know and extend those results. Once I established that I could do everything that I used to do competently and correctly, I’m ready to extend things. This post will benefit from that capability.

One of the big things that Claude does is do this with ease and speed. I estimate that in 2 days of intermittent use of Claude, I was able to produce what would have taken me a month in Mathematica and writing my own code to produce. It also created figures and effortlessly produced a report that I could look at to check the quality of the work and its correctness. The key to doing this work is writing good prompts at the start. Next, the output needs to be carefully examined and verified. The prompt and the work from AI can be wrong. It makes mistakes, and your job is to find them.

In what follows, the results were produced in this manner.

You can get the interface tracking by using an anti-upwind method. The key is to do this in a way that does not go unstable. You take the anti-upwind flux and make sure it is nonlinearly stable with a limiter. TVD theory gives you that theory. The result is quite similar to SLIC.

You can keep this simple by using the so-called Superbee limiter. This is a good, but not great solution. You can amp up the performance by making the limiter sharper, including the CFL number in the form. It is still second-order when flow is resolved, but very sharp for interfaces.

“…a steeper representation of discontinuities, particularly contact discontinuities.” – Phillip Colella and Paul R. Woodward

Starting with Colella and Woodward’s PPM method, this can be improved. A way to do this is to use an anti-upwind high-order method. For example, you find an anti-upwind method using a minmod limiter. This is what we used in the IP method. The key is knowing that the more diffusive a method is in the upwinding, the more anti-diffusive it will be. This is what you do for PPM, ENO, and WENO methods to great effect. PPM also had a heuristic contact-material interface detection. It’s all that simple.

I’ve included some figures to demonstrate what these methods look like geometrically and in performance. Enjoy. This is a really simple and relatively bulletproof method. Going to multi-dimensions and avoiding pathological results is important. Anti-dissipation is playing with fire and must be done with extreme care.

“It can scarcely be denied that the supreme goal of all theory is to make the irreducible basic elements as simple and as few as possible without having to surrender the adequate representation of a single datum of experience.” – Albert Einstein

References

Rider, William J., and Douglas B. Kothe. “Reconstructing volume tracking.” Journal of computational physics 141, no. 2 (1998): 112-152.

Scardovelli, Ruben, and Stéphane Zaleski. “Direct numerical simulation of free-surface and interfacial flow.” Annual review of fluid mechanics 31, no. 1 (1999): 567-603.

Noh, William F., and Paul Woodward. “SLIC (simple line interface calculation).” In Proceedings of the fifth international conference on numerical methods in fluid dynamics June 28–July 2, 1976 Twente University, Enschede, pp. 330-340. Berlin, Heidelberg: Springer Berlin Heidelberg, 2005.

Harten, Amiram. “The artificial compression method for computation of shocks and contact discontinuities.” Maths Comput 32 (1978): 363.

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

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

Shu, Chi-Wang. “High order ENO and WENO schemes for computational fluid dynamics.” In High-order methods for computational physics, pp. 439-582. Berlin, Heidelberg: Springer Berlin Heidelberg, 1999.