Computation of the van Keken Problem with the Volume-of-Fluid Interface Tracking Method#

This section is a co-production of Jonathan Robey and E. G. Puckett. One can also model the van Keken problem with the volume-of-fluid (VOF) interface tracking algorithm in . In fact, this problem is particularly well-suited to being computed with the VOF method, since it consists of two distinct, immiscible fluids and interface tracking algorithms are specifically designed not to allow the two fluids to mix. In particular, assuming the computation is sufficiently well-resolved, the fluids will not mix at sub-grid scales over the entire duration of the computation. However, note that this implies that all computations of the van Keken problem made with the VOF method must necessarily be with discontinuous initial conditions. Finally, since one is often interested in a high resolution image of the shape of the interface between the two fluids, one can use the output of the VOF method to examine the interface within individual cells or within regions consisting of groups of cells.

Another advantage the VOF method has over modeling moving interface problems with compositional fields is that for problems in which the interface occupies a relatively small part of the computational domain, all of the computational work in the VOF method is done in cells that lie in a neighborhood of the interface, rather than in the entire computational domain.


Fig. 130 Evolution of the root mean square velocity as a function of time for computations of the van Keken problem made with the VOF interface tracking algorithm with five different global mesh refinements. Since the VOF initial conditions are discontinuous, the above results should really be compared with the computations with discontinuous initial conditions on the left in Figure [fig:vk-2]. However, the above results also compare extremely favorably with the computations with smoothed, continuous initial conditions for the compositional field on the right in Figure [fig:vk-2]. As in Figure [fig:vk-2], 5 global refinements correspond to a 32 \times 32 mesh and 10 global refinements correspond to a 1024 \times 1024 mesh.#


Fig. 131 The results of two computations of the van Keken problem made with the VOF interface tracking algorithm overlaid upon each other at t_\text{end}=2000. This visualization shows the reconstructed boundary between the two materials at the final time t_\text{end} as computed on a uniform grid with 7 and 8 levels of refinement. The boundaries between the materials are displayed as contours of the fields \tilde{\psi}^7(t_\text{end}) (black) and \tilde{\psi}^8,(t_\text{end}) (bright green), which are generated by the visualization postprocessor. The contours for the reconstructed material boundaries are superimposed on a color gradient visualization of the material composition for the computation with 8 levels of refinement in order to make the regions with each fluid type more evident. Compare with the fourth image on the right in Figure [fig:vk-1].#

As noted above, when the interface is discontinuous, the van Keken problem is a version of the Rayleigh-Taylor problem, which is unstable to perturbations of all wavelengths1 (e.g. see (Chandrasekhar 1961)). Therefore, it is extremely sensitive to the initial conditions. In order to address this sensitivity, we do not use the default approach of computing the initial material volume fractions using a composition quadrature. Instead we compute the initial volume fractions using a signed distance function \(\phi\) as follows (Robey 2019; Robey and Puckett 2019).

First we create the function \(\phi\), which has the following two properties:

  1. it is positive in the region that contains one of the fluids, which we will refer to as fluid 1, and negative in the complement of this region, which we will refer to as fluid 2, and 2) at each point in the domain the magnitude \(| \phi |\) of \(\phi\) is the distance to the boundary between the two fluids or materials. In the computations shown here, we use an approximation \(\tilde{\phi}\) to \(\phi\) such that the difference between \(\tilde{\phi}\) and \(\phi\) is small enough for the purposes of making the computations high-quality. The use of an approximation as opposed to the function itself is due to obtaining an appropriate function for almost all nontrivial (not a line or a circle) boundaries is extremely difficult. In this case, because the boundary is close to horizontal, we use the vertical distance to the boundary based on the argument that the gradient will not differ sufficiently from \(1\) to induce errors in the initialization computation. The primary advantage of choosing this particular initialization algorithm is that it allows us to more accurately reproduce the initial condition on a sub-grid scale than would otherwise be possible on the coarser grid on which we compute the time evolution of the interface.

The relevant sections of the parameter file for this type of initialization of the VOF method appears immediately above. In particular, the combination of Number of initialization samples with the level set initialization type indicates that our initialization will consist of dividing each grid cell into \(16 \times 16\) subcells and the distance to the given initial interface \(f(\mathbf x)\), provided in Function expression, is computed in each of the 256 subcells. We then use this information to compute a piecewise linear interface approximation to \(f(\mathbf x)\). The volume fraction in each subcell is then found in the manner described in (Robey 2019; Robey and Puckett 2019). This initialization procedure provides a much finer and thus, more accurate, initial condition than the standard VOF initialization procedure described above.

While the visualization configuration in a typical parameter file is sufficient for most purposes, when using the VOF method one has the ability to see the division between the fluids reconstructed by the VOF algorithm in each cell. This is accomplished by plotting the zero contour of a field \(\tilde\psi\) that is generated to be \(0\) on the reconstructed interface, positive in the region with fluid 1, and negative in the region with fluid 2. However \(\tilde{\psi}\) does not satisfy the requirement that the magnitude is equal to the distance to the interface as would be required for the signed distance function \(\phi\). The modifications to the parameter file that are necessary in order to draw the reconstructed boundary as a contour are shown immediately below. The full configuration file for this version of the benchmark problem can be found at cookbooks/van-keken-vof/van-keken-vof.prm.

