Solving Partial Differential Equations with Finite Elements

Introduction

The aim of this tutorial is to give an introductory overview of the finite element method (FEM) as it is implemented in NDSolve. The notebook introduces finite element method concepts for solving partial differential equations (PDEs). First, typical work flows are discussed. The setup of regions, boundary conditions, and equations is followed by the solution of the PDE with NDSolve. The visualization and animation of the solution is then introduced, and some theoretical aspects of the finite element method are presented. Classical PDEs such as the Poisson and Heat equations are discussed. Coupled PDEs are also introduced with an example from structural mechanics.

Why Finite Elements?

Explicit closed-form solutions for partial differential equations (PDEs) are rarely available. The finite element method (FEM) is a technique to solve partial differential equations numerically. It is important for at least two reasons. First, the FEM is able to solve PDEs on almost any arbitrarily shaped region. Second, the method is well suited for use on a large class of PDEs. While it is almost always possible to conceive better methods for a specific PDE on a specific region, the finite element method performs quite well for a large class of PDEs. In summary, the finite element method is important since it can deal with:

  • arbitrarily shaped regions
  • a large class of partial differential equations

What Is Needed for a Finite Element Analysis

To solve partial differential equations with the finite element method, three components are needed:

  • a discrete representation of a region, i.e. a mesh
  • a partial differential equation
  • boundary conditions that link the equation with the region

This section deals with partial differential equations and their boundary conditions. Finite element meshes can be generated with ToElementMesh. A treatise on the generation of finite element meshes can be found in "Element Mesh Generation".

The Scope of the Finite Element Method as Implemented in NDSolve

The current version of the implementation of the finite element method supports the following features:

  • stationary and transient PDEs in 1, 2, and 3 dimensions
  • single and coupled PDEs with derivatives up to second order in space and arbitrary order in time
  • linear PDEs with variable coefficients
  • linear variable coefficient Dirichlet conditions, generalized Neumann (Robin) boundary values, and a mixture of these boundary conditions are possible
  • automatic static mesh generation with curved boundaries
  • first- and second-order meshes with manual refinement options
  • intercept the solution process and access low-level data at any time

The following sections provide an overview of this functionality.

For convenient use of the finite element functionality, load the finite element package.

Load the finite element package.
In[1]:=
Click for copyable input

Regions

To perform a finite element analysis, a region over which the partial differential equation is to be solved needs to be defined.

Rectangular regions may be specified using for each of the spatial independent variables.

Solve a PDE and visualize the solution over a rectangular region.
In[3]:=
Click for copyable input
Out[4]=

Regions of arbitrary shape may be specified using the notation where Ω is a region so that RegionQ[Ω] gives True.

One way to specify a non-rectangular region is by using Boolean predicates. A predicate is a function that either returns True or False. Since the term predicate has many meanings in mathematical literature, the terminology Boolean predicate is used to emphasize that the predicate returns a Boolean. To describe regions, Boolean predicates can be used as in functions like RegionPlot and RegionPlot3D, for example.

Region setup and visualization.
In[5]:=
Click for copyable input
Out[6]=

This tutorial concentrates on solving partial differential equations with the finite element method, without emphasis on the creation of regions and meshes. More detailed information on this topic can be found in "Element Mesh Generation".

Classical Partial Differential Equations

The Coefficient Form of Partial Differential Equations

What types of equations can be solved with the finite element method as implemented in NDSolve? Consider a single linear partial differential equation in :

The PDE is defined in . Here is the dependent variable for which a solution is sought. The coefficients , , , and are scalars; , , and are vectors; and is an × matrix.

What follows are some well-known PDEs and their corresponding coefficients. To illustrate the generality of (1), the components that are relevant to a specific equation are black, while the non-relevant components are gray.

The Laplace equation simply contains a diffusive term:

To model Poisson's equation, only a small modification is needed; add a load term :

Helmholtz's equation adds a reaction term :

Convection-diffusion-reaction type equations are another common class of PDEs. Compared to the previous examples, these have an additional convection term :

