NDSolve Options for Finite Elements

Overview

If you have never used the finite element method implemented in the Wolfram Language, this tutorial is probably not a good starting place. To get an overview of the finite element method, a first reading should be the tutorial Solving Partial Differential Equations with Finite Elements. This tutorial explains how to fine-tune details of the finite element method.

Differential equations come in many forms, such as ordinary differential equations, partial differential equations, nonlinear differential equations and time-dependent differential equations, to name a few. Many algorithms exist to solve these various kinds of differential equations and have various advantages and disadvantages that make them more or less suitable to solve a particular problem. What complicates matters further is that many of these algorithms have tuning parameters that make solving a particular differential equation faster, less memory intensive or even possible at all. In the Wolfram Language, the NDSolve family of functions collects some of these algorithms and provides an interface to them. Based on a symbolic analysis of the differential equation given, an automatic algorithm can take place. This works well in many cases, as Wolfram Research puts a lot of effort into the equation analysis and automatic algorithm choice to solve a particular differential equation. What is also done is to automatically find reasonable default values for tuning parameters.

While in general the automatic algorithm selection works well, there are scenarios where it is desirable to choose a different algorithm or modify tuning parameters.

The purpose of this tutorial is to give an overview of what algorithms and algorithm options are available for the finite element method implemented in NDSolve. The tutorial also explains which algorithms are automatically chosen and how to choose different algorithms. For each of the methods available, tuning parameters are presented and discussed or links to other sections of the finite element method documentation with more examples are given.

This tutorial speaks a lot about NDSolve. Everything mentioned here also holds true for the sister function NDSolveValue. The difference between NDSolve and NDSolveValue is in the form in which they return the computed solutions. While NDSolve will return a list of rules with InterpolatingFunction objects, NDSolveValue will directly return InterpolatingFunction objects.

Evaluating some parts of this tutorial may use longer run times and larger amounts of memory than usual for typical documentation examples.

NDSolve Options

In this section, general options of NDSolve are presented and put into perspective with the finite element method.

What follows is a list of the options and what they mean with regard to the finite element method.

AccuracyGoal

If the AccuracyGoal option is specified, its value will be propagated to all algorithms NDSolve uses. It is thus possible to have different values of AccuracyGoal for different algorithms. For example, a different AccuracyGoal for time integration and mesh generation may be used when specified via MeshOptions.

In the current version, the finite element method does not make use of automatic adaptive mesh generation. Nonuniform mesh cell sizes can be achieved through the specification of a MeshRefinementFunction.

Compiled

The Compiled option is used for various numerical functions and specifies whether the expressions they work with should automatically be compiled. The option does not have an effect for finite element computations, since the finite element method is available in machine precision and the method always automatically compiles PDE coefficients and boundary conditions as needed. For more information, see the section Efficient Evaluation of PDE Coefficients of the Finite Element Method Usage Tips tutorial.

DependentVariable

DependentVariable has the same meaning as outside of the finite element method context.

DiscreteVariable

DiscreteVariable has the same meaning as outside of the finite element method context.

EvaluationMonitor

For time-dependent PDEs, EvaluationMonitor has the same meaning as outside of the finite element method context. For stationary PDEs, the option has no meaning. If an EvaluationMonitor for a nonlinear stationary PDE is needed, then the EvaluationMonitor needs to be specified through options for FindRoot. This is explained in the section on FindRootOptions.

An example of how to use a monitor for monitoring the progress during a time integration for a PDE that uses the finite element method is given here.

ExtrapolationHandler

Interpolation functions in the Wolfram Language will extrapolate if they are queried outside of their domain. In a finite element analysis, this is not wanted, since the extrapolation will almost certainly not give the physically correct value. The interpolation functions NDSolve returns from a finite element analysis will return Indeterminate as an extrapolation value. This behavior can be switched off.

The following PDE is going to be solved over the solution domain.

The resulting interpolation function can be queried at points that do not belong to the solution domain. This will result in a warning message and Indeterminate as an extrapolated value.

To change the behavior, you can have the interpolation function return an extrapolation value and quiet the warning message.

A query point outside the solution domain now returns an extrapolation value and no warning message.

The "DefaultExtrapolationHandler" system option can be set to change the default behavior.

