Numerical methods

One
distinguishing component of the present SciDAC project is the
assessment and development of numerical methods suitable for problems
with interactions between turbulence, shock waves, and material
interfaces. Computational turbulence is a field where the subtle
details of the numerical method can have a dramatic impact on the
accuracy of the results -- one may even draw the wrong scientific
conclusions if one is not careful. The literature is filled with
numerical methods that perform well on shock waves or turbulence in
isolation, but much less is known about how these methods treat shocks
and turbulence simultaneously.

The
computational challenge of predicting shock/turbulence interactions
stems from the fundamentally different physics at play. Shock waves
are thin regions wherein flow properties change rapidly over a
distance roughly equal to the molecular mean free path; hence, they
are essentially strong discontinuities in the flow field. Turbulence,
on the other hand, is a chaotic phenomenon with broadband spatial and
temporal scales of motion. Most shock-capturing methods rely on strong
numerical dissipation to artificially smooth the discontinuity, such
that it can be resolved on the computational grid. However, the
artificial dissipation necessary for capturing shocks has a
deleterious effect on turbulence. An additional problem is the fact
that shock-capturing schemes are typically based on one-dimensional
Riemann solutions that are not strictly valid in multiple dimensions.
This can lead to anisotropy errors and grid-seeded perturbations.
Other complications arising from upwinding, flux limiting, operator
splitting, etc..., can seriously degrade performance and generate
significant errors.

The first
milestone of the present project is a comparison between a
set of promising methods on a carefully chosen suite of benchmark
problems. We have divided this item into two sub-items, by first
considering single-fluid cases only. The chosen
single-fluid problems represent the key physics with increasing
complexity:

1. Taylor-Green vortex: tests the accuracy and numerical stability on
broadband turbulence without shocks.

2. Spherical Noh implosion: tests the capability to achieve full shock
compression, correctly predict the shock speed, and maintain spherical
symmetry on an infinite Mach number shock without turbulence.

3. 1D Shu-Osher problem: tests the capability of predicting a
shock/entropy interaction, specifically to avoid damping of the
post-shock entropy fluctuations.

4. 2D shock/entropy/vorticity interaction: extends the 1D Shu-Osher
problem to two dimensions, with the associated induced unsteady
shock-movement.

5. High-Mach number decaying turbulence problem: includes broadband
vortical, acoustic, and entropy motions, as well as turbulence-induced
weak shocks.

We point out
that this benchmark study was undertaken in the spirit
of learning about the strengths and weaknesses of different
numerical approaches, with the obvious hope of devising future
improvements. Therefore all problems were computed by four different
numerical methods (i.e., codes) without *ad hoc* "tuning" for
each problem. It is important to realize that most methods can be "tuned" to give accurate results for the first four benchmark problems
-- but in the process sacrificing accuracy (or even stability) on
the remaining problems. No such tuning was performed here.
The four codes all represent different computational approaches, but
every one could be considered state-of-the-art within its context. A
brief summary of each code is given here:

1. The *Miranda *code (Cook, Phys. Fluids 2007}) uses very
high-order operators (8th-10th order) to approximate the
Navier-Stokes equations, with augmented artificial fluid
properties (like viscosity) to capture shocks and other
discontinuities. The guiding philosophy behind this method is to
first regularize the equations (through the artificial properties)
and then to solve the regularized equations using a very
high-order method. The code has been run on up to 65,536
processors on the BlueGene/L machine.

2. The Hybrid code (Larsson *et al.*, CTR Annual Research Briefs
2007)
is based on the guiding idea that turbulence and shock waves are
fundamentally different phenomena and thus should be treated by
different numerical methods. Hence the *Hybrid* method relies
on a sensor to identify regions of shock waves, and then applies a
shock-capturing (WENO) scheme in these regions with a high-order central
difference scheme elsewhere. This hybridization ensures minimal
numerical dissipation in the method. The code has been run on up
to 4,096 processors on the Franklin Cray XT-4 machine.