We made a number of computations of the van Keken problem with the VOF method in order to compare the wall clock times with computations using a DG compositional field. We ran both on the same cluster at global refinements 5–8 using one node with four CPUs and refinements 9 and 10 using two nodes with 16 CPUs. Our results are shown in Table 1. In all of the computations shown in Table 1 we used a CFL number of \(\sigma=0.5\). Due to the change in the CFL number from \(\sigma = 1.0\) in Table [tab:runtime-table] to \(\sigma = 0.5\) in Table 1 and the difference between HPC clusters on which the computational results shown in the two tables were made, we can’t make a direct quantitative comparison between the data in Tables [tab:runtime-table] and 1.

However, we can compare the required run time for a VOF computation to that for a DG computation. We note that the use of the VOF advection algorithm significantly reduces the required computation time in all cases, frequently requiring less than half the time required by the DG compositional field.

We now examine the RMS velocity data shown in Figure 1. Other than for the case of 5 levels of uniform global refinement, the curves for the RMS velocities for \(6\), \(7\), \(8\), \(9\) and \(10\) levels of refinement in Figure 1 are nearly indistinguishable.

Upon examining the solution at the final time, we note that the general structure of the solution shown in Figure 2 matches the form and the general structure found in other versions of this benchmark such as the fourth image on the right in Figure [fig:vk-1]. We also note that the differences in the shape of the interface based on a single refinement as shown in Figure 2 are minor, although still slightly visible. This is to be expected as refinement is a perturbation of the initial condition at a smaller wave length.

```{table} Comparison of runtimes for the van Keken problem with VOF and a DG compositional field, in which the initial conditions for DG smoothed are as described in section [paragraph:van-keken compositional fields] above. The times shown are for the full computation, ending at \(t_\text{end} = 2000\) with a CFL number of \(\sigma=0.5\) in both cases. All of these computations were made with ASPECT version 2.2.0-pre (master, commit ef542ecc2) in release mode on the Peloton2 cluster at U.C. Davis. We note that the change in the CFL number \(\sigma\) and the differing choice of cluster makes a direct quantitative comparison between this table and Table [tab:runtime-table] invalid due to too many confounding factors.



Global Refinement

Number of Processors





1.33 minutes

2.57 minutes



8.51 minutes

19.5 minutes



1.15 hours

2.49 hours



8.53 hours

19.6 hours



16.30 hours

2.72 days



5.17 days

>6.00 days

The consistency of the results shown here differs noticeably from the behavior
of the problem with discontinuous initial conditions when computed with the
FEM and DG advection algorithms. One possible reason for these differences is
the specialized initialization procedure used for the volume of fluid method,
which permits a much more consistent initialization by reducing the variation
in the initial condition when the initial mesh is refined.

To study this feature of our algorithm and the sensitivity of the problem to
the precise initial condition, we vary the size of the initial interface
perturbation and examine the sensitivity of the final results to a small
change in the initial conditions. Specifically, we vary the amplitude $a$ of
the cosine function in the initial conditions, as shown below.

``` prmfile

Fig. 132 Computations of the van Keken problem made with the VOF interface tracking algorithm showing the evolution of the RMS velocity as a function of time for small changes in the amplitude a of the cosine function in the initial condition at 7 levels of refinement. Compare to Figures [fig:vk-6] and 1.#

In these computations we vary the value of \(a\) from its usual value of \(a = 0.02\) to \(5\% = 0.001\) below its usual value to \(5\%\) above its usual value in increments of \(0.01\). In other words, we compare the values for \(a =0.019\), \(0.020\), and \(0.021\). Upon examination of Figure 3, we see a visible variation in the location of the second peak, although the overall shape of the curve remains consistent with the curves in Figure 1. The size of this variation in the initial conditions cannot be expected to be reproduced using the standard compositional quadrature initialization procedure for VOF unless the cell size is on the scale of the change in the value of \(a\); i.e., \(h \, \lessapprox \, \Delta a = 0.001\). We also note that the smoothing parameter which would produce a \(10^{-3}\leq C \leq 1 - 10^{-3}\) band on the order of the same size as the amplitude variation shown here, would be approximately \(2.8 \cdot 10^{-4}\). This perturbation is much smaller than any of the changes in width of the smoothed regions in the computations shown in Figure [fig:vk-6]. In summary, these results demonstrate the sensitivity of the discontinuous version of the van Keken problem to even extremely small variations in the initial conditions.

Chandrasekhar, S. 1961. Hydrodynamic and Hydromagnetic Stability. New York: Dover.

Robey, Jonathan M. 2019. “On the Design, Implementation, and Use of a Volume-of-Fluid Interface Tracking Algorithm for Modeling Convection and Other Processes in the Earth’s Mantle.” Applied Mathematics, University of California, Davis.

Robey, Jonathan M., and Elbridge Gerry Puckett. 2019. “Implementation of a Volume-of-Fluid Method in a Finite Element Code with Applications to Thermochemical Convection in a Density Stratified Fluid in the Earth’s Mantle.” Computers & Fluids 190 (2): 217–53.

1 This is true whether the two fluids have the same viscosity or different viscosities.

            compositional fields