# Finite Element Programming

Introduction | Nonlinear PDEs |

Finite Element Data within NDSolve | Transient PDEs |

Passing Finite Element Options to NDSolve | Coupled PDEs |

A Workflow Overview | Large-Scale FEM Analysis |

Stationary PDEs |

## 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 NDSolve`FEM` package provides a lower-level interface that gives extensive control for each part of the solution process.

In[1]:= |

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

- 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 NDSolve`ProcessEquations, NDSolve`Iterate and NDSolve`ProcessSolutions and a data object NDSolve`StateData 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 NDSolve`StateData. The following section explores how to access the finite element data from the NDSolve`StateData object. The finite element data is stored in an NDSolve`FEM`FiniteElementData object and accessed through the "FiniteElementData" property.

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

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

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 NDSolve`Iterate finds the solution of a system of equations with LinearSolve.

In[5]:= |

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

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

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

The same is possible for transient PDEs.

Options for ToBoundaryMesh and ToElementMesh can be passed to NDSolve by specifying "MeshOptions", as shown above. Additionally, options to InitializePDECoefficients and PDESolve and their sub-options can be passed to NDSolve via "InitializePDECoefficientsOptions" and "PDESolveOptions", respectively. Also, the integration order and interpolation order can be specified via options.

## 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 NDSolve`FEM`FiniteElementData object is created through the functions introduced in the following.

The solution is found in three stages:

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

InitializePDECoefficients | initialize partial differential equation coefficients |

InitializeBoundaryConditions | initialize boundary conditions |

InitializePDEMethodData | initialize partial differential equation method data |

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.

DiscretizePDE | discretize initialized partial differential equations |

DiscretizeBoundaryConditions | discretize initialized boundary conditions |

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

DeployBoundaryConditions | deploy discretized boundary conditions into discretized partial differential equations |

LinearSolve | solve linear systems of equations |

Once the solution is found, InterpolatingFunction objects can be constructed.

ProcessPDESolutions | generates InterpolatingFunction objects |

As an alternative for stages two and three, the function PDESolve can be made use of.

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 being solved by a numerical method such as the finite element method, three components are needed:

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

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

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

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

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 NDSolve`ProcessEquations. 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 NDSolve`ProcessEquations creates. With this it is possible to use NDSolve`ProcessEquations 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 ({x,y}) variables and the region (Ω) needs to be given. The same is true for the finite element functions. The variable and solution data is 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.

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 "Space" 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 Ω.

##### 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 preceding PDE with the PDE to be modeled.

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

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.

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.

##### InitializeBoundaryConditions

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

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

##### InitializePDEMethodData

To initialize the finite element method data, which is needed in subsequent discretization 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 and an ElementMesh is generated from the NumericalRegion.

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.

All mesh generation options as discussed in the "Element Mesh Generation" tutorial can be specified.

During the method initialization, an ElementMesh object will be generated and stored in the NumericalRegion object.

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.

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

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.

In[37]:= |

##### 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, NDSolve`ProcessEquations can be utilized for this purpose.

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

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

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.

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

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

The "StiffnessMatrix" 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 "StiffnessMatrix" 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.

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

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.

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.

### Post-Processing

As a post-processing step, an interpolation function can be created from the solution obtained from LinearSolve. For this, the solution data is updated with the solution found and ProcessPDESolutions then constructs the InterpolatingFunction.

In[63]:= |

ProcessPDESolutions has the advantage that it also works well for multiple dependent variables. In the specific case at hand, an alternative to ProcessPDESolutions is to create the interpolation function directly with ElementMeshInterpolation. Here the first argument is the mesh and the second argument are the values at the mesh nodes.

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

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

It is also possible to extract the ElementMesh from the interpolating function.

## Nonlinear PDEs

In the following section, the solution of nonlinear stationary partial differential equations will be discussed. This will be done at several levels. First, a top-level example is given, and subsequently the section will drill deeper into how this works and is implemented.

### Top-Level Example

To illustrate the solution process of a nonlinear partial differential equation, a minimal surface over a region with a particular boundary is to be found. The minimal surface can be found by solving the following nonlinear equation:

The coefficient is nonlinear because it contains the dependent variable ; specifically, it contains the gradient of the dependent variable . To formulate nonlinear PDEs, any of the coefficients in

that are detailed in equation 1 can depend on both and and partial derivatives of the first order like . They cannot, however, depend on higher than first-order derivatives such as .

In general, nonlinear equations that make use of the finite element method as a solution method can be set up in the typical way that is common for NDSolve. Since NDSolve does not hold its arguments, the input equations will be evaluated before they reach NDSolve. To avoid premature evaluation, PDEs can be specified in Inactive form. Nonlinear equations that are not wrapped in Inactive may produce unwanted and unnecessary higher-order derivatives that NDSolve then cannot solve, especially if derivatives of the dependent variable are present in the diffusion coefficient , as is the case here. In such a case, the PDE has to be inactivated. Inactivating a PDE is the most general way of specifying PDEs, and if in doubt, this should be the method chosen.

In[2]:= |

In[4]:= |

Very similar to FindRoot, NDSolve needs an initial seeding for each dependent variable when solving stationary nonlinear partial differential equations with the finite element method. The initial seeding can be specified with InitialSeeding.

If an initial seeding is not specified, it is taken to be zero.

If the equation at hand is not in inactive form, NDSolve will not be able to solve the PDE.

For more information on inactive PDEs, refer to the section about Formal Partial Differential Equations.

In the sections that follow, the general internal solution process that NDSolve takes for nonlinear equations like the preceding are presented. The same stages as in the linear case are present in the nonlinear case.

### Initialization Stage

The initialization stage for nonlinear equations is mostly equivalent to the linear case. The only difference in the nonlinear case is that, like for FindRoot, a starting seed needs to be specified. This is done by setting the solution data component "DependentVariables" to the initial seed.

In[8]:= |

In[10]:= |

The initial seed can be a numeric value or any expression based on the independent variables, and in this case. For highly nonlinear problems, one approach is to solve a less nonlinear problem and use the resulting InterpolatingFunction as a initial seed for the highly nonlinear problem. If no initial seed is specified, is taken as the initial seed.

In[11]:= |

### Discretization and Solution Stage

Solving nonlinear equations is an iterative process. At the lowest level, the iteration process is done by FindRoot. The first step is to linearize the PDE coefficients. The details of this process follow. The initial seed is then used to evaluate the linearized PDE. This creates a linear system of equations that can be solved with LinearSolve. The result of that process is fed back into the PDE, and the updated linearized PDE is discretized with the new approximation for the solution. This process is repeated until, ideally, the solution converges.

The function PDESolve orchestrates the call of FindRoot, takes care of the linearization and the repeated discretization and returns a solution data object that contains the solution or, if the algorithm failed to converge, the best approximation to the solution PDESolve could find.

In[14]:= |

The solution that is found is stored in the solution data object returned.

### Post-Processing

An InterpolatingFunction object can be created from the solution data.

### The Linearization Process

FindRoot by default uses the Newton–Raphson method to solve nonlinear equations. The way the method works is that it starts with an assumed solution, the initial seed, and tries to improve on it by using the linearized form of the equations.

To understand the linearization process better, it is helpful to review the derivation of the Newton method for a scalar function. Let be a nonlinear function of and the unknown solution such that

Now take an arbitrary seed and subtract that from the unknown solution :

which results in a residual . After rearranging the equation to

Taylor expand around and obtain:

Ignoring higher-order terms you arrive at:

Assuming that the derivative of is not zero, you can get to an iteration process:

Generalizing this for an iteration count , you obtain:

We look at an example of two algebraic equations.

In[17]:= |

In[19]:= |

In[20]:= |

This procedure is continued until a solution is found. As a next step, the iterative procedure is wrapped in a function that computes a step given a vector. While the norm of two consecutive steps is above a threshold, a new next step is computed.

In[23]:= |

One thing to note in the given example is that both and the Jacobian have constant parts and parts in them that depend on . When solving a nonlinear PDE with the finite element method, it makes sense to compute constant parts once prior to the iterative process and to add them to the nonconstant parts that are computed in each step.

