Finite Element Programming

Introduction

NDSolve provides a high-level, one-step interface for solving partial differential equations with the finite element method. However, you may want to control the steps of the solution process with more detail. The package provides a lower-level interface that gives extensive control for each part of the solution process.

To use the finite element functions, the package needs to be loaded.
In[1]:=
Click for copyable input

The low-level functions in the package may be used for a variety of purposes:

  • to better understand what NDSolve does internally and how NDSolve finds solutions
  • to better understand what options are available, what their usage is, and when they are beneficial
  • to intercept the solution process at various stages and provide access to intermediate data
  • to enable development of specific, finite-element-based solvers, not only to solve PDEs but also other areas of numerics
  • to use NDSolve as an equation preprocessor

Finite Element Data within NDSolve

The NDSolve framework has lower-level functions , , and and a data object that provide some access to the solution process for all types of differential equations. The general use of these functions is described in "Components and Data Structures".

Consistent with the rest of NDSolve, the finite element method has its data accessible through . The following section explores how to access the finite element data from the object. The finite element data is stored in an object and accessed through the property.

Set up the object.
In[3]:=
Click for copyable input
Out[3]=

When there is no time dependence in the problem, as in this case, the display form of the object will indicate this by displaying . For transient equations, the time interval over which it is to be integrated will be shown.

The finite element data can be extracted from the state object.

Get the finite element data.
In[4]:=
Click for copyable input
Out[4]=

All data objects generated during the solution process are documented, for example, FiniteElementData, as are their constructor functions. This tutorial will explain the usage of these functions.

The display form for a FiniteElementData object shows the degrees of freedom of the discretized partial differential equation (PDE). The degrees of freedom are typically the number of nodes in the mesh times the number of dependent variables. The degrees of freedom indicate the number of rows and columns of the matrices that are to be assembled.

For a steady state problem, invoking finds the solution of a system of equations with LinearSolve.

Compute the system solution.
In[5]:=
Click for copyable input

The solution is then stored in the finite element data object.

Extract the solution values.
In[6]:=
Click for copyable input
Out[6]//Short=

The solution is represented as an ×1 matrix, where is the degrees of freedom.

Extract the solution as an InterpolatingFunction and plot the result.
In[7]:=
Click for copyable input
Out[7]=
In[8]:=
Click for copyable input
Out[8]=

Passing Finite Element Options to NDSolve

In the subsequent sections, the finite element functions internal to NDSolve will be discussed. Some of these functions have options. It is always possible to pass these options to NDSolve without making use of internal finite element functions. These are given in the Method option to NDSolve.

NDSolve[,Method{"PDEDiscretization"{"FiniteElement", femopts}}]use the finite element method with options femopts for a steady-state problem
NDSolve[,Method{"PDEDiscretization"{"MethodOfLines", "SpatialDiscretization"->{"FiniteElement", femopts}}}]use the finite element method with options femopts for the spatial discretization of a transient problem

Specifying finite element method options.

Call NDSolve with finite element options specified.
In[9]:=
Click for copyable input
Out[9]=

The same is possible for transient PDEs.

Solve a transient PDE with NDSolve and finite element options specified.
In[10]:=
Click for copyable input
Out[10]=

A Workflow Overview

The following sections show how NDSolve solves a finite element model step by step. NDSolve essentially uses the same functions as seen in these sections. Furthermore, the content of the object is created through the functions introduced below.

The solution is found in three stages:

  • initialization,
  • discretization
  • solving

During the initialization stage, the finite element method data is set up and stored in the FEMMethodData object. The PDE and boundary conditions are analyzed and classified into different components, and these results are stored in the PDECoefficientData and BoundaryConditionData objects, respectively.

InitializePDECoefficientsinitialize partial differential equation coefficients
InitializeBoundaryConditionsinitialize boundary conditions
InitializePDEMethodDatainitialize partial differential equation method data

Initialization functions.

In the second stage, the PDE and boundary conditions are discretized and stored in DiscretizedPDEData and DiscretizedBoundaryConditionData. "Discretized" essentially means that the continuous PDE and boundary conditions are approximated by discrete versions represented by so-called system matrices.

DiscretizePDEdiscretize initialized partial differential equations
DiscretizeBoundaryConditionsdiscretize initialized boundary conditions

Discretization functions.

The last step is then to merge the discretized PDE with the discretized boundary conditions (DeployBoundaryConditions) and to use LinearSolve to solve the equations.

DeployBoundaryConditionsdeploy discretized boundary conditions into discretized partial differential equations
LinearSolvesolve linear systems of equations

Solution stage functions.

Each of these stages will be addressed. Before that, however, a model equation is set up.

The Partial Differential Equation Problem Setup

To make a PDE susceptible to be solved by a numerical method such as the finite element method, three components are needed:

  • a partial differential equation (PDE)
  • a region
  • boundary conditions

The following setup will be used as a model PDE in the subsequent explanations.

Specify a partial differential equation operator.
In[11]:=
Click for copyable input
Out[11]=
Create a region.
In[12]:=
Click for copyable input
Out[12]=
Set up boundary conditions at the left- and right-hand sides of the domain.
In[13]:=
Click for copyable input
Out[13]=
In[14]:=
Click for copyable input
Out[14]=

