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.
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.
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. 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:
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.
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.
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.
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.
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.
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.
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:
Look at an example of two algebraic equations.
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.
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.
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.
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.
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.
The last step is to find the solution to the linearized PDE. For this an affine covariant Newton method is used.
An InterpolatingFunction object can be created from the solution data.
The previous outlined procedure is limited in that it does not work if a PeriodicBoundaryCondition is present.
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.
Next, the initial conditions are set up and NDSolve will time integrate the discretized system of equations.
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.
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.
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 t0 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.
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.
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.
Previous versions of NDSolve had the Jacobian specified as "Jacobian{Automatic,Sparsesparsity}", where "sparsity" was the sparsity pattern of the damping matrix. From Version 12.2 moving onwards, this has been superseded by the possibility to specify a sparsity pattern for both the stiffness and the damping matrix. While the old Jacobian setup can still be used, it will result in a slower time integration, and thus the switch to specifying both the Jacobians for the stiffness and damping matrix is highly encouraged.
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.
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.
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
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 depends on 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.
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.
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.
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.
For more information on the physics of a phase change, please refer to the Heat Transfer PDE Modeling tutorial that also shows how phase change can be modeled at the NDSolve level.
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.
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.
The initialization stage is done in the same manner as shown in the previous section.
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.
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
To compute the sparsity pattern, one makes use of a helper function to construct the stiffness and damping matrices at a specific point in time. You need to create stiffness and damping matrix sparsity patterns that are valid throughout the entire integration time. Thus, extract the sparsity patterns of the matrices constructed at an arbitrary point in time, making sure that the possible time-dependent coefficients are active at that point in time. The exact nonzero time used, however, is not important, since at this point you are only concerned about the pattern of the stiffness and damping matrices.
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:
Both the domain of the integral and the differential equation are from 0 to 5.
The initialization stage is done in the same manner as shown in the previous section.
Initialize all coefficients except the integral coefficient.
Start, as previously, by computing the stationary components.
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.
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.
The sparsity patterns are computed 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 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 pressure of one unit is applied; a NeumannValue enforces a boundary load by specifying a pressure on the beam in the negative y direction.
The initialization has the same components as in previous sections. First, the variable and solution data is set up.
For coupled PDEs, more than one dependent variable needs to be given.
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 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.
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.
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.
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 pressure that is going to be applied in the y direction by specifying a pressure through a NeumannValue.
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 
.
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.
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.
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.
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.
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.
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.
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.