The linearization process for the nonlinear PDE in the coefficient form given in equation 2 follows the same steps. Start from equation 3 and set and to obtain a time-independent coefficient form:

For convenience, rearrange the equation to

and rename the coefficients to

This reformulates the time-independent coefficient form to

To perform the linearization, now proceed in the same manner as for the scalar case above. and are functions of , the unknown solution, such that

Now take an arbitrary seed and subtract that from the unknown solution :

which results in a residual . After rearranging the equation to

Expand with a Taylor series around and obtain:

Ignoring higher-order terms, one arrives at:

For and , the derivatives with respect to and are computed:

The actual solution of the linearized PDE makes use of the coefficient form given in equation 4 with the following coefficients for an arbitrary number of dependent variables:

where are integers. where is the number of dependent variables. For scalar equations in , .

After the explanation of Newton's root-finding method and how it applies to partial differential equations, it is instructive to continue to work through the example of solving

at a deeper level, showing how the call to PDESolve presented in the section Discretization and Solution Stage works. In a first step, the linearized coefficients are constructed.

In[31]:= |

The partial derivative coefficients for in equation 7 can be computed by using D.

The partial derivative coefficients for in equation 8 can be computed by using D.

The partial derivative coefficients for in equation 9 can be computed by using D.

The partial derivative coefficients for in equation 10 can be computed by using D.

These coefficients are conveniently computed with a utility function LinearizePDECoefficients from the finite element context.

### Solving the Linearized PDE

Once the PDE is linearized, the coefficients are split. This is done because the coefficients are needed in different parts of the root-finding algorithm. The coefficient form

All terms on the left will be assembled into the stiffness matrix. The stiffness matrix is the Jacobian for the PDE. The stiffness matrix is analog to the term from the preceding scalar example. Since is a tangent to , the stiffness matrix is also called the tangent stiffness matrix. The terms on the right of the PDE correspond to . The splitting of the PDECoefficientData is done with function SplitPDECoefficients.

The coefficients and are now in the variable linLoadPDEC while the rest is in linStiffnessPDEC. There are no damping and mass matrix coefficients for this particular problem, so linDampingPDEC and linMassPDEC are empty.

In[47]:= |

In[49]:= |

In[51]:= |

The function femRHS given here also takes care of generalized Neumann values even though none are present in the example under consideration. In fact, the function given here is very close to the actual implementation.

Note that in function femRHS only the load vector is made use of. This is because linLoadPDEC, generated with the SplitPDECoefficients command, only contains the load coefficient and the load derivative coefficient .

Once the variables nonlinear and nonlinearLoad are no longer needed, they are set to Null to free the memory they point to.

In[52]:= |

Again, the function femJacobian is very close to the actual implementation and variables no longer needed are set to Null.

For the construction of the Jacobian, the linearized stiffness components are used. These include the diffusion coefficient , the conservative convection coefficient , the convection coefficient and the reaction coefficient . All of these are stored in the linStiffnessPDEC PDECoefficientData generated with the SplitPDECoefficients command.

It is noteworthy that in the implementation, all coefficients that are linear will only be assembled once. In both the functions femRHS and femJacobian, only the nonlinear coefficients will be computed and added to the previously computed linear components. This avoids unnecessary computations of linear coefficients that remain unchanged in each iteration. Last, the nonlinear and linear boundary conditions are deployed.

In[53]:= |

The last step is to find the solution to the linearized PDE. For this an affine covariant Newton method is used.

In[54]:= |

In[55]:= |

An InterpolatingFunction object can be created from the solution data.

## 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 Conditions—Introduction

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.

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.

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.

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

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

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.

In[13]:= |

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

In[15]:= |

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

In[77]:= |

### Transient PDE with Stationary Coefficients and Stationary Boundary Conditions

Consider the following model PDE

subject to the boundary conditions Null at Null and Null at Null. 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.

In[80]:= |

For more information about the relation between BoundaryMeshRegion objects and ElementMesh objects, please consult the "Element Mesh Generation" tutorial.

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.

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.

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 only needs to be done once, and the solvers are used for each time step, it is often advantageous to use node reordering for transient problems.

