In the previous episodes we solved a few equations with Elmer. We did some choices when we setup the solver parameters. What those parameters do, and how should we set them? This is perhaps the trickiest part in FEM (beside making the mesh right). In this episode we will step back and look at those solver options more closely. This post is really not meant to be an exhaustive explanation. For that, refer to the Elmer documentation.

Linear System and Nonlinear System

Every time we set the Elmer Solver Settings, we typically do it through the Solver control window of ElmerGUI. ElmerGUI is very convenient, especially to start getting our hands dirty with FEM, but it is essentially a GUI layer around various utilities Elmer provides. One of the most important operations it does is to prepare a solver input file (sif) for us. As we will see later in the series, not all Elmer models are supported by the GUI, so we will need to write at least parts of the sif by ourselves at some point in the future. In any case, what Solver settings we should choose for our given problems?

Turns out this question is actually quite hard to answer.

Keep your Equation in Mind

We should not think of Elmer, or any other solver for the matter, as a replacement for our brain. Rather, we should try to understand our problem as best as possible, so to setup the solver in the best possible way. Turns out that this is not always easy, and maybe we will need some trial and error of different solver options, but figuring out the details of our problem will give us a jump-start. With Elmer most issues you might face are about not setting up the solver properly. These, in my experience, are the most common issues after mesh related issues.

Your Equation

Look at the equation you aim to solve. For example, so far we dealt with Linear Elasticity and Helmholtz equations, which I am going to report below as formulated in the Elmer models (notation might be slightly different with respect the last episodes, but it will be clarified here).

Linear Elasticity Equation

$$\rho\frac{\partial^{2}\mathbf{d}}{\partial t^{2}}-\nabla\cdot\tau=\mathbf{f}$$

Where $\rho$ is the volumetric mass density of the body we are modelling, $\mathbf{d}$ is the displacement (from rest) field (a 3D vector), $\tau$ is the stress tensor (more details here) and $\mathbf{f}$ is a body force field, a force value we apply to each point of the body we are modelling (it is a 3D vector and each point of the body can have a different force applied to it).

To make things clearer, let’s expand the divergence:

$$\rho\frac{\partial^{2}\mathbf{d}}{\partial t^{2}}-\sum_{i=x,y,z}\sum_{k=x,y,z}\frac{\partial\tau_{i,k}}{\partial k}=\mathbf{f}$$

Where $\tau_{i,k}$ are the various components of the stress tensor.

Helmholtz Equation

$$\left(k^{2}-jkD\right)P+\nabla^{2}P=0$$

Where $k$ is the wave number, $D$ is an optional dumping factor, $j$ is the imaginary unit and $P$ is the spatial part of the acoustic pressure field. To make things clearer let’s expand the Laplace operator:

$$\left(k^{2}-jkD\right)P+\frac{\partial^{2}P}{\partial x^{2}}+\frac{\partial^{2}P}{\partial y^{2}}+\frac{\partial^{2}P}{\partial z^{2}}=0$$

See here for more information.

First Question: are they Linear?

To understand whether the PDE is linear, we need to look at what happens to our unknown field. Being a PDE, our unknown field will undergo differentiation. If all the terms of the PDE are linear, then the PDE is linear. A linear term is a term in which a derivative is used as is, or eventually multiplied by a number or any function different from the unknown field. For example, if $F$ was an unknown scalar field (but it is the same for vector fields) the terms below would be all examples of linear terms:

$$\frac{\partial F}{\partial x},\space 4\frac{\partial^{7}F}{\partial z^{7}},\space x^{4}\frac{\partial^{3}F}{\partial y^{3}},\space g\left(x,y,z\right)F$$

While the terms below are nonlinear:

$$F\frac{\partial F}{\partial x},\space \ln\left(\frac{\partial^{5}F}{\partial x^{5}}\right),\space \left(\frac{\partial^{3}F}{\partial y^{3}}\right)^{6},\space \frac{g\left(x,y,z\right)}{F}$$

Let’s look at the Helmholtz equation first. $P$ is our unknown field, and clearly there are only linear terms in the Helmholtz equation. So, the equation is linear.

Let’s look now at the linear elasticity equation instead. The equation itself is, well, linear, as the derivatives of the unknown field are linear terms. However, there might be a catch! The stress tensor $\tau$ depends on the strain, and the strain in turns depend on our unknown field $\mathbf{d}$ and its derivatives. In linear materials, or if we use a linearised formulation, then the dependency is still linear, and in the double sum above we will have linear terms of mixed derivatives of the unknown field. So, the linear elasticity equation is linear as long as we treat the material as linear. Things simplify even more if the material is isotropic (same elasticity properties along every dimension) as in that case we can figure out the stress tensor out of two numbers: Poisson ratio and Young modulus.

Note that we could say a similar thing about acoustics: the acoustic pressure will disturb the material properties. However, the Helmholtz equation is derived from the linear wave equation which, in turn, is derived in the approximation of adiabatic wave propagation, so these effects are not accounted for. We will see how things change for more advanced acoustic equations.

So, be sure you understand the material parameters too, as they might depend non-linearly on the unknown field. In our previous studies we specified isotropic linear parameters for linear elasticity, so our equations were linear.

The Linear System Parameters

