FEATool  v1.8
Finite Element Analysis Toolbox
 All Files Functions Pages
Solvers

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

The default option in FEATool is to use the best linear solvers available in Matlab and Octave, which currently is the Umfpack and SuiteSparse direct linear solvers. These solvers are very robust for a wide range of problems, but are not as memory efficient or fast as specialized and dedicated solvers can be. Due to this FEATool also features GUI and CLI interfaces to specialized external solvers, such as the computational fluid dynamics (CFD) solver OpenFOAM, and the general FEA and PDE solver FEniCS. These interfaces allows one to set up and define multiphysics problems in Matlab and Octave with the FEATool GUI and CLI interfaces, export and solve them with the external solvers, and automatically import the solutions for further (post-)processing. The OpenFOAM and FEniCS external solvers are discussed in the following sections.


Solve Mode

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.


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.

gui_solve_50.png

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


OpenFOAM CFD Solver

OpenFOAM, is very popular, flexible, and efficient finite volume FVM based computational fluid dynamics (CFD) solver for both incompressible and compressible laminar and turbulent flow regimes. The stable and tested OpenFOAM 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 high performance CFD simulations directly in Matlab and Octave, typically yielding a magnitude or more speedup compared to the built-in solvers. Additionally, iterative and multigrid solvers, such as included in OpenFOAM, uses much less memory than direct solvers allowing for much larger simulations.


OpenFOAM Usage

To use OpenFOAM to solve Navier-Stokes fluid 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 and re-imports the computed solution into FEATool. 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. Settings option allows for manually editing the OpenFOAM control dictionary and case files and overriding the default settings.

openfoam_dlg_50.png

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

  • Time scheme - Selects between Steady State, Euler, Crank-Nicolson, Backward, and Local Euler time discretization schemes.
  • Time step - Specifies the time step size.
  • End time - Sets the simulation maximum end time.
  • Stopping criteria - Specifies the stopping criteria for steady state simulations (corresponding to the residualControl parameter).
  • Turbulence model - Selects between the pre-defined turbulence models Laminar and k-epsilon.


Command Line (CLI) Usage

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

FEA = OPENFOAM( FEA, VARARGIN ) Export, solves, or imports the solved problem described in the FEATool FEA problem struct using the OpenFOAM CFD solver. Accepts the following property/value pairs, moreover also accepts the openfoam_data property/value pairs to set the OpenFOAM dicts during export.

Input          Value/{Default}        Description
-----------------------------------------------------------------------------------
modes          check, export,         Command mode(s) to call
               solve, import
data           default                Default OpenFOAM data and parameter dict
application    string {simpleFoam}    OpenFOAM application binary to run
ddtScheme      string {steadyState}   Time stepping scheme
tolres         scalar/vector {1e-5}   Stopping criteria for residuals (Simple)
startTime      scalar {0.0}           Simulation start time
endTime        scalar {100}           Simulation end time
deltaT         scalar {1.0}           Time step size
maxDeltaT      scalar {0.1}           Maximum time step size
maxCo          scalar {0.5}           Maximum Courant number
writeInterval  scalar {100}           Solution output write interval
purgeWrite     scalar {0}             Specified number of output solutions
transportModel string {Newtonian}     Transport model
simulationType string {laminar}       Simulation type laminar/RAS/LES
RASModel       string {kEpsilon}      RAS turbulence model
LESModel       string {Smagorinsky}   LES turbulence model
nu             scalar {1.0}           Kinematic viscosity (constant)
interp         true                   Interpolate solution to grid points
casedir        default                OpenFOAM case directory
foamdir        default                OpenFOAM installation directory
logfname       featool-openfoam       OpenFOAM log/output filename
logfid         1                      Log file/message output file handle
tolres         1e-5                   Stopping criteria for initial residuals
clear          boolean {true}         Clear output and log files


Installation

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

For Microsoft Windows systems it is generally recommended to install and use the pre-compiled blueCFD-core (2017) binaries from blueCAPE http://bluecfd.com/Core. For Ubuntu Linux OpenFOAM can be installed with the apt-get package manager

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 call the OpenFOAM solver binaries simpleFoam and pimpleFoam through the use of the openfoam.sh and openfoam.bat shell script which may be modified to suit the users need (which can be found in the openfoam folder).


Limitations

The following limitations currently apply to the OpenFOAM solver.

  • Single Navier-Stokes physics mode
  • Single geometry object with one subdomain
  • Constant density and viscosity
  • Constant initial conditions
  • Constant velocity/pressure boundary conditions


Advanced Usage

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 or Octave.

For large scale simulations then it is instead recommended to first export a FEATool model with the openfoam export command. 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).


Tutorial

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


FEniCS General FEA Solver

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


GUI Usage

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

fenics_dlg_50.png

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.


Command Line Usage

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.


Notes

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.


Installation

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 required dependencies such as the Python programming language interpreter)

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