Feel++ Coefficient Form PDE Toolboxes in Python
1. Introduction
The toolbox solves a coefficient form partial differential equation (PDE) using the finite element method. In coefficient form, the PDE is expressed as a linear combination of differential operators with coefficients that may vary spatially. The PDE is discretized using a finite element method, where the solution is approximated by a linear combination of basis functions defined on a mesh. The coefficients of the differential operators are also discretized on the mesh, in time and depend on parameters. The resulting system of linear equations is then solved to obtain the numerical solution.
3. System of PDEs
A lot of PDE(s) can be writen in a generic form, and depends mainly on the definition of coefficients. The generic form that we use is describe by the next equation, find \(u: \Omega \subset \mathbb{R}^d \longrightarrow \mathbb{R}^n\) with \(d=2,3\) and \(n=1\) (\(u\) is a scalar field) or \(n=d\)(\(u\) is a vector field) such that
We call this generic form by Coefficient Form PDE and the coefficients are

\(d\) : damping or mass coefficient

\(c\) : diffusion coefficient

\(\alpha\) : conservative flux convection coefficient

\(\gamma\) : conservative flux source term

\(\beta\) : convection coefficient

\(a\) : absorption or reaction coefficient

\(f\) : source term
All the coefficient may depend on the unknown \(u\) and on the space variable \(x\), time \(t\), parameters \(\mu\) and other unknowns \(u_1, \dots, u_N\).
parameters \(\mu\) may depend on the unknown \(u\) and on the space variable \(x\), time \(t\) and other unknowns \(u_1, \dots, u_N\). 
Many problems are multiphysics (i.e. several unknowns) and the generic form can be extended naturally into a system of generic PDE. The system of \(N\) equations reads then: For \(i=1\dots N\), find \(u_i: \Omega \subset \mathbb{R}^d \longrightarrow \mathbb{R}^{n_i}\) with \(d=2,3\) and \(n_i=1\) (\(u_i\) is a scalar field) or \(n_i=d\)(\(u_i\) is a vector field) such that
\[\begin{eqnarray*}
d_i \frac{\partial u_i}{\partial t} + \nabla \cdot \left( c_i \nabla u_i  \alpha_i u_i + \gamma_i \right) + \beta_i \cdot \nabla u_i + a_i u_i = f_i \quad \text{ in } \Omega
\end{eqnarray*}\]
The coefficients of each equation can be defined by an expression that depends of the current unknown or unknowns of other equations. 
4. Coefficients
We need also to respect some constraint on the coefficient shape as described in the next table.
Coefficient  Shape if Scalar Unknown  Shape if Vectorial Unknown 

