13. The good, the bad and the ugly
Here I share some thoughts on our solver design.
I consider minimalism to be a good thing when doing a coding project and since I did my best to achieve this feature, I do believe our solver is a minimalist design for the functionality put on display in the gallery. It nicely served to highlight the concepts of the steps typically also tackled by more useful solvers. Further, I think it helps to better appreciate alternative methods. The user interface is also relatively good as I do think the codes of the Gallery are quire intuitive.
Another good thing is that we have written our solver in the C programming language. This is a proper language for high-performance computing and may be the only thing that you could take from this solver if you are to code a better one.
Our solver is time inefficient. The discrete approximations are not very accurate and therefore require relatively small discretization elements (time steps and mesh-element sizes) to achieve an acceptable solution. To make things worse, our implementation of the Poisson-problem solver is highly inefficient and the number of required sweeps for a given spatial problem actually increases with the resolution, making the solver not well behaved. For this reason, virtually none of the choices we have made in this regard are used in more useful solver designs.
We have already encountered the problem where we cannot run our code with a large number of grid points. Memory allocation for the solution fields on the stack is obviously not the way to go.
In an earlier chapter, I wrote about the modular design of this solver for the purpose of future improvement. However, if almost every thing needs improving, a complete rewrite makes more sense.
What do other solvers do different?
As mentioned in the introduction, there is no consensus on how to solve flow problems most efficiently and numerical methods are often developed with a specific type of flow problems in mind. Aside the from coding details, I will list some solver-design choice alternatives which could have some relevant depending on the application.
Functions and features
Solvers typically advertise their merits by listing the type of flow problems they can handle. Examples include; multi-phase flows, flow in complex geometries, coupling to other solvers (multi physics), etc.
Method of Lines
The method of lines, where the mesh only discretizes the spatial dimensions, is often used in any solver that is not specifically investigating alternatives.
The finite difference method is not particularly unpopular. However, the finite volume method (FVM) is often preferred for solving problems described by a flux-divergence evolution equation. the FVM is formulated such that the approximate solution to a conservation equation inherits this conservation property exacty. A second prominent alternative is the finite element method. Apart from the mathematical rigor and the grid-node layout flexibility, I cant think of many fundamental advantages for it in this context.
Virtually all finite difference codes use a staggered grid, where the velocity components are defined at the faces in between grid cells. In practice, this allows to easily do an exact projection.
Our Cartesian grid is excellent for meshing the square domain. However, it would be nice if we could focus the cells a bit towards the regions of interest. This could be achieved with stretching and squeezing the mesh, but if you are doing many computations in complex geometries where only certain parts of the domain require a high resolution, an unstructured mesh would be more flexible. Noting that a quadtree structure may represent the best of both worlds.
Forward Euler time advancement
The forward Euler method is prone to instabilities and requires many small time steps for an accurate solution. It can thus become a computationally expensive option. Popular time-advancement schemes (both ex and implicit) can typically be casted in a Runge-Kutta formulation. Special attention is often given to the formulations that require least storage in memory. For some reason, linear multistep methods are also still around.
First-order accurate upwind advection
Although our first-order method has some excellent theoretical properties, the first order accuracy introduces large (diffusive) errors. There really exists a full zoo of advection-scheme alternatives. I like to highlight so-called compact numerical schemes as they can have excellent dispersion properties for a given degree of accuracy.
2nd-order accurate diffusion
This was quite an OK choice, although its explicit stability is a huge issue for viscous problems. For this reason, there are quite a few solvers that employ a so-called implicit + explicit (IMEX) scheme, where only the viscous term is treated implicitly.
Iterative Poisson solver
Iterative Poisson solving is often synonymous for also using multigrid acceleration because it does not really work well without it (as we have seen). Further, the sweeping design can be altered in many ways. Most notably, under relaxation may help to improve the convergence rate. Apart from the iterative method, on periodic grids, a (fast) Fourier transform would be a better option. This is often referred to as a spectral method, and we could actually have done all our differencing approximations in spectral space. Finally, (direct and indirect) matrix solver are also a popular choice to tackle linear problems.
Chorin’s projection method
Chorin’s projection method is quite popular but other operator splitting methods exist.