In the past months, we have made significant advances in putting the NUMO prototype, which delivered results reported in this blog so far, back into the GNuMe framework. GNuMe stands for [G]alerkin [Nu]merical [M]odeling [E]nvironment and is a framework that hosts both NUMO and NUMA – Non-hydrostatic Unified Model of the Atmosphere. This merger allows NUMO to leverage a lot of  NUMA technology (both current and future, as the codes share the same infrastructure) and more rapid advancement in NUMO development. The benefits go both ways, as the development of NUMO brings in new features to NUMA as well. Now that we are happy with how the code infrastructure looks like, we can go back to evolving NUMO into a truly unstructured non-hydrostatic model.

Ever since we have obtained the first lock-exchange results, two issued bothered me. Firstly, we experienced some unphysical instabilities and mesh imprinting close to domain boundaries, particularly in the domain corners. Those instabilities eventually led to simulation crash at later times, way past the usual end time of the simulation when the fronts reach the wall. Regardless, we definitely cannot allow this in a production code. Secondly, we compared our results to those available in the literature, but there was quite a bit of spread between reported Froude (Fr) numbers among different researchers. Recently I have revisited those simulations, and I believe I found answers to both issues.

Since we were getting mostly correct physical behavior of the flow, I suspected  that the instabilities are simply due to not enough stabilization of the temperature field. Following Hiester et al. (2011), we used kinematic viscosity ν=1e-6 m^2/s for velocity, and no thermal diffusion. Given no numerical diffusion in the high-order continuous Galerkin method used for those simulations, we always require some stabilization, either through artificial diffusion, filtering, or other methods. I was using Boyd-Vandeven filters, but this was clearly not enough. I decided to introduce some physical diffusion to see if this will stabilize the simulation. Since I was locked in to kinematic viscosity of water at 20C, I decided to use matching thermal diffusion of κ=0.148e-6 m^2/s. The result was a very stable lock-exchange simulation, which I run till 500s with no trace of growing instabilities. Here is the video of the flow:

This led to a further investigation of the lock-exchange results. As mentioned before, Hiester et al. (2011) report the Grashoff number is Gr=1.25e6, with kinematic viscosity ν=1e-6 m^2/s and no thermal diffusivity (κ=0 m^2/s), which leads to Schmidt number Sc =  ν/κ = ∞. A similar set-up is found in the simulation by Fringer et al. (2006), but the reported Froude number is significantly different than Hiester et al. (see the table here). The authors, however, compare the result with the DNS simulation by Härtel et al. (2000), which provides yet another Fr value, and note the possible cause of this discrepancy. The simulation of Härtel et al. was made (apart from significantly higher resolution) with Sc=0.71, which with ν=1e-6 m^2/s gives the thermal diffusivity of κ=1.41e-6 m^2/s. Fringer et al. point out that even though they do not explicitly use thermal diffusivity, the numerical diffusion of the method (first-order upwinding) makes the Sc number lower than this of Härtel et al., and therefore their fronts are slower.

This led me to try two simulations with different values of Sc number, one with Sc=0.71 like in Härtel et al., and the other with Sc=6.74, which is the simulation shown in the video above. I have compiled the results in a table below, taking into account Sc numbers (all simulations with Gr=1.25e6).

It turns out that NUMO results match very well the DNS simulation by Härtel et al. when using matching Schmidt numbers, and the simulation of Hiester et al. when using high Sc. Even though Hiester is not using any thermal diffusivity explicitly, it is hard to evaluate the numerical diffusion of their model. The front speed in Fringer et al. is clearly slower than in other simulations, indicating much lower Schmidt number due to high numerical diffusion of the model.

This result gives me another bit of confidence in NUMO results, as it agrees with the DNS simulation very well, while also matching the result of high resolution adaptive finite element code of Hiester et al. (2011) by using only physical diffusivity.

References

Fringer, O.B., Gerritsen, M. & Street, R.L., 2006. An unstructured-grid, finite-volume, nonhydrostatic, parallel coastal ocean simulator. Ocean Modelling, 14(3–4), pp.139–173.

Härtel, C., Meiburg, E. & Necker, F., 2000. Analysis and direct numerical simulation of the flow at a gravity-current head. Part 1. Flow topology and front speed for slip and no-slip boundaries. Journal of Fluid Mechanics, 418, pp.189–212.

Hiester, H.R., Piggott, M.D. & Allison, P.A., 2011. The impact of mesh adaptivity on the gravity current front speed in a two-dimensional lock-exchange. Ocean Modelling, 38(1–2), pp.1–21. Available at: http://dx.doi.org/10.1016/j.ocemod.2011.01.003.

Once we have verified the correctness of the implementation of CG/DG machinery, it is time to move on to more physical-like benchmarks. This time we are not able to rely on exact analytic solutions as we benchmark NUMO on a lock-exchange case, which is a standard test for non-hydrostatic ocean codes. An additional benefit of this set-up is that it can be studied experimentally. But what is the lock-exchange case exactly?

In a box of 0.8m x 0.1m x 0.1m we put two fluids side by side, one is warm and light and the other one cold and heavy. The temperature difference between them is set to 1C. To compare this set-up with other results from the literature, we use a Grasshoff number Gr=1.25e+06, which, with water viscosity of 1e-06, implies the thermal conductivity of 1e-03. We assume that there is no thermal diffusion present. We let the two fluids interact, with the difference in buoyancy of the fluids driving the flow. The movie below shows the simulation result. Note that even though we show only a 2D result, this is a full 3D simulation and we are showing only a middle plane of the solution.

The domain is bounded by a no-slip wall on the bottom, and all the other walls have free-slip condition prescribed. To compare this result with the literature, we track the velocity of the no-slip and free-slip fronts and report it as non-dimensional Froude number.

The figure shows the Froude number plotted against the distance that each front traveled. We compare this result with the simulation by Hiester et al. and see that both results match very closely, even though the simulation of Hiester is much better resolved (dx=0.00025m, while we use an effective resolution of 0.001m, but also a higher order method). Close to the initial position both fronts go through a phase of acceleration, and at x=0.4m they hit the side-walls. There is, however, a period of relative steady velocity between x=0.2m and x=0.3m, so we average the Froude number there and compare with other literature results.

As expected, the NUMO result matches closely recent Hiester et al. publication. However, NUMO is also very close to Hartel et al. very-high resolution DNS. This result builds our confidence in the correctness of the buoyancy effects in our code, as well as the performance of temperature (and salinity, but not featured here) transport.

So what is next? This result is an important stepping stone to a full-fledged non-hydrostatic ocean model. We are ready to move towards more realistic, non-box geometries and idealized ice/ocean interaction tests in ISOMIP suite. To be able to call ourselves a real ocean model we still need the free-surface implementation, which is currently in the works.