More examples can be found on the reference page for ElementMeshInterpolation.

InitialSeeding

An initial seeding may be needed for nonlinear stationary PDEs. InitialSeeding provides a mechanism to specify an initial guess for Newton's method to find a solution to the nonlinear PDE. The reference page of InitialSeeding has more details.

InterpolationOrder

The only value for the option InterpolationOrder besides Automatic with respect to NDSolve is InterpolationOrderAll. Typically, NDSolve stores the results from every timestep during a time integration. When InterpolationOrderAll is specified, then additional information is stored in the InterpolatingFunction and a different interpolation method may be used that allows for a more accurate interstep interpolation. This is sometimes referred to as dense output. The advantage is that the inter step accuracy is very high; however, since more information needs to be stored, the size of the InterpolatingFunction increases. An example can be seen on the reference page for InterpolationOrder.

For the stationary case, this option does not have an effect. The interpolation order of the InterpolatingFunction is then derived from the finite element mesh order.

MaxStepFraction

For time-dependent PDEs, MaxStepFraction has the same meaning as outside of the finite element method context. For stationary PDEs, the option has no meaning.

MaxSteps

For time-dependent PDEs, MaxSteps has the same meaning as outside of the finite element method context. For stationary PDEs, the option has no meaning.

MaxStepSize

For time-dependent PDEs, MaxStepSize has the same meaning as outside of the finite element method context. For stationary PDEs, the option has no meaning.

Method

The usage of the Method option of NDSolve is the main purpose of this document and will be shown in the following.

NormFunction

The option NormFunction has no effect with the finite element method.

PrecisionGoal

If a PrecisionGoal is specified, its value will be propagated to all algorithms NDSolve uses. It is thus possible to have different values of PrecisionGoal for different algorithms. For example, a different PrecisionGoal for time integration and mesh generation may be used when specified via MeshOptions.

In the current version, the finite element method does not make use of automatic adaptive mesh generation. Nonuniform mesh cell sizes can be achieved through the specification of a MeshRefinementFunction.

StartingStepSize

For time-dependent PDEs, StartingStepSize has the same meaning as in the no finite element case. For stationary PDEs, the option has no meaning.

StepMonitor

For time-dependent PDEs, StepMonitor has the same meaning as outside of the finite element method context. For stationary PDEs, the option has no meaning. If a StepMonitor for a nonlinear stationary PDE is needed, then the StepMonitor needs to be specified through options for FindRoot. This is explained in the section on FindRootOptions.

An example of how to use a monitor for monitoring the progress during a time integration for PDE that uses the finite element method is given here.

WorkingPrecision

WorkingPrecision has the same meaning as outside of the finite element method context.

The Method Option for Solution Stages

The finite element method (FEM) is a method to discretize the spatial components of partial differential equations. The FEM is a sister method of the tensor product grid method (TPG), which has the same purpose as the FEM but uses a different method to compute a result and has different advantages and disadvantages. The tensor product grid method is typically more accurate and solves a differential equation faster, but it cannot be used on arbitrary regions and for time-independent (stationary) equations, for which the finite element method can be used. The tensor product grid method is discussed in more detail in the context of the Method of Lines tutorial.

The choice of algorithms used by NDSolve can be influenced by setting method options in NDSolve. To understand how options for NDSolve are structured, it is instructive to know that NDSolve before Version 10 was mostly for solving time-dependent differential equations. Some time-independent problems could be solved by using methods like the shooting method, but in essence, the algorithms in NDSolve were for time integration of differential equations. In Version 10, the finite element method was added, and with it support to solve stationary partial differential equations.

NDSolve typically solves differential equations by going through several different stages, depending on the type of equation. With Method{s₁m₁,s₂m₂,…}, stage s_i is handled by method m_i. The actual stages used and their order are determined by NDSolve, based on the problem to solve.

Possible solution stages are:

All spatial discretization options for the finite element method are set up through the "PDEDiscretization" stage. The "BoundaryValues" stage is the only irrelevant option for the finite element method. All remaining solution stages are only relevant for time-dependent PDEs.

While a general overview of the methods and options NDSolve accepts is given in the reference page under the Details and Options and the Options section and also with more detail in the Advanced Numerical Differential Equation Solving in the Wolfram Language tutorial, this notebook explains the structure of the NDSolve options and the interaction with the finite element method.

