FEATool Multiphysics
v1.11 Finite Element Analysis Toolbox |

Solvers

This chapter describes the solver options available with FEATool Multiphysics. Solvers are used to compute solutions to the involved discretized equations.

In addition to the built-in multiphysics solvers FEATool also features integrated GUI and CLI interfaces to specialized external solvers, the computational fluid dynamics (CFD) solver OpenFOAM and the general FEM and PDE solvers FEniCS and Firedrake. These interfaces allows one to set up and define multiphysics problems in MATLAB with the FEATool GUI and CLI interfaces, export and solve them with the external solvers. Afterwards, the solutions can automatically be imported into the FEATool GUI for further post-processing. Moreover, this multi-simulation solver approach allows for several solvers and numerical methods to be used for the same problem, which is useful to compare and validate results and solutions. The OpenFOAM and FEniCS external solvers are discussed in the following sections.

After 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 *Solve* 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 by using the **Solver Call Command...** option, allowing for custom and external solvers to be used.

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 OpenFOAM external CFD solver.
- Solves general finite element problems with the FEniCS external FEA 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 the built-in **Stationary** (solvestat), **Time-Dependent** / instationary (solvetime), or **Eigenvalue** (solveeig) solver to be selected. The built-in solvers are monolithic, meaning that all equations and 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 already **Computed 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 *Linear solver* used to solve the sparse linear systems can be selected from the built-in **backslash** operator (Umfpack), **mumps**, or the iterative **gmres** and **bicgstab** solvers. The default built-in option in FEATool is the direct solver Umfpack/SuiteSparse. Alternatively, the direct solver mumps can be up to 20% faster for selected problems. These solvers are very robust and fast for a large range of problems, but are not as memory efficient as iterative solvers such as *gmres* or *bicgstab* with ILU(k) preconditioning. Users can also user their own linear solvers by overriding the solvelin function.

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

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 full, standard row-sum lumped, diagonal lumped, or consistent HRZ-lumped **Mass matrix type** representations.

In the *Eigenvalue Settings* frame one can control the **Number of eigenvalues** to be computed, as well as the **Type of eigenvalues** to find.

Fluid-structure interactions (FSI) problems feature a dedicated solver which handles the physics coupling and allows for movement of the grid and domains by using a solid extension mesh motion technique with Jacobian-based stiffening [1]. The following solver settings are available for FSI problems.

The FSI solver supports time-dependent problems with either fully-implicit or semi-implicit discretizations (with explicit geometry treatment). The *Initial Condition* section allows for either initial conditions from a pre-defined constant or expression, or interpolated from a previously computed solution.

The *Time Dependent Settings* section allows for prescribing the **Initial time**, **Time step** size, and **Final time**. Specific **Output times** can also be supplied as a vector of specific times to save and output (default [t_{0}:t_{step}:t_{end}]). Moreover, the order and accuracy of the employed time discretization schemes can also be selected (1st to 5th for the fluid domain using a BDF discretization, and Newmark scheme for the solid domain).

[1] A. Quarteroni, A. Manzoni, F. Negri. Reduced Basis Methods for Partial Differential Equations. An Introduction, Springer, 2016.

OpenFOAM, is very popular, flexible, and efficient finite volume method (FVM) based computational fluid dynamics (CFD) solver. OpenFOAM features a wide selection of specialized CFD solvers for both incompressible and compressible, laminar and turbulent, iso-thermal and non-isothermal flow regimes. The stable and tested code base has been used and validated in many academic and commercial CFD projects.

The *OpenFOAM-FEATool* CFD solver integration makes it easy to perform CFD simulations directly in MATLAB, typically yielding a magnitude or more speedup compared to the built-in solvers. Additionally, advanced iterative and multi-grid linear parallel solvers, such as included in OpenFOAM, uses much less memory than direct solvers allowing for significantly larger simulations.

To use OpenFOAM to solve Navier-Stokes, Swirl, or Compressible Euler flow problems, press the corresponding *Solve Mode* toolbar button labeled **OpenFOAM**, instead of the default solve button. This will open the OpenFOAM solver settings and control dialog box.

