CFD Solver Benchmark and Comparison - OpenFOAM, FEniCS and FEATool

Multi-Simulation and External Solvers in FEATool

With the external solver and multi-simulation interfaces, FEATool Multiphysics makes it easy to use multiple solvers which can be suitable for different classes of problems. Using more than one solver can also be very helpful in validating physics models and simulations, as one can compare results directly for several solution methods.

With FEATool, direct knowledge about how the external solvers actually work is not required, as model and mesh conversion and export, solver control, solution interpolation and re-import are all handled conveniently and automatically by the external solver interface. This saves a significant amount of time as models only need to be set up and defined once with the easy FEATool GUI or MATLAB CLI interface (instead of once for each solver).

In the following, the built-in FEATool multiphysics solver, dedicated OpenFOAM computational fluid dynamics (CFD) solver, and general FEniCS finite element analysis (FEA) solver will be compared for the same benchmark test case.

Flow Around Cylinder Benchmark Test Case

A classic laminar and incompressible flow around cylinder CFD test case is used to evaluate the efficiency of the solution methods [1]. Very accurate reference solutions for drag, lift, and pressure difference have been established for this benchmark [2]. The test case is described in more detail in the flow around a cylinder tutorial example which can found in the FEATool User’s Guide. A video tutorial is also available showing step by step how to set up the model in the FEATool GUI

Although the benchmark simulations can be performed just as well in the GUI, it is usually more convenient to write a MATLAB m-script file, which then can be run automatically (and repeated for reproducibility). FEATool m-script model files can either be written completely from scratch, or conveniently exported from the GUI using the Save As M-Script Model… option in the File menu (this will translate each GUI action to a corresponding MATLAB function call and export it to an editable plain text file).

The computational grid employed in this study is the same structured quadrilateral mesh used in reference [2]. The cylinder flow benchmark grid is constructed using the structured grid primitives and joined together. Up to four successive uniform refinements (five grid levels) are used to observe the convergence behavior with respect to the mesh size. The table below shows the grid statistics as well as the total number of unknowns (degrees of freedom) for the solvers.

FEATool Multiphysics Solver

FEATool features dedicated solvers for both stationary and time dependent coupled multiphysics problems. PDEs are are discretized with the finite element method (FEM) as a single fully coupled system. This approach means that all equations and unknowns are solved together in one step. In contrast to traditional splitting and segregated solution approaches where equations are decoupled and solved sequentially, a fully coupled (monolithic) solver typically features better convergence and robustness, at the cost of system size and increased memory requirements. The ensuing sparse linear systems are in turn solved with the direct solver built into MATLAB, which currently defaults to SuiteSparse’s UMFPACK.

In the present benchmark the Navier-Stokes equations are discretized with the Q₂P_-1 finite element shape function space for quadrilateral cells (second order velocity and first order discontinuous pressure approximation). To call the FEATool stationary solver the solvestat function is used

fea.sol.u = solvestat( fea, 'jac', jac );

where fea is a FEATool finite element data struct, and jac a user defined analytic Newton Jacobian (which is faster than using a numerical Jacobian). Fixed-point Picard iterations will be used for non-linear problems if the jac argument is omitted. The resulting solution vector is stored in the fea.sol.u field.

The resulting drag c_D, and lift c_L coefficients, and pressure difference between front and rear of the cylinder, Δp, can be seen in the table below. The required CPU time in seconds is denoted by t_solve and the number of non-linear iterations it. As can be seen the solutions converge nicely towards the reference values as the mesh density increases.

OpenFOAM CFD Solver

OpenFOAM is the leading leading free, open source software for computational fluid dynamics (CFD) [3]. Featuring dedicated solvers and support for many types of flow regimes such as incompressible and compressible, turbulent, non-isothermal, and multiphase flows, OpenFOAM is a very versatile flow solver package. To discretize the equations OpenFOAM uses a cell centered finite volume method (FVM), with operator splitting and iterative linear solvers. For steady state calculations, such as featured in this benchmark test case, the Simple scheme is employed with pseudo time stepping.

OpenFOAM is a quite complicated and complex software code, designed for the command line, and usually takes significant time and effort to learn how to use properly. To make this easier, FEATool includes dedicated functionality automate solving fluid mechanics problems with OpenFOAM. Instead of using the solvestat function to call the FEATool solver as above, one can instead use the openfoam function call