For some options, shortcuts or shorter versions that come from still-supported legacy option specifications can be given. This notebook uses the full option structure. This has the advantage that the structure of the options is much clearer.

How to Tell What Method Has Been Used

Sometimes, it is of interest to know if the finite element method has been used as a solution method to solve a particular differential equation. The easiest method to do so is by inspecting the returned InterpolatingFunction. If the function contains an ElementMesh, then the finite element method was used. If no ElementMesh is present, then some other method has been used.

If the default algorithm selection of NDSolve did not choose the finite element method, the use of the FEM can be forced by specifying a Method option.

If no solution can been found by NDSolve, then there is typically an error message. Error messages that involve the finite element method typically start with NDSolve::femXYZ.

In common message cases, following the link that becomes visible once one clicks on the three red dots in front of the message icon leads to more information about the message and how to circumvent the message. Making use of the messaging system in the language is explained in the workflow Understand Error Messages.

What Triggers the Use of the Finite Element Method

NDSolve is designed in such a way that it chooses the method to solve the equation at hand automatically. This means that any single one of the following triggers will invoke the the use of the finite element method:

• Specifying a region with Element[x,dom]
• dom is a region (RegionQ[dom] is True)
• dom is an ElementMesh object
• dom is a NumericalRegion object

• The PDE is elliptic (stationary) and the embedding dimension is greater than 1

• The PDE is specified using Inactive

• The boundary conditions are specified using DirichletCondition, NeumannValue or PeriodicBoundaryCondition

• When too few initial conditions are given for a time-dependent PDE problem, the problem is solved as a stationary problem

If one desires to use the alternative tensor product grid method, then these triggers need to be eliminated from the call to NDSolve.

Solve a PDE specified in a way that automatically uses the finite element method:
Solve the PDE specified in a way that uses the tensor product grid method:

The finite element method as implemented is currently available in real- and complex-valued machine precision, so arbitrary-precision solutions to PDEs with the finite element method are currently not possible.

Finite Element Method Options for Stationary Partial Differential Equations

The following box shows the general structure for specifying options and suboptions for the finite element method in the stationary, that is, the time-independent case.

Use the finite element method to discretize a PDE.

The usage of Method{"FiniteElement"} is a shorthand for Method{"PDEDiscretization""FiniteElement"}. Be careful with this shorthand though, as for a time-dependent PDE, this shorthand will force the solution as a time-independent PDE.

The finite element method as implemented in NDSolve separates the solution process into smaller stages. The sequence of stages and the working of internal FEM functions are documented in the Finite Element Programming tutorial.

Some of the internal functions have their own set of options that can be specified within NDSolve by setting up suboptions for the finite element method.

Forcing the usage of the finite element method can lead to unexpected results if not done properly in the case of time-dependent PDEs. There is a subtle difference in forcing the usage of the finite element method for the entire PDE or if only the spatial component should be discretized with the finite element method. Often this second method is what is wanted.

The following box summarizes the low-level finite element functions available that have options.

Finite element method low-level functions with options.

The next sections discuss how the low-level finite element function options can be set up as suboptions from NDSolve and give links to relevant sections of the documentation that show more detailed usage of these options.

InitializePDECoefficients

Once the PDE is parsed, the suitability of the extracted PDE coefficients is verified. This includes tests that, for example, the coefficients are given in the proper dimensions or that proper dimensions can be deduced from the given input. The coefficients can be autocompiled. This functionality is provided by the function InitializePDECoefficients. All options that InitializePDECoefficients takes can also be specified on the NDSolve level.

InitializePDECoefficientOptions for the finite element method.

What follows is an example on how to specify options for InitializePDECoefficients. The verification process evaluates the coefficient at a point inside the simulation domain. Should this evaluation point be, for example, at a singularity of the coefficient, then the coefficient is rejected.

Set up a coefficient that has a singularity at :

For a one-dimensional mesh with bounds from 0 to 1, the verification coordinate happens to be at .

The coefficient is singular at the verification coordinate for this region:

Change the verification data such that the test coordinate is not at the singularity of :

More information can be found in the options section of InitializePDECoefficients.

