FlexPDE uses estimates of the modeling error to control mesh refinement and timestep size. This note describes the methods used and the interpretation of the reports.

Spatial Error

The Galerkin Finite Element method uses integrals of the PDE's to form the discretized equations at the mesh nodes.

Each nodal equation requires that the weighted integral of the associated equation over the mesh cells surrounding the node be satisfied within a convergence tolerance. In FlexPDE this tolerance is taken to be a relative error of (ERRLIM * OVERSHOOT) in the norm of the solution vector.

In a regular hexagonal 2D mesh, for example, the Galerkin method requires that each hexagonal set of six triangular mesh cells must produce a weighted integral residual of zero.

This method at no point imposes any conditions on the integral over a single mesh cell, and conceivably on could have cancelling errors in adjacent cells.

In FlexPDE, we choose to use the individual cell integrals as a measure of the mesh quality. If the aggregate (eg 6-cell) integral is correct but the individual cells show large error, then the mesh must be refined.

The fundamental system which is solved by FlexPDE can be indicated as R=G(U)=0, where R is the residual and G(U) is the Galerkin integral of the PDE for variable U. If the residual over an individual cell is R, we can write J*dU=R, where J is the Jacobian matrix of derivatives of the Galerkin integral with respect to the nodal values, and dU is the error in U which produces the residual R.

J is of course the coupling matrix which is solved to produce the solution U. We don't want to completely repeat the solution process just to get an error estimate, so we use D, the diagonal components of J, to produce the error estimate dU=Inv(D)*R.

The "RMS Error" reported by FlexPDE in the Status Panel is just the root-mean-square average of dU/range(U) over the cells of the problem, while the reported "MAX Error" is the largest error dU/range(U) seen in any cell.

Mesh cells for which dU/range(U) > ERRLIM are split in the mesh refinement pass.

Notice that the error measure is not a guarantee that the computed solution is "accurate" to within the stated error, that is, that the computed solution differs from the "true" solution by no more than the stated error. The error estimate is a local measure of how much variation of the solution would produce the computed error in the cell integral. Deviations from the "true" solution might accumulate over the domain of the problem, or they might cancel in neighboring regions.

Temporal Error

In time dependent problems, an estimate must also be formed of the error in integrating the equations in time.

FlexPDE integrates equations in time using a second-order implicit Backward Difference Formula (Gear method).

In order to measure temporal error, FlexPDE stores an additional timestep of values previous to the three points of the quadratic solution, and fits a cubic in time to the sequence at each node. The size of the cubic term implies the error in the quadratic solution, and is used to either increase or decrease the timestep in order to keep the RMS temporal error within the range specified by ERRLIM.

The three-point integration method requires an independent method to create data for the initial interval. FlexPDE uses a comparison of one-step and two-step trapezoidal rule integration to adapt the initial timestep to a range that produces acceptable error.

The temporal error estimate is not currently reported on the status panel.

FlexPDE Error Controls

There are several SELECT controls that can be used to alter the behavior of FlexPDE in regard to error measures.

The basic control is ERRLIM, which specifies the desired relative error in the solution variables, and controls both spatial and temporal measures. Smaller ERRLIM causes more mesh subdivision and smaller timesteps. Larger ERRLIM allows cruder meshes and, in principle, larger timesteps. However, a large ERRLIM can allow oscillations to develop, ultimately causing severe timestep cuts and a slower overall execution. It is rarely advisable to use an ERRLIM value larger than the default 0.001.

XERRLIM and TERRLIM are analogous to ERRLIM, but refer specifically to the spatial and temporal controls, allowing separate control of the two processes.

PREFER_SPEED and PREFER_STABILITY are two alternative strategies that can be requested.

PREFER_SPEED limits the Newton-Raphson iteration in time-dependent problems to one iteration per timestep, on the theory that the timestep control will keep the step size in the range where this is a valid and stable procedure. It also selects use of the L4 norm (mean fourth power) average of temporal errors over the problem domain. This is the most computationally efficient method, but in strongly nonlinear problems, it can become unstable.

PREFER_STABILITY allows up to three Newton-Raphson iterations per timestep in nonlinear time-dependent problems, and uses the L16 norm (mean 16th power) average of temporal errors over the problem domain. More Newton iterations makes the method more expensive, but more stable in severely nonlinear problems. The higher power norm allows detection of spatially localized activity, before oscillations can get out of hand.

The two components of the "PREFER..." controls can be set independently by the NRUPDATE and TNORM controls.