Finite Element Method Usage Tips

Introduction

The aim of this tutorial is to point out possible issues when using the finite element method with NDSolve and offer best practices to avoid potential issues.

Load the finite element package.
In[1]:=
Click for copyable input

Monitoring Progress of Time Integration of Transient Partial Differential Equations

Time integration of partial differential equations can take time. This is especially true if the boundary conditions or coefficients are time dependent or if the PDE region at hand is three dimensional. An easy way to monitor progress of an integration of a transient PDE is by using an EvaluationMonitor.

Using the EvaluationMonitor within NDSolve works as in other cases. Here is an example wave equation with a time-dependent Dirichlet condition.

Specify an ellipsoidal region with a circular hole in one focal point.
In[1]:=
Click for copyable input
Create and visualize a mesh of the region.
In[4]:=
Click for copyable input
Out[5]=
Specify a transient Dirichlet boundary condition on the circular hole that is active until time .
In[6]:=
Click for copyable input
Set up the notebook's status area for displaying the current integration time.
In[8]:=
Click for copyable input
Time-integrate the PDE with the status displayed in the status area.
In[10]:=
Click for copyable input
Out[10]=

Once the initialization of the PDE is done, the status area will display the current time during the time integration.

Cleanup of the status area.
In[11]:=
Click for copyable input

Solving Memory-Intensive PDEs

The amount of memory needed by NDSolve to solve a particular problem with the finite element method is proportional to the number of equations and the space dimension in which you are operating. For coupled systems in 3D, the amount of memory provided by the underlying hardware can be a barrier. The following section identifies easy-to-use options to alleviate this potential barrier.

During the finite element analysis there are two key memory bottlenecks. The first is the creation of the discretization and the second is the solution (LinearSolve) of the system of equations. NDSolve provides options to both mesh generation and the LinearSolve step that have an effect on the memory requirement during discretization and solving. As a rough estimate, the number of coordinates in the internal generated mesh multiplied by the spatial dimension gives the degrees of freedom (number of rows and columns) of the discretized system matrix to be generated.

As an example, consider a stress operator from elastostatics. The stress operator is a system of three coupled PDEs in 3D space. Since the main focus of this section is to create a large system of equations, no additional information about the PDE itself is presented.

Define the stress operator, a Young's modulus, and a Poisson ratio.
In[14]:=
Click for copyable input

A series of measurements illustrates how various options influence the solution time and the amount of memory needed to find a solution. For this purpose, a helper function is implemented. The function receives as input any number of NDSolve options and returns the time and the amount of memory in MB needed to solve the PDE. It is also instructive to get a feeling for the size of the system matrices NDSolve internally establishes. The size of the system matrices is its degrees of freedom, or in other words, the number of rows and columns in the system matrices.

Create a helper function that calls NDSolve with options. The measured time, memory needed, and degrees of freedom to find the solution are printed.
In[15]:=
Click for copyable input

The default settings in NDSolve are set up such that the solution is accurate and found fairly efficiently. To make this example applicable for a larger class of computers, tetrahedron elements are used. The use of these results in smaller system matrices.

Calling NDSolve with settings that enforce tetrahedron elements in the mesh.
In[16]:=
Click for copyable input

With the default setting, NDSolve will call LinearSolve with no special options. This will route to a direct solver that is very accurate but has a high requirement of available memory.

Visualize the deformed structure.
In[17]:=
Click for copyable input
Out[19]=

One way to reduce the memory consumption is to have NDSolve generate a mesh with fewer elements.

Instruct NDSolve to use a smaller mesh.
In[20]:=
Click for copyable input

This has reduced both solution time and memory requirements at the price of accuracy. How many elements a mesh needs to represent an accurate solution is problem dependent. A good method is to start with a coarse mesh and refine until the solutions stop differing substantially.