Mesh Generation Options

Once the PDE coefficients and boundary conditions are initialized, the finite element mesh is generated. The mesh is a fundamental part of the finite element method. The quality of the mesh has a large influence on the accuracy of the solution of the equation. At the same time, using a mesh with too many elements will unnecessarily prolong and strain the solution process. All options that ToBoundaryMesh, ToElementMesh and ElementMesh take can also be specified on the NDSolve level. The mesh generation options are specified as a suboption of the "FiniteElement" method.

MeshOptions for the finite element method.

As a specific example of how to specify any option for the generation of an ElementMesh, the specification of a MaxCellMeasure is shown. All other mesh generation options are specified in the same manner.

Fine-tuning options for mesh generation are documented in the respective reference pages under the Details and Options and the Options section of ToBoundaryMesh, ToElementMesh and ElementMesh. Also, there is extensive documentation on finite element mesh generation, so mesh generation will not be further discussed here.

PDESolveOptions

Once the PDE coefficients, boundary conditions and mesh are set up, the PDE is discretized. This is the process where the continuous PDE is converted into a discrete system of equations. This system of equations needs to be solved, which is possible after the boundary conditions are deployed. The function that takes care of the discretization, boundary condition deployment and solution of the equations is PDESolve. All options that PDESolve takes can also be specified on the NDSolve level.

PDESolveOptions for the finite element method.

Two important suboptions for the PDE discretization are "LinearSolver" for all linear problems and additionally, "FindRootOptions" for nonlinear problems.

LinearSolver

Once the coefficients, boundary conditions and method data are initialized, PDESolve will assemble a possibly large sparse system of equations of the form . Here is the sparse stiffness matrix, is the right-hand side, the load vector, and is the solution sought. The details creating the sparse system of equations are covered in Finite Element Programming tutorial.

By default, the system of equations will be solved with LinearSolve. This solution process can be time and memory consuming. Generally speaking, linear solvers are either direct or iterative. Direct solvers are based on variations of Gaussian elimination and only depend weakly on the stiffness matrix , which makes direct solvers very robust. On the other hand, iterative solvers typically use less memory and can be faster with an appropriate PDE-dependent preconditioner but are less robust. The finite element method uses the efficient direct PARDISO solver as the default linear solver.

Options for LinearSolve are specified in the following manner:

"LinearSolver" options for the finite element method.

The "LinearSolver" options allow you to specify options for LinearSolve or to replace the system LinearSolve function with a custom linear solver. The defaults are given in the following.

Specify options for LinearSolve or a custom solver.

A discussion about solving large-scale finite element applications and how LinearSolve can be optimized for this is given in the Finite Element Method Usage Tips tutorial.

Direct linear solver

If a solution to a discretized PDE exists, a direct solver will find the solution up to machine precision. Direct solvers just work and this is a major advantage. The default linear solver for the finite element method is the "Pardiso" solver. This direct solver is optimized for speed and memory efficiency.

For coupled PDEs with complicated 3D geometries, a direct solver can require more than the computer's available RAM memory. Before trying to make use of an iterative solver, the default "Pardiso" solver allows for some portions of the solution process to be stored outside of the RAM memory. This process is called out-of-core solution. The usage of this option is shown in the section on solving large-scale equations in the Finite Element Method Usage Tips tutorial.

If the direct solver still runs out of memory, another option is to try to make use of an iterative solver.

Iterative linear solver

Iterative solvers are very tricky to use. It is only recommended to make use of iterative solvers for finite element applications when the direct solver needs more RAM memory than available. This mostly comes up in 3D models with coupled PDEs.

For demonstration purposes, however, a single PDE model in 3D is used.

Using the direct solution as a reference, you can estimate the relative error of the iterative solution. Note that the direct solution also has errors from discretization, but these are negligible compared to the iterative solution quality.

The iterative solver needs less memory to solve the PDE than the direct solver. The first difference to note is that iterative solvers need a minimal tolerance the solution should satisfy. A setting of Tolerance10^-3 is a good starting value. In this case, the iterative solver is also faster than the direct solver.