The PDEs considered so far are stationary, i.e. they have no time dependence. The heat equation adds time dependence to the Poisson equation. It has the following form:

Similarly, the wave equation is given as:

Equation (1) provides the components for modeling a range of different phenomena, since it provides spatial derivatives up to order 2.

Poissons Equation with Dirichlet Conditions

The Poisson equation is to be solved over a region with Dirichlet boundary conditions. First, a region needs to be defined where the equation will be solved. Then the equation and boundary conditions are defined. Finally, the equation is solved over the region.

One way to specify a region is by using Boolean predicates.

Region setup and visualization.
In[7]:=
Click for copyable input
Out[8]=

Next, define a Poisson operator.

PDE operator setup.
In[9]:=
Click for copyable input

The third step is to set up the boundary conditions. The first Dirichlet boundary condition enforces the dependent variable to a value of 0 wherever the evaluation of yields True on the boundary of the region . In this case, this is the upper-left part of the boundary. The second Dirichlet boundary condition specifies the value of the dependent variable whenever yields True on the boundary of region .

Dirichlet boundary conditions setup.
In[10]:=
Click for copyable input

In the solution step, NDSolve needs the PDE, the boundary conditions, and the region.

Solve the PDE.
In[11]:=
Click for copyable input
Out[11]=

Internally, the region is converted to a finite element mesh and a finite element method is used to solve the equation automatically.

To visualize the solution, a contour plot can be used.

Visualize the solution.
In[12]:=
Click for copyable input
Out[12]=

Partial Differential Equations and Boundary Conditions

Two primary types of boundary conditions are used. Dirichlet boundary conditions prescribe a condition on the dependent variable of value or at a part of the boundary of :

On the other hand, generalized Neumann boundary values, also known as Robin values, specify a value . The value prescribes a flux over the outward normal on some part of the boundary:

The finite element method is based on the weak form of the differential equation. This form is obtained by taking equation (1), multiplying it by a so-called test function , and integrating over the region :

Integration by parts gives:

This process is done internally. The integrand in the boundary integral in (11) can be replaced with (9) and yields:

specifies the value of the boundary integral integrand (11) in the weak form and thus the name NeumannValue. Other boundary conditions are conceivable, but currently not implemented.

In most cases, Dirichlet boundary conditions need not be associated with a particular equation and can be specified independent of the equation. They are specified as equations. Generalized Neumann values, on the other hand, are specified by giving a value, since the equation satisfied is implicit in the value. Neumann values are mathematically tied to the PDE.

For practical reasons, the NeumannValue is part of the equation. Consider a system of coupled PDEs. One would like to be able to unambiguously specify any given Neumann value to any given single PDE of that system of PDEs. It is not possible to derive (unambiguously) from the Neumann value with which PDE equation the value should be associated. Making a NeumannValue part of a PDE equation solves this problem. By making the NeumannValue part of the equation, the association that Neumann values belong to with the PDE equation becomes unambiguous.

{op[u[x1,]] ΓN,ΓD}partial differential equation with differential operator op, Neumann boundary values , and Dirichlet boundary condition equations
{op1[u[x1,],v[x1,],] ΓN1,op2[u[x1,],v[x1,],] ΓN2,,ΓD}system of coupled partial differential equation operators , Neumann boundary values , and Dirichlet boundary condition equations

Specifying partial differential equations with boundary conditions.

Both DirichletCondition and NeumannValue also need as a second argument a predicate describing the location on the boundary where the conditions/values are to be applied.

A partial differential equation typically needs at least one Dirichlet boundary condition on some part of the region to be uniquely solvable. For this reason, Dirichlet boundary conditions are also called essential boundary conditions.

The use of generalized Neumann values arise from the boundary integral from the weak form, and so they are also commonly referred to as natural conditions.

A pure Neumann boundary equation is a generalized Neumann boundary equation (9) where :

If is also not specified, then a Neumann zero boundary value is automatically implied ():