Visualize the difference of the computed results at a cross section.
In[21]:=
Click for copyable input
Out[22]=
Visualize the difference of node distances as an Image3D rescaled to the range of 0 to 0.01.
In[23]:=
Click for copyable input
Out[25]=

A way to further reduce the amount of memory needed for computing the solution is to use a mesh that is only first-order accurate.

Instruct NDSolve to use a first-order-accurate mesh.
In[26]:=
Click for copyable input
Visualize the difference of the computed results at a cross section.
In[27]:=
Click for copyable input
Out[28]=
Visualize the difference of node distances as an Image3D rescaled to the range of 0 to 0.01.
In[29]:=
Click for copyable input
Out[31]=

This has reduced both solution time and memory requirements at the price of accuracy. Making use of a first-order mesh results in much smaller system matrices, which can be solved more efficiently. Using a first-order mesh may significantly reduce accuracy.

The previous suggestions worked by reducing the size of the system matrices internally generated. A different approach is to instruct NDSolve to use a different solver method for LinearSolve. All options that LinearSolve take can be directly given to NDSolve.

Instruct NDSolve to use the "Pardiso" solver.
In[32]:=
Click for copyable input
Visualize the difference of the computed results at a cross section.
In[33]:=
Click for copyable input
Out[34]=
Visualize the difference of node distances as an Image3D rescaled to the range of 0 to 0.01.
In[35]:=
Click for copyable input
Out[37]=

While this takes longer to solve, it uses substantially less memory while at the same time retaining the accuracy, compared to the default solution. The reason that this method is not the default is that the default LinearSolve method is guaranteed to converge if the discretized PDE has a solution. This benevolent property is not true for all solvers implemented in LinearSolve.

The options discussed above can also be combined. A more in-depth treatise of options that can be given to NDSolve can be found in the finite element tutorials.

Extrapolation of Solution Domains

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.

This shows the solution domain.
In[38]:=
Click for copyable input
Out[40]=

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

In[41]:=
Click for copyable input
This solves the equation.
In[42]:=
Click for copyable input

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.

In[43]:=
Click for copyable input
Out[43]=

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

This solves the equation but gives extrapolation values for query points outside the solution domain.
In[44]:=
Click for copyable input

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

In[45]:=
Click for copyable input
Out[45]=

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

In[46]:=
Click for copyable input
Out[46]=

Stabilization of Convection-Dominated Equations

Consider a convection-diffusion type of equation.

In a situation where the convective component becomes large, the solution of the PDE may become unstable, as shown in the following example.

Specify a convection-diffusion equation operator.
In[47]:=
Click for copyable input
Out[48]=
Specify boundary conditions.
In[49]:=
Click for copyable input
Specify a mesh.
In[50]:=
Click for copyable input
Out[50]=
Solve a convection-diffusion equation with a large convective component.
In[51]:=
Click for copyable input
Out[51]=

The message indicates that the solution may not be stable.

Visualize the solution of a convection-diffusion equation with a large convective component.
In[52]:=
Click for copyable input
Out[52]=

The result is disastrous. There are a few things that can be done in such situations. First, however, investigate how the Peclet number is computed.

Compute the minimal diameter of the mesh elements.
In[53]:=
Click for copyable input
Out[53]=
Compute the Peclet number from the convection and diffusion coefficients and the minimal element diameter.
In[54]:=
Click for copyable input
Out[54]=

One way to improve the stability of the solution is to add so-called artificial diffusion. This artificial diffusion is then added to the diffusive term in order to stabilize the PDE.

Estimate the amount of artificial diffusion needed.
In[55]:=
Click for copyable input
Out[57]=
Verify that the addition of the artificial diffusion would return a stable result.
In[58]:=
Click for copyable input
Out[58]=

Note that the newly computed Peclet number is smaller than the mesh order.

Solve a convection-diffusion equation with a large convective component.
In[59]:=
Click for copyable input
Visualize the solution of a convection-diffusion equation with a large convective component.
In[60]:=
Click for copyable input
Out[60]=
Visualize a cross section of the solution.
In[61]:=
Click for copyable input
Out[61]=

This is already much better. The disadvantage with this method is that the PDE is modified. This may not always be a good approach. As an alternative, the mesh can be modified. For more information about mesh generation, please consult "Element Mesh Generation".

Generate and visualize a mesh with a refined boundary.
In[62]:=
Click for copyable input
Out[63]=
Solve a convection-diffusion equation with a large convective component on a predefined mesh.
In[64]:=
Click for copyable input
Visualize the solution of a convection-diffusion equation with a large convective component.
In[65]:=
Click for copyable input
Out[65]=
Visualize a cross section of the solution.
In[66]:=
Click for copyable input
Out[66]=

The solution is improved compared to the previous example, at the cost of some additional computational time. In principal, it would be sufficient to only refine the right and top boundaries of the mesh. This would then decrease the computational time.

The Coefficient Form of a PDE

In this example, we are going to demonstrate the importance of the coefficient form of PDEs for the finite element method. The coefficient form for a single dependent variable PDE is given in equation (1).

The PDE is defined in . Here is the dependent variable for which a solution is sought. The coefficients , , , and are scalars; , , and are vectors; and is an × matrix.

The key point of this example is that this equation (and coupled versions of it) is in essence the only equation the finite element method can solve. To emphasize this point, consider the simpler equation (2).

We set up the equation and its coefficients. As a spatial domain for the PDE, we use , and the time integration will be set up to stop at . We set the diffusion coefficient to 1, which will internally get transformed to the identity matrix. For the damping coefficient , we use a piecewise continuous function that has the value of 1 for values of , and 2 otherwise.

Specify the coefficient and set the equation (3) up.
In[28]:=
Click for copyable input
Solve the equation.
In[30]:=
Click for copyable input
Out[30]=
Plot the solution at .
In[31]:=
Click for copyable input
Out[31]=

One could be tempted to reformulate equation (4) in the following manner:

Set up equation (5).
In[32]:=
Click for copyable input
Solve the equation.
In[33]:=
Click for copyable input
Out[33]=
Plot the solution at .
In[34]:=
Click for copyable input
Out[34]=
Plot the difference of solution 1 and solution 2.
In[35]:=
Click for copyable input
Out[35]=

Clearly, solution 1 and solution 2 differ. To investigate what happens, we write a little helper function that extracts the actually parsed PDE coefficients from NDSolve's data structures. For more information on the internal workings of the finite element method, please consult the Finite Element Programming tutorial.

Specify a helper function to extract PDE coefficients.
In[36]:=
Click for copyable input
Inspect the PDE coefficients from equation (6).
In[37]:=
Click for copyable input
Out[37]=

The coefficients for equation (7) look as expected.

Inspect the PDE coefficients from equation (8).
In[39]:=
Click for copyable input
Out[39]=

The coefficients for equation (9) come as a surprise. The coefficient got pulled into the diffusion coefficient . Why is that? The key is to understand that the finite element method can only solve this type of equation.

Note that there is no coefficient in front of the term . To get equations of the type to work, is reset to and is adjusted to get rid of the derivative caused by . Here is an example.

In[59]:=
Click for copyable input
Out[61]=

Following the same procedure for the equation (10), we get:

In[81]:=
Click for copyable input
Out[81]=

This is the effectively same as explicitly specifying . Mathematically, because of the discontinuity, the coefficient β is like a Dirac delta function, which would be omitted in the time integration anyway.

Click for copyable input
Inspect the PDE coefficients.
In[82]:=
Click for copyable input
Out[82]=

Formal Partial Differential Equations

Formal partial differential equations are equations that have Inactive components. This section explains when PDE need to be inactive. To start with an example, a plane stress operator is defined. For Young's modulus and Poisson's ratio, values for steel are used.

Define a plane stress operator.
In[67]:=
Click for copyable input

The first eigenvalue and eigenfunction of a beam held fixed at the left with length 10 meters and height 0.7 meters are computed.

Compute the eigenfunctions for a beam.
In[69]:=
Click for copyable input
Compute resonant frequency from the eigenvalue.
In[71]:=
Click for copyable input
Out[71]=

As a means of comparison, the first resonance frequency utilizing EulerBernoulli beam theory is computed analytically.

In[72]:=
Click for copyable input
Out[73]=

Compare the difference in results; about half the difference is attributable to the different theoretical approaches.

In[74]:=
Click for copyable input
Out[74]=

This is expected.

Next, the Active version of the PDE is inspected and the resulting resonance frequency is computed.

In[75]:=
Click for copyable input
Out[77]=

The above result is wrong. What happened? In short, the non-inactive version cannot account for unsymmetric coefficients. Looking at the original inactive PDE of the four coefficient matrices, two, the off-diagonal ones, are not symmetric.

In[78]:=
Click for copyable input
Out[78]=

In other words, even before NDEigensystem (and other related functions like NDSolve) can parse the PDE in its intended form, it will have evaluated.

In[79]:=
Click for copyable input
Out[79]=
In[80]:=
Click for copyable input
Out[80]=
In[81]:=
Click for copyable input
Out[81]=

From the evaluated PDE, there is no way to recover the antisymmetric coefficients. For this reason Inactive has been introduced. For antisymmetric coefficients, the Inactive version of the PDE is a must.

To inspect the coefficients that have been parsed, NDSolve`ProcessEquations can be used.

In[82]:=
Click for copyable input
Out[83]=

NeumannValue and Formal Partial Differential Equations

This section explains the interaction between formal partial differential equations and Neumann values. In the stationary case, a partial differential equation can be expressed as follows:

On the other hand, generalized Neumann boundary values specify a value . The value prescribes a flux over the outward normal on some part of the boundary:

The derivation of the Neumann boundary value is outlined in Partial Differential Equations and Boundary Conditions. It is important to note that the coefficients , , and in equation (11) are not specified directly by the use of NeumannValue. The coefficients , , and are specified in equation (12). The use of NeumannValue will replace with the value of on the boundary.

Now, the usage of Div and Grad alone is not enough to fully specify equation (13). Consider the specification of

and its corresponding Neumann value:

Consider a case where is given as , is the identity matrix and is a zero vector. Note that the evaluation of the divergence results in

where is and is the divergence of , in this example. Equations (14) and (15) are equivalent.

The evaluation poses a problem, however. Even though the equations are equivalent, the Neumann specifications are not. The corresponding Neumann value for equation (16) has implicitly changed to

To illustrate this further, a helper function is created that parses a PDE and extracts the PDE coefficient data from the finite element data, as explained in "Finite Element Programming".

Create a helper function that generates the NDSolve finite element data and extracts the PDE coefficient data.
In[84]:=
Click for copyable input
Set up the PDE and extract the PDE coefficient data.
In[85]:=
Click for copyable input
Out[87]=

PDECoefficientData is a data structure that contains the parsed PDE coefficients.

The diffusion coefficient is a matrix with .
In[88]:=
Click for copyable input
Out[88]=
The convection coefficient is .
In[89]:=
Click for copyable input
Out[89]=
The reaction coefficient has the value of .
In[90]:=
Click for copyable input
Out[90]=
No conservative convection coefficient is present.
In[91]:=
Click for copyable input
Out[91]=

During evaluation, the coefficient split into a and term. The proper mechanism to input differential equations where coefficients are evaluated prematurely is by either Inactivate or directly through the use of Inactive.

Inactivate a PDE.
In[93]:=
Click for copyable input
Out[93]=
Extract the PDE coefficient data.
In[94]:=
Click for copyable input
Out[94]=
The diffusion coefficient is .
In[95]:=
Click for copyable input
Out[95]=
The conservative convection coefficient is .
In[96]:=
Click for copyable input
Out[96]=
No convection coefficient is present.
In[97]:=
Click for copyable input
Out[97]=
No reaction coefficient is present.
In[98]:=
Click for copyable input
Out[98]=

If a γ term is present one might use the following inactive PDE:

Inactivate a PDE with a term.
In[99]:=
Click for copyable input
Out[100]=
Extract the PDE coefficient data.
In[101]:=
Click for copyable input
Out[101]=
The load derivative coefficient is .
In[102]:=
Click for copyable input
Out[102]=

The following example illustrates the difference in meaning of an active and inactive PDE with respect to Neumann values, from a practical perspective. The example shows a diffusion confined in a region subject to a biasing potential .

Specify a rectangular region with a hole inside.
In[103]:=
Click for copyable input
Specify a biasing potential and its gradient .
In[104]:=
Click for copyable input

Dirichlet boundary conditions are set up at the top and bottom () of the region with a potential of 1 and 0, respectively.

Specify Dirichlet boundary conditions at .
In[106]:=
Click for copyable input

At the remaining two sides () natural boundary conditions (Neumann zero) apply. This means zero flux across these two edges.

The PDE to be modeled is given by:

and its corresponding Neumann value at is

Specify the inactive PDE.
In[107]:=
Click for copyable input
Out[107]=
Solve the inactive PDE.
In[115]:=
Click for copyable input
Inspect the PDE coefficients.
In[116]:=
Click for copyable input
Out[117]=

First the satisfaction of the Neumann values at is inspected.

Plotting the normal derivative of the inactive operator over the edges with the Neumann zero value implicitly specified. The normal derivative is in the direction of .
In[118]:=
Click for copyable input
Out[118]=

With zero flux boundary conditions at (no flux across these two edges) and no source terms, the flux into the box at should be the same as the flux out of the box at .

Plotting the normal derivative of the inactive operator over the edges with the Dirichlet condition. The normal derivative is in the direction of .
In[119]:=
Click for copyable input
Out[119]=

The flux into the box is roughly the same as the flux out of the box. As a comparison, the activated PDE is solved.

Solve the activated PDE.
In[120]:=
Click for copyable input
Investigate the PDE coefficients.
In[121]:=
Click for copyable input
Out[122]=

The coefficient is parsed as a convection and a reaction coefficient.

Investigate the convection and reaction coefficients.
In[123]:=
Click for copyable input
Out[123]=

In the activated case, the boundary condition computed is .

Plotting the normal derivative of the activated operator over the edges with the Dirichlet condition. The normal derivative is in the direction of .
In[124]:=
Click for copyable input
Out[124]=

The Relation between NeumannValue and Boundary Derivatives

The following example explains the relation of a NeumannValue and boundary derivatives, and is set up such that the parameters are easy to modify. First, the region is specified.

Specify a region.
In[14]:=
Click for copyable input
Out[15]=

Here is the PDE to be solved.

The PDE has one dependent variable . Note that since this is a one-dimensional system, the diffusion coefficient is a 1×1 matrix, which can be specified as a number.

Set up the PDE operator.
In[4]:=
Click for copyable input
Out[5]=

Dirichlet boundary conditions specify a value for the dependent variable .

Neumann values are specified by the following form.

Next, the boundary conditions are set up.

Specify boundary conditions for the equation.
In[6]:=
Click for copyable input
In[8]:=
Click for copyable input
Solve the system of equations.
In[11]:=
Click for copyable input
Out[11]=
Compute the analytical solution using DSolveValue.
In[12]:=
Click for copyable input
Out[12]=

Note how in the DSolveValue formulation of the equation, the generalized Neumann boundary equation has to be given explicitly, and how the coefficients match those in the PDE.

To compare the analytical result to the numerical result, make a plot of the difference between the two.

Plot the difference between the analytical and numerical solutions.
In[13]:=
Click for copyable input
Out[13]=

Ordering of Dependent Variable Names

In some cases, when analyzing coupled systems of PDEs, the finite element method as implemented in NDSolve may return what seem to be unexpected results. This section explains when this happens and presents solutions to avoid the problem. In the last section, a more in-depth explanation is presented.

First a system of equations with appropriate boundary conditions and a region are set up.

Specify a region.
In[125]:=
Click for copyable input
Specify a system of three equations.
In[126]:=
Click for copyable input
Specify accompanying boundary conditions.
In[127]:=
Click for copyable input
Specify the PDE system.
In[128]:=
Click for copyable input

The dependent variables for this example are u, v and w. The system is solved and a visualization of the dependent variable w is given.

Solve the system and visualize the dependent variable w.
In[129]:=
Click for copyable input
Out[130]=

In a second experiment, the dependent variable name w is changed to p and the system of equations is solved yet again.

Solve the system with dependent w replaced with p and visualize a dependent variable p.
In[131]:=
Click for copyable input
Out[132]=

The result is unexpected. Before explaining in detail what happens, a way to enforce the expected behavior is presented.

NDSolve offers a finite element method option to specify an interpolation order that should be used for the dependent variables during the discretization of the PDE. The specification of the interpolation order will implicitly enforce the same ordering of the equations and the dependent variables.

Specify the interpolation order of the dependent variables and plot the solution of p.
In[135]:=
Click for copyable input
Out[136]=

Generally, specifying the same interpolation order for each dependent variable is fine. In this specific case, the equations are fluid flow equations (Stoke's flow), and stabilization can be enforced if the pressure variable's Interpolation order is one less than the interpolation order for the velocities u and v.

To investigate this behavior further, an NDSolve state object is created.

Create an NDSolve state object.
In[137]:=
Click for copyable input

Note that the ordering of the dependent variables in the FEMMethodData in the FiniteElementData is {p,u,v}.

Inspect the dependent variables in the state object.
In[138]:=
Click for copyable input
Out[138]=

Inspecting the diffusion and convection coefficients of the PDECoefficientData reveals the structure of the coupled PDE coefficients.

Display the coupled PDEs' diffusion coefficients.
In[139]:=
Click for copyable input
Out[139]//MatrixForm=
Display the coupled PDEs' convection coefficients.
In[140]:=
Click for copyable input
Out[140]//MatrixForm=

When the finite element interpolation order method option is set, the coefficient structure is adjusted to the canonical order of the dependent variables.

Create a NDSolve state object.
In[141]:=
Click for copyable input
Display the coupled PDEs' diffusion coefficients.
In[142]:=
Click for copyable input
Out[142]//MatrixForm=
Display the coupled PDEs' convection coefficients.
In[143]:=
Click for copyable input
Out[143]//MatrixForm=

In fact, reordering the equations' input to match the canonical order of the dependent variables also gives the expected result.

Solve the coupled PDE with different equation ordering.
In[144]:=
Click for copyable input
Out[145]=

A remaining question is, why can this reordering not happen automatically?

The inherent challenge is that the boundary conditions specified can only be applied to diagonal components of the coefficients. The reason is that the way the boundary conditions are entered causes the information about how a specific boundary condition relates to a specific equation in the PDE to be lost.

Inspect the internal data of the parsed boundary conditions.
In[146]:=
Click for copyable input
Out[147]=

Here each sublist has two parts. The first part indicates the row and column in which the second partthe boundary conditionwill be applied during the discretization.

Another option is to associate the Dirichlet boundary conditions with the respective equation to which they are related.

Create a system of coupled PDEs that have Dirichlet boundary conditions associated with them.
In[148]:=
Click for copyable input
Solve the coupled PDE and visualize part of the solution.
In[149]:=
Click for copyable input
Out[150]=

In summary, the best option is to specify the interpolation order.

Specify the interpolation order of the dependent variables and plot the solution of p.
In[151]:=
Click for copyable input
Out[152]=

Related Tutorials