3. The *ADPDIS3D* code (Yee and Sjogreen, J. Comput. Phys. 2007)
combines the best of non-dissipative high order spatial schemes with
shock-capturing methods via an efficient containment of numerical
dissipation. These low dissipative high order methods include limiting
and filtering with flow sensors. Flow sensors are used as an adaptive
procedure to analyze the computed flow data and indicate the amount,
location and type (e.g., spectral filter, compact filter or nonlinear
(shock) filter) of built-in numerical dissipation that can be
eliminated or further reduced. Such a filter method consists of two
steps, a full time step using a spatially high order non-dissipative
base scheme, followed by an adaptive multistep filter consisting of
the products of wavelet-based flow sensors and linear and nonlinear
numerical dissipations. The idea to control numerical dissipation is
very general and can be used in conjunction with high order spectral,
spectral element, finite volume and finite difference compact or
non-compact central spatial base schemes. Any shock-capturing scheme
can be used as nonlinear dissipation, (e.g., the dissipative portion
of the TVD, MUSCL, or (W)ENO methods), usually with flux limiters. The
linear filter can be the standard spectral or compact filter, or the
product of a high-order linear dissipation and an appropriate flow
sensor. By design, the flow sensors, spatial base schemes, and linear
and nonlinear dissipation models are stand-alone modules. Therefore, a
whole class of low dissipative high-order filter schemes can be easily
derived. With this family of adaptive filter schemes, a wide spectrum
of flow regimes can be simulated using the same code. Central spatial
base schemes of order up to 14 and the dissipative portion of
shock-capturing schemes of order up to 9 have been implemented.

4. The final code differs from the others by using shock-fitting,
rather than shock-capturing, to handle the shock waves.
In other words, the exact shock-jump relations are imposed at the
shock, which is therefore treated exactly.
The method yields superior accuracy for geometrically simple shocks,
but has yet to be applied to truly complex shock-structures.

Figure 1. Entropy profiles after
shock-turbulence interaction. Results show that: 1) the Miranda
code performs extremely well on density disturbances; 2) standard
WENO is highly dissipative, which kills the entropy fluctuation;
3) the Hybrid approach significantly improves on this, despite
using the identical WENO scheme a the shock. |
Figure 2. Temperature rms-profiles in
decaying isotropic turbulence. Results show that: 1) The Miranda
code excessively damps thermodynamic fluctuations; 2) standard
WENO and the Hybrid method are less dissipative. |

The benchmark comparison has brought several points to the
surface. Most important is that the comprehensive nature of the
problem suite has unmercifully brought out some of the weaknesses in
every method, and one could argue that some of these findings were
made possible only through the comprehensive suite. Figure 1 shows the profiles of entropy
after the shock-interaction in the Shu-Osher problem. The standard
WENO method is severely dissipative, and the *Hybrid*
philosophy of minimizing the dissipation by switching schemes is a
vast improvement. The *Miranda* method is highly accurate for
this problem, being capable of essentially capturing the full
entropy fluctuation even on this coarse (200 points) grid.
Figure 2 shows the evolution of
temperature fluctuations in decaying isotropic turbulence on coarse
(64^3) grids. Here the results are the opposite: this shows how
the behavior of numerical methods on computational turbulence with
shock waves is complex. The under-prediction by the *Miranda*
method is most likely caused by the artificial bulk viscosity;
improvements to the method that alleviate this issue have been and
continue to be investigated.

In addition to the sample giving above, some other findings from the
comprehensive assessment include:

1. While high-order WENO schemes are state-of-the-art in terms of
handling pure shock waves, they add significant numerical
dissipation on broadband turbulence and decrease the well-resolved
range of scales by a factor of two or more.

2. The solution-adaptive philosophy in the *Hybrid* method
drastically reduces the numerical dissipation, leading to much
more accurate treatment of broadband turbulence.

3. More work is needed on reliable shock sensors, that work in a
broad range of problems.

4. The *Miranda* treatment of shocks through artificial
fluid properties works well on well-resolved motions, but needs
improvement on broadband turbulence.