Note that when the Neumann value is zero, the term drops out in (12). This is typically a case where there is no flux or stress at the boundary. When no boundary conditions are explicitly specified for a part of the region boundary , a Neumann zero value is assumed.

Poissons Equation with Neumann Values

Consider the solution of Poisson's equation with generalized Neumann boundary values.

Use that same region and the PDE operator.
In[13]:=
Click for copyable input

Poisson's equation is given as:

The corresponding Neumann boundary equation is:

Note how the normal derivative of the Neumann type equation relates to the PDE. In the Neumann boundary equation, you specify the values of and on the right-hand side. The left-hand side is implicitly specified though the PDE. If the coefficient in the PDE is given as an identity matrix, then the coefficient on the left-hand side of the Neumann equation is that same coefficient .

The Neumann boundary value is specified though the values and . Typically when , one speaks of a Neumann boundary value, and in the case , one speaks of a generalized Neumann or Robin boundary value. It is common to use the term Neumann boundary condition for both the generalized Neumann and Neumann boundary value.

Specify a generalized Neumann value and Dirichlet boundary condition.
In[15]:=
Click for copyable input
In[16]:=
Click for copyable input

It is important to observe that the constants in the term of the Neumann boundary equation exactly match the constants in from the weak form of the PDE (11). Because these terms are required to match, NeumannValue only requires the user to enter the right-hand side of the boundary equation. For example, to enter the boundary equation , it suffices to enter NeumannValue[1,]. The correct expression for the left-hand side of the boundary equation (i.e. ) is implied by the term in the weak form of PDE (12). This close relation of the Neumann boundary value and the PDE is further accommodated by making the NeumannValue part of the equation.

Solve the PDE.
In[17]:=
Click for copyable input
Out[17]=
Plot the solution.
In[18]:=
Click for copyable input
Out[18]=

In the previous examples, not all boundaries had boundary conditions explicitly attributed; for those parts of the boundaries, the Neumann zero condition is used.

Formal Partial Differential Equations

In the preceding example, the PDE is set up by making use of Laplacian. The exact same PDE can be formulated with Div and Grad. To do so, use the coefficient in . Here is an × matrix where is the number of dimensions in .

In[2]:=
Click for copyable input
Out[2]=

To prevent evaluation, Inactive can be used.

Set up the PDE operator.
In[3]:=
Click for copyable input
Out[4]=

The formal version of the PDE allows you to easily modify the coefficient matrix c. The class of PDEs NDSolve can solve with finite element method, as currently implemented, do not require Inactive operators.

Also, note that while it is possible to specify a coefficient as a function, this deprives the internal code of the ability to optimize the coefficient integration. In this specific case, the coefficient is constant and with the optimization the coefficient is only integrated once. A black box function would need continuous evaluations over the region.

Partial Differential Equations with Variable Coefficients

The PDE coefficients may depend on space and time. Here indicates the temporal variable and the spatial variables:

The example of a shielded micro strip illustrates the effect of variable coefficients. It is an electrostatics application that can be modeled with a Poisson equation where the coefficient is the relative permittivity and depends on space:

We construct a geometry where the boundary separates two materials. In each material, has a different value.

First, a boundary mesh is created. For more information about creating element meshes, please consult ToBoundaryMesh.

Create a boundary mesh.
In[27]:=
Click for copyable input
Out[28]=
Visualize the boundary mesh.
In[29]:=
Click for copyable input
Out[29]=

The conversion of the boundary mesh to an element mesh is done with ToElementMesh.

Convert the boundary mesh to a full mesh and visualize the resulting mesh.
In[30]:=
Click for copyable input
Out[31]=

Next, the PDE is set up. The load term of the PDE is the charge density divided by the vacuum permittivity :

The lower part of the geometry is silicon, while the upper part is vacuum. Silicon has a relative permittivity of about 11.7, while vacuum has a value of 1.

In[32]:=
Click for copyable input
Out[32]=

The charge density and vacuum permittivity are given as follows.

In[33]:=
Click for copyable input
Set up the partial differential operator.
In[35]:=
Click for copyable input
Out[35]=