The complete list of iterative solvers available and their options can be found on the Linear Solve reference page. For the finite element method, probably BiCGStab and GMRES are the important algorithms, since they work on general matrices. Using the iterative solvers without a preconditioner is not recommended, as that will result in a poor convergence rate and thus a long wait for the solution to converge.

The convergence rate of the iterative solver depends on the condition number of the stiffness matrix . The lower the condition number, the better. To actually inspect the condition number, the techniques to extract the system matrices presented in the Finite Element Programming tutorial need to be used. The condition number is influenced by various factors, such as the geometry or the scale of the PDE coefficients. To make iterative solvers work, the condition number needs to be reduced as much as possible.

A higher-quality mesh can result in a lower condition number. The quality of a mesh can be improved by specifying the option "MeshQualityGoal"1 during the mesh generation process. For more details on this option, please refer to the section ElementMesh quality in the ElementMesh Generation tutorial.

Another way to proceed is to not reduce the condition number stiffness matrix but to provide LinearSolve with a matrix that is similar but has a lower condition number. This can be achieved by applying a preconditioner matrix to the stiffness matrix . Say and are preconditioner matrices which, when applied to the system , decrease the condition number of the resulting system of equations . The goal is to generate an with lower condition number than . By applying a preconditioner, you can solve one of these equations:

• for a left preconditioner

• and for a right preconditioner

• and for both-side preconditioners

The quest is to find a good preconditioner. The best preconditioner is . But since computing is effectively done during the direct solution process, and that process uses too much memory, an approximation to needs to be found. Currently, a good preconditioner for all kinds of stiffness matrices is not known.

However, an ILU preconditioner is often used as a general-purpose preconditioner. Here ILU stands for incomplete LU (decomposition). Three variants if the ILU preconditioners are implemented:

	ILUO	incomplete LU with no fill in
	ILUT	threshold ILU
	ILUTP	threshold ILU with column permutation

The ILU0 preconditioner is the simplest one among these three preconditioners. Its main advantage is that it does not introduce any additional nonzero elements (= fill in) into the ILU0 matrix. So the original matrix and the ILU0 matrix have the same sparsity structure.

The functioning of the ILUT preconditioner is more complex than that of the ILU0 preconditioner. ILUT uses two parameters: drop tolerance and fill in. It applies a drop strategy for each number in a row as it is computed. If the absolute value of the number is less than the drop tolerance, then that element in the preconditioner matrix is set to zero. Additionally, only a fixed number of the largest entries are kept at all; the rest are set to zero and thus dropped.

Use the ILUT preconditioner with a drop tolerance and fill-in specification:
Estimate the difference in relative solution error:

In this case, the solution quality could be improved with roughly the same computational time and memory requirements.

One useful iterative solver option is the "StartingVector". What can be done is to use a direct solver with a mesh that has fewer elements or uses a mesh order of 1, and use the result of this simulation as a starting point for a simulation on a refined mesh. This procedure is explained in the solving large-scale problems section of the Finite Element Method Usage Tips tutorial.

Making use of an iterative solver for nonlinear PDEs is disadvantageous for two reasons. First, a direct solver can store its decomposition and reuse that in the nonlinear solution process, which is explained in the FindRoot method tutorial on the affine covariant Newton method. Using an iterative solver that cannot store a decomposition will most likely be slower. A second reason why solving a nonlinear PDE with an iterative solver can be problematic is that when the solution process fails, it is much harder to tell the cause of the failure. Did the solution process fail because the Newton method was not able to find a solution, say because the PDE does not have a solution in the first place, or did the solution process fail because of a poor quality result during an iterative solve step? If an iterative solver has to be used for a nonlinear PDE, then a higher tolerance of Tolerance10^-6 may be necessary.

In general, it is difficult to give general advice on the use of iterative solvers, since their performance depends on the PDE and geometric model used.

Custom linear solver

It is possible to replace the system function call to LinearSolve with a custom linear solver.

The solution process for nonlinear PDEs also makes use of a linear solver. The exact details of how the nonlinear solver works are explained in the finite element programming notebook in the section about nonlinear PDEs. A difference compared to the linear case is, however, that a factorization of the system matrix is done and then applied to different right-hand sides. To reflect this need, the preceding wrapper code is extended to allow for a one-argument form like LinearSolve. In the nonlinear PDE case, the custom linear solve needs to return a function like a LinearSolveFunction that can be applied to a right-hand side during the solution process.

