V&V » RieCG: Rayleigh-Taylor unstable configuration

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 test case is to assess time dependent fluid motion in the presence of Rayleigh-Taylor unstable conditions, i.e., opposing density and pressure gradients. The derivation of this test problem is given in [1].

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(x_i) & = \rho_0 - \left( \beta_1 x_1^2 + \beta_2 x_2^2 + \beta_3 x_3^2 \right) \\ u_i(x_i,t) & = f(t) g(x_i) = \cos(k \pi t) \left( \begin{array}{ccc} x_3 \sin (\pi x_1) \\ x_3 \cos (\pi x_2) \\ -\frac{1}{2} \pi x_3^2 \left[ \cos (\pi x_1) - \sin(\pi x_2) \right] \end{array} \right) \\ p(x_i) & = p_0 + \alpha \left( \beta_1 x_1^2 + \beta_2 x_2^2 + \beta_3 x_3^2 \right) \\ E(x_i,t) & = \frac{p(x_i)} {\rho \left( \gamma -1 \right) }+ \frac{1}{2} f^2(t) g_j g_j \end{split} \]

and the source terms

\[ \begin{split} S_\rho & = u_i \frac{\partial \rho}{\partial x_i} \\ S_{u,i} & = \rho g_i \frac{\partial f}{\partial t} + f g_i S_{\rho} + \rho f^2 g_j \frac{\partial g_i}{\partial x_j} + \frac{\partial p}{\partial x_i} \\ S_E & = \rho g_i g_i f \frac{\partial f}{\partial t} + \left[ \frac{p}{\rho \left( \gamma -1 \right) } + \frac{1}{2} f^2 g_i g_i \right] S_{\rho} \\ & \quad + \rho f g_i \left[ \frac{1}{\rho \left( \gamma -1 \right)} \frac{\partial p}{\partial x_i} - \frac{p}{\rho^2 \left(\gamma -1 \right)} \frac{\partial \rho}{\partial x_i} + f^2 g_j \frac{\partial g_j}{\partial x_i} \right] + f g_i \frac{\partial p}{\partial x_i} \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 is time-dependent. The simulation was run for a single time unit computing the numerical errors during time stepping and the errors at t=1 are used to assess the convergence rate. A time step of 0.001 was used for the coarsest mesh and was successively halved for the finer meshes.

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

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

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 non-stationary Rayleigh-Taylor"

term = 3.0
ttyi = 10

solver = "riecg"

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

part = "rcb"

problem = {
  name = "rayleigh_taylor",
  alpha = 1.0,
  beta = { 1.0, 1.0, 1.0 },
  p0 = 1.0,
  r0 = 1.0,
  kappa = 1.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 rayleigh_taylor.q

Visualization

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

paraview out.e-s.0.32.0

The figures below depict the velocity field at the surface at different simulation times, showing the reversal of the velocity components in time, indicating the spatial and temporal dynamics of the velocity field.

Image Image Image Image Image Image Image Image Image

The figures below show the density, pressure, energy, and velocity on the surface at the final simulation time.

Image Image Image Image

Numerical errors and convergence rate

The table below shows the numerical L1 errors of the flow quantities computed. 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 total number of points of the mesh, $\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. 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. The following table summarizes the L1 errors in all flow variables at t=1.

MeshL1(r)L1(u)L1(v)L1(w)L1(e)
750K4.72e-38.50e-48.84e-41.05e-32.31e-2
6M1.35e-32.85e-42.80e-43.53e-46.62e-3
48M3.73e-48.72e-58.24e-51.03e-41.82e-3
p1.861.711.761.781.86

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 approaches the expected value of 2.0.

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.