The potential is set to 0 wherever . A second potential is set to 1000 in the inner strip.

Set up the boundary conditions.
In[36]:=
Click for copyable input
Solve the PDE.
In[37]:=
Click for copyable input
Out[37]=
Visualize the result.
In[38]:=
Click for copyable input
Out[38]=

Heat Equation

The heat equation is a time-dependent Poisson equation

where the dependent variable depends on the spatial coordinates and time .

Both Dirichlet boundary conditions and Neumann boundary values may also depend on time.

The overall procedure to solve PDEs remains the same: a region needs to be specified and a PDE with boundary conditions needs to be set up.

Set up the spatial region.
In[39]:=
Click for copyable input

The dependent variable now depends on time and space.

Set up the PDE operator, being sure to include the time as an argument of .
In[40]:=
Click for copyable input

Also, the boundary conditions may depend on space and time.

Set up the boundary conditions.
In[41]:=
Click for copyable input
In[42]:=
Click for copyable input

To solve first-order time-dependent problems, NDSolve needs an initial condition and a time region.

Solve the PDE.
In[43]:=
Click for copyable input
Make an animation of the solution.
In[44]:=
Click for copyable input
Out[44]=

Note that it is more efficient to convert the mesh once and use the discretized mesh during the plotting instead of rediscretizing the region for every animation frame.

Wave Equation

The wave equation is a second-order time-dependent equation:

Define circular region .
In[45]:=
Click for copyable input
Define the wave equation PDE operator.
In[46]:=
Click for copyable input
Dirichlet boundary conditions are set to a fixed value of on all boundaries.
In[47]:=
Click for copyable input

The wave equation is a second-order time-dependent PDE, and initial conditions and the derivative of the initial condition must be specified.

Initial conditions for the wave equation.
In[48]:=
Click for copyable input
In[49]:=
Click for copyable input
Solve the PDE and have NDSolve find initial values.
In[50]:=
Click for copyable input

To make an animation of the resulting InterpolatingFunction, use Plot3D and Animate:

Make an animation of the solution.
In[51]:=
Click for copyable input
Out[51]=

Systems of Partial Differential Equations

The Coefficient Form of Systems of Partial Differential Equations

Consider a system of two equations with dependent variables and :

For a region in , each is an × matrix, and each of , , and are vectors.

In this case, the generalized Neumann equations are

One-Dimensional Coupled Partial Differential Equation

The following example is set up such that the parameters are easy to modify. First, the region is specified.

Specify a region.
In[52]:=
Click for copyable input
Out[53]=

Here is the PDE to be solved:

The PDE has two dependent variables and . The coupling is between two diffusion operators. Note that since this is a one-dimensional system, each coefficient is a 1×1matrix, which is equivalent to a number.

Set up the PDE operators.
In[54]:=
Click for copyable input

Each equation has its own set of boundary conditions. In the most general form, Dirichlet boundary conditions are now systems of conditions:

The current implementation does not support cross-coupling in Dirichlet boundary conditions. The terms and cannot be given.

Neumann values specify a system of the following form:

Next, the boundary conditions are set up.

Specify boundary conditions for the first equation.
In[57]:=
Click for copyable input
In[59]:=
Click for copyable input
Specify boundary conditions for the second equation.
In[61]:=
Click for copyable input
In[63]:=
Click for copyable input
Solve the system of equations.
In[65]:=
Click for copyable input

Here is the analytical solution of the system of PDEs computed with DSolve.

The analytical solution.
In[66]:=
Click for copyable input
Out[66]=

Note how in the DSolve formulation of the system of equations the generalized Neumann boundary equations have to be given explicitly, and how their coefficients match those in the PDE.

To compare the analytical result to the numeric result, make a plot of the difference between the two.

Plot the difference between the analytical and numerical solutions.
In[67]:=
Click for copyable input
Out[67]=

This example also makes apparent why the NeumannValue is part of the equation. The equations for the Neumann value are given by:

Consider the boundary value :

In[68]:=
Click for copyable input
Out[68]=

There is no way to know from the Neumann value with which equation it should be associated. Making the Neumann value part of the equation disambiguates with which equation the Neumann value is associated.

Structural Mechanics

In structural mechanics, the deformation of objects under load is computed. Structural mechanics provide good examples of coupled PDEs. As an example region, consider a beam five units long and one unit high.

The PDE that models the physics is a plane stress relation and is valid for thin objects. Here, "thin" means thin relative to the other dimensions of the object. The plane stress formulation depends on two dependent variables and the independent variables. The two dependent variables give the component-wise contribution of the deformed object. As material data, the Young's modulus and the Poisson ratio need to be specified.

Define a plane stress operator.
In[69]:=
Click for copyable input

Furthermore, a boundary load may be specified. In this case, the boundary load applies at the right edge of the beam and a force of unit is applied in the negative direction.

On the left-hand side of the region, both dependent variables are held fixed with 0 displacement.

Set up the boundary conditions so that the beam is held fixed at the left.
In[71]:=
Click for copyable input
Solve the equation.
In[72]:=
Click for copyable input

In order to show the deformed beam, an element mesh must be created for the beam at rest, which is deformed by the interpolating functions given.

Visualize the boundary of the beam at rest versus the deformed mesh of the beam.
In[73]:=
Click for copyable input
Out[74]=

More information about the usage of can be found in "Element Mesh Visualization".

The following contour plots show the deformation of the beam. The dependent represents the deformation in the direction. It can be seen that the top-right part of the beam experiences a shift in the positive direction. The lower-right part is moved in the negative direction.

In[75]:=
Click for copyable input
Out[75]=

The contour plot of the dependent represents the deformation in the direction. The right-hand side of the beam undergoes a deformation in the negative direction.

In[76]:=
Click for copyable input
Out[76]=

Note that the color scale of the two contour plots is not the same.

For completeness, a PDE that models the physics a plane strain relation is presented. The plane strain formulation is, in contrast to the plane stress formulation, valid for thick objects. Here, "thick" means thick relative to the other dimensions of the object. The plane strain formulation depends on two dependent variables and the independent variables. The two dependent variables give the component-wise contribution of the deformed object. As material data, the Young's modulus and the Poisson ratio need to be specified.

Define a plane strain operator.
In[77]:=
Click for copyable input

The boundary conditions this time are prescribed displacements on both sides. On the left-hand side of the region, both dependent variables are held fixed with 0 displacement. On the right-hand side, a displacement of one unit in the direction of negative axes is given and the displacement on the direction is held fixed.

Set up the boundary conditions so that the beam is held fixed at the left and the right.
In[79]:=
Click for copyable input
Solve the equation.
In[80]:=
Click for copyable input

In order to show the deformed beam, an element mesh must be created for the beam at rest, which is deformed by the interpolating functions given.

Visualize the boundary of the beam at rest versus the deformed mesh of the beam.
In[81]:=
Click for copyable input
Out[82]=

The following contour plots show the additive deformation of the beam. The dependent represents the deformation in the direction and the dependent represents the deformation in the direction.

In[83]:=
Click for copyable input
Out[83]=

Utilities for Solving Partial Differential Equations

The next sections contain miscellaneous utilities for solving partial differential equations.

Monitoring Progress of Time Integration of Transient Partial Differential Equations

Time integration of partial differential equations can take time. This is especially true if the boundary conditions or coefficients are time dependent or if the PDE region at hand is three-dimensional. An easy way to monitor progress of an integration of a transient PDE is by using an EvaluationMonitor.

Using the EvaluationMonitor within NDSolve works as in other cases. Here is an example wave equation with a time-dependent Dirichlet condition.