In this case, as there is no ImplicitRegion specification, the ElementMesh object may be directly used to create the NumericalRegion object.

Populate the solution data. The "Time" component is the initial time t_{0} where the integration starts. Setting the "DependentVariables" and "TemporalDerivatives" components sets the initial condition and derivative of the initial condition.

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

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.

In[95]:= |

In[96]:= |

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.

In[101]:= |

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.

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.

In[107]:= |

A single interpolating function for time and the spatial variables can be constructed. For this, the time steps at which the transient interpolation has been set up need to be extracted.

Next, the interpolation values are set up. For each time step, a list of values at the nodes is required.

In[110]:= |

Next, the actual interpolation function is constructed.

### 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, 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 ×1. Next, establish a second, much smaller system of equations with degrees of freedom:

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

In[113]:= |

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

### 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 is 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 .

The initialized coefficients store the coefficients specified.

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

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

There are no transient matrix components in this example.

The "Stationary" 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.

In the above example, the string argument "Stationary" 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.

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.

In[135]:= |

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 preceding 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 cannot 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.

In[136]:= |

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.

The visualization is done in the same manner as before.

## 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 in the following.

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

In[2]:= |

In the plane stress operator, Y 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.

The boundary conditions are as follows: at the left-hand side, the beam is held fixed in both the x and y directions. 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 y direction.

In[7]:= |

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

In[10]:= |

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

In[11]:= |

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

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.

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.

In[15]:= |

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

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.

The initialized PDE coefficients are transformed into a discrete form.

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

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

The system of equations can then be solved with LinearSolve.

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.

### A Swinging Beam—Transient Coupled PDEs

For a transient coupled PDE, not much more is needed. The following example shows how to deflect the beam and have 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 is represented as follows:

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

The transient coefficients also have the mass matrix coefficients set.

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

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.

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.

In[36]:= |

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.

In[38]:= |

In[44]:= |

### A Swinging and Dynamically Loaded Beam

In this section, the previous example will be taken one step further by adding a time-dependent body load and a time-dependent boundary condition.

To do so, the variable and solution data need to have the temporal variable set.

In[46]:= |

In[47]:= |

For the new PDE model, first add a term acting on the body load to model expansion and contraction of the beam in the x direction with a sinusoidal motion.

Secondly, add a time-dependent boundary force that is going to be applied in the y direction.

In[50]:= |

Initialize the PDE coefficients and boundary conditions.

Precompute the discretized stationary PDE coefficients and boundary conditions.

Now write a helper function that will get the current time, , and and values. These are combined in single vectors and . The current time needs to be updated in the solution data. The stationary system matrices are extracted from the discretized PDE. It is important to realize the these are not recomputed, as they do not change in each time step. Next, the transient components of the PDE and boundary conditions are discretized. The transient and stationary system matrices are added together, and the Rayleigh damping is recomputed. Following that, both stationary and transient boundary conditions are deployed. The function returns the discretized version of .

In[55]:= |

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.

In[56]:= |

For the specification of the sparsity pattern, use the stationary discretization of the PDE. This is reasonable, since the transient component is only in the load vector. In other words, the following is only correct if there are no time-dependent components in the mass and stiffness matrix.

In[58]:= |

In[61]:= |

In[64]:= |

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

In[2]:= |

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

- Precompute geometric data (InitializePDEMethodData).

- Discretize the PDE (DiscretizePDE).

- Solve the system of equations (LinearSolve).

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.

In[3]:= |

In[4]:= |

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 "MeshOrder" 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 "MeshOrder" also reduces the quality of the solution that can be expected, as has been discussed in "Solving Partial Differential Equations with Finite Elements".

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 "MeshElementBlocks" 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 "MeshElementBlocks" too large will slow down the assembly process.

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 NDSolve`ProcessEquations is used solely for equation processing, the geometric preprocessing can be switched off.

In[14]:= |

In[15]:= |

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.

In[23]:= |

Next, there is a very subtle point to consider. The system matrices are stored in the DiscretizedPDEData data object. Those have been extracted in the preceding 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.

In[27]:= |