To solve the PDE over the region with the given boundary conditions, NDSolve can be used directly.

Use NDSolve to solve the PDE.
In[15]:=
Click for copyable input
Out[15]=

This is a typical workflow that uses the one-step capability of NDSolve.

To compare with using the functions to solve step by step, start by extracting and inspecting the data used by NDSolve.

Extract the object.
In[16]:=
Click for copyable input
Out[16]=
Get the finite element data.
In[17]:=
Click for copyable input
Out[17]=
Inspect the properties.
In[18]:=
Click for copyable input
Out[18]=

The PDECoefficientData has been created by a call to InitializePDECoefficients. The BoundaryConditionData has been created by a call to InitializeBoundaryConditions, and the FEMMethodData has been created by a call to InitializePDEMethodData. These data objects hold data for subsequent computations. The creation of these data objects is discussed in the next section.

Stationary PDEs

Initialization Stage

The initialization stage creates all data objects that are created during a run of . This section will illustrate that process. It is possible to skip this section and continue with the discretization stage and make use of the initialized data structures creates. With this it is possible to use as an equation preprocessor, for example, for a new numerical discretization method.

Variable and Solution Data

For NDSolve to find a solution, the specification of the dependent (u) and independent () variables and the region (Ω) needs to be given. The same is true for the finite element functions. The variable and solution data are stored in lists that make a specific data structure corresponding to different components. These data lists may be generated using NDSolve`VariableData and NDSolve`SolutionData. Elements of the data lists can be filled in or modified using NDSolve`SetSolutionDataComponent.

Create the variable data with dependent and independent variable names.
In[19]:=
Click for copyable input
Out[19]=

Accompanying the variable data is a solution data structure. The solution data structure is essentially the numerical incarnation of the variable data. For example, the numerical version of the component of the solution data is the NumericalRegion that will also contain the ElementMesh. "Element Mesh Generation" gives more information about the creation of the mesh from the region Ω.

Specify a NumericalRegion.
In[20]:=
Click for copyable input
Out[20]=
Create an ElementMesh from the region.
In[21]:=
Click for copyable input
Out[21]=

The numerical region now has an ElementMesh.

Inspect the ElementMesh created from the region.
In[22]:=
Click for copyable input
Out[22]=
Visualize the element mesh.
In[23]:=
Click for copyable input
Out[23]=
Create the solution data with the component set.
In[24]:=
Click for copyable input
Out[24]=
InitializeFEMData

To initialize finite element method data, which is needed in subsequent steps, InitializePDEMethodData is used. Currently, the only discretization method available in this framework is the finite element method. Thus, by default, InitializePDEMethodData generates a FEMMethodData object. During the initialization, a method-specific setup is performed; for example, among other things, the interpolation and integration order of the method are set up.

Initialize the finite element data with the variable and solution data.
In[25]:=
Click for copyable input
Out[25]=

The initialized PDE method data is stored in the FEMMethodData object. The display form of FEMMethodData shows degrees of freedom, the interpolation order for each dependent variable, and integration order.

Query degrees of freedom, element, and integration order.
In[26]:=
Click for copyable input
Out[26]=

The degrees of freedom of a system of equations result from a combination of the nodes in the element mesh, the number of dependent variables, and the interpolation order of the dependent variables if multiple dependent variables are given. When one dependent variable is given, the degrees of freedom correspond to the nodes in the ElementMesh object.

For one dependent variable, the degrees of freedom correspond to the number of coordinates in the mesh.
In[27]:=
Click for copyable input
Out[27]=

The is the order of accuracy used to integrate the finite element operators.

Increase the integration order.
In[28]:=
Click for copyable input
Out[28]=

For a first-order element mesh, the default integration order is to use second-order accurate integration points, while for a second-order element mesh, the default is to use fourth-order accurate integration points. The maximum supported integration order is fifth order.

Reducing the integration order will result in a faster solution, but potentially a less-accurate solution. Increasing the integration order will result in a longer time to solve the finite element model.

Inspect that the state object created through generates the same FEMMethodData object.
In[29]:=
Click for copyable input
In[30]:=
Click for copyable input
Out[30]=

After the method data is set, the PDE is initialized.

InitializePDECoefficients

Consider the following partial differential equation in one dependent variable :

.

That PDE is defined in . The coefficients , , , and are scalars. , , and are vectors of length and is an × matrix. The input coefficients to InitializePDECoefficients are found by correspondence of the above PDE with the PDE to be modeled.

The coefficients of the model PDE.
In[31]:=
Click for copyable input
Out[31]=
In[32]:=
Click for copyable input
Out[32]=

All coefficients that are not explicitly specified are assumed to be zero.

Initialize the partial differential equation coefficients.
In[33]:=
Click for copyable input
Out[33]=

The initialized coefficients are stored in the PDECoefficientData data structure. The display form shows the system size, i.e. the number of dependent variables (1), and the operator size, i.e the space dimension (2) of the initialized coefficients.

Extract the system size and the spatial dimension of the initialized coefficients.
In[34]:=
Click for copyable input
Out[34]=
Extract the raw coefficients.
In[35]:=
Click for copyable input
Out[35]=

During the initialization, the coefficients are classified into different categories of coefficients, such as stationary coefficients and transient coefficients. This will be discussed later in more detail.

In[36]:=
Click for copyable input
Out[36]=
InitializeBoundaryConditions

Similar to the partial differential equation initialization, the boundary conditions need to be initialized.

Initialize the boundary conditions.
In[37]:=
Click for copyable input
Out[37]=

Boundary condition coefficients are classified into categories similar to the PDE coefficients.

In[38]:=
Click for copyable input
Out[38]=
Extracting the Initialized Finite Element Data from the NDSolve`StateData