Specify an ellipsoidal region with a circular hole in one focal point.
In[84]:=
Click for copyable input
Create and visualize a mesh of the region.
In[87]:=
Click for copyable input
Out[88]=
Specify a transient Dirichlet boundary condition on the circular hole that is active until time .
In[89]:=
Click for copyable input
Set up the notebook's status area for displaying the current integration time.
In[91]:=
Click for copyable input
Time-integrate the PDE with the status displayed in the status area.
In[93]:=
Click for copyable input
Out[93]=

Once the initialization of the PDE is done, the status area will display the current time during the time integration.

Cleanup of the status area.
In[94]:=
Click for copyable input

Solving Memory-Intensive PDEs

The amount of memory needed by NDSolve to solve a particular problem with the finite element method is proportional to the number of equations and the space dimension operated in. For coupled systems in 3D, the amount of memory provided by the underlying hardware can be a barrier. The following section identifies easy-to-use options to alleviate this potential barrier.

During the finite element analysis there are two key memory bottlenecks. The first is the creation of the discretization and the second is the solution (LinearSolve) of the system of equations. NDSolve provides options to both mesh generation and the LinearSolve step that have an effect on the memory requirement during discretization and solving. As a rough estimate, the number of coordinates in the internal generated mesh multiplied by the spatial dimension gives the degrees of freedom (number of rows and columns) of the discretized system matrix to be generated.

As an example, consider the NavierCauchy equation from elastostatics. The NavierCauchy equation is a system of three coupled PDEs in 3D space. Since the main focus of this section is to create a large system of equations, no additional information about the PDE itself is presented.

Define the NavierCauchy operator.
In[2]:=
Click for copyable input
Define the Lamé parameters.
In[3]:=
Click for copyable input

A series of measurements illustrates how various options influence the solution time and the amount of memory needed to find a solution. For this purpose, a helper function is implemented. The function receives as input any number of NDSolve options and returns the time and the amount of memory in MB needed to solve the PDE. It is also instructive to get a feeling for the size of the system matrices NDSolve internally establishes. The size of the system matrices are its degrees of freedom, or in other words, the number of rows and columns in the system matrices.

Create a helper function that calls NDSolve with options. The measured time, memory needed, and degrees of freedom to find the solution are printed.
In[4]:=
Click for copyable input
Click for copyable input

The default settings in NDSolve are set up such that the solution is accurate and found fairly efficiently.

Calling NDSolve with default settings.

With the default setting, NDSolve will call LinearSolve with no special options. This will route to a direct solver that is very accurate but has a high requirement of available memory.

Visualize the deformed structure.
In[6]:=
Click for copyable input
Out[8]=

One way to reduce the memory consumption is to have NDSolve generate a mesh with fewer elements.

Instruct NDSolve to use a smaller mesh.

This has reduced both solution time and memory requirements at the price of accuracy. How many elements a mesh needs to represent an accurate solution is problem dependent. A good method is to start with a coarse mesh and refine until the solutions stop to differ substantially.

Visualize the difference of the computed results at a cross section.
In[10]:=
Click for copyable input
Out[11]=
Visualize the difference of node distances as an Image3D rescaled to the range of 0 to 0.05.
In[12]:=
Click for copyable input
Out[19]=

A further way to reduce the amount of memory needed for computing the solution is to use a mesh that is only first-order accurate.

Instruct NDSolve to use a first-order accurate mesh.
Visualize the difference of the computed results at a cross section.
In[21]:=
Click for copyable input
Out[22]=
Visualize the difference of node distances as an Image3D rescaled to the range of 0 to 0.05.
In[23]:=
Click for copyable input
Out[25]=

This has reduced both solution time and memory requirements at the price of accuracy. Making use of a first-order mesh results in much smaller system matrices, which can be solved more efficiently. Using a first-order mesh may significantly reduce accuracy.

The previous suggestions worked by reducing the size of the system matrices internally generated. A different approach is to instruct NDSolve to use a different solver method for LinearSolve. All options that LinearSolve take can be directly given to NDSolve.

Instruct NDSolve to use the solver.
Visualize the difference of the computed results at a cross section.
In[27]:=
Click for copyable input
Out[28]=
Visualize the difference of node distances as an Image3D rescaled to the range of 0 to 0.05.
In[29]:=
Click for copyable input
Out[31]=