\(d\) 
scalar 
scalar 
\(c\) 
scalar or matricx 
scalar or matrix 
\(\alpha\) 
vectorial 
scalar or matrix 
\(\gamma\) 
vectorial 
matrix 
\(\beta\) 
vectorial 
vectorial 
\(a\) 
scalar 
scalar 
\(f\) 
scalar 
vectorial 
6. Boundary Conditions
Here are supported boundary conditions
6.1. Dirichlet
Dirichlet boundary condition is a type of boundary condition commonly used in partial differential equations. It specifies the value of the solution at the boundary of the domain. In other words, it prescribes the behavior of the solution at the boundary.
The condition may depend on the space variable \(x\), time \(t\), parameters \(\mu\) and other unknowns \(u_1, \dots, u_N\) than the current one.
For example, in a heat transfer problem, a Dirichlet boundary condition may specify the temperature at the boundary of the domain. In a fluid flow problem, a Dirichlet boundary condition may specify the velocity or pressure at the boundary.
The Dirichlet boundary condition is essential in determining a unique solution to a PDE problem. Without it, the solution would be underdetermined, and there would be an infinite number of solutions that satisfy the PDE.
The shape of the Dirichlet condition is the same as the unknown shape.
The user provides the expression for \((g_i)_{i=1\dots N}\) in the .json
file only if Dirichlet conditions are used.
6.2. Neumann
Neumann boundary condition is another type of boundary condition commonly used in partial differential equations. It specifies the normal derivative of the solution at the boundary of the domain. In other words, it prescribes the flux of the solution across the boundary.
For example, in a heat transfer problem, a Neumann boundary condition may specify the heat flux at the boundary of the domain. In a fluid flow problem, a Neumann boundary condition may specify the normal stress or shear stress at the boundary.
The Neumann boundary condition is also essential in determining a unique solution to a PDE problem. It provides additional information about the behavior of the solution at the boundary, which complements the Dirichlet boundary condition. Together, the Dirichlet and Neumann boundary conditions form a complete set of boundary conditions that fully specify the PDE problem.
The Neumann conditions may depend on the space variable \(x\), time \(t\), parameters \(\mu\) and the unknowns \(u_1, \dots, u_N\).
The shape of the Neumann condition is the same as the unknown shape.
The user provides the expression for \((g_i)_{i=1\dots N}\) in the .json
file only if Neumann conditions are used.
6.3. Robin
Robin boundary condition is a type of boundary condition that combines both Dirichlet and Neumann boundary conditions. It specifies a linear combination of the solution and its normal derivative at the boundary of the domain. In other words, it prescribes both the value and the flux of the solution at the boundary.
For example, in a heat transfer problem, a Robin boundary condition may specify a heat transfer coefficient that relates the temperature difference between the boundary and the surrounding medium to the heat flux at the boundary. In a fluid flow problem, a Robin boundary condition may specify a slip coefficient that relates the velocity difference between the boundary and the surrounding medium to the shear stress at the boundary.
The Robin boundary condition is useful in modeling situations where the boundary is in contact with a medium that has a different thermal or mechanical behavior than the domain. It provides a more realistic and accurate description of the physical problem than using only Dirichlet or Neumann boundary conditions.
The Robin conditions may depend on the space variable \(x\), time \(t\), parameters \(\mu\) and the unknowns \(u_1, \dots, u_N\) other than the current one.
The shape of the Robin condition is the same as the unknown shape.
The user provides the expression for \((\eta_i)_{i=1\dots N}\) and \((\zeta_i)_{i=1\dots N}\) in the .json
file only if Robin conditions are used.
7. Finite Element Approximation
The following table summarizes the supported finite element approximation for each type of unknown. Here’s an example of an asciidoc table with a list of finite elements supported by the Coefficient form PDE toolbox:
Finite Element  polynomial order  Description 