As an alternative to the manual initialization of the finite element method data, PDE coefficients, and boundary conditions, can be utilized for this purpose.

Extract the finite element data from the .
In[39]:=
Click for copyable input
Out[40]=
Inspect that the objects created are the same.
In[41]:=
Click for copyable input
Out[41]=
In[42]:=
Click for copyable input
Out[42]=

In some cases, details like arithmetic order or generated unique variable names may be different, but otherwise the expressions are the same.

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

Discretization Stage

During the discretization stage, the continuous PDE and boundary conditions are approximated with discrete variants. This process is at the heart of the partial differential equation solving process.

DiscretizePDE

The initialized PDE coefficients are transformed into a discrete form.

Compute the discretized partial differential equation.
In[46]:=
Click for copyable input
Out[46]=

The display form of the DiscretizedPDEData shows the first dimension of the system matrices.

Extract all system matrices.
In[47]:=
Click for copyable input
Out[47]=
Visualize the stiffness matrix.
In[48]:=
Click for copyable input
Out[48]=
Extract the system matrices separately.
In[49]:=
Click for copyable input
Out[49]=

During the matrix assembly process, access to the finite elements is possible. Note that the form of the finite elements is such that the matrix assembly process is efficient.

Extract the computed finite elements during the matrix assembly.
In[50]:=
Click for copyable input
Out[51]//Short=
DiscretizeBoundaryCondition

DiscretizeBoundaryConditions computes a discretized version of the boundary conditions given. In a later step, the discretized boundary conditions are then deployed into the actual system matrices.

Discretize the initialized boundary conditions.
In[52]:=
Click for copyable input
Out[52]=

The DirichletCondition and NeumannValue boundary specifications are represented differently. The contributions from NeumannValue are to the system load matrix and the system stiffness matrix.

Components from the generalized Neumann boundary value.
In[53]:=
Click for copyable input
Out[53]=

The component of the discretized boundary condition has nonzero contributions if the NeumannValue has a dependency on a dependent variable; otherwise, it is zero. In other words, if the Neumann value is a generalized Neumann value, then there is a contribution to the component.

The DirichletCondition representation is given through a matrix that contains positions. These positions specify which degrees of freedom should have the condition applied. A second matrix is returned that contains the values that are to be applied.

Components from the DirichletCondition boundary condition.
In[54]:=
Click for copyable input
Out[54]=

Solution Stage

DeployBoundaryConditions

To have the discretized boundary conditions take effect, they need to be deployed into the system matrices. To save memory, this is done in place; no new matrices are generated.

Deploy the boundary conditions in place.
In[55]:=
Click for copyable input

DeployBoundaryConditions has attribute HoldFirst, so that the system matrices may be modified in place without making copies first. This means, however, that the input to DeployBoundaryConditions needs to be given in terms of symbols that have matrices as values rather than actual matrices. If the unmodified system matrices are needed later, they can be re-extracted from the DiscretizedPDEData. In other words, if the unmodified system matrices are no longer needed and computer memory consumption is crucial, then the DiscretizedPDEData needs to be cleared, e.g. with ClearAll.

During the deployment of the boundary conditions, the contributions from NeumannValue are added to the stiffness and load matrix. Then the system matrices are modified such that the DirichletCondition is satisfied.

After the boundary condition deployment, the system matrices have changed; the system matrices typically contain fewer elements.
In[56]:=
Click for copyable input
Out[56]=

Note that generally only after the deployment of the boundary conditions is the system of equations solvable.

Equation Solving

The system of equations can then be solved with LinearSolve.

Solving the system of equations.
In[57]:=
Click for copyable input
Out[57]//Short=

Post Process

An interpolation function based on an ElementMesh can also be created. In this case, the second argument to is the values at the mesh nodes.

Create the interpolation function with values at the nodes.
In[58]:=
Click for copyable input
Out[58]=
Visualize the interpolation function.
In[59]:=
Click for copyable input
Out[59]=

Typically, an interpolation function extrapolates if queried outside of its domain. This behavior can be changed.

Evaluate outside the domain using extrapolation.
In[60]:=
Click for copyable input
Out[60]=
Switch off automatic extrapolation in interpolation functions by having it return Indeterminate.
In[61]:=
Click for copyable input
Out[61]=
In[62]:=
Click for copyable input
Out[62]=

Transient PDEs

Transient PDEs are time-dependent PDEs. In the simplest case, a transient PDE has a time derivative of the dependent variable as a part of the PDE. In this case, an additional system matrix is created that represents the discretized time derivative. The time integration of the system matrices can then be performed, for example, by NDSolve. Other than that the workflow is much the same as for a stationary PDE. More involved examples have coefficients and/or boundary conditions that are dependent. A later section will show how to time integrate those.