While this takes longer to solve, it uses substantially less memory while at the same time retaining the accuracy compared to the default solution. The reason that this method is not the default is that the default LinearSolve method is guaranteed to converge if the discretized PDE has a solution. This benevolent property is not true for all solvers implemented in LinearSolve.

The options discussed above can also be combined. A more in-depth treatise of options that can be given to NDSolve can be found in the finite element tutorials.

Extrapolation of Solution Domains

The interpolation functions NDSolve returns will extrapolate if they are queried outside of their domain. In a finite element analysis, this may not be wanted since the extrapolation will almost certainly not give the physically correct value. This behavior can be switched off.

This shows the solution domain.
In[95]:=
Click for copyable input
Out[97]=

The following PDE is going to be solved over the solution domain.

In[98]:=
Click for copyable input
This solves the equation.
In[99]:=
Click for copyable input

The resulting interpolation function can be queried at points that do not belong to the solution domain. This will result in a warning message and an extrapolated value.

In[100]:=
Click for copyable input
Out[100]=

To change the behavior, you can have the interpolation function return Indeterminate and quiet the then-unnecessary warning message.

This solves the equation but gives Indeterminate for query points outside the solution domain.
In[101]:=
Click for copyable input

A query point outside the solution domain now returns Indeterminate and no warning message.

In[102]:=
Click for copyable input
Out[102]=

Stabilization of Convection Dominated Equations

Consider a convection-diffusion type of equation

In a situation where the convective component becomes large, the solution of the PDE may become unstable, as shown in the following example.

Specify a convection-diffusion equation operator.
In[103]:=
Click for copyable input
Out[104]=
Specify boundary conditions.
In[105]:=
Click for copyable input
Solve a convection-diffusion equation with a large convective component.
In[106]:=
Click for copyable input
Out[106]=

The message indicates that the solution may not be stable.

Visualize the solution of a convection-diffusion equation with a large convective component.
In[107]:=
Click for copyable input
Out[107]=

The result is disastrous. There are a few things that can be done in such situations. First, however, we investigate how the Peclet number is computed.

Extract the mesh from the interpolating function.
In[108]:=
Click for copyable input
Compute the minimal diameter of the mesh elements.
In[109]:=
Click for copyable input
Out[109]=
Compute the Peclet number from the convection and diffusion coefficients and the minimal element diameter.
In[110]:=
Click for copyable input
Out[110]=

One way to improve the stability of the solution is to add so-called artificial diffusion. This artificial diffusion is then added to the diffusive term in order to stabilize the PDE.

Estimate the amount of artificial diffusion needed.
In[111]:=
Click for copyable input
Out[113]=
Verify that the addition of the artificial diffusion would return a stable result.
In[114]:=
Click for copyable input
Out[114]=

Note that the newly computed Peclet number is smaller than the mesh order.

Solve a convection-diffusion equation with a large convective component.
In[115]:=
Click for copyable input
Visualize the solution of a convection-diffusion equation with a large convective component.
In[116]:=
Click for copyable input
Out[116]=
Visualize a cross section of the solution.
In[117]:=
Click for copyable input
Out[117]=

This is already much better. The disadvantage with this method is that the PDE is modified. This may not always be a good approach. As an alternative, the mesh can be modified. For more information about mesh generation, please consult the "Element Mesh Generation" tutorial.

Generate and visualize a mesh with a refined boundary.
In[118]:=
Click for copyable input
Out[119]=
Solve a convection-diffusion equation with a large convective component on a predefined mesh.
In[120]:=
Click for copyable input
Visualize the solution of a convection-diffusion equation with a large convective component.
In[121]:=
Click for copyable input
Out[121]=
Visualize a cross section of the solution.
In[122]:=
Click for copyable input
Out[122]=

The solution is improved compared to the previous example, at the cost of some additional computational time. In principal, it would be sufficient to only refine the right and top boundary of the mesh. This would then decrease the computational time.