Project Topics
Get Complete Project Material File(s) Now! »
Basilisk – State of the art: October 2015
When I started my PhD on 01/10/2015, the Basilisk tracer advection solver was based on the work of Bell, Colella and Glaz (BCG) (Bell et al. (1989)) and it still remains the work-horse solver for Basilisk. The BCG solver is a full Navier{Stokes solver, with formulations for the convective terms, diusive terms and distinct projection methods with time marching schemes. Given the context of this chapter I will be focussing on the convective term discretization with temporal scheme of the BCG solver in the next section, while the rest of the solver will be discussed in detail in chapter 4 on Navier{Stokes equations.
Basilisk is a powerful adaptive cartesian grid solver for PDEs, and the way adaptivity is implementedin Basilisk is through the use of wavelet function. This method will be discussed in detail in section 2.10 in this chapter. The error estimations and the general order of the numerical method, which was considered state of the art for Basilisk at the time of the start of my PhD was a O(2) method based on bilinear prolongation functions.
We intend to build higher-order methods for advection on Basilisk. For this we will be looking at WENO schemes for advection, implement the higher-order temporal schemes as introduced in section 2.3 and build a higher-order method to implement adaptivity.
Advection scheme by Bell, Colella and Glaz
Bell et al. (1989) introduced a second-order projection method for Navier Stokes equations in their 1988 paper, which improved upon the original projection method of Chorin (1967). The BCG algorithm, similar to Chorin’s method, rst solves the diusion-convection equation to predict an intermediate velocity eld, which is subsequently projected on a space of divergence-free vector elds. Where the two methods dier, is that in the BCG method the diusion-convection step and the projection step are coupled in a unique way to achieve a second-order temporal discretization which is not the case for Chorin’s original algorithm, which has a rst-order temporal accuracy.
Convection term
The convection term is mathematically expressed as (U r)U n+1=2 its formulation is based on the work of Colella (1990) and Vanleer (1983). The method involved here is dierent from the conventional upwind dierencing methods, in that, these unsplit second-order Godunov methods couple the spatial and the temporal discretization by making sure that the information is propagated only along the characteristics, leading to robust schemes with higher-order discretizations and improved phase error. There are four steps in the process, namely Reconstruction, characteristic extrapolation, Riemann problem and ux computation.
Reconstruction: This step involves the computation of limited slope proles in each cell. The gradient is computed using a slope limiter method as described by equations 2.14 & 2.15. The limiting algorithm is constructed to prevent the centered slopes from introducing new maxima and minima in the velocity eld. A new parameter, named determines the exact limiting scheme. If the value is equal to 1 we get the minmod limiter, which is the most dissipative scheme, while with a value of 2 we get the superbee limiter which is the least dissipative. d1 = (ui ui1).
Weighted Essentially Non Oscillatory (WENO) scheme
WENO schemes are based on ENO (Essentially Non Oscillatory) schemes, which were rst introduced by Harten et al. (1987). The key idea behind the ENO scheme is to use a higher-order ux reconstruction using a broad stencil in regions of smooth solutions while breaking down the higher-order stencil into multiple lower-order stencils, and then to use only the smoothest stencil among several candidate stencils to approximate the ux in regions of solution discontinuities. ENO schemes are uniformly higher-order accurate right up to the regions of solution discontinuities (shocks / hydraulic jumps) and are very robust to use. However there are some drawbacks as well. Firstly, wherever the solution or its derivative approaches a zero value, even round-o perturbations near that region can cause completely dierent choices of stencils. Secondly, ENO schemes require heavy usage of logical statements while choosing the appropriate stencil, which results in poor performance while implementation on vector supercomputers.
The WENO schemes of Liu et al. (1994) overcomes these drawbacks while maintaining the robustness of the ENO schemes. The methodology behind WENO diers from the methodology behind the ENO scheme in the way it approximates the numerical ux. While the ENO scheme uses only one of the candidate stencils in regions of solution discontinuities, the WENO scheme uses a convex combination of all the candidate stencil contributions. Each of the candidate stencil is assigned a weight based on a stencil smoothness indicator and this determines the eective contribution of the candidate stencil to the nal ux approximation. The weights are assigned in such a way that in regions of smooth solutions the eective reconstruction is a higher-order one, while in regions near the discontinuities the stencils containing the discontinuity are assigned a nearly-zero weight. WENO schemes remove the dependence on the logical stencil choosing statements making them run twice as fast compared to ENO codes on vector supercomputers. Also, the WENO schemes are also insensitive to round-o errors, unlike the ENO schemes. The method of Liu et al. (1994) was a third-order accurate nite volume method. Later Jiang and Shu (1996) derived the framework for implementing an arbitrary order WENO scheme, whose fth-order implementation remains the most popularly applied version of the WENO scheme. Next, based on the work of Jiang and Shu (1996), a step by step guide to implement the classical fth-order WENO scheme is presented (which we will be using for our application). For simplicity, we will start with a one-dimensional case and subsequently extend the application to higher-dimension cases. We look at the initial value problem given by eq. 2.1, the semi-discrete form of which can be expressed using eq. 2.24.
Implementing WENO-5 stencils on Basilisk
The complete stencil for a 1D – WENO-5 scheme includes data from seven points. As can be seen from the left and the right side formulations, the interpolation on the xi1=2 face requires volume averages from – fUi3;Ui2;Ui1;Ui;U1+1;U1+2g cells, similarly the interpolation on the xi+1=2 face requires volume averages from – fUi2;Ui1;Ui;Ui+1;U1+2;U1+3g cells. Hence, to discretize the convection term, we need access to seven volume averages or for three neighbors on each side for each cell. However, the way basilisk implements its stencils in the adaptive quadtree, a programmer gets access to only two neighbors on each side of a cell, or ve cells in total. (More details here – Basilisk-Stencils).
To have a work-around for this problem, we use an indirect approach of computing a O(2) cell-centered gradient eld, with the formulation given in eq. 2.30. After applying boundaryconditions on this gradient eld we now have access to a 7-point stencil, where the third neighbor on each side can be referenced by using the eq. 2.31.
Test case { Advection in a 1D Domain
In this section, we take two dierent one-dimensional test cases to observe the performance of the developed WENO schemes in Basilisk. Both cases involve the passive advection of a tracer eld under the in uence of a given velocity eld. The results and performance of the WENO methodolgy are compared and contrasted with the available second-order, Bell, Collella and Glaz scheme.
Passive advection of a 1D smooth tracer eld
Our rst test case involves the passive transport of a one dimensional sinusoidal eld given by the initial condition Tracer(x; 0) = sin(2x) over a periodic domain x 2 [0:5; 0:5], under the in uence of a uniform velocity eld given by ux = 0:1 . The time marching is carried out over one time period viz. t ! 0 to 10, and the advected solution is compared with the initial condition. The maximum errors are noted for a range of simulations carried out over various grid resolutions using both the BCG scheme and the new WENO scheme. The runs are repeated over two CFL numbers.
Test case – Passive advection in a 2D domain
In this section we will take a look at two dierent two-dimensional test cases to observe the performance of the WENO schemes in Basilisk. Both cases involve the passive advection of a tracer eld under the in uence of a given velocity eld. The results and performance of the WENO methodology are compared and contrasted with the available second-order Bell, Collella and Glaz scheme.
Periodic tracer in a uniform velocity eld
The 2D domain chosen is given by : (x; y) 2 [0:5 : 0:5;0:5 : 0:5], where the periodic tracer function has the formulation given in eq. 2.48. The velocities are initialized with constant values as shown in eq. 2.49. The time marching is carried out over one entire time period and the error is estimated by comparing the initial and the nal tracer data. We use the WENO non-limited scheme along with RK4 time marching. The simulations are run for two distinct CFL numbers. Tr(x; y; t = 0) = sin(2x)sin(2y) (2.48) u = 0:1^x + 0:1^.
Table of contents :
1 Introduction
1.1 Background
1.2 Computational uid dynamics
1.3 Chorin’s projection method
1.4 Error of a numerical scheme
1.5 Basilisk
1.6 Grids
1.6.1 Basilisk adaptivity
1.7 Fields
1.8 Boundary Conditions
1.9 State of the art for Basilisk at the start of my PhD
1.10 Outline of this PhD Thesis
2 Advection Solver
2.1 Context
2.2 Godunov’s scheme and the Riemann problem
2.3 Temporal Scheme – Runge{Kutta schemes
2.4 Basilisk – State of the art: October 2015
2.5 Advection scheme by Bell, Colella and Glaz
2.5.1 Convection term
2.6 Weighted Essentially Non Oscillatory (WENO) scheme
2.6.1 Left and right side reconstruction
2.6.2 Fifth-order WENO { Formulations
2.6.3 Implementing WENO-5 stencils on Basilisk
2.6.4 Literature review
2.7 Test case { Advection in a 1D Domain
2.7.1 Passive advection of a 1D smooth tracer eld
2.7.2 Passive advection of a 1D discontinuous tracer eld
2.8 WENO in 2D and 3D cases
2.8.1 Transverse sweeps – Gaussian quadratures
2.9 Test case – Passive advection in a 2D domain
2.9.1 Periodic tracer in a uniform velocity eld
2.9.2 Compact tracer in a solid body rotation
2.10 Multi-resolution analysis
2.10.1 Wavelet Transform – Lifting Algorithm
2.10.2 Wavelet transform and Basilisk adaptivity
2.10.3 Restriction operator
2.10.4 Prolongation operator
2.10.5 Fifth-order prolongation
2.10.6 Testing the order of the prolongation operator
2.11 Advection of a tracer under rotation and stretching
2.11.1 Uniform grid computations
2.11.2 Adaptive grid computation for the tracer advection problem
2.12 Conclusion
3 Poisson{Helmholtz Solver
3.1 Context
3.2 State of the Art: October 2015
3.3 Numerical Algorithm { Poisson Solver
3.3.1 Iterative Methods
3.3.2 Multigrid Methods
3.3.3 Discretization Scheme { Second-order solver
3.3.4 Discretization Scheme { Fourth-order solver
3.3.5 Higher dimension cases
3.3.6 Boundary Conditions
3.4 Results for the 9-point stencil
3.4.1 Uniform grid { Direct problem
3.4.2 Uniform grid { Inverse problem
3.4.3 Non-uniform grid { Direct problem
3.5 Convergence studies on adaptive grids
3.6 Conclusion
3.6.1 Applications of the Poisson{Helmholtz solver
4 Navier{Stokes Solver
4.1 Governing equations
4.2 Literature survey
4.3 Navier{Stokes solver by Bell, Colela and Glaz
4.3.1 Temporal discretization
4.3.2 Projection Algorithm
4.3.3 Viscous dissipation terms
4.4 Higher-order method for Navier{Stokes equations
4.4.1 Time-marching schemes
4.4.2 Convection term – WENO interpolation and Riemann Solver
4.4.3 Projection Algorithm
4.4.4 Viscous dissipation term
4.4.5 Test case: higher-order semi-implicit viscosity solver
4.5 Taylor{Green Vortex
4.6 Taylor{Green vortex with uniform background ow
4.7 Taylor{Green vortex with viscosity
4.8 Conclusion & Future scope
5 Explicit Saint-Venant Schemes
5.1 Background
5.2 Basilisk O(2) Saint-Venant solver
5.2.1 First-order well balanced method
5.2.2 Second-order well balanced method
5.2.3 Riemann Solver
5.2.4 Time marching scheme { predictor-corrector algorithm
5.3 Test case – Linear surface gravity wave
5.3.1 WENO-based explicit Saint-Venant solver
5.4 Conclusion
6 Conclusion & Perspectives
