V&V » KozCG: Sod's shocktube

This example uses KozCG in Inciter to compute the Sod problem [1], widely used to test numerical methods to model discontinuous solutions. The problem consists of a pipe in which two different states of a gas is initially separated by a diaphragm. When the diaphragm is raptured at t=0, a wave structure is created and propagated in time. While the solution is 1D, we modeled the process in 3D. The available analytic solution enables verification of the accuracy of the numerical method and the correctness of the software implementation.

Equations solved

The system of equation solved is

\[ \begin{split} \frac{\partial U}{\partial t} + \frac{\partial F_j}{\partial x_j} = 0,\qquad U = \begin{Bmatrix} \rho \\ \rho u_i \\ \rho E \end{Bmatrix}, \enskip F_j = \begin{Bmatrix} \rho u_j \\ \rho u_iu_j + p\delta_{ij} \\ u_j(\rho E + p) \end{Bmatrix} \end{split} \]

detailed at the KozCG page describing its numerical method.

Problem setup

The problem is modeled as a 3D tube with a circular cross section of unit length with a diameter of 0.05, with the specific heat ratio as 1.4, and the following initial conditions

\[ \begin{split} u_i & = 0.0 \\ \rho(x_1) & = \left\{ \begin{array}{cc} 1.0 & x_1 < 0.5 \\ 0.125 & x_1 \ge 0.5 \end{array} \right. \\ p(x_1) & = \left\{ \begin{array}{cc} 1.0 & x_1 < 0.5 \\ 0.1 & x_1 \ge 0.5 \end{array} \right. \end{split} \]

where $x_1$ is the distance along the length of the tube. Dirichlet boundary conditions were applied on all flow variables sampling their initial values at the end-faces of the domain and symmetry (free-slip) conditions were set along the surface of the tube. The numerical solution was advanced using the CFL number 0.5 using four different meshes whose properties are list in the following table:

MeshPointsTetrahedrah
07,29434,1810.008393
151,794273,4480.004138
2389,1392,187,5840.002060
33,014,27717,500,6720.001028

Here h is the average edge length along the length of the tube.

Code revision to reproduce

To reproduce the results below, use code revision fadd412 and the control file below.

Control file

-- vim: filetype=lua:

print "Sod shocktube"

term = 0.2
ttyi = 10
cfl = 0.5

solver = "kozcg"
fctclip = true
fctsys = { 1, 2, 5 }

part = "rcb"

problem = {
  name = "sod"
}

mat = { spec_heat_ratio = 1.4 }

bc_sym = {
  sideset = { 2 }
}

bc_dir = {
 { 1, 1, 1, 1, 1, 1 },
 { 3, 1, 1, 1, 1, 1 }
}

fieldout = {
  iter = 10000
}

diag = {
  iter = 1,
  format = "scientific"
}

Run using mesh 1 on 4 CPUs

./charmrun +p4 Main/inciter -i tube01.exo -c sod.q

Visualization

ParaView can be used for interactive visualization of the numerically computed fields as

paraview out.e-s.0.4.0

The surface mesh colored by the fluid density at t=0 and t=0.2 on mesh 1 are depicted below.

Image Image

The numerical solution is also extracted at t=0.2 along the center of the tube. The density, velocity, total energy, and pressure are depicted below for the successively finer meshes together with the analytical solution.

Image Image Image Image

References

  1. J.R. Kamm, Enhanced verification test suite for phyics simulation codes, LA-14379, SAND2008-7813, 2008.