Transient PDE with Stationary Coefficients and Stationary Boundary ConditionsIntroduction

The first example is a heat equation in 1D. Consider the following model PDE

in the spatial region 0 to 1 and a time domain from 0 to 1. Boundary and initial conditions are 0 everywhere.

Use NDSolve directly to find a solution to the PDE.
In[2]:=
Click for copyable input
Out[2]=

In this case, the method option specifies that the spatial variable should be discretized with the finite element method and that time integration is done with the method of lines. This approach is called semi-discretization. In contrast, the full discretization would discretize both the temporal variable and the spatial variables.

Plot the numerical solution.
In[3]:=
Click for copyable input
Out[3]=

The process of finding the solution to the PDE is similar to the stationary case. The main reason is that neither the PDE coefficients nor the boundary conditions depend on time. To find the solution of the transient system, the approach of semi-discretizing the PDE is taken. The PDE is discretized in space () with the finite element method. The result of that FEM discretization is then a system of ordinary differential equations in the temporal variable (). NDSolve will be used to time integrate that system of ordinary differential equations.

First, the variable data is created and populated.

Create the variable data with dependent and independent variable names.
In[4]:=
Click for copyable input
Out[4]=

Next, the solution data is set up with an initial time and a mesh.

Specify a NumericalRegion.
In[5]:=
Click for copyable input
Out[5]=
Create an ElementMesh from the region.
In[6]:=
Click for copyable input
Out[6]=
Create the solution data with the and component set.
In[7]:=
Click for copyable input
Out[7]=

The PDE coefficients, boundary conditions, and method data are initialized.

Initialize the finite element data with the variable and solution data.
In[8]:=
Click for copyable input
Out[8]=
Initialize the partial differential equation coefficients.
In[9]:=
Click for copyable input
Out[9]=
Initialize the boundary conditions.
In[10]:=
Click for copyable input
Out[10]=

The discretization is performed next. Even though a transient system is considered, the discretization has only stationary components since neither the PDE coefficients nor the boundary conditions have a dependency on time.

Compute the discretized partial differential equation.
In[11]:=
Click for copyable input
Out[11]=
Discretize the initialized boundary conditions.
In[12]:=
Click for copyable input
Out[12]=
Extract all system matrices.
In[13]:=
Click for copyable input
Deploy the boundary conditions in place.
In[14]:=
Click for copyable input

Next, the initial conditions are set up and NDSolve will time integrate the discretized system of equations.

Set up initial conditions based on the boundary conditions.
In[15]:=
Click for copyable input
Time integrate the system of equations with NDSolve.
In[17]:=
Click for copyable input
Out[17]=

The NDSolve options are chosen such that they match the automatically selected options from the call to NDSolve above where NDSolve controlled all aspects of the solution process.

Set up a function that, given a time , constructs and memorizes an interpolating function.
In[18]:=
Click for copyable input
Visualize the difference between the automatic and manual solutions.
In[20]:=
Click for copyable input
Out[20]=

Transient PDE with Stationary Coefficients and Stationary Boundary Conditions

Consider the following model PDE

subject to the boundary conditions at and at . For more information about markers, please consult the "Element Mesh Generation" tutorial.

The PDE has a first-order time derivative, but other than that the PDE coefficient and the boundary conditions are not time dependent. One strategy to solve the PDE is to first discretize the PDE and boundary conditions and then the resulting set of ODEs can be time integrated with NDSolve as a time integrator. This process is how NDSolve internally solves transient PDEs and will be shown in the next section.

For this example, a boundary mesh with markers has been predefined. The geometry models a Scanning Electrochemical Microscopy (SECM) probe. This specific probe is a heptode as it has seven electrodes embedded in a circular glass structure.

Use a predefined boundary mesh.
In[2]:=
Click for copyable input
Convert the BoundaryMeshRegion to a boundary ElementMesh.
In[3]:=
Click for copyable input
Out[3]=

The predefined boundary mesh region has integer markers set to the boundary faces and points. Those markers can be visualized and used as positions where boundary conditions are applied during the solution process.

Visualize the point marker according to their grouping.
In[4]:=
Click for copyable input
Out[4]=

Since the boundary conditions only involve Dirichlet conditions, only the point elements and markers are of interest. Dirichlet conditions always operate on point elements. For completeness, the boundary mesh also has markers for the boundary elements implemented. These were needed to specify Neumann type boundary values, which always operate on boundary elements.

In[5]:=
Click for copyable input
Out[5]=
Create a full first-order mesh with optimized node connectivity.
In[6]:=
Click for copyable input
Out[6]=

Strictly speaking, the node reordering is optional. When the node reordering method option is set to True, the bandwidth the system matrices have is minimized. The lower bandwidth can be beneficial for iterative solvers. Since the node reordering need only be done once, and the solvers are used for each time step, it is often advantageous to use node reordering for transient problems.

Set up variable data.
In[7]:=
Click for copyable input
Out[7]=
Set up solution data.
In[8]:=
Click for copyable input
Out[8]=

Populate the solution data. The component is the initial time where the integration starts. Setting the and components sets the initial condition and derivative of the initial condition.