Now that we know whether our equation is linear or nonlinear, we can understand better what the various parameters of our solver mean.

Linear PDEs are turned by FEM into linear systems of algebraic equations that are formulated as matrix equations. These matrices are big, as they are assembled from the various equations at each single node of the mesh.

These matrix equations can be solved directly. We did so for linear elasticity. This is a very good method, as it will solve the matrix equation exactly within numerical precision (see here), but direct methods do not scale well for very big problems (very fine meshes and/or high mesh order) as they will require massive amounts of memory and CPU time. If your problem is small enough, end even better, if it is 2D or 1D, I would recommend you start off with the Umfpack direct method.

To work around this limitation of direct methods, iterative methods can be used, where the linear system is solved by iterative algorithms. This is where things become complicated.

Iterative methods bring in convergence. So we ideally want two things:

• Converge in a small number of iterations.
• Compute each iteration in a reasonably small time.

But what convergence means? We touched on this in the FEM in a Nutshell episode, and we will say something more here. In this case, it means that the solution is not varying a lot between consecutive iterations, hence we are reaching the final correct solution of the linear system. We can specify the convergence tolerance for Elmer, and I would suggest to start from the default value of $10^{-10}$. Additionally, we can specify the maximum number of iterations. After the number of iterations is reached, Elmer will stop, but be aware that our solution will most likely be inaccurate if we did not reach convergence, for this reason we should tick the Abort if the solution did not converge box (if we are using ElmerGUI).

There are many different iterative methods that one could use, each one best suited for a certain type of linear system. I would recommend to start with BiCGStabl, with a BiCGStabl order of $2$. This method should give you a smoother convergence in most cases.

However, there is one more thing: preconditioning. Preconditioning is an operation carried on the matrices to aid convergence rate. The best preconditioning has to be searched by trial and error. I normally start with ILUT and a default ILUT tolerance of $10^{-3}$.

Finally, Elmer provides Multigrid methods, but I am not familiar with those.

For more information, please refer to the Elmer Solver Manual, especially this chapter. As it is not easy to nail down what the matrices for complex geometries and/or complex (multi)physics will look like, trial end error becomes important. For big problems you will use iterative or multigrid methods.

As suggested in the The FEM Pipeline episode, I would recommend that you use one or more simpler and smaller simplified problems to monitor both convergence and accuracy: quickly iterate among your simplified problems with different solution methods and preconditioning, until you are satisfied you minimised both the number of iterations needed for convergence and the time each iteration needs, all while ensuring good accuracy of your model. Solver convergence information will be provided by the Elmer log, while accuracy of the FEM solution to its analytical counterpart (or, eventually, measurement of real systems) has to be assessed with the most appropriate method, which varies by problem. Be aware that not all iterative methods bring smooth convergence.

This is clearly quite hard to do in practice, and the problems we will solve in the future will allow us to touch this development cycle with hand.

The Nonlinear System Parameters

Well, equations can be nonlinear, as we seen. So, how can the maths of linear systems deal with them? Well, the nonlinear system can be linearised and solved iteratively. These are the so called nonlinear iterations. Most solvers for nonlinear equations have tunable linearisation methods, so we will see most of this stuff later. For our linear equations we could (and perhaps should) set the Max. iterations to $1$, as we essentially do not need them. We did so for linear elasticity, but I left the default value for our spherical source because I forgot about it on purpose so that I could explain a few things here.

Elmer reports the nonlinear convergence plot during the solver execution, which is shown above for our Pulsating Sphere problem. We can see that, in fact, we did only one nonlinear iteration before exiting, as the solution did not change in the slightest between the first solution and the first nonlinear iteration, as expected for linear problems.

We will touch more on this when we get to nonlinear PDEs.

Solver Specific Options

Well, these are solver specific, so they change between solvers. There is very little general stuff to say about this. It is best to refer to the Elmer Models Manual to check what we can do. With the Solver Specific Options we can change what we solve for (for example, setup Eigen Analysis for linear elasticity, or set special numerical algorithms which can aid accuracy and convergence as well).

Conclusion

This post was more about raising awareness about how nontrivial it is to solve problems with FEM. The best possible way to ensure that we know how to solve a problem properly is do a few studies with benchmark solutions, as we are doing in this first episodes, and compare them to known analytical solutions. To summarise, the steps to setup our study properly are as follows:

• Read the Elmer Models Manual and figure out whether you need or not the Nonlinear iterations. Be careful that non-linearity can be hidden into the materials dependency on the unknown fields.
• Set a few benchmark problems against known analytical solutions (such as our spherical source or room model) to test various Direct (if the problem doesn’t have too many nodes), Iterative or Multigrid solvers, with different preconditioning, as well as different solver specific settings. We want:
  • Convergence in a reasonably small number of iterations.
  • Reasonably fast iterations.
  • High accuracy of our solution.
• Once we have figured out what solver parameters give us the best results, consistently among different benchmarks, then we can take the same parameters to more complex problems, knowing they will work alright (solver specific settings are probably going to be tightly related to particular problems, though).

So, now that we are aware of this, we will discuss our solver settings choices a little more depth in the next episodes.

License Information

This work is licensed under a Creative Commons Attribution 4.0 International License.

• Introductory
• Physics
• Elmer