The *OpenFOAM* solver dialog box allows one to automatically use *OpenFOAM* to solve CFD problems. The **Solve** button automatically export, solves, (interpolates) and re-imports a computed solution into FEATool. **Stop** halts the solution process and plots the intermediate solution, while **Cancel** terminates and discards the current solution (Please note that it can take some time for the OpenFOAM solver register a halt event and stop). **Export** allows for the OpenFOAM dictionary files for external manual processing and editing. If the **Close automatically** checkbox is marked, the OpenFOAM dialog box will be automatically closed and FEATool switch to Postprocessing mode after the solution has been computed.

When OpenFOAM solver is running one can switch between the *Log* and *Convergence* tabs to see and monitor the real-time solver output log and convergence plots, respectively.

The following OpenFOAM solver settings can be controlled though the dialog box

**Time discretization scheme**- Select between*Steady State*,*Euler*,*Crank-Nicolson*,*Backward*, and*Local Euler*time discretization schemes.**Initial Conditions**may be specified as either constant or subdomain expressions for the flow variables, a potential flow solution, or use an existing solution.**Time step**- Specifies the time step size.**End time**- Prescribes the simulation maximum end time.**Stopping criteria**- Specifies the stopping criteria for steady state simulations (corresponding to the*residualControl*parameter).

**Turbulence Modeling** allows for *laminar*, *k-epsilon*, and *k-omega* turbulence models to be used with standard wall functions. If not known and prescribed directly, the turbulent quantities *k*, *epsilon*, and *omega* are estimated as

\[ \begin{aligned} & k = 3/2\cdot (|\mathbf{u}|\cdot I_{turb})^2 \\ & \epsilon = C_{\mu}^{3/4}\cdot k^{3/2}/l_{turb} \\ & \omega = C_{\mu}^{-1/4}\cdot k^{1/2}/l_{turb} \end{aligned} \]

where \(C_{\mu} = 0.09\) is a turbulence modeling constant, \(I_{turb}\) the turbulent intensity (1% - 10% recommended), and \(l_{turb}\) a length scale for the turbulent eddies (default 8% of the boundary length). Press **Recompute** after changing the *Intensity* and/or *Length* scale to recompute the turbulent inlet quantities.

Note that the OpenFOAM solver binaries are not included with the FEATool distribution and must be installed (and compiled) separately.

The OpenFOAM MATLAB solver integration has been verified with OpenFOAM version 5. For Microsoft Windows systems it is recommended to install and use the pre-compiled blueCFD-core (2017) binaries from blueCAPE/blueCFD.

For Ubuntu based Linux and Windows Subsystem for Linux system, OpenFOAM can be installed by opening a terminal shell and running the following commands (which automatically also installs any required dependencies)

sudo add-apt-repository http://dl.openfoam.org/ubuntu sudo sh -c "wget -O - http://dl.openfoam.org/gpg.key | apt-key add -" sudo apt-get update sudo apt-get -y install openfoam5

For other systems, please follow the installation recommendations from the official OpenFOAM website http://openfoam.org.

The FEATool openfoam function will attempt to perform system calls the OpenFOAM solver binaries *simpleFoam*, *pimpleFoam*, *rhoCentralFoam*, *potentialFoam*, and *collapseEdges*. It is therefore necessary that the binaries are installed and properly set up so they can be called from system script files (bash scripts on Linux and MacOS and bat/vbs on Windows).

The following limitations currently apply to the OpenFOAM solver.

- Single Navier-Stokes, Swirl Flow, or Compressible Euler physics mode
- Single geometry object with one subdomain
- Isothermal with constant density and viscosity

The openfoam function can be used instead of *solvestat* and *solvetime* to solve CFD problems with OpenFOAM on the MATLAB command line (CLI). In the following is an example of laminar steady flow in a channel solved with OpenFOAM

fea.sdim = {'x','y'}; fea.grid = rectgrid(50,10,[0,2.5;0,0.5]); fea = addphys(fea,@navierstokes); fea.phys.ns.eqn.coef{1,end} = { 1 }; fea.phys.ns.eqn.coef{2,end} = { 1e-3 }; fea.phys.ns.bdr.sel(4) = 2; fea.phys.ns.bdr.coef{2,end}{1,4} = '2/3*0.3'; fea.phys.ns.bdr.sel(2) = 4; fea = parsephys(fea); fea = parseprob(fea); % Alternative to calling: fea.sol.u = solvestat( fea ); fea.sol.u = openfoam( fea ); postplot(fea,'surfexpr','sqrt(u^2+v^2)')