Initialize the method data.
In[9]:=
Click for copyable input
Out[9]=
Initialize the PDE coefficients.
In[10]:=
Click for copyable input
Out[10]=
Initialize the PDE boundary conditions.
In[11]:=
Click for copyable input
Out[11]=
Discretize the system.
In[12]:=
Click for copyable input
Out[12]=
In[13]:=
Click for copyable input
Out[13]=
Extract all system matrices and deploy the boundary conditions.
In[14]:=
Click for copyable input
Out[14]=
In[15]:=
Click for copyable input

Note that the first transient matrix, the damping matrix, is populated. This will be needed during the time integration later.

Visualize the stiffness matrix.
In[16]:=
Click for copyable input
Out[16]=

Note the structure of the stiffness matrix. Using "NodeReordering"->False during the mesh generation process will result in much less structured system matrices.

Before computing the transient solution, inspect the stationary solution.

Find the stationary solution.
In[17]:=
Click for copyable input
Alternatively, find the stationary solution using a different solver.
In[18]:=
Click for copyable input
Visualize the solution as a surface plot.
In[19]:=
Click for copyable input
Out[19]=
Construct an interpolation function that returns Indeterminate if queried outside of its domain and does not warn.
In[20]:=
Click for copyable input
Out[20]=
ContourPlot through the direction.
Visualize the solution as a 3D contour plot within the model.
In[22]:=
Click for copyable input
Out[22]=

The PDE discretization resulted in a system of ODEs. The time integration of the ODEs can be done with NDSolve. For this end, initial conditions that match the boundary conditions need to be set up.

Set up initial conditions based on the boundary conditions.
In[23]:=
Click for copyable input
Visualize the initial conditions.
In[25]:=
Click for copyable input
Out[25]=

To speed up the time integration of NDSolve, it is possible to specify the sparsity pattern of the Jacobian through an option. This is not necessary when the system is solved as a whole though NDSolve. It is necessary here because NDSolve does not try to decipher the structure for vector variables. The sparsity pattern of the Jacobian is available through the pattern of the damping matrix.

Set up the sparsity pattern for the Jacobian based on the damping matrix.
In[26]:=
Click for copyable input
Out[26]=
Find the solution to the transient PDE with NDSolve as a time integrator.
In[27]:=
Click for copyable input
Out[27]=
Out[28]=

Note that using an EvaluationMonitor slows the solution process down a little bit but allows live monitoring of the solution process.

The solution shows how the transient system reaches the steady state computed above.

Visualize the reaching of steady state at various time steps.
In[29]:=
Click for copyable input
In[30]:=
Click for copyable input
Out[30]=

Model Order Reduction of Transient PDEs with Stationary Coefficients and Stationary Boundary Conditions

As a scenario for the transient case where access to the system matrices is important, we consider model order reduction. The question model order reduction is addressing is the following, Considering a transient system of equations, is there an equivalent system of equations that is smaller than the original while preserving most of the original properties?

Consider a discretized system of equations with degrees of freedom of the form

The damping matrix and the stiffness matrix are of size and the load matrix is of size . We will now establish a second, much smaller system of equations with degrees of freedom:

with .

Because the degrees of freedom of the second, reduced, system are much smaller than the original system , the time integration will be much more efficient. The following algorithm establishes the projection vector that maps between and :

Here is of size , is of size , and is an approximation discrepancy.

As a rough explanation, the projection vector is found as follows: Given a nonzero starting vector and a matrix , the algorithm finds orthonormal such that

The model order reduction projection code.
In[32]:=
Click for copyable input

The procedure returns the projection vector that has the dimensions of the reduced system times the original system.

Find an approximate system of equations that has only degrees of freedom.
In[33]:=
Click for copyable input
Out[33]=
Extract the reduced system matrices and inspect the dimensions.
In[34]:=
Click for copyable input
Out[35]=
Efficiently time integrate the reduced system.
In[36]:=
Click for copyable input
Out[38]=
Out[39]=
Compare the norm of the full system and the reduced system.
In[40]:=
Click for copyable input
Out[42]=

Transient PDEs with Transient Coefficients

The next step is to consider transient PDEs, where one or more of the coefficients is time dependent. Consider the following model PDE with a time-dependent load ,

subject to the same stationary boundary conditions as in the previous example. In this example, use a load function of

Also, variable and solution data are the same as in the previous example. The FEM method data can also be reused.

The PDE Initialization now has an additional component, the load coefficient, which it sets to the time variable .

Initialize the PDE coefficients.
In[92]:=
Click for copyable input
Out[92]=

The initialized coefficients store the coefficients specified.

Extract the stored coefficients.
In[93]:=
Click for copyable input
Out[93]=

In order to perform efficient time integration, the coefficients are classified into stationary and transient coefficients. Stationary coefficients are evaluated only once, much like in the previous example. Those coefficients that have time-dependent components need to be reevaluated from time step to time step and added to the stationary discretization.

The Initialization of the PDE coefficients internally performs a classification of all coefficients. This classification can be inspected and made use of. Each classification, e.g. the transient coefficients, is represented by two lists containing positions. The first list contains scalar coefficients while the second contains vector and matrix coefficients. To illustrate further, consider the general PDE form

