According to the tutorial of FEniCS (software that solves PDEs using the Finite Elements method), solving the Poisson equation in 2D is programming the “Hello World” of PDEs.
The mathematical aspects
We have to consider that the Poisson equation involves these two conditions:
-∇² u (x) = f(x), x in Ω (01)
u(x) = uD (x), x on δΩ (02)
In (01), f(x) is a known function, Ω is the domain of the solution, ∇² is the Laplace operator (Δ), and we have to solve u(x). In 2D, we can represent -∇² u (x) as: -(δ²u/δx²) – (δ²u/δy²) = f(x,y).
In (02), ∂Ω is the boundary of Ω, and u=uD represents the boundary condition.
FEniCS uses the Finite Elements method to solve PDEs, then it is compulsory to express the PDE (01) in a variational form. In order to get the variational form, we can multiply the PDE by a function called v(test function) with a premise that v=0 on δΩ, then we can integrate it over the domain Ω, having this:
(03)
Analysing the right side of (03), we can do integration by parts to finally get a result such:
(04)
Both v(TestlFunction) and u(TrialFunction) are going to be defined in the same V domain that contains a mesh with the number of vertex selected to discrete the domain. Since the functions located in the domain V are continuous, a technique used by FEniCS is to discrete to make it computational. Then, FEniCS only considers a set of functions called V_h which are also included in V. Later, it is also convenient to assign the left side of (04) to the variable L, and the right side of the equation to the variable a, as follows:
(05)
The author in the tutorial indicates that he has constructed a problem with a known analytical solution so that it can be easily checked the correctness in computational terms. The following values have been selected in this matter:
f(x,y)=−6,
uD(x,y)=1+x²+2y²,
Ω=[0,1]×[0,1].
The computational aspects
Let’s do a review line by line throughout this “hello world” FEniCS implementation.
from fenics import *
Firstly, we are importing all the libraries of FEnicCS to be used in the program. To install FEniCS, I used Python3, Dolfin and podman on Fedora 32. Other tools can be found here.
mesh = UnitSquareMesh(8, 8)
We are setting a mesh of a shape of square of dimension 1(Unit square mesh). It will be divided by 8 parts in 2D, 8 on x axis, and 8 on y.
V = FunctionSpace(mesh, 'P', 1)
The previous mesh defined, is used in a finite element function space called V. The Lagrange method of order 1 is applied to solve the equation with the Poisson restrictions explained before. A note aside, is the use of order greater than 1, can generate more errors.
u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)
The u_D function is a known input which defines the boundary conditions. In this case, the formula implies a polynomial order 2. x[0] denotes the value of x, and x[1] denotes y.
def boundary(x, on_boundary):
return on_boundary
The on_boundary function has previously been coded on FEniCS to validate if a point belongs to the current boundary of the domain, or not.
bc = DirichletBC(V, u_D, boundary)
The boundary condition (bc) is stated in the space V, to be applied on boundary of the function u_D. By using the type of boundary condition “Dirichlet”, it will specify the values that a solution needs to take along the boundary in the domain V.
u = TrialFunction(V)
It was mentioned that the u function is the PDE we are looking for. In order to find the solution, it was also explained that we need the variational formula. In this case, “u” is the new problem to solve and it will be called as TrialFunction in a TrialSpace V.
v = TestFunction(V)
It was taking into consideration a v function that will have the value of zero in the border. This v function is going to be called TestFunction in the V TestSpace. It will be computed against the trial function u to satisfy the variational problem.
f = Constant(-6.0)
The f function, is a function that is known to solve the EDP. In this case, is a constant = -6.
a = dot(grad(u), grad(v))*dx
The left side of the variational formula is going to be assigned to a. It will compute the trial function u versus all the v test functions in V.
L = f*v*dx
The right side of the variational formula considers the function f on all the test functions v.
u = Function(V)
Since u is the Function that solves the EDP, we need to prepare this variable to receive the final solution in the domain V.
solve(a == L, u, bc)
The operation to solve u requires that, a is equal to L, as well as the calculation of bc. The solution will be assigned to u.
plot(u)
The solution u can be plot and could be visualized using tools such as Paraview or Visit.
plot(mesh)
We can also display the mesh to the see the highness of the PDE in regards its domain.
vtkfile = File('poisson.pvd')
We are going to save our solution in the vtkfile variable, under the name of poisson.pvd.
vtkfile << u
Again, u is the final answer (PDE) that will be stored in the vtkfile.
error_L2 = errornorm(u_D, u, 'L2')
After finding u, we need to know how far we are in comparison to the correct solution. It is known that the U_D is the solution, then we have to calculate the error with the error_L2 formula which utilizes spaces of Gilbert to find the norm in u.
vertex_values_u_D = u_D.compute_vertex_values(mesh)
The error is also computed in the vertexes of the mesh, in this case, the ones in the border.
vertex_values_u = u.compute_vertex_values(mesh)
This line compute the error in the vertexes of the mesh inside the border of the PDE.
import numpy as np
We are going to import the library numpy with the alias np to use max and abs functions:
error_max = np.max(np.abs(vertex_values_u_D - vertex_values_u))
This calculates the maximum error found in the solution including border and content.
print('error_L2 =', error_L2)
This line only prints the error found in the calculation of error_L2.
print('error_max =', error_max)
This line prints the maximum error found in the PDE.
The simulation
Paraview shows both, the domain and the PDE as follows:
You can find more PDEs tested from the tutorial in my github. Thanks to Manuel Merino!
Sólo para Mí! 