Pch<k> 
0,1,2 
Piecewise continuous scalar functions of arbitrary degree \(k\) 
Pchv<k> 
0,1,2 
Piecewise continuous vectorial functions of arbitrary degree \(k\) 
RT<k> 
k=0 
RaviartThomas element of degree \(k\) 
This table lists various finite elements supported by Coefficient form PDE, along with a brief description of each element.
Feel++ supports a wider range of finite elements, including piecewise arbitrary order polynomials, as well as mixed finite elements such as RaviartThomas and BrezziDouglasMarini elements or Nedelec’s first families. 
8. Time scheme
8.1. Backward Differences Formula
The backward difference formula scheme is a numerical method for approximating the derivative of a function. It is commonly used in numerical analysis and is defined as:
\(f'(x_n) ≈ \frac{1}{\Delta t} \left(\alpha f(x_n) + \beta f(x_{n1}) + \gamma f(x_{n2}) + \dots\right)\)
where \(\Delta t\) is the time step size, \(x_n\) is the point at which the derivative is approximated, and \(\alpha\), \(\beta\), \(\gamma\), etc. are coefficients that depend on the order of the scheme.
For example, the firstorder backward difference formula scheme is:
\(f'(x_n) ≈ (f(x_n)  f(x_{n1})) / \Delta t\)
whereas the secondorder backward difference formula scheme is:
\(f'(x_n) ≈ (3f(x_n)  4f(x_{n1}) + f(x_{n2})) / (2 \Delta t)\)
and so on. The backward difference formula scheme is useful when the function is only known at discrete points and its derivative needs to be approximated at those points.
8.2. Theta scheme
The theta scheme is a numerical method for solving partial differential equations. It is a finite difference scheme that combines the forward and backward difference formulas to obtain a weighted average of the solutions at the current and previous time steps. The theta scheme is defined as:
where \(u_i^n\) is the numerical solution at the \(i\)th spatial point and \(n\)th time step, \(\Delta t\) and \(\Delta x\) are the time and spatial step sizes, and \(\theta\) is a parameter that determines the weighting between the current and previous time steps.
When \(\theta=0\), the scheme reduces to the backward difference formula, whereas when \(\theta=1\), it reduces to the forward difference formula.
For \(\theta=0.5\), the scheme is known as the CrankNicolson scheme, which is a popular choice due to its stability and accuracy. The theta scheme is widely used in numerical simulations of heat transfer, fluid flow, and other physical phenomena.
9. Stabilized finite element methods
Stabilized finite element methods are a class of numerical methods used to solve partial differential equations (PDEs). These methods are designed to overcome the limitations of traditional finite element methods, which can suffer from numerical instabilities and inaccuracies when applied to certain types of PDEs, such as those with convectiondominated or highly oscillatory solutions.
Stabilized finite element methods introduce additional terms to the weak form of the PDE, which act as stabilizers to improve the accuracy and stability of the numerical solution. These terms are typically chosen to balance the effects of convection, diffusion, and reaction in the PDE, and to ensure that the numerical solution satisfies certain physical and mathematical constraints.
There are several types of stabilized finite element methods, including streamline diffusion, PetrovGalerkin, and leastsquares methods. Each method has its own strengths and weaknesses, and the choice of method depends on the specific problem being solved.
Stabilized finite element methods have been successfully applied to a wide range of PDE problems, including fluid dynamics, heat transfer, and structural mechanics. They have also been used in the development of advanced simulation tools for engineering and scientific applications.
Coefficient form PDE toolbox provides the possibility to use stabilized finite element methods (GaLS and SUPG) for equations such as the advection diffusion reaction equation.
9.1. GaLS
The Galerkin least squares formulation is a stabilized finite element method used to solve partial differential equations. The variational formulation of the Galerkin least squares method is given by:
Find $u \in V$ such that
where \(V\) is a finite element space, \(a\) is the bilinear form defined by, for \((u,v) \in V \times V)\)
and \(l\) is the linear form defined by
Here, \(\epsilon\), \(\beta\), and \(\gamma\) are positive constants that control the balance between diffusion, convection, and reaction in the PDE, and \(\alpha\) is a positive constant that controls the strength of the stabilization term on the boundary. The Galerkin least squares method introduces additional terms to the bilinear form that act as stabilizers to improve the accuracy and stability of the numerical solution. These terms are chosen to minimize the residual of the PDE in a least squares sense, and are typically expressed in terms of the gradient of the solution and its higherorder derivatives. The Galerkin least squares method has been shown to be effective in solving a wide range of PDE problems, including convectiondominated and highly oscillatory problems.
where \(\tau\) is a positive constant that controls the strength of the stabilization term. The first term in the additional terms is a diffusion term that penalizes the gradient of the solution, while the second term is a convection term that penalizes the solution itself. The third term is a reaction term that ensures that the numerical solution satisfies certain physical and mathematical constraints. The Galerkin least squares method with these additional terms has been shown to be effective in improving the accuracy and stability of the numerical solution for a wide range of PDE problems.
9.2. SUPG
The Streamline Upwind Petrov Galerkin (SUPG) method is a stabilized finite element method used to solve partial differential equations. The additional terms in the SUPG method are given by:
where \(\tau\) is a positive constant that controls the strength of the stabilization term, \(\beta\) is a vector field that represents the direction and magnitude of the convection term in the PDE, and $\Delta$ is the Laplace operator. The SUPG method introduces an additional term that penalizes the gradient of the solution in the direction of the convection term, which improves the accuracy and stability of the numerical solution for convectiondominated problems. The SUPG method has been shown to be effective in solving a wide range of PDE problems, including fluid dynamics and heat transfer.
9.3. Coefficient form PDE toolbox
Given a cfpdes
equation named myeq
, SUPG and GaLS can be used as stabilisation methods.
To enable them use, in the commandline or .cfg
file, the option cfpdes.eq.stabilitsation=1
and define the stabilisation type cfpdes.eq.stabilitsation.type=gls #supg#unusualgls #gls
type  


default option 



Examples are available here