The coefficients , , , and are scalars; , , and are vectors; and is an × matrix. The first of the two lists returned represents the scalar coefficients (, , , and ) and the second represents the vector and matrix coefficients (, , , and ).

Extract the position of the transient initialized PDE coefficients.
In[94]:=
Click for copyable input
Out[94]=

There is a single transient component present in the model PDE.

Extract the transient scalar coefficients from all coefficients.
In[95]:=
Click for copyable input
Out[95]=

There are no transient matrix components in this example.

Extract the transient matrix coefficients from all coefficients.
In[96]:=
Click for copyable input
Out[96]=

The coefficient represents all coefficients that are stationary, including those that are not specified and hence are zero.

This classification into stationary and transient components allows for a separate computation of the stationary and transient discretization of the PDE in an efficient manner.

Discretize and store the stationary components of the PDE.
In[97]:=
Click for copyable input
Out[98]=

In the above example, the string argument is optional, as this is the default. Note that the load vector of the discretized stationary PDE coefficients is empty; the load vector has only transient components.

Discretize the transient components of the system at a specific point in time.
In[99]:=
Click for copyable input
Out[101]=

The load vector of the discretized transient PDE coefficients is populated. The time has been substituted with the numerical value 3. Also note that no other system matrices are populated, as they do not have transient coefficients.

The rationale behind the classification of PDE coefficients is that stationary coefficients can be discretized once, stored, and then reused in the time integration loop. This means that only coefficients that exhibit a time variable dependence need to be discretized during each step. This saves a lot of computational effort. Also note that since important geometric information has been precomputed and stored in the FEMMethodData object, the computational time for each transient component is kept to a minimum.

Reuse the discretized PDE boundary conditions.
In[102]:=
Click for copyable input
Out[102]=
A function to compute the system matrices and re-deploy the stationary boundary conditions.
In[103]:=
Click for copyable input

The computation of the right-hand side makes use of the stationary coefficients of the PDE that have been computed already. In this example, the damping matrix has no transient components and the factorization can be cached. If the damping matrix has time-dependent components, the above computation of the right-hand side needs to be modified, since then during each time step the damping matrix is updated and the LinearSolve step can not be cached; the factorization then needs to be performed for each time step.

Since there are no transient boundary conditions, the initial conditions are set up in the same manner as before.

Set up initial conditions based on the boundary conditions.
In[104]:=
Click for copyable input

The sparsity pattern is set again derived from the damping matrix. In this case, the boundary conditions have not been deployed into the damping matrix.

Set up the sparsity pattern for the Jacobian based on the damping matrix.
In[106]:=
Click for copyable input
Out[106]=
In[107]:=
Click for copyable input
Out[107]=
Find the solution to the transient PDE with NDSolve as a time integrator.
In[108]:=
Click for copyable input
Out[108]=
Out[109]=

The visualization is done in the same manner as before.

Visualize the reaching of steady state at various time steps.
In[111]:=
Click for copyable input
Out[112]=

Coupled PDEs

This section explains how coupled partial differential equations are solved. In principle, the workflow is much the same as for single PDEs. The differences will be pointed out below.

Deformation of a Beam under Load

As an example, a structural mechanics problem, the deformation of a beam, is considered. The deformation can be modeled as a coupled PDE in two dependent variables.

A plane stress formulation as a coupled PDE.
In[113]:=
Click for copyable input

In the plane stress operator, is Young's modulus and ν is Poisson's ratio. Those are material properties.

The beam to be deformed is of 5 units in length and 1 unit in height.

Create a mesh region with an explicit bounding box and a visualization of the mesh.
In[115]:=
Click for copyable input
Out[117]=

The boundary conditions are as follows: at the left-hand side, the beam is held fixed in both and direction. On the right-hand boundary, a downward force of one unit is applied; a NeumannValue enforces a boundary load on the beam in the negative direction.

Set up boundary conditions at the left- and right-hand side of the domain.
In[118]:=
Click for copyable input

The initialization has the same components as in previous sections. First, the variable and solution data are set up.

Populate the variable data with two dependent variables and and independent variable names.
In[121]:=
Click for copyable input

For coupled PDEs, more than one dependent variable needs to be given.

Set the component of the solution data.
In[122]:=
Click for copyable input

The initialization of the method data is done in the same manner as for a single dependent variable.

Initialize the finite element data with the variable and solution data.
In[123]:=
Click for copyable input
Out[123]=

The display form of FEMMethodData shows the degrees of freedom, the interpolation order for each dependent variable, and integration order.

The degrees of freedom of a system of equations result from a combination of the nodes in the element mesh, the number of dependent variables, and the interpolation order of the dependent variables.

The degrees of freedom are related to the number of dependent variables and the nodes in the mesh.
In[124]:=
Click for copyable input
Out[124]=
Out[125]=

In the coupled PDE case, the coefficients become a little more involved. A coupled PDE is in essence a conglomeration of single PDEs; however, each component is present for each dependent variable. Consider the following partial differential equation in two dependent variables , :

The coupled PDE is defined in . The coefficients and are scalars in each PDE. , , and are vectors of length and is an × matrix in each PDE. The input coefficients to InitializePDECoefficients are found by correspondence of the above coupled PDE with the coupled PDE to be modeled.

