FEATool
v1.7.1 Finite Element Analysis Toolbox |

Solvers

This chapter describes the solver options available with FEATool. Solvers are used to compute solutions to the involved finite element FEM discretizations.

After a geometry, grid, and equation and boundary conditions have been defined one can change to solve mode by pressing the mode button.

The **Solve** menu similarly allows for switching to **Solve Mode** and also opening the **Solver Settings...** dialog box. The menu also features an option to **Get Initial Solution** which computes the expressions defined in the *Initial Condition* fields of the equation settings dialog box and returns this as the solution. Moreover, the default solver call command can be edited and changed through the **Solver Call Command...** option thus allowing for custom and external solvers to be employed.

In solve mode the available toolbar buttons and actions are the following.

- Calls the default solver as defined by the
*solver settings*with initial conditions prescribed in the equation settings. After the problem has been solved FEATool switches to*postprocessing*mode and plots the solution according to the chosen plot types. - Re-solves the problem using the last solution as initial value. Note that for this to work the geometry, grid, physics modes, and finite element shape functions must all remain the same.
- Solves fluid flow problems with the FeatFlow external CFD solver.
- Solves general finite element problems with the FEniCS external FEM solver.
- Opens the
*Solver Settings*dialog box.

The **Solver Settings** dialog box shown below and can be opened by pressing the corresponding toolbar button. The solver settings allow advanced control over the default built-in FEATool solver.

Firstly, the **Solver Type** option allows for either stationary monolithic solver solvestat, or the time dependent instationary solver solvetime to be chosen. Monolithic means that all dependent variables are discretized and solved together in a large coupled system, rather than iteratively solving a segregated set of smaller decoupled systems.

**Initial Condition** controls whether the initial values should be computed from the equation expressions in the initial value fields, or taken from an old solution in which case the solution can be chosen from the corresponding list box. Note that for this to work the geometry, grid, physics modes, and finite element shape functions must all remain the same.

The **Time Dependent Settings** controls time-stepping parameters when using the time-dependent solver. The **Time step** selects a fixed time step size for single step schemes, or the average time step size for multi-step schemes. Time dependent simulations will terminate either when the accumulated time has reached the specified **Simulation time**, or the normed changes in the dependent solution variables are smaller than the **Time stopping criteria** indicating a stationary state. The **Time stepping scheme** allows choosing between the single-step first order Backward Euler and second order Crank-Nicolson schemes, as well as the second order (multi-step) Fractional-step-theta scheme. Lastly, one can choose between using standard row-sum lumped and full **Mass matrix** representations.

The non-linear solver is used then dependent variables are detected in the equation coefficients, which is the typical for multiphysics problems. Per default standard fixed point iterations are used while Newton iterations can be used if an analytical Jacobian is supplied by the user (see the Navier-Stokes fluid flow models for examples of this).The **Non-Linear Solver Settings** allows control over the **Maximum number of non-linear iterations**, and the **Non-linear relaxation parameter** which is the fraction of the new to old solution to weigh in the next non-linear iteration. The non-linear solver will terminate and stop when either of the **Defect stopping criteria** or **Solution changes stopping criteria** are met which can be chosen to be measured both absolute or relatively. Stopping criteria for the non-linear solver

In the **General Settings** frame one there is a setting to control the numerical integration rule and order used in the matrix assembly operations (it is necessary to increase this for higher order >2 shape functions).

FEniCS is a flexible and comprehensive finite element FEM and partial differential equation PDE modeling and simulation toolkit with Python and C++ interfaces along with many integrated solvers. As both FEATool and FEniCS discretize equations employing a weak finite element formulation it is quite straightforward to translate FEATool syntax and convert it to FEniCS python scripts. Similar to what has been done with the FeatFlow CFD solver the FEATool-FEniCS integration allows for easy conversion, exporting, solving, and importing FEATool Multiphysics models to FEniCS directly from the GUI, as well as the Matlab and Octave command line interfaces.

