2D DNS Godunov-type solver or why I love Julia HPC?

Main features:

Solver: Density-based unstructured FVM
Spacial discretization: MUSCL 2nd order
Limiters: MinMod
Time discetization: Euler first-order (TVD RK2)
Gradients: Combined Green-Gauss noded based and LSQ
Inviscid fluxes: AUSM+ (Roe flux splitting)
Parallelization: multi-threading
GPU acceleration: inviscid fluxes (too slowly yet)

GPU calcs are still very slow! Challenges:

- copy/paste memory CPU-GPU operations
- x64 double precision math for inviscid Flux calculations
- most consumer Nvidia cards are dedicated to x32
- even Nvidia V100 and A100 types lose its performance from x32-> x64
Hybrid CPU-GPU solver performance is 50% slowly to CPU solver

Summary and outlook

The prime goal is to develop the in-house innovate DNS Riemann-based solver for advanced computation of the turbulent separated flows, heat transfer and combustion. The original version has been implemented using C++ environment, however, the recent development continues based on the super fast Julia language accompanied with the state-of-the-art hybrid CPU-GPU parallelism.

Numerical excercises (II): 2D Taylor-Green vortex for the diffusion term

The Taylor–Green vortex is an unsteady flow of a decaying vortex, which has an exact closed form solution of the incompressible Navier–Stokes equations. Here, implementation of the Diffusion Term approximation is tested within the compressible DNS framework. The box domain of size [-2π → 2π,-2π → 2π] is considered. The Reynolds and Mach numbers are set as Re = 4000 and M = 0.1, respectevily. Properties of air (γ = 1.4) are set as: P = 100000 [Pa], ρ = 1 [kg/m3] and μ = 1e-2 [Pa*s]. Pressure contours and decay of the total kinetic energy up to t* = t*U/L = 10 are provided. The results obtained using the triamgular-cell grid (b) are quite consitent with the analytical solution (a).


Numerical excercises (I): GPUx32 for the inviscid fluxes

Here, the prominent test case by Woodward Colella (1984) is used to compare solutions obtained by the convention solver (CPUx64) and hybrid solver, when inviscid fluxes are computed using GPUx32. Time evolution of the L2 norm of numerical diffusion displays the increasing divergence between them. Meanwhile, contours plotted at t = 4 are quite similar yet.

Test#6> Interaction of an oblique shock with a spatially-developing shear layer

In this study the interaction of an oblique shock with a 2d-developing shear layer (a) is analyzed by means of DNS. The purpose of this numerical experiment is to enhance flow mixing by the incident shoch. Here, the inlet vorticity thickness δ is 1.44e-3 [m] and the convective Mach number, Mc = 0.48. The oblique shock wave angle Β is 33 ℃. White noise perturbations (2.5%) are added to the inlet velocity to trigger flow transition. The numerical schrilen visualization of the instantaneous flow field, obtained on the 1200x480 triangular (b) and 1800x640 quad (c) grids, displays local features of the shear layer and can be further assessed by considering the Richtmyer-Meshkov instability.

References: Boukharfane, R., Ferrer, P.J.M., Arnaud Mura, A., Giovangigli, V., On the role of bulk viscosity in compressible reactive shear layer developments, ICCFD10-2018-0086

Test#5> Supersonic flow over a circular cylinder at M = 3.5

Numerical and experimental studies indicate that the supersonic flow over a stationary solid obstacles at high Re-numbers, becomes unstable with large asymmetric recirculation zone and vortex shedding in downstream. The flow visualization, based on the Numerical Schlieren (M = 3.5), shows many interesting features both in the wake and in the region between the detached bow-shock and the obstacle, which can be a subject for the further research to validate the numerical method.

References: Chaudhuri, A., Hadjadj, A., Chinnayya, A., On the use of immersed boundary methods for shock/obstacle interactions, J. Comput. Phys, 230, 1731-1748 (2011)


Test#4> Schardin’s problem

A vertical plane shock wave is striking a two-dimensional wedge. Both, numerical interferogram (left) and Schlieren (right) pictures highlight the complex shock–vortex dynamics of the Schardin’s problem. The incident shock wave runs from left to right at M = 1.3. An unsteady flow is developed around the obstacle accompanied by the complex wave pattern, including a cylindrical shock, reflected off the wedge surface, and a shock wave and the Mach stem, which propagate normal to the wedge surfaces, joining the incident and reflected shock waves at a triple point.

References: Schardin, H., High frequency cinematography in the shock tube, J. Photo. Sci., 5(2), 17-19 (1957)


Test#3> A single- and double-Mach reflection problem

This video shows numerical investigation to replicate the the single and double oblique shock-wave reflection configurations studied in the experiments of Henderson and Gray (1981). Here, the Euler equations are solved.

References: Henderson, L.F., Gray, P.M., Experiments on the diffraction of strong blast waves, Proc. R. Soc. Lond. A, 377(1771), 363-378 (1981)


Test#2> A Mach 3 Wind Tunnel With a Step

For more than fifty years supersonic Mach 3 flow over forward step is a canonical benchmark to test different computation schemes. The image demonstrates the performance of the solver to reproduce this challenging benchmark on the baseline 80x20 grid. The corner of the step is a singularity point, which can produce significant artificial numerical errors depend on applied boundary conditions and numerical schemes used. Here, solution was obtained on the grid with the smoothed corner. No any special treatment of the singularity point was used for both cases. One can see clearly that the solver demonstrates the adequate performance to reproduce all flow features.

References: Woodward, P, Colella, P., The numerical simulation of two-dimensional fluid flow with strong shocks, J. Comput. Phys., 54, 115-173 (1984)

Test#1> 2D shock reflection problem

In this problem, the incident shock angle is φ = 29˚and the free-stream Mach number and M = 2.9. The computational domain is set to [0 4.1] and [0 1] with a uniform triangular grid with size of 0.05. The appropriate analytical boundary conditions (supersonic inflows) are applied along the left and top boundaries of the domain, while the right and lower boundaries are treated as supersonic outflow and rigid solid surface, respectively.

References: Yee, H.C, Warming, R.F., Harten, A., Implicit Total Variation Diminishing (TVD) Schemes for Steady-State Calculations, NASA TM 84342 (1983)