The coefficients of the model PDE.
In[126]:=
Click for copyable input

In the coupled case, there is a diffusion coefficient matrix for each dependent variable in each PDE.

Initialize the partial differential equation coefficients.
In[127]:=
Click for copyable input
Out[127]=

The display form of PDECoefficientData shows the system size, i.e. the number of dependent variables (two) and the operator size, i.e the space dimension (two) of the initialized coefficients.

Next, the boundary conditions are set up. Each of the dependent variables has a list of boundary conditions. This associates the boundary conditions with specific parts of the system of PDEs. The first set of boundary conditions is associated with the first PDE of the coupled system, the second set of boundary conditions is associated with the second set PDE of the coupled system, and so on. This allows you to uniquely associate boundary conditions with parts of PDEs.

Initialize the boundary conditions.
In[128]:=
Click for copyable input
Out[128]=

The initialized PDE coefficients are transformed into a discrete form.

Compute the discretized partial differential equation.
In[129]:=
Click for copyable input
Out[129]=

The display form of the DiscretizedPDEData shows the first dimension of the system matrices.

Extract the system matrices.
In[130]:=
Click for copyable input
Out[130]=
Visualize the stiffness matrix.
In[131]:=
Click for copyable input
Out[131]=

The visualization of the stiffness matrix shows that the structure of the underling PDE coefficients is a system of two dependent variables. The degrees of freedom are split between the two dependent variables.

Compute how the incidents are split between the two dependent variables.
In[132]:=
Click for copyable input
Out[132]=
Discretize the initialized boundary conditions.
In[133]:=
Click for copyable input
Out[133]=
Deploy the boundary conditions in place.
In[134]:=
Click for copyable input

The system of equations can then be solved with LinearSolve.

Solve the system of equations.
In[135]:=
Click for copyable input
Out[135]//Short=

Interpolation functions for each dependent variable can be created. Since the solution form LinearSolve is one solution vector containing both solutions to and , the solution vector needs to be split up.

Create the interpolation function with values at the nodes.
In[136]:=
Click for copyable input
Out[136]=
Out[137]=
Visualize the interpolation function.
In[138]:=
Click for copyable input
Out[138]=
Out[139]=
Visualize the deformed structure in red.
In[140]:=
Click for copyable input
Out[140]=
In[141]:=
Click for copyable input
Out[141]=

A Swinging BeamTransient Coupled PDEs

For a transient coupled PDE, not much more is needed. The following example shows how to deflect the beam and have a Rayleigh damping dampen the motion out. In Rayleigh damping, the damping matrix is a linear combination of the mass matrix and the stiffness matrix. The structure of the system matrices for the two dependent variable coupled transient PDEs are represented as follows:

where the damping system matrix is a combination of the mass matrix and the stiffness matrix .

The additional coefficients of the transient model PDE.
In[142]:=
Click for copyable input
Out[142]=

The transient coefficients also have the mass matrix coefficients set.

Initialize the transient partial differential equation coefficients.
In[143]:=
Click for copyable input
Out[143]=
Compute the discretized partial differential equation.
In[144]:=
Click for copyable input
Out[144]=

The display form of the DiscretizedPDEData shows the first dimension of the system matrices.

Extract the system matrices.
In[145]:=
Click for copyable input
Out[145]=

Note that the damping matrix sparse array is empty. This is because the mass matrix coefficient has been set and the default damping matrix coefficient used here is zero.

Visualize the mass matrix.
In[146]:=
Click for copyable input
Out[146]=

The mass matrix shows the structure of the PDE coefficients where the block diagonals are populated and the off diagonal blocks are zero.

The Rayleigh damping is computed as a linear combination of the mass and stiffness matrices.

Compute the Rayleigh damping matrix.
In[147]:=
Click for copyable input
Deploy the boundary conditions.
In[148]:=
Click for copyable input

The initial condition and the derivative of the initial condition are set to zero. This is reasonable since the Dirichlet boundary conditions are all zero and the beam starts from rest.

Set up the initial conditions.
In[149]:=
Click for copyable input
Set up the sparsity pattern for the Jacobian based on the damping matrix.
In[151]:=
Click for copyable input
Out[151]=
In[152]:=
Click for copyable input
Out[152]=
Find the solution to the transient PDE with NDSolve as a time integrator.
In[157]:=
Click for copyable input
Out[157]=
Out[158]=
Visualize the swinging beam.
In[159]:=
Click for copyable input
In[160]:=
Click for copyable input
Out[160]=

Large-Scale FEM Analysis

The finite element method, as implemented in NDSolve, has been optimized for speed. This optimization comes at a cost. It is common in numerics to be able to "trade" speed versus memory consumption. This means that a faster implementation will need more memory to solve a problem at hand. It is, however, possible to tune the implemented finite element method through various means to use less memory, at the cost of some additional time that is needed to find the solution to the PDE. Some suggestions have been presented in "Solving Partial Differential Equations with Finite Elements". The following sections discuss a range of additional options to decrease the amount of memory needed during a finite element solution. The suggestions range from general advice to specific options.

For this section to be illustrative, it is good to start with a fresh kernel.