FEniCS is aimed at supporting and solving general systems of PDEs, such as found in engineering and coupled multiphysics problems. As both FEATool and FEniCS discretize the equations in the same way the solutions they produce should be virtually identical. One of the many advantages of using FEniCS then, in addition to now supporting the Python scripting language as well as Matlab and Octave, is that FEniCS solvers have support for both distributed (MPI) and shared memory (OpenMP) parallel execution allowing for larger models to be solved faster (FEniCS has been tested on problem sizes up to 10^{8} degrees of freedom run on 512 CPUs in parallel, where Matlab and Octave in contrast is limited to the serial sparse linear solvers, the Umfpack direct solver and a selection of iterative ones.

A button labeled **FEniCS** is present in the *Solve Mode* toolbar of the FEATool GUI. Pressing this button will open a dialog box where the current model equations and parameters has been translated to a FEniCS python script. The script can there be inspected and edited, as well as changing the output file name and system shell command to run it (which is *bash* by default).

The FEniCS solver dialog box also features **Solve**, **Auto**, **Export**, and **Import** buttons. If the **Auto** toggle button is engaged the **Solve** button will automatically try to **Export**, **Solve**, and **Import** the computed solution. With the **Auto** toggle disengaged the **Solve** button only execute the FEniCS python script in a bash shell if present (not export or import). The **Export** option saves the current grid in Dolfin XML format and exports the FEniCS simulation script and **Import** will import an existing solution if it matches with the FEATool problem definition.

Additional options include the base *File Name* and overriding the default bash *Solve Command*.

The installed fenics Matlab function can also be used on the command line (CLI) to manually perform the export, solve, and import actions. In the following is an example to set up a simple heat transfer model on a unit circle with a unit heat source term, and *T=0* fixed temperature on all the boundaries. In the FEATool Matlab m-script language the model can look like the following

fea.sdim = {'x' 'y'}; fea.grid = quad2tri(circgrid()); fea = addphys(fea, @heattransfer); fea.phys.ht.eqn.coef{6,end} = {1}; % Heat source term. fea.phys.ht.bdr.sel(:) = 1; % Zero temperature boundary conditions. fea = parsephys(fea); % Parse physics mode and fea problem struct. fea = parseprob(fea); fea.sol.u = solvestat(fea); postplot(fea, 'surfexpr', 'T')

By running the command **fenics(fea, 'mode', 'export')** a grid file *featool-fenics-mesh.xml* and FEniCS simulation script file *featool-fenics.py* will automatically be generated in the current directory

from fenics import * # Mesh and subdomains. mesh = Mesh("featool_fenics_mesh.xml") # Finite element and function spaces. E0 = FiniteElement("P", mesh.ufl_cell(), 1) V = FunctionSpace(mesh, E0) T = Function(V) T_t = TestFunction(V) # Model constants and expressions. rho = Constant(1) cp = Constant(1) k = Constant(1) u = Constant(0) v = Constant(0) q = Constant(1) # Bilinear forms. a = ( (rho*cp*u)*T*T_t.dx(0) + (k)*T.dx(0)*T_t.dx(0) + \ (rho*cp*v)*T*T_t.dx(1) + (k)*T.dx(1)*T_t.dx(1) )*dx # Linear forms. f = q*T_t*dx # Boundary conditions. dbc0 = DirichletBC(V, Constant(0), 0) dbc1 = DirichletBC(V, Constant(0), 1) dbc2 = DirichletBC(V, Constant(0), 2) dbc3 = DirichletBC(V, Constant(0), 3) dbc = [dbc0, dbc1, dbc2, dbc3] # Initial conditions. assign(T, interpolate( Constant(0), V)) # Solve. solve( a - f == 0, T, dbc ) # Output. import numpy as np T_h = T.compute_vertex_values(mesh) np.savetxt("featool_fenics_sol.txt", np.column_stack((T_h))) # Postprocessing. plot(T, title = "T") interactive()

The generated FEniCS python FEM script is longer and somewhat more verbose since all the FEATool physics mode defaults must be explicitly expressed.

Calling **fenics(fea, 'mode', 'solve')** will attempt to solve the problem if python and FEniCS is installed and set up correctly. Alternatively, the script can be run by itself by a valid FEniCS installation.

The solution process will generate the file *featool-fenics-sol.txt* containing the nodal solution(s) which can be imported back into FEATool with **fea = fenics(fea, 'mode', 'import')** after which it can be postprocessed and visualized with the usual FEATool and Matlab functions.

The FEniCS-FEATool integration and problem file export should work for general multiphysics problems using both the GUI and command line. However, a subset of FEATool functionality is currently not supported (the known ones are listed here)

- Non-linear Dirichlet fixed value boundary conditions are not supported. The exported FEniCS C++ boundary condition expressions do not currently allow for dependent variables / solution unknowns.
- As FEniCS does not feature a built-in time dependent solver. Solvers for instationary problems must currently be implemented manually.
- Switch and logical expressions (for example
*x>1 & y<=0*) must be manually rewritten to conform with the C++ expression syntax. - Point source terms are not supported by the FEniCS non-linear solution form (which is used by default since it also handles linear problems).
- Although FEATool and FEniCS uses the same finite element FEM basis functions the internal ordering is different, thus the FEniCS solution will be exported in the grid points (
*P1*space). Thus even if higher-order FEM spaces are used solution accuracy will be lost during FEATool import. This issue might not be that important for visualization but when calculating quantities, boundary and subdomain integrals, and expression evaluations it would currently be more accurate to perform that in FEniCS and python. - FEniCS does not currently support quadrilateral or hexahedral grid cells and must be converted to triangles and tetrahedra, respectively. Grid conversion of these types of grids can be performed directly in the GUI or with the quad2tri and hex2tet commands. Typically this will not be an issue since FEATool and the automatic mesh generator gridgen by default creates simplex grids (line segments, triangles, and tetrahedra).
- In addition to supporting FEniCS, the exported python simulation scripts should also be compatible with the Firedrake project solver which also uses the FEniCS Unified Form Language (UFL) for problem definitions.

Note that the FEATool distribution does not include a Python interpreter or FEniCS itself which must be installed separately. The FEniCS homepage provides instructions how to install FEniCS on *nix systems and pre-configured Docker Linux images.

For systems running *Windows 10* FEniCS can be installed with the Ubuntu Bash Windows Subsystem for Linux by simply opening a Windows Bash shell and running the FEniCS on Ubuntu commands (which automatically also installs Python if required)

sudo add-apt-repository ppa:fenics-packages/fenics sudo apt-get update sudo apt-get install --no-install-recommends fenics sudo apt-get dist-upgrade

To allow FEniCS and python plotting and visualization on Windows one must also install an X window server such as Xming (note that plotting of the solution variables is commented and disabled in the FEATool-FEniCS scripts by default).

FEATool uses the default built-in linear solvers available in Matlab and Octave which currently is the direct Umfpack and SuiteSparse solver suites. These are very robust for a wide range of problems but not as efficient either in memory or speed as a specialized dedicated solver can be. Starting with FeatFlow, FEATool provides functionality for using more efficient specialized external solvers.

FeatFlow, is very fast and efficient FEM based computational fluid dynamics (CFD) solver for the incompressible Navier-Stokes equations. The FeatFlow solver makes it easy to perform high performance CFD simulations with FEATool Multiphysics directly in Matlab and Octave.

FeatFlow is a finite element CFD code based on using an efficient FEM discretization (Rannacher-Turek non-conforming ansatz functions) together with a geometric multigrid approach for solving linear systems [1],[2]. This results in a very fast and computationally efficient solver typically yielding a magnitude or more of speedup compared to the built-in direct solver in Matlab (UMFPACK). Additionally, an iterative solver, such as in FeatFlow, uses much less memory than a direct one allowing for much larger simulations. Being a very stable and tested CFD code, FeatFlow has been used and validated in many commercial CFD projects, and studies have shown FeatFlow to be significantly more efficient than both Ansys CFX and OpenFOAM with respect to the total CPU time required to achieve a target accuracy [3].

The video below is a tutorial showing how to set up and solve a backwards facing step CFD problem with the FEATool GUI. For this backwards facing step benchmark example the FeatFlow solver is about 50 times faster compared to using the default built-in solver.

FeatFlow strictly employs quadrilateral grids in two dimensions. In the FEATool GUI one can use the **Grid > Convert Grid Cells** menu option to change an unstructured triangular simplex grid into quadrilateral one. On the Matlab and Octave command lines one can also use the tri2quad, gridgen_quad, and grid primitive functions such as rectgrid, circgrid etc. to generate quadrilateral grids. Alternatively, one can also import a pre-made external quadrilateral grid in any of the import format that are supported in FEATool, such as for example from GiD and Gmsh.

In contrast to most CFD and physics simulation codes it is in FeatFlow desirable to start with a very coarse grid. FeatFlow will then internally uniformly refine this grid a prescribed number of times to generate the multigrid level hierarchies, and also adapt boundaries to the geometric boundary parametrization. In this way one can achieve optimal conditions for the geometric multigrid solver. This also means that the output from FeatFlow will correspond to a much finer grid and have a higher quality solution.

FeatFlow can be used instead of the usual solve button **=** to call the external FeatFlow solver for Navier-Stokes physics mode models. The FeatFlow solver control panel can be accessed in *Solve Mode* of the FEATool GUI by pressing the button labeled **FeatFlow**.

The *FeatFlow* solver dialog box allows one to automatically use *FeatFlow* to solve CFD problems. If the **Auto** toggle button is engaged the **Solve** button will automatically try to **Export**, **Solve**, and **Import** the computed solution. With the **Auto** toggle disengaged the **Solve** button only triggers the corresponding FeatFlow CC2D or CC3D solver (not export or import). The **Export** export button manually exports the current grid and model definition to FeatFlow data files. **Import** manually reads a computed FeatFlow solution and imports it into the FEATool GUI (refining the initial grid as necessary determined by the FeatFlow *Grid refinement level* setting).

The FeatFlow solver settings are described as follows

**Grid refinement level**- Number of grid refinements of the imported grid to compute on. FeatFlow uses a geometric multigrid solver approach where the initial coarse input grid is uniformly refined several times to generate the finer multigrid levels. Thus it is often desirable to start with a quite coarse input grid.**Artificial stabilization**- Select between streamline diffusion and upwinding for artificial stabilization of convective terms.**Stabilization parameter**- Stabilization tuning parameter.**Solver type**- Select either*Stationary (monolithic)*or*Time-Dependent*problem and solver type.**Time stepping scheme**- Choose between implicit 2nd order*Fractional-step-theta*,*Crank-Nicolson*, and 1st order*Backward Euler*time stepping schemes.**Time step**- Sets the macro time step size (internally FeatFlow takes three sub-time steps for each macro time step).**Simulation time**- Maximum/final simulation time.**Time stopping criteria**- Stationary limit to stop a time dependent simulation.**Full settings**- Access all FeatFlow CC2D and CC3D solver parameters (see the FeatFlow documentation for a description of these) [1].

The featflow function can also be used on the command line as follows:

FEA = FEATFLOW( FEA, VARARGIN ) Export, solves, or imports the solved problem described in the FEA finite element struct using the FeatFlow CFD solver. Accepts the following property/value pairs

Input Value/{Default} Description ----- ---------------------------- --------------------------- mode export/solve/import Command mode to call data default Default solver parameters solver '' FeatFlow solver to use fdir featflow Processing directory fname featool-featflow FeatFlow output filename wbash C:\Windows\System32\bash.exe Windows bash path clear boolean {true} Clear output and log files

The following limitations apply to the FeatFlow solver.

- Single geometry object with quadrilateral/hexahedral grid
- The solvers are statically compiled with a 240 million double precision size memory limit
- Single Navier-Stokes physics mode
- Constant density and viscosities
- Zero initial conditions
- Prescribed velocity, wall, and neutral Neumann outflow boundary conditions
- Prescribed velocity boundary condition must be linear expressions and may not include dependent variables (velocities or pressure)

The FeatFlow solvers *cc2d* and *cc3d* included with FEATool are capable of performing large scale simulations. Although technically possible, for memory and stability reasons it is not advised to do this from within Matlab or Octave. For large scale simulations then it is instead recommended to manually launch the FeatFlow *cc2d* / *cc3d* binaries from the system command line. To do this first make sure the proper FeatFlow *data* and *grid* files have been exported. The *data* files *cc2d.dat* / *cc3d.dat*, *bdry.dat*, and *expr.dat* should be located in the *featool/featflow/data* subdirectory, and the grid files *featool-featflow.tri* / *featool-featflow.prm* in the *featool/featflow/adc* directory. The *ccXd* binaries can then simply be launched as, for example

cc2d

Convergence output will be printed to the terminal screen as well as mirrored in the log file *featflow/data/ff.out*. As memory FeatFlow is statically allocated, four different binaries may be supplied for each solver with a postfix *1* through *4*. If the grid and memory requirements are too large for one binary one can try using the one with a higher number. Result output files can be found in the *featool/featflow/output* subdirectory for the pressure, velocities, cell centers, and cell edge/face midpoints. Moreover, the ascii input data files can be tuned accordingly to provide for example binary, GMV, and AVS output data (see the FeatFlow documentation for a full explanation of the solver input parameters).

An example of setting up and solving a flow over a backwards facing step problem using the FeatFlow external solver is given in the tutorial section Fluid Dynamics - Flow Over a Backwards Facing Step.

[1] S. Turek, C. Becker, FeatFlow, Finite element software for the incompressible Navier-Stokes equations, User Manual, Release 1.1, Heidelberg, 1998.

[2] S. Turek, Efficient Solvers for Incompressible Flow Problems: An Algorithmic and Computational Approach, Series: Lecture Notes in Computational Science and Engineering ,Volume 6, Springer-Verlag, 1999.

[3] E. Bayraktar, O. Mierka, S. Turek, Benchmark Computations of 3D Laminar Flow Around a Cylinder with CFX, OpenFOAM and FeatFlow, International Journal of Computational Science and Engineering, 7, 3, 253-266, 2012.