With the custom linear solver code, you have direct access to the system matrices right before the solution process.

In order to pass options to the custom linear solver, these options need to be set up.

Options that are now specified behind the custom linear solver specification will be passed on to the custom linear solver.

FindRootOptions

Options for FindRoot are only applicable for nonlinear time-independent partial differential equations.

FindRootOptions for the finite element method.

The default FindRoot method for the finite element method is the affine covariant Newton method that has options and use cases documented in the notebook The Affine Covariant Newton Method. To provide an overview here of how options for FindRoot can be set, a nonlinear PDE is solved.

As a specific example, monitors are used to monitor step, function and Jacobian evaluation. All other options mentioned as explained in the notebook The Affine Covariant Newton Method can be specified in the same manner.

The number of function evaluations, steps and Jacobian evaluations is important to be able to monitor, as that will show if other options given have an effect on the efficiency with which nonlinear PDEs can be solved.

When the discretization used is the finite element method, there are no known cases where it is beneficial to use a different method for FindRoot than the default affine covariant Newton method.

As seen in this case, the number of Jacobian evaluations is larger with the standard Newton method compared to the affine covariant Newton method. Because the evaluation of the Jacobian matrix is very expensive for the finite element method, the affine covariant Newton method tries to minimize the number of Jacobian evaluations.

All root-finding algorithms based on Newton's method need a starting value. The default for the finite element method is 0. This may not always be the correct choice. The starting value may be modified using the InitialSeeding option of NDSolve. The scope section has an example showing how the solution of a linear version of a differential equation can be used solve a nonlinear version of the same equation.

For more fine-tuning options available and examples of their use, refer to the method's documentation given in The Affine Covariant Newton Method. The theory and implementation of the linearization used in the finite element method is discussed in the section Nonlinear PDEs of the Finite Element Method Programming tutorial.

Remaining "FiniteElement" Method Options

AccuracyGoal

Specify a different AccuracyGoal for the finite element method algorithms than for the rest of NDSolve.

BoundaryTolerance

Setting the boundary tolerance to None will enable inconsistent DirichletCondition detection:

In the preceding case, the boundary conditions are inconsistent at two coordinates.

The default "BoundaryTolerance" is Automatic and does not warn about inconsistent boundary conditions but will eliminate duplicate boundary conditions that may come up during the inconsistent boundary condition analysis. Setting "BoundaryTolerance" to Infinity will not perform any inconsistency check.

IntegrationOrder

IntegrationOrder is a submethod option for InitializePDEMethodData. The IntegrationOrder is the order of accuracy used to integrate the finite element operators. In the case where an ElementMesh with "MeshOrder" 2 is given, the Automatic setting chooses an integration order of 4. For an ElementMesh with "MeshOrder" 1, an integration order of 2 is chosen. The maximum order is 5.

Sometimes, it is of interest to know if the finite element method has been used as a solution method to solve a particular differential equation. The easiest method to do so is by inspecting the returned InterpolatingFunction. If the function contains an ElementMesh, then the finite element method was used. If no ElementMesh is present, then some other method has been used.

InterpolationOrder

With an ElementMesh having "MeshOrder" , the InterpolationOrder can be at most and "InterpolationOrder"Automatic defaults to the order . Multiple dependent variables may have a different interpolation order, but at least one interpolation order needs to be set to the maximum mesh order . Having multiple different interpolation orders is equivalent to specifying what other finite element software know as P2/P1 or Q2/Q1 elements. This can be useful for fluid dynamics simulations and is shown in more detail in the section Fluid Flow in the Solving Partial Differential Equations with Finite Elements tutorial.

LinearSolveMethod

This is a legacy option and should not be used any longer. Instead, specify options for LinearSolve through PDESolveOptions.

PrecisionGoal

Specify a different PrecisionGoal for the finite element method algorithms than for the rest of NDSolve.

PrecomputeGeometryData

NDSolve precomputes and stores data that is based on the geometry. Switching this off is never beneficial on the NDSolve level. Only low-level finite element programming functions may benefit from not precomputing finite element data. PrecomputeGeometryData is a submethod option for InitializePDEMethodData.