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 NDSolve`FEM` 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.

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

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.

Set up the NDSolve`StateData object.

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.

Get the finite element data.

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.

Compute the system solution.

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

Extract the solution values.

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

Extract the solution as an InterpolatingFunction and plot the result.

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.

The same is possible for transient PDEs.

Solve a transient PDE with NDSolve and finite element options specified.

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. The usage and purpose of the internal finite element method function options are documented on the respective reference pages and in a dedicated tutorial for Finite Element Method 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:

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:

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.

Create the variable data with dependent and independent variable names.

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

Specify a NumericalRegion.
Create the solution data with the "Space" component set.
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.

The coefficients of the model PDE.

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

Initialize the partial differential equation coefficients.

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.
Extract the raw 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.

Initialize the boundary conditions.

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.

Initialize the finite element data with the variable and solution data.

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, interpolation order and integration order.

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

Initialize the finite element data with the variable and solution data, with options for the mesh generation.
Query degrees of freedom, interpolation order and integration order.

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

Extract the ElementMesh object from the NumericalRegion.

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.

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

Increase the integration order.

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 NDSolve`ProcessEquations generates the same FEMMethodData object.
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.

Extract the finite element data from the NDSolve`StateData.
Inspect that the objects created are the same.

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.

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

Extract all system matrices.
Visualize the stiffness matrix.
Extract the system matrices separately.

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

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.

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.

Components from the DirichletCondition boundary condition.

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.

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.

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.

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.

Update the solution data sd with the solution.
Create an InterpolatingFunction object.

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.

Create the interpolation function with values at the nodes.
Visualize the interpolation function.

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

Evaluate outside the domain using extrapolation.
Switch off automatic extrapolation in interpolation functions by having it return Indeterminate.

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

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

Extract the ElementMesh from an 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.

Set up a nonlinear equation.
Set up a region and boundary conditions on the region.

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.

Solve the nonlinear equation with an initial seed specified.

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

Visualize the solution.

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

Solving the activated nonlinear equation results in a message.

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.

Set up the numerical region and variable data.
Set up the solution data with an initial seed of .

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.

Specify the method data, boundary conditions and PDE coefficients.

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.

Use PDESolve to solve a nonlinear PDE.

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.

Visualize the solution.

The Linearization Process

FindRoot by default uses the NewtonRaphson 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:

with .

Look at an example of two algebraic equations.

.

Set up the nonlinear function .
Compute the Jacobian .
Set up a function to evaluate the Jacobian.
Specify an initial seeding .
Compute the first Newton step from the initial seeding and save the result in .
Compute the second Newton step from and save the result in .

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.

Set up a function to compute the Newton steps.
Solve the nonlinear problem.
Check that the computed result solves the nonlinear equations to the requested accuracy.

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

and

.

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.

Extract dependent and independent variables.
Set up variables for the number of dependent and independent variables.
Construct a list of the dependent variables.
Construct a list of the derivatives of the dependent variables.
Extract the coefficients from the PDECoefficientData object, including any nonlinearity.
Construct from equation 5.
Construct from equation 6.

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

Inspect the linearized PDE diffusion coefficients .

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

Inspect the linearized PDE conservative convection coefficients .

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

Inspect the linearized PDE convection coefficients .

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

Inspect the linearized PDE reaction coefficients .

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

Compute the linearized PDE coefficients.
Inspect the linearized diffusion coefficients.
Inspect the linearized load-derivative coefficients.

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

is regrouped such that

.

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.

Split the linearized PDE coefficients.

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.

Evaluate the initial seed on the mesh.
Discretize the linear components of the linearized PDE.
Discretize the stationary boundary conditions.
Specify a right-hand-side function.

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.

Specify a function to compute the Jacobian.

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.

Deploy the boundary conditions into the initial seeding.

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

Use FindRoot to solve the nonlinear PDE.
Update the solution data.

An InterpolatingFunction object can be created from the solution data.

Using ProcessPDESolutions to create an InterpolatingFunction object.
Visualize the solution.

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

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.

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

Specify a NumericalRegion.
Create the solution data with the "Space" and "Time" component set.

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

Initialize the partial differential equation coefficients.
Initialize the boundary conditions.
Initialize the finite element data with the variable and solution data.
Extract the ElementMesh from the NumericalRegion.

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.
Discretize the initialized boundary conditions.
Extract all system matrices.
Deploy the boundary conditions in place.

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.
Time integrate the system of equations with NDSolve.

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.

Set up a function that given a time , constructs and memorizes an interpolating function.
Visualize the difference between the automatic and manual solutions.

Transient PDE with Stationary Coefficients and Stationary Boundary Conditions

Consider the following model PDE

subject to the boundary conditions u[t,x,y,z]0 at ElementMarker1 and u[t,x,y,z]1 at ElementMarker2. 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.
Convert the BoundaryMeshRegion to a boundary ElementMesh.

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.

Visualize the point markers according to their grouping.

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.

Create a full first-order mesh with optimized node connectivity.

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.

Set up variable data.

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

Set up solution data.

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

Initialize the PDE coefficients.
Initialize the PDE boundary conditions.
Initialize the method data.
Discretize the system.
Extract all system matrices and deploy the boundary conditions.

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

Visualize the stiffness matrix.

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.
Alternatively, find the stationary solution using a different solver.
Visualize the solution as a surface plot.
Construct an interpolation function that returns Indeterminate if queried outside of its domain and does not warn.
ContourPlot through the x direction.
Visualize the solution as a 3D contour plot within the model.

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.
Visualize the initial conditions.

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.
Find the solution to the transient PDE with NDSolve as a time integrator.

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.

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.

Extract the times at which the interpolating function was constructed.

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

Extract the solution values at each time step.

Next, the actual interpolation function is constructed.

Construct the single interpolating function.
Inspect the single interpolating function at a point in space for all time steps.

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:

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.

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.
Extract the reduced system matrices and inspect the dimensions.
Efficiently time integrate the reduced system.
Compare the norm of the full system and the reduced 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 depends on the time variable .

Initialize the PDE coefficients.

The initialized coefficients store the coefficients specified.

Extract the stored coefficients.

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.

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

Extract the transient scalar coefficients from all coefficients.

There are no transient matrix components in this example.

Extract the transient matrix coefficients from all coefficients.

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.

Discretize and store the stationary components of the PDE.

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.

Discretize the transient components of the system at a specific point in time.

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.
A function to compute the system matrices and redeploy the stationary boundary conditions.

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.

Set up initial conditions based on the boundary conditions.

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.
Find the solution to the transient PDE with NDSolve as a time integrator.

The visualization is done in the same manner as before.

Visualize the reaching of steady state at various time steps.

Transient PDEs with Nonlinear Transient Coefficients

This section considers transient PDEs, where one or more of the coefficients are nonlinear. Consider the following model PDE with a nonlinear damping matrix and a nonlinear diffusion coefficient ,

A nonlinear heat equation is used as a particular example. Consider

with temperature , density , temperature-dependent heat capacity and temperature-dependent conductivity . In order to keep things simple here, and focus as much as possible on the FEM programming aspect the phase change is modeled in a 1D domain and parameters are chosen for simplicity.

The preceding PDE will be used to model phase changes. A common way to model a change from one phase to the next is to make use of a smooth transition function.

Create a smoothed step function.
Visualize the transition of smoothed step function from phase 1 to phase 2.
Set up nonlinear parameters.
Specify a 1D domain.

On the left boundary there is a constant inflow of energy and on the right boundary a fixed temperature is set up. Over time, the initial phase of the material will change from the left to the right end of the domain.

Set up the boundary conditions.

The initialization stage is done in the same manner as shown in the previous section.

Set up the numerical region, the variable and solution data.
Initialize the PDE coefficient, boundary conditions and the method data.

Note that the nonlinear coefficients and can be directly used.

To set up the equation for NDSolve one also proceeds in much the same way as in the previous section. The idea is to extend what in the previous section was called , to also compute the nonlinear coefficients and include the nonlinear damping matrix. This will put all PDE coefficients in a single function called discretizePDEResidual to compute the residual for the DAE solver:

Start, as previously, by computing the stationary components.

Compute the stationary components of the PDE coefficients and boundary conditions.

The discretization function gets the time , the dependent variable and the derivative of the dependent variable as arguments. Then the time-dependent PDE coefficients and boundary conditions are computed and added to the stationary components. After that, the nonlinear components are computed and also added. The function then returns the discretized version of

Create a function to compute the discretized PDE in each time step.
Evaluate the initial conditions such that the values are set to everywhere.

To compute the sparsity pattern, one makes use of a helper function to construct the damping matrix at a specific point in time. You need to create a damping matrix sparsity pattern that is valid throughout the entire integration time. Thus, extract the sparsity pattern of the damping matrix constructed at an arbitrary point in time but that has its time-dependent coefficients active. The exact nonzero time used, however, is not important, since at this point you are only concerned about the pattern of the damping matrix. This is not an issue in this specific case, since the PDE modeled here does not have purely time-dependent coefficients in the damping matrix.

A helper function to construct the damping matrix at a specific point in time.
Compute the sparsity pattern at a time 10^-6.
Solve the equations.
Create an InterpolatingFunction from the time steps.
Visualize the solution.

Transient PDEs with Integral Coefficients

The method presented in the section Transient PDEs with Nonlinear Transient Coefficients can be extended to solve transient integral partial differential equations. The idea is to construct an ElementMeshInterpolation in each time step and make use of NIntegrate to integrate over the current interpolating function. This approach is not limited to integration: in fact, any operation that gives a number can be used on the interpolating function.

As an example, consider the following transient integro-differential equation:

Expand the equation:

Both the domain of the integral and the differential equation are from 0 to 5.

Specify a 1D domain.
Specify a DirichletCondition on the left and right ends of the domain.

The initialization stage is done in the same manner as shown in the previous section.

Set up the numerical region, the variable data and solution data.

Initialize all coefficients except the integral coefficient.

Initialize the PDE coefficient, boundary conditions and the method data.

Start, as previously, by computing the stationary components.

Compute the stationary components of the PDE coefficients and boundary conditions.

In the discretization function named discretizePDEResdual, load the linear components, add the transient and nonlinear components, then construct the current interpolating function. Next, numerically integrate over the constructed function and construct the missing PDE coefficients that contain the integral. Discretize the missing coefficient and add the result to the system matrices. The function then returns the discretized version of:

In order to create an ElementMeshInterpolation in each time step, extract the mesh from the NumericalRegion.

Extract the ElementMesh.
Create a function to compute the discretized PDE in each time step.

Strictly speaking, the transient and nonlinear component discretization of the integral coefficients is not necessary in this case because the integral coefficient is linear. The code includes computation for a transient and nonlinear integral coefficient for generality.

Evaluate the initial conditions such that the values are set to everywhere.

The sparsity pattern is computed in the same manner as before.

A helper function to construct the damping matrix at a specific point in time.
Compute the sparsity pattern at a time 10-6.
Solve the equations.
Create an InterpolatingFunction from the time steps.
Visualize the solution.

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.

A plane stress formulation as a coupled PDE.

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.

Create a mesh region with an explicit bounding box and a visualization of the mesh.

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 pressure of one unit is applied; a NeumannValue enforces a boundary load by specifying a pressure on the beam in the negative y direction.

Set up boundary conditions at the left- and right-hand side of the domain.

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

Populate the variable data with two dependent variables u and v and independent variable names.

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

Set the "Space" component of the solution data.

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.

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 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 the coupled case, there is a diffusion coefficient matrix for each dependent variable in each PDE.

Initialize the partial differential equation coefficients.

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.

The initialized PDE coefficients are transformed into a discrete form.

Compute the discretized partial differential equation.

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

Extract the system matrices.
Visualize the stiffness matrix.

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.

Compute how the incidents are split between the two dependent variables.
Discretize the initialized boundary conditions.
Deploy the boundary conditions in place.

The system of equations can then be solved with LinearSolve.

Solve the system of equations.

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.
Visualize the interpolation function.
Visualize the deformed structure in red.

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 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 additional coefficients of the transient model PDE.

The transient coefficients also have the mass matrix coefficients set.

Initialize the transient partial differential equation coefficients.
Compute the discretized partial differential equation.

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

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

Visualize the mass matrix.

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.
Deploy the boundary conditions.

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.
Set up the sparsity pattern for the Jacobian based on the damping matrix.
Find the solution to the transient PDE with NDSolve as a time integrator.
Visualize the swinging beam.

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.

Set the temporal variable in the variable data.
Set an initial time in the solution data.
Initialize the method data.

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.

An additional time-dependent body load for the transient model PDE.

Secondly, add a time-dependent boundary pressure that is going to be applied in the y direction by specifying a pressure through a NeumannValue.

A time-dependent boundary load.

Initialize the PDE coefficients and boundary conditions.

Initialize the transient partial differential equation coefficients.
Initialize the boundary conditions.

Precompute the discretized stationary PDE coefficients and boundary conditions.

Compute the discretized stationary partial differential equation.
Compute the discretized stationary boundary conditions.

The equation is rearranged:

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 .

A function to compute the transient system matrices and deploy boundary conditions.

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 initial conditions based on the boundary conditions.

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.

Set up the sparsity pattern for the Jacobian based on the damping matrix.
Find the solution to the transient PDE with NDSolve as a time integrator.
Visualize the dynamically loaded swinging beam.

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.

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