The model parameters are taken from the ex_navierstokes1 example script model. Most of the *ex_navierstokes** models in the examples directory features a *solve* flag which can be used to enable the OpenFOAM solver instead of the default solver. The OpenFOAM function can also be embedded in custom user-defined m-scripts, which can use all standard MATLAB functions and toolboxes.

The OpenFOAM solvers are capable of performing large scale simulations, and although technically possible to use FEATool to do this, for memory and stability reasons it is not advised to do this from within MATLAB.

For large scale simulations it is recommended to first export a FEATool model with the openfoam export command or *Export* button in the *OpenFOAM Settings* dialog box. This will generate corresponding OpenFOAM case files. Then one can manually launch the OpenFOAM solver the system command line, and possibly manually import the solution when the solution has finished (consult the OpenFOAM documentation for a full explanation how to use OpenFOAM as a stand alone solver).

An example of setting up and solving a turbulent flow over a backwards facing step problem using the OpenFOAM external solver is given in the tutorial section of the Turbulent Flow Over a Backwards Facing Step (OpenFOAM) example.

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 OpenFOAM 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 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, 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 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')

Similar to the solvetime and solvetime functions, the FEniCS solver can be called with a fenics function call

fea = fenics( fea );

This automatically converts the FEATool problem defined in the *fea* struct and calls FEniCS to compute a solution to the problem. The returned solution will be assigned to the *fea.sol.u* field.

If required the solver *modes* or phases can be called separately (available modes are *check*, *export*, *solve*, and *import*). For example, for the *export* mode

fenics( fea, 'modes', 'export' )

generates a Dolfin XML grid file *featool-fenics-mesh.xml* and FEniCS simulation script file *featool-fenics.py* as well as control shell files to call the FEniCS solver. For example, the FEniCS *Python* FEM script for the model above is shown below

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)))

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

The solution process will generate the file *featool-fenics-sol.txt* containing the nodal solution(s) which can be manually imported back into FEATool with

fea = fenics( fea, 'modes', 'import' );

which it can be postprocessed and visualized with the usual FEATool and MATLAB functions.

For reference, the following command line arguments are supported by the fenics solver function

Input Value/{Default} Description ----------------------------------------------------------------------------------- modes data, export, solve, import Command mode(s) to call (default all) data default FEniCS input data file to use fname featool-fenics FEniCS base filename root fdir string {pwd} Directory to write mesh and data files order scalar {2} Integration order used in c-expressions scmd system default Custom system solve command string ccmd string FEniCS installation bash check command python -c "import dolfin;print(dolfin.dolfin_version())" wbash C:\Windows\System32\bash.exe Windows bash executable path clear boolean {true} Clear output and log files pause boolean {false} Pause external log terminal fid scalar {1} File identifier for output ([]=no output)

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)

- 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. - 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).
- 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).

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 (Docker is currently not supported by FEATool).

To use the FEniCS-FEATool solver integration on Microsoft Windows systems requires *Windows 10* or later, with the Windows Subsystem for Linux installed.

For Ubuntu based Linux and Windows Subsystem for Linux system, FEniCS can be installed by opening a terminal shell and running the following commands (which automatically also installs required dependencies such as the *Python* programming language interpreter)

sudo apt-get install software-properties-common sudo add-apt-repository ppa:fenics-packages/fenics sudo apt-get update sudo apt-get install --no-install-recommends fenics

Although not necessary for FEATool-FEniCS functionality, to allow FEniCS and plotting and visualization with *Python* on Windows it is also necessary to install an X window server such as Xming (note that plotting and output of solution variables is disabled in the FEATool-FEniCS script export by default).

The following functions can be used on the command line for to solve *fea* problems.

Function | Description |
---|---|

fenics | Solve problem with the external FEniCS solver |

openfoam | Solve problem with the external OpenFOAM solver |

solveeig | Solve eigenvalue/frequency problem |

solvestat | Solve stationary problem |

solvetime | Solve time-dependent/instationary problem |