V&V » RieCG: Vortical flow

This example uses RieCG in Inciter to compute a problem whose analytic solution is known, therefore it can be used to verify the accuracy of the numerical method and the correctness of the software implementation.

The purpose of this problem is to test velocity errors generated by spatial operators in the presence of 3D vorticity and in particular the superposition of planar and vortical flows, analogous to vorticity stretching. The derivation of this test problem, using the method of manufactured solutions, is given in [1], preceded by a simplified version published in [2]. The combination of vorticity and velocity gradients is a fundamental source of kinetic energy in flows that transition to shear-driven turbulence, and hence an important ingredient of computing many practical engineering and atmospheric flows.

Equations solved

The system of equation solved is

\[ \begin{split} \frac{\partial U}{\partial t} + \frac{\partial F_j}{\partial x_j} = S,\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}, \enskip S = \begin{Bmatrix} S_\rho \\ S_{u,i} \\ S_E \end{Bmatrix} \end{split} \]

detailed at the RieCG page describing its numerical method, with the exact solution

\[ \begin{split} \rho & = 1 \\ u_i(x_i) & = \left( \begin{array}{ccc} \alpha x_1 - \beta x_2 \\ \beta x_1 + \alpha x_2 \\ - 2 \alpha x_3 \end{array} \right) \\ p(x_3) & = p_0 - 2 \rho \alpha^2 x_3^2 \\ e & = p(x_3) \rho \left( \gamma -1 \right) \end{split} \]

and the source terms

\[ \begin{split} S_\rho & = 0 \\ S_{u,i} & = \left( \begin{array}{c} \rho \left( \alpha^2 - \beta^2 \right) x_1 - 2 \rho \alpha \beta x_2 \\ \rho \left( \alpha^2 - \beta^2 \right) x_2 + 2 \rho \alpha \beta x_1 \\ 0 \end{array} \right) \\ S_E & = u_i S_{u,i} + \frac{8 \rho \alpha^3 x_3^2}{\gamma - 1} \end{split} \]

For more details, see Papers.

Problem setup

The simulation domain is a cube centered around the point {0,0,0}. The initial conditions are sampled from the analytic solution at t=0. We set Dirichlet boundary conditions on the sides of the cube, sampling the analytic solution. The numerical solution does not depend on time and approaches steady state due to the source term which ensures equilibrium in time. As the numerical solution approaches a stationary state, the numerical errors in the flow variables converge to stationary values, determined by the combination of spatial and temporal errors which are measured and assessed.

Code revision to reproduce

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

Control file

-- vim: filetype=lua:

print "Euler equations computing vortical flow"

term = 3.0
ttyi = 10

solver = "riecg"

dt = 0.002      -- 750K
--dt = 0.001      --   6M
--dt = 0.0005      --  48M

part = "rcb"

problem = {
  name = "vortical_flow",
  alpha = 1.0,
  kappa = 1.0,
  p0 = 10.0,
}

mat = { spec_heat_ratio = 5/3 }

bc_dir = {
  { 1, 1, 1, 1, 1, 1 },
  { 2, 1, 1, 1, 1, 1 },
  { 3, 1, 1, 1, 1, 1 },
  { 4, 1, 1, 1, 1, 1 },
  { 5, 1, 1, 1, 1, 1 },
  { 6, 1, 1, 1, 1, 1 }
}

fieldout = {
  iter = 1000
}

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

Run using the 750K mesh on 32 CPUs

./charmrun +p32 Main/inciter -i unitcube_750K.exo -c vortical_flow.q

Visualization

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

paraview out.e-s.0.32.0

The initial and final velocity, pressure, and total energy are depicted below.

Image Image Image Image

Numerical errors

To estimate the order of convergence, the numerical solution is computed using three different meshes, whose properties of are tabulated below.

MeshPointsTetrahedrah
750K132,651750,0000.02
6M1,030,3016,000,0000.01
48M8,120,60148,000,0000.005

Here h is the average edge length. The next figure shows the time evolution of the L1 error in fluid density for the three different meshes. The L1 error is computed as

\[ \left| \left| \varepsilon \right| \right|_1 = \frac{\sum_{v=1}^n V^v \left| \hat{U}^v - U^{v} \right|} {\sum_{v=1}^n V^v} \]

where n is the number of points, $\hat{U}^v$ and $U^v$ are the exact and computed solutions at mesh point v, and $V^v$ denotes the volume associated to mesh point v.

Image

Gnuplot commands to reproduce the above plot:

set terminal png size 800,600 enhanced font "Helvetica,16"
set output "vortical_flow_L1r_riecg.png"
set xlabel "time"
set ylabel "L1 error"
set grid
plot [] [:9e-5] "diag.750K" u 2:20, "diag.6M" u 2:20, "diag.48M" u 2:20

Convergence rate

The following table summarizes the L1 errors in all flow variables at the stationary state.

MeshL1(r)L1(u)L1(v)L1(w)L(e)
750K6.87e-56.83e-54.15e-57.92e-51.41e-3
6M1.91e-51.62e-59.67e-62.00e-53.75e-4
48M5.02e-63.93e-62.23e-64.97e-69.89e-5
p1.922.052.122.011.92

Here r, u, v, w, and e are the fluid density, velocity in the x, y, z directions, and the internal energy, respectively. The order of convergence is estimated as

\[ p = \frac{ \log \left| \left| \varepsilon \right| \right|_1^{m+1} - \log \left| \left| \varepsilon \right| \right|_1^m}{ \log h^{m+1} - \log h^m} \]

where m is a mesh index. As expected, the asymptotic convergence of the computed numerical solution in all variables approaches 2.

References

[1] J. Waltz and T.R. Canfield and N.R. Morgan and L.D. Risinger and J.G. Wohlbier, Manufactured solutions for the three-dimensional Euler equations with relevance to Inertial Confinement Fusion, Journal of Computational Physics, 267, 196–209, 2014.
[2] J. Waltz, T.R. Canfield, N.R. Morgan, L.D. Risinger, J.G. Wohlbier, Verification of a three-dimensional unstructured finite element method using analytic and manufactured solutions, Computers & Fluids, 81: 57-67, 2013.