fea = openfoam( fea );

This automatically converts incompressible Navier-Stokes flow problems to OpenFOAM case files (including extruding 2D meshes to 3D), solves them with OpenFOAM, and imports the resulting solutions back into FEATool. (Note that the actual OpenFOAM and external solver binaries have to be installed separately, as they are not included with the FEATool distribution.)

The resulting benchmark quantities are presented in the table below (here it designates the number of pseudo time steps required to reach steady state). The results produced by OpenFOAM are not as accurate as for the FEATool solver. This can possibly be attributed to the cell centered constant discretization which can not yield the same accuracy as a higher order one.

FEniCS FEA Solver

FEniCS is designed to quickly translate scientific models into efficient finite element code and be a computing platform for solving partial differential equations (PDEs). The FEniCS Project finite element analysis (FEA) solver is, similar to FEATool Multiphysics, designed to be able to solve general systems of coupled equations.

FEniCS features its own PDE syntax and language to which FEATool models can be directly converted and exported. This can be accomplished with help of the fenics interface FEATool function as

fenics( fea, 'mode', 'export' );
fenics( fea, 'mode', 'solve'  );
fea = fenics( fea, 'mode', 'import' );

where the mode argument specifies which solver phase or action to take. (This will in a future release be simplified to a single call as for the other solvers.)

In contrast to FEATool, FEniCS currently only supports simplex mesh cell shapes (triangles in 2D and tetrahedra in 3D). The structured benchmark grid is therefore subdivided into triangles for the FEniCS simulations. As for FEM shape functions FEniCS here employs the P₂P₁ Taylor-Hood mixed finite element space for the Navier-Stokes equations. As for sparse linear solver the default (direct PETSc) solver is used.

From the results shown in the table below it can be seen that convergence results similar to those for FEATool is achieved. This is also expected as FEniCS per default also employs a fully coupled monolithic approach (although operator splitting is supported but has to be implemented manually). As the FEniCS results are interpolated to the grid points before re-import into MATLAB and FEATool, some accuracy in the velocity space is also expected to have been lost. (In critical applications FEniCS post-processing might better be performed in FEniCS/Python as it is not quite clear how to export higher order unknowns, or how the degrees of freedom are internally reordered.)

Solver Comparison

The resulting CPU timings as a function of accuracy (cost vs. accuracy) are shown in the figures below. With respect to drag and lift coefficients the FEATool Multiphysics solver performs significantly better than both OpenFOAM and FEniCS (although technically FEniCS should eventually be able to achieve the same accuracy). For the pressure drop measure, FEniCS and FEATool was roughly equal while the OpenFOAM results were about two orders less accurate for the same computational cost. As for convergence FEATool featured better convergence for the drag coefficient, while lift and pressure was linear as for FEniCS. OpenFOAM seemed in particular to have issues with the finest grid, as well as the lift coefficient (which is the most difficult quantity to compute accurately).

To summarize one can say that for this steady laminar flow test case (Re = 20) fully coupled monolithic solvers, such as used in FEATool and FEniCS, seem to be the best choice. Here OpenFOAM computed overall less accurate and uncertain results for a given CPU effort. In particular OpenFOAM showed weak scaling with respect to mesh density (increasing number of iterations to converge). The benefit of OpenFOAM it that due to the splitting approach and iterative solvers, less memory is required and larger mesh sizes can be used. In the future it might be interesting to extend the study to instationary test cases which might more be in favor of OpenFOAM and operator splitting (which can be implemented as custom solver approaches in both FEATool and FEniCS).

CFD Benchmark Download

The CFD benchmarking script used here has been written as a convenient and self contained FEATool MATLAB m-script file and can be downloaded from the link below

References

[1] G. Nabh, On higher order methods for the stationary incompressible Navier-Stokes equations, PhD Thesis, 1998, University of Heidelberg. Preprint 42 / 98, 1998.

[2] V. John, G. Matthies, Higher-order finite element discretizations in a benchmark problem for incompressible flows, International Journal for Numerical Methods in Fluids, 2001, 37:8:885-903.

[3] The OpenFOAM homepage, The OpenFOAM Foundation, London, United Kingdom, 2018.