Generally, the Wolfram Language stores references to previously computed results. This mechanism can result in an accumulation of data that, if given some thought, is not needed for future computations. The mechanism that controls the history length can be set to not remember previous output.

Set the history length to 0 to not store previous output.
In[2]:=
Click for copyable input

The general procedure in the finite element discretization is as follows:

  • Mesh generation.

During the finite element analysis, there are two key memory bottlenecks. The first is the creation of the stiffness matrix (DiscretizePDE), and the second is the solution (LinearSolve) of the system of equations. Options to both mesh generation and geometric data precomputation have an effect on the memory requirement during discretization and solving. Furthermore, by default, the discretization step itself is performed in two steps. First all elements for a specific system matrix are computed, and then the matrix is assembled from those elements. An additional distinction is to be made if a given memory-saving action taken will merely need more time for the computation, or if that specific action also has an effect on the quality of the solution.

The options, which can also be passed to NDSolve directly, are best discussed with an example. Typical examples that require a great deal of memory are 3D application examples and/or involve many coupled PDEs.

Define the NavierCauchy operator.
In[3]:=
Click for copyable input
Define the Lamé parameters.
In[4]:=
Click for copyable input
Boundary conditions.
In[5]:=
Click for copyable input

The first crucial step is the element mesh generation. There are various aspects that influence the memory needed to solve the PDE over the mesh. The has a large effect on how much memory is needed. For example, a first-order tetrahedron has four nodes, whereas a second-order tetrahedron has 10 nodes. Reducing the also reduces the quality of the solution that can be expected, as has been discussed in "Solving Partial Differential Equations with Finite Elements".

Create the variable data with dependent and independent variable names.
In[9]:=
Click for copyable input
Out[9]=
Create a boundary mesh.
In[10]:=
Click for copyable input
Out[10]=

During the matrix assembly in the upcoming discretization stage, all elements from a mesh are first computed and then assembled. This means that at one point in time the elements and the system matrix are in memory. This can be avoided by setting the option to a positive integer. Then the mesh elements are split into a number of blocks. Once an element mesh block is computed, it is then assembled into the system matrix. Then the block is discarded and the next block is assembled into the same system matrix, until all blocks have been computed and assembled. With this scheme, the maximum memory consumption is the fully assembled system matrix and the size of the largest element mesh block and its computed elements. Setting too large will slow down the assembly process.

Create a mesh with 10 mesh element blocks and measure the memory usage in megabytes.
In[11]:=
Click for copyable input
Out[11]=
Inspect the number of mesh elements blocks in the mesh.
In[12]:=
Click for copyable input
Out[12]=
Inspect the size of the mesh.
In[13]:=
Click for copyable input
Out[13]=
Create the solution data with the component set.
In[14]:=
Click for copyable input
Out[14]=

During equation processing, a call to InitializePDEMethodData is made. During this call, geometric data that is needed for the finite element method is precomputed. If is used solely for equation processing, the geometric preprocessing can be switched off.

Initialize the finite element data with the variable and solution data.
In[15]:=
Click for copyable input
Extract the raw diffusion coefficients and delete the state data.
In[16]:=
Click for copyable input
Initialize the partial differential equation coefficients.
In[18]:=
Click for copyable input
Out[18]=
Initialize the boundary conditions.
In[19]:=
Click for copyable input
Out[19]=

During the discretization step, the precomputed geometry data is needed to compute the contribution of each element to the system matrices. Note that for transient systems, switching off the geometric precomputation will require that this data is recomputed in each time step. This will further increase the time needed to solve the system, but for very large systems this may be an option to obtain a solution at all.

Initialize the method data.
In[20]:=
Click for copyable input
Out[20]=
Investigate the current memory usage.
In[21]:=
Click for copyable input
Out[21]=
Compute the discretized partial differential equation.
In[22]:=
Click for copyable input
Out[22]=
Investigate the current memory usage.
In[23]:=
Click for copyable input
Out[23]=
Extract the load vector.
In[24]:=
Click for copyable input
Extract the stiffness matrix.
In[25]:=
Click for copyable input
Out[25]=
Investigate the memory needed to store the stiffness matrix.
In[26]:=
Click for copyable input
Out[26]=
Discretize the initialized boundary conditions.
In[27]:=
Click for copyable input
Out[27]=

Next, there is a very subtle point to consider. The system matrices are stored in the DiscretizedPDEData data object. Those have been extracted above and stored into the symbols load and stiffness. The deployment of the boundary conditions will modify these system matrices. But since they are stored in the DiscretizedPDEData (it is possible to extract the unmodified system matrices once more), a duplicate version will be created. To avoid this and in order to save memory, the DiscretizedPDEData object can be deleted such that only one reference to the system matrices exists.

Clear unnecessary data.
In[28]:=
Click for copyable input
Investigate the current memory usage.
In[30]:=
Click for copyable input
Out[30]=
Deploy the boundary conditions in place.
In[31]:=
Click for copyable input
Investigate the current memory usage.
In[32]:=
Click for copyable input
Out[32]=
Solving the system of equations.
In[33]:=
Click for copyable input
Out[33]=

MaxMemoryUsed in this case does not report all memory used, since a library is called over which the kernel does not have memory control.