import dolfinx
print(f'DOLFINx version: {dolfinx.__version__}')
DOLFINx version: 0.9.0
\(~\)
Dedalus solves differential equations using spectral methods. It’s open-source, written in Python, and MPI-parallelized
\(~\)
FEniCS is a popular open-source computing platform for solving partial differential equations (PDEs) with the finite element method (FEM). FEniCS enables users to quickly translate scientific models into efficient finite element code. With the high-level Python and C++ interfaces to FEniCS, it is easy to get started, but FEniCS offers also powerful capabilities for more experienced programmers. FEniCS runs on a multitude of platforms ranging from laptops to high-performance computers
The latest stable release of FEniCSx is version 0.9, which was released in October 2024. The easiest way to start using FEniCSx
on MacOS and other systems is to install it using conda
:
$ conda create -n fenicsx
$ conda activate fenicsx
$ conda install -c conda-forge fenics-dolfinx mpich pyvista petsc petsc4py
$ conda install ipykernel
$ python -m ipykernel install --user --name fenicsx --display-name "FEniCSx"
\(~\)
import dolfinx
print(f'DOLFINx version: {dolfinx.__version__}')
DOLFINx version: 0.9.0
\(~\)
The FEniCS
Project is a research and software initiative focused on developing mathematical methods and software for solving partial differential equations (PDEs). Its goal is to provide intuitive, efficient, and flexible tools for scientific computing. Launched in 2003, the project is the result of collaboration among researchers from universities and research institutes worldwide. For the latest updates and more information, visit the FEniCS Project
The latest version of the FEniCS
project, FEniCSx
, consists of several building blocks, namely DOLFINx
, UFL
, FFCx
, and Basix
. We will now go through the main objectives of each of these building blocks
DOLFINx
is the high performance C++ backend of FEniCSx
, where structures such as meshes, function spaces and functions are implemented. Additionally, DOLFINx
also contains compute intensive functions such as finite element assembly and mesh refinement algorithms. It also provides an interface to linear algebra solvers and data-structures, such as PETSc
UFL
is a high-level form language for describing variational formulations with a high-level mathematical syntaxFFCx
is the form compiler of FEniCSx
; given variational formulations written with UFL
, it generates efficient C codeBasix
is the finite element backend of FEniCSx
, responsible for generating finite element basis functionsThe goal of this tutorial is to demonstrate how to apply the finite element to solve PDEs using FEniCS
. Through a series of examples, we will demonstrate how to:
Solve linear PDEs (such as the Poisson equation)
Solve time-dependent PDEs (such as the heat equation)
Solve non-linear PDEs
Solve systems of time-dependent non-linear PDEs
Important topics include: how to set boundary conditions of various types (Dirichlet, Neumann, Robin), how to create meshes, how to define variable coefficients, how to interact with linear and non-linear solvers, and how to post-process and visualize solutions
Authors: Hans Petter Langtangen, Anders Logg, Jørgen S. Dokken
The goal of this section is to solve one of the most basic PDEs, the Poisson equation, with a few lines of code in FEniCSx
. We start by introducing some fundamental FEniCSx
objects, such as Function
, functionspace
, TrialFunction
and TestFunction
, and learn how to write a basic PDE solver. This will include:
How to formulate a mathematical variational problem
How to apply boundary conditions
How to solve the discrete linear system
How to visualize the solution
The Poisson equation is the following boundary-value problem
\[\begin{aligned} -\nabla^2 u(\mathbf{x}) &= f(\mathbf{x})&&\mathbf{x} \in \Omega\\ u(\mathbf{x}) &= u_D(\mathbf{x})&& \mathbf{x} \in \partial\Omega \end{aligned}\]Here, \(u=u(\mathbf{x})\) is the unknown function, \(f=f(\mathbf{x})\) a prescribed function, the Laplace operator (often written as \(\Delta\)), \(\Omega\) the spatial domain, and \(\partial\Omega\) the boundary of \(\Omega\). The Poisson problem, including both the PDE, \(-\nabla^2 u = f\), and the boundary condition, \(u=u_D\) on \(\partial\Omega\), is an example of a boundary-value problem, which must be precisely stated before we can start solving it numerically with FEniCSx
\[-\frac{\partial^2 u}{\partial x^2} - \frac{\partial^2 u}{\partial y^2} = f(x,y)\]
The unknown \(u\) is now a function of two variables, \(u=u(x,y)\), defined over the two-dimensional domain \(\Omega\)
Solving a boundary value problem in FEniCSx
consists of the following steps:
FEniCSx
As we have already covered step 1, we shall now cover steps 2-4
FEniCSx
is based on the finite element method, which is a general and efficient mathematical technique for the numerical solution of PDEs. The starting point for finite element methods is a PDE expressed in variational form
The basic recipe for turning a PDE into a variational problem is:
The function \(v\) which multiplies the PDE is called a test function. The unknown function \(u\) that is to be approximated is referred to as a trial function. The terms trial and test functions are used in FEniCSx
too. The test and trial functions belong to certain function spaces that specify the properties of the functions
In the present case, we multiply the Poisson equation by a test function \(v\) and integrate over \(\Omega\):
\[\int_\Omega (-\nabla^2 u) v~\mathrm{d} x = \int_\Omega f v~\mathrm{d} x\]
Here \(\mathrm{d} x\) denotes the differential element for integration over the domain \(\Omega\). We will later let \(\mathrm{d} s\) denote the differential element for integration over \(\partial\Omega\), the boundary of \(\Omega\)
A rule of thumb when deriving variational formulations is that one tries to keep the order of derivatives of \(u\) and \(v\) as small as possible. Here, we have a second-order differential of \(u\), which can be transformed to a first derivative by employing the technique of integration by parts. The formula reads
\[-\int_\Omega (\nabla^2 u)v~\mathrm{d}x = \int_\Omega\nabla u\cdot\nabla v~\mathrm{d}x- \int_{\partial\Omega}\frac{\partial u}{\partial n}v~\mathrm{d}s\]
where \(\frac{\partial u}{\partial n}=\nabla u \cdot \vec{n}\) is the derivative of \(u\) in the outward normal direction \(\vec{n}\) on the boundary
Another feature of variational formulations is that the test function \(v\) is required to vanish on the parts of the boundary where the solution \(u\) is known
In the present problem, this means that \(v\) is \(0\) on the whole boundary \(\partial\Omega\). Thus, the second term in the integration by parts formula vanishes, and we have that
\[\int_\Omega \nabla u \cdot \nabla v~\mathrm{d} x = \int_\Omega f v~\mathrm{d} x\]
\(~\)
If we require that this equation holds for all test functions \(v\) in some suitable space \(\hat{V}\), the so-called test space, we obtain a well-defined mathematical problem that uniquely determines the solution \(u\) which lies in some function space \(V\). Note that \(V\) does not have to be the same space as \(\hat{V}\). We call the space \(V\) the trial space. We refer to the equation above as the weak form/variational form of the original boundary-value problem. We now properly state our variational problem: Find \(u\in V\) such that
\[\int_\Omega \nabla u \cdot \nabla v~\mathrm{d} x = \int_\Omega f v~\mathrm{d} x\qquad \forall v \in \hat{V}\]
For the present problem, the trial and test spaces \(V\) and \(\hat{V}\) are defined as
\[\begin{aligned} V&=\{v\in H^1(\Omega) \,\vert\, v=u_D&\text{on } \partial \Omega \}\\ \hat{V}&=\{v\in H^1(\Omega) \,\vert\, v=0 &\text{on } \partial \Omega \} \end{aligned}\]In short, \(H^1(\Omega)\) is the Sobolev space containing functions \(v\) such that \(v^2\) and \(\vert \nabla v \vert ^2\) have finite integrals over \(\Omega\). The solution of the underlying PDE must lie in a function space where the derivatives are also continuous, but the Sobolev space \(H^1(\Omega)\) allows functions with discontinuous derivatives. This weaker continuity requirement in our weak formulation (caused by the integration by parts) is of great importance when it comes to constructing the finite element function space. In particular, it allows the use of piecewise polynomial function spaces. This means that the function spaces are constructed by stitching together polynomial functions on simple domains such as intervals, triangles, quadrilaterals, tetrahedra and hexahedra
The variational problem is a continuous problem: it defines the solution \(u\) in the infinite-dimensional function space \(V\). The finite element method for the Poisson equation finds an approximate solution of the variational problem by replacing the infinite-dimensional function spaces \(V\) and \(\hat{V}\) by discrete (finite dimensional) trial and test spaces \(V_h\subset V\) and \(\hat{V}_h \subset \hat{V}\). The discrete variational problem reads: Find \(u_h\in V_h\) such that
\[ \begin{aligned} \int_\Omega \nabla u_h \cdot \nabla v~\mathrm{d} x &= \int_\Omega fv~\mathrm{d} x && \forall v \in \hat{V}_h. \end{aligned} \]
This variational problem, together with suitable definitions of \(V_h\) and \(\hat{V}_h\), uniquely define our approximate numerical solution of the Poisson equation. Note that the boundary condition is encoded as part of the test and trial spaces. This might seem complicated at first glance, but means that the finite element variational problem and the continuous variational problem look the same
We will introduce the following notation for variational problems: Find \(u\in V\) such that
\[\begin{aligned} a(u,v)&=L(v)&& \forall v \in \hat{V} \end{aligned}\]
For the Poisson equation, we have:
\[\begin{aligned} a(u,v) &= \int_{\Omega} \nabla u \cdot \nabla v~\mathrm{d} x,\\ L(v) &= \int_{\Omega} fv~\mathrm{d} x \end{aligned}\]
In the literature \(a(u,v)\) is known as the bilinear form and \(L(v)\) as a linear form. For every linear problem, we will identify all terms with the unknown \(u\) and collect them in \(a(u,v)\), and collect all terms with only known functions in \(L(v)\).
To solve a linear PDE in FEniCSx
, such as the Poisson equation, a user thus needs to perform two steps:
In this section, you will learn:
DOLFINx
Up to this point, we’ve looked at the Poisson problem in very general terms: the domain \(\Omega\), the boundary condition \(u_D\), and the right-hand side \(f\) were all left unspecified. To actually solve something, we now need to pick concrete choices for \(\Omega\), \(u_D\), and \(f\)
A good strategy is to set up the problem in a way that we already know the exact solution. That way, we can easily check whether our numerical solution is correct. Polynomials of low degree are usually the best choice here, because continuous Galerkin finite element spaces of degree \(r\) can reproduce any polynomial of degree \(r\) exactly
To test our solver, we’ll construct a problem where we already know the exact solution. This approach is known as the method of manufactured solutions. The idea is simple:
Step 1: Choosing the exact solution
Let’s take a quadratic function in 2D:
\[ u_e(x,y) = 1 + x^2 + 2y^2 \]
Step 2: Computing the source term
If we insert \(u_e\) into the Poisson equation, we obtain
\[f(x,y) = -6, \;\; u_D(x,y) = u_e(x,y) = 1 + x^2 + 2y^2\]
Notice that this holds regardless of the domain shape, as long as we prescribe \(u_e\) on the boundary
Step 3: Choosing the domain
For simplicity, let’s use the unit square:
\[\Omega = [0,1] \times [0,1]\]
Step 4: Summary
This small example illustrates a very powerful strategy:
This workflow is at the heart of the method of manufactured solutions, and it provides a simple but rigorous way to validate our solver
Generating simple meshes
The next step is to define the discrete domain, called the mesh. We do this using one of FEniCSx
’s built-in mesh generators
Here, we create a unit square mesh spanning \([0,1]\times[0,1]\). The cells of the mesh can be either triangles or quadrilaterals:
import numpy as np
from mpi4py import MPI
from dolfinx import mesh
= mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral) domain
Notice that we need to provide an MPI communicator. This determines how the program behaves in parallel:
MPI.COMM_WORLD
, a single mesh is created and distributed across the number of processors we specify. For example, to run the program on two processors, we can use:$ mpirun -n 2 python tutorial_01.py
MPI.COMM_SELF
, each processor will create its own independent copy of the mesh. This can be useful when running many small problems in parallel, for example when sweeping over different parametersDefining the finite element function space
Once the mesh is created, the next step is to define the finite element function space \(V\). DOLFINx
supports a wide variety of finite element spaces of arbitrary order. For a full overview, see the list of Supported elements in DOLFINx
When creating a function space, we need to specify: 1. The mesh on which the space is defined 2. The element family (e.g., Lagrange, Raviart–Thomas, etc.) 3. The polynomial degree of the element
In DOLFINx
, this can be done by passing a tuple of the form ("family", degree)
, as shown below:
from dolfinx.fem import functionspace
= functionspace(domain, ("Lagrange", 1)) V
The next step is to create a function that stores the Dirichlet boundary condition. We then use interpolation to fill this function with the prescribed values
from dolfinx import fem
= fem.Function(V)
uD lambda x: 1 +x[0]**2 +2 *x[1]**2) uD.interpolate(
Now that we have defined the boundary data (which, in this case, is also the exact solution of our finite element problem), we need to enforce it on the boundary of the mesh.
The first step is to identify which parts of the mesh correspond to the outer boundary. In DOLFINx
, the boundary is represented by facets (that is, the line segments making up the outer edges in 2D or surfaces in 3D).
We can extract the indices of all exterior facets using:
= domain.topology.dim
tdim = tdim -1
fdim
domain.topology.create_connectivity(fdim, tdim)= mesh.exterior_facet_indices(domain.topology) boundary_facets
This gives us the set of facets lying on the boundary of our discrete domain. Once we know where the boundary is, we can proceed to apply the Dirichlet boundary conditions to all degrees of freedom (DoFs) located on these facets
For our current problem, we are using the “Lagrange” degree-1 function space. In this case, the degrees of freedom (DoFs) are located at the vertices of each cell, so every boundary facet contains exactly two DoFs
To identify the local indices of these boundary DoFs, we use dolfinx.fem.locate_dofs_topological
. This function takes three arguments:
Once we have the boundary DoFs, we can create the Dirichlet boundary condition as follows:
= fem.locate_dofs_topological(V, fdim, boundary_facets)
boundary_dofs = fem.dirichletbc(uD, boundary_dofs) bc
Defining the trial and test function
In mathematics, we usually distinguish between the trial space \(V\) and the test space \(\hat{V}\). For the present problem, the only difference between the two would be the treatment of boundary conditions
In FEniCSx
, however, boundary conditions are not built directly into the function space. This means we can simply use the same space for both the trial and test functions
To express the variational formulation, we make use of the Unified Form Language (UFL)
import ufl
= ufl.TrialFunction(V)
u = ufl.TestFunction(V) v
Defining the source term
Since the source term is constant throughout the domain, we can represent it using dolfinx.fem.Constant
:
from dolfinx import default_scalar_type
= fem.Constant(domain, default_scalar_type(-6)) f
While we could simply define the source term as f = -6
, this has a limitation: if we later want to change its value, we would need to redefine the entire variational problem. By using dolfinx.fem.Constant
, we can easily update the value during the simulation, for example with f.value = 5
Another advantage is performance: declaring f
as a constant allows the compiler to optimize the variational formulation, leading to faster assembly of the resulting linear system
Defining the variational problem
Now that we have defined all the components of our variational problem, we can write down the weak formulation:
= ufl.dot(ufl.grad(u), ufl.grad(v)) *ufl.dx
a = f *v *ufl.dx L
Notice how closely the Python syntax mirrors the mathematical expressions:
\[a(u,v) = \int_{\Omega} \nabla u \cdot \nabla v \,\mathrm{d}x, \;\; L(v) = \int_{\Omega} f v \,\mathrm{d}x\]
Here, ufl.dx
represents integration over the domain \(\Omega\), i.e. over all cells of the mesh. This illustrates one of the major strengths of FEniCSx
: variational formulations can be written in Python in a way that almost directly matches their mathematical form, making it both natural and convenient to specify and solve complex PDE problems
Expressing inner products
The inner product
\[\int_\Omega \nabla u \cdot \nabla v \,\mathrm{d}x\]
can be expressed in different ways in UFL
. In our example, we wrote it as: ufl.dot(ufl.grad(u), ufl.grad(v)) *ufl.dx
. In UFL, the dot operator performs a contraction: it sums over the last index of the first argument and the first index of the second argument. Since both \(\nabla u\) and \(\nabla v\) are rank-1 tensors (vectors), this reduces to a simple dot product.
For higher-rank tensors, such as matrices (rank-2 tensors), the appropriate operator is ufl.inner, which computes the full Frobenius inner product. For vectors, however, ufl.dot
and ufl.inner
are equivalent
Forming and solving the linear system
Having defined the finite element variational problem and boundary conditions, we can now create a dolfinx.fem.petsc.LinearProblem
. This class provides a convenient interface for solving
the variational problem:
Find \(u_h\in V\) such that \(a(u_h, v)=L(v), \;\; \forall v \in \hat{V}\)
In this example, We will use PETSc
as the linear algebra backend, together with a direct solver (LU factorization)
For more details on Krylov solvers and preconditioners, see the PETSc-documentation. Note that PETSc
is not a required dependency of DOLFINx
, so we explicitly import the DOLFINx
wrapper to interface with PETSc
. Finally, to ensure that the solver options passed to the LinearProblem
apply only to this specific KSP solver, we assign a unique option prefix
from dolfinx.fem.petsc import LinearProblem
= LinearProblem(
problem =[bc],
a, L, bcs={
petsc_options# Direct solver using LU factorization
"ksp_type": "preonly",
"pc_type": "lu"
}
)
# Solve the system
= problem.solve()
uh
# Optionally, view solver information
#ksp = problem.solver
#ksp.view()
Using problem.solve()
, we solve the linear system of equations and return a dolfinx.fem.Function
containing the solution
Computing the error
Finally, we want to compute the error to check the accuracy of the solution. We do this by comparing the finite element solution uh
with the exact solution. We do this by interpolating the exact solution into the the \(P_2\)-function space
= fem.functionspace(domain, ("Lagrange", 2))
V2 = fem.Function(V2)
uex lambda x: 1 +x[0]**2 +2 *x[1]**2) uex.interpolate(
We compute the error in two different ways. First, we compute the \(L^2\) norm of the error, defined by
\[E=\sqrt{\int_\Omega (u_D-u_h)^2 \,\mathrm{d} x}\]
We use UFL
to express the \(L^2\) error, and use dolfinx.fem.assemble_scalar
to compute the scalar value. In DOLFINx
, assemble_scalar
only assembles over the cells on the local process. This means that if we use 2 processes to solve our problem, we need to gather the solution to one. We can do this with the MPI.allreduce
function
= fem.form(ufl.inner(uh -uex, uh -uex) *ufl.dx)
L2_error = fem.assemble_scalar(L2_error)
error_local = np.sqrt(domain.comm.allreduce(error_local, op=MPI.SUM)) error_L2
Secondly, we compute the maximum error at any degree of freedom (dof). The finite element solution uh
can be expressed as a linear combination of the basis functions \(\phi_j\) spanning the space \(V\):
\[u = \sum_{j=1}^N U_j \phi_j\]
When we call problem.solve()
, we obtain all coefficients \(U_1\), \(\dots\), \(U_N\). These coefficients are the degrees of freedom (dofs). We can access the dofs by retrieving the underlying vector from uh
However, note that a second-order function space contains more dofs than a first-order space, so the corresponding arrays cannot be compared directly. Fortunately, since we have already interpolated the exact solution into the first-order space when defining the boundary condition, we can compare the maximum values at the dofs of the approximation space
= np.max(np.abs(uD.x.array -uh.x.array))
error_max
# Only print the error on one process
if domain.comm.rank == 0:
print(f"Error_L2 : {error_L2:.2e}")
print(f"Error_max : {error_max:.2e}")
Error_L2 : 8.24e-03
Error_max : 5.33e-15
Plotting the mesh using pyvista
We will visualize the mesh using pyvista, a Python interface to the VTK toolkit. To begin, We convert the mesh into a format compatible with pyvista
using the function dolfinx.plot.vtk_mesh
. The first step is to create an unstructured grid that pyvista
can work with
You can check the current plotting backend with:
import pyvista
from dolfinx import plot
domain.topology.create_connectivity(tdim, tdim)= plot.vtk_mesh(domain, tdim)
topology, cell_types, geometry = pyvista.UnstructuredGrid(topology, cell_types, geometry) grid
PyVista
supports several backends, each with its own advantages and limitations. For more information and installation instructions, see the pyvista documentation
We can now use the pyvista.Plotter
to visualize the mesh. We will show it both as a 2D and as a 3D warped representation.
In the jupyter notebook, we use the default setting pyvista.OFF_SCREEN=False
, which will renders the plots directly within the notebook
from pathlib import Path
= Path("fenicsx/fundamentals")
results_folder =True, parents=True)
results_folder.mkdir(exist_ok
= pyvista.Plotter(off_screen=True)
plotter =True)
plotter.add_mesh(grid, show_edges
plotter.add_axes()
plotter.view_xy()
# if not pyvista.OFF_SCREEN:
# plotter.show()
# HTML 저장
"fenicsx/fundamentals/unit_square_mesh.html") plotter.export_html(
Plotting a function using pyvista
We want to plot the solution uh
. Since the function space used to define uh
is disconnected from the one used to define the mesh, we first create a mesh based on the DOF coordinates of the function space V
. We then use dolfinx.plot.vtk_mesh
, passing the function space as input, to generate a mesh whose geometry is based on these DOF coordinates
= plot.vtk_mesh(V) u_topology, u_cell_types, u_geometry
Next, we create the pyvista.UnstructuredGrid
and add the DOF-values to the mesh
= pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry)
u_grid "u"] = uh.x.array.real
u_grid.point_data["u")
u_grid.set_active_scalars(
= pyvista.Plotter(off_screen=True)
u_plotter
u_plotter.add_mesh(
u_grid, =True,
show_edges={
scalar_bar_args"title": "u",
"fmt": "%.1f",
"color": "black",
"label_font_size": 12,
# "vertical": True,
"n_labels": 10,
},
)
u_plotter.add_axes()
u_plotter.view_xy()
# if not pyvista.OFF_SCREEN:
# u_plotter.show()
# HTML 저장
"fenicsx/fundamentals/poisson_solution_2D.html") u_plotter.export_html(
We can also warp the mesh by scalar to make use of the 3D plotting
= u_grid.warp_by_scalar()
warped
= pyvista.Plotter(off_screen=True)
u_plotter2
u_plotter2.add_mesh(
warped, =True,
show_edges={
scalar_bar_args"title": "u",
"fmt": "%.1f",
"color": "black",
"label_font_size": 12,
#"vertical": True,
"n_labels": 10,
},=True)
show_scalar_bar
u_plotter2.add_axes()
#if not pyvista.OFF_SCREEN:
# u_plotter2.show()
# HTML 저장
"fenicsx/fundamentals/poisson_solution_3D.html") u_plotter2.export_html(
External post-processing
For post-processing outside of Python, it is recommended to save the solution to a file using either dolfinx.io.VTXWriter
or dolfinx.io.XDMFFile
, and then visualize it in ParaView. This approach is particularly useful for 3D visualization
from dolfinx import io
= results_folder / "poisson"
filename
with io.VTXWriter(domain.comm, filename.with_suffix(".bp"), [uh]) as vtx:
0.0)
vtx.write(
with io.XDMFFile(domain.comm, filename.with_suffix(".xdmf"), "w") as xdmf:
xdmf.write_mesh(domain) xdmf.write_function(uh)
Authors: Hans Petter Langtangen, Anders Logg, Jørgen S. Dokken
In the first FEniCSx
example, we solved a simple problem that allowed us to easily verify the implementation. In this section, we shift our focus to a more physically relevant problem, one that produces solutions with a more interesting structure
Specifically, we aim to compute the deflection \(D(x,y)\) of a two-dimensional circular membrane of radius \(R\), subjected to a load distribution \(p(x,y)\). The governing PDE model is
\[ \begin{aligned} -T \nabla^2D&=p \quad\text{in }\; \Omega=\{(x,y)\,\vert\, x^2+y^2\leq R^2 \} \end{aligned} \]
Here, \(T\) is the tension in the membrane (constant), and \(p\) is the external pressure load. The boundary of the membrane has no deflection. This implies that \(D=0\) is the boundary condition. We model a localized load as a Gaussian function:
\[ \begin{aligned} p(x,y)&=\frac{A}{2\pi\sigma}e^{-\frac{1}{2}\left[\left(\frac{x-x_0}{\sigma}\right)^2 +\left(\frac{y-y_0}{\sigma}\right)^2\right]} \end{aligned} \]
The parameter \(A\) is the amplitude of the pressure, \((x_0, y_0)\) the location of the maximum point of the load, and \(\sigma\) the “width” of \(p\). We will take the center \((x_0,y_0)\) to be \((0,R_0)\) for some \(0<R_0<R\) Then we have
\[ \begin{aligned} p(x,y)&=\frac{A}{2\pi\sigma}e^{-\frac{1}{2}\left[\left(\frac{x}{\sigma}\right)^2 +\left(\frac{y-R_0}{\sigma}\right)^2\right]} \end{aligned} \]
There are many physical parameters in this problem, and we can benefit from grouping them by means of scaling. Let us introduce dimensionless coordinates \(\bar{x}=\frac{x}{R}\), \(\bar{y}=\frac{y}{R}\), and a dimensionless deflection \(w=\frac{D}{D_e}\), where \(D_e\) is a characteristic size of the deflection. Introducing \(\bar{R}_0=\frac{R_0}{R}\), we obtain
\[ \begin{aligned} -\frac{\partial^2 w}{\partial \bar{x}^2} -\frac{\partial^2 w}{\partial \bar{y}^2} &=\frac{R^2A}{2\pi\sigma TD_e}e^{-\frac{R^2}{2\sigma^2}\left[\bar{x}^2+(\bar{y}-\bar{R}_0)^2\right]}\\ &=\alpha e^{-\beta^2 \left[\bar{x}^2+(\bar{y}-\bar{R}_0)^2\right]} \end{aligned} \]
for \(\bar{x}^2+\bar{y}^2<1\), where \(\alpha = \frac{R^2A}{2\pi\sigma TD_e}\) and \(\beta=\frac{R}{\sqrt{2}\sigma}\)
With an appropriate scaling, \(w\) and its derivatives are of order unity. Thus, the left-hand side of the scaled PDE is also of order unity, while the right hand side is characterized by the parameter \(\alpha\). This suggests selecting \(\alpha\) to be on the order of unity; in this particular case, we choose \(\alpha=4\). (It is also possible to derive the analytical solution in scaled coordinates and verify that the maximum deflection is \(D_e\) when \(\alpha = 4\), thereby determining \(D_e\))
With \(D_e=\frac{R^2A}{8\pi\sigma T}\) and omitting the bars, the scaled problem becomes
\[ \begin{aligned} -\nabla^2 w = 4e^{-\beta^2 \left[x^2+(y-R_0)^2\right]} \end{aligned} \]
to be solved over the unit disc, with \(w=0\) on the boundary
In this formulation, only two parameters remain: the dimensionless width of the pressure distribution, \(\beta\), and the location of the pressure peak, \(R_0 \in [0,1]\). As \(\beta \to 0\), the solution approaches the special case \(w = 1 - x^2 - y^2\)
Finally, given a computed scaled solution \(w\), the physical deflection is recovered as
\[ \begin{aligned} D=\frac{AR^2}{8\pi\sigma T}w \end{aligned} \]
Author: Jørgen S. Dokken
In this section, we will solve the deflection of the membrane problem. After finishing this section, you should be able to:
GMSH
Python API and load it into DOLFINx
ufl.SpatialCoordinate
to create a spatially varying functionufl.Expression
into an appropriate function spacedolfinx.Function
at any pointCreating the mesh
To create the computational geometry, we use the Python-API of GMSH. We start by importing the gmsh-module and initializing it
# $ conda install -c conda-forge python-gmsh
import gmsh
if not gmsh.isInitialized():
gmsh.initialize()
The next step is to set up the membrane and start the computations using the Gmsh CAD kernel, which will generate the necessary data structures behind the scenes. When calling addDisk
, the first three arguments give the \(x\), \(y\), and \(z\) coordinates of the circle’s center, while the last two specify the radii in the \(x\)- and \(y\)-directions
= 1.0 # radius of the membrane
R = gmsh.model.occ.addDisk(0.0, 0.0, 0.0, R, R)
membrane
# Synchronize the CAD kernel with the model
gmsh.model.occ.synchronize()
Next, we turn the membrane into a physical surface, so that Gmsh
recognizes it when generating the mesh. Since a surface is a two-dimensional entity, we pass 2
as the first argument, the membrane’s entity tag as the second argument, and the physical tag as the last argument. In a later demo, we will explain exactly when and why this physical tag is important
= 2
gdim = 1
physical_tag
# Remove any existing physical groups with the same (dim, tag)
for dim, tag in gmsh.model.getPhysicalGroups():
if dim == gdim and tag == physical_tag:
gmsh.model.removePhysicalGroups([(dim, tag)])
# Now safely add the new physical group
gmsh.model.addPhysicalGroup(gdim, [membrane], physical_tag)
1
Finally, we generate the two-dimensional mesh. We set a uniform mesh size by modifying the GMSH
options
"Mesh.CharacteristicLengthMin", 0.05)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.05)
gmsh.option.setNumber( gmsh.model.mesh.generate(gdim)
Info : Meshing 1D...
Info : Meshing curve 1 (Ellipse)
Info : Done meshing 1D (Wall 0.000585125s, CPU 0.000627s)
Info : Meshing 2D...
Info : Meshing surface 1 (Plane, Frontal-Delaunay)
Info : Done meshing 2D (Wall 0.0578568s, CPU 0.085236s)
Info : 1552 nodes 3103 elements
Interfacing with GMSH in DOLFINx
We will import the Gmsh-generated mesh directly into DOLFINx
using the dolfinx.io.gmshio
interface. In this example, we did not specify which process created the gmsh
model, so a copy of the model has been created on each MPI process. However, we want to work with a single mesh distributed across all processes. To achieve this, we take the model generated on rank 0 of MPI.COMM_WORLD
, and distribute it to all available ranks
This import will also provide two mesh tags: one for cells marked by physical groups and one for facets marked by physical groups. Since we did not add any physical groups of dimension gdim -1
, facet_markers
object will contain no entities
from mpi4py import MPI
from dolfinx.io import gmshio
from dolfinx.fem.petsc import LinearProblem
= 0
gmsh_model_rank = MPI.COMM_WORLD
mesh_comm = gmshio.model_to_mesh(
domain, cell_markers, facet_markers
gmsh.model,
mesh_comm,
gmsh_model_rank, =gdim
gdim )
We define the function space as in the previous tutorial
from dolfinx import fem
= fem.functionspace(domain, ("Lagrange", 1)) V
import pyvista
from dolfinx.plot import vtk_mesh
# Extract topology from mesh and create pyvista mesh
= vtk_mesh(V)
topology, cell_types, x = pyvista.UnstructuredGrid(topology, cell_types, x)
grid
= pyvista.Plotter(off_screen=True)
plotter =True)
plotter.add_mesh(grid, show_edges
plotter.add_axes()
plotter.view_xy()
# if not pyvista.OFF_SCREEN:
# plotter.show()
# HTML 저장
"fenicsx/fundamentals/membrane_mesh.html") plotter.export_html(
Defining a spatially varying load
The right-hand side pressure function is defined using ufl.SpatialCoordinate
together with two constants: one for \(\beta\) and one for \(R_0\)
from dolfinx import default_scalar_type
import ufl
= ufl.SpatialCoordinate(domain)
x
= fem.Constant(domain, default_scalar_type(12))
beta = fem.Constant(domain, default_scalar_type(0.3))
R0
= 4 *ufl.exp(-beta**2 *(x[0]**2 +(x[1] -R0)**2)) p
Interpolation of a ufl
expression
Since we defined the load p
as a spatially varying function, we want to interpolate it into a suitable function space for visualization. For this, we use dolfinx.Expression
, which can take any ufl
expression along with a set of points on the reference element. In practice, we use the interpolation points of the target space. Because p
varies rapidly in space, we choose a high-order function space to represent it
= fem.functionspace(domain, ("Lagrange", 5))
Q = fem.Expression(p, Q.element.interpolation_points())
expr = fem.Function(Q)
pressure pressure.interpolate(expr)
We next plot the load on the domain
= pyvista.UnstructuredGrid(*vtk_mesh(Q))
p_grid "p"] = pressure.x.array.real
p_grid.point_data[
= p_grid.warp_by_scalar("p", factor=0.5)
warped_p "p")
warped_p.set_active_scalars(
= pyvista.Plotter(off_screen=True)
load_plotter
load_plotter.add_mesh(
warped_p,=True,
show_edges=True
show_scalar_bar
)
load_plotter.add_axes()
load_plotter.view_xy()
# if not pyvista.OFF_SCREEN:
# load_plotter.show()
# HTML 저장
"fenicsx/fundamentals/membrane_load.html") load_plotter.export_html(
Create a Dirichlet boundary condition using geometrical conditions
The next step is to define the homogeneous boundary condition. Unlike in the first tutorial, here we use dolfinx.fem.locate_dofs_geometrical
to identify the degrees of freedom on the boundary. Since our domain is the unit circle, these degrees of freedom lie at coordinates \((x, y)\) satisfying \(\sqrt{x^2 + y^2} = 1\)
import numpy as np
def on_boundary(x):
return np.isclose(np.sqrt(x[0]**2 +x[1]**2), 1)
= fem.locate_dofs_geometrical(V, on_boundary) boundary_dofs
Since our Dirichlet condition is homogeneous (u=0
on the entire boundary), we can create the boundary condition with dolfinx.fem.dirichletbc
by specifying a constant value, the boundary degrees of freedom and the function space on which to apply it
= fem.dirichletbc(default_scalar_type(0), boundary_dofs, V) bc
Defining and solving the variational problem
The variational problem is the same as in our first Poisson problem, where f
is replaced by p
= ufl.TrialFunction(V)
u = ufl.TestFunction(V)
v
= ufl.dot(ufl.grad(u), ufl.grad(v)) *ufl.dx
a = p *v *ufl.dx
L
= LinearProblem(
problem
a,
L, =[bc],
bcs={"ksp_type": "preonly", "pc_type": "lu"}
petsc_options
)= problem.solve() uh
We plot the deflection \(u_h\) over the domain \(\Omega\)
# Set deflection values and add it to plotter
"u"] = uh.x.array
grid.point_data[= grid.warp_by_scalar("u", factor=25)
warped
= pyvista.Plotter(off_screen=True)
u_plotter
u_plotter.add_mesh(
warped, =True,
show_edges=True,
show_scalar_bar="u")
scalars
# if not pyvista.OFF_SCREEN:
# u_plotter.show()
# HTML 저장
"fenicsx/fundamentals/membrane_u.html") u_plotter.export_html(
Plotting along a line in the domain
Another useful way to compare the deflection and the load is by plotting them along the line \(x=0\). To do this, we simply define a set of points along the \(y\)-axis and then evaluate the finite element functions \(u\) and \(p\) at those points
= 0.001 # Avoid hitting the outside of the domain
tol = np.linspace(-1 +tol, 1 -tol, 101)
y
= np.zeros((3, 101))
points 1] = y
points[= []
u_values = [] p_values
Since a finite element function is expressed as a linear combination of all its degrees of freedom,
\[u_h(x) = \sum_{i=1}^N c_i \, \phi_i(x)\]
where \(c_i\) are the coefficients of \(u_h\) and \(\phi_i\) are the basis functions, we can, in principle, evaluate the solution at any point in \(\Omega\)
However, because a mesh usually contains a large number of degrees of freedom (i.e. \(N\) is large), it would be inefficient to evaluate every basis function at each point. Instead, we first determine which mesh cell the point \(x\) belongs to. This can be done efficiently using a bounding box tree, which allows for a fast recursive search through the mesh cells
from dolfinx import geometry
= geometry.bb_tree(domain, domain.topology.dim) bb_tree
We can now compute which cells each point collides with by using dolfinx.geometry.compute_collisions_points
. This function returns, for every input point, a list of cells whose bounding boxes overlap with that point. Since different points may correspond to a different number of cells, the results are stored in a dolfinx.cpp.graph.AdjacencyList_int32
. The cells associated with the i-th point can be accessed with links(i)
Because a cell’s bounding box covers more of \(\mathbb{R}^n\) than the cell itself, we must check whether the point actually lies in the cell. This is done with dolfinx.geometry.select_colliding_cells
, which measures the exact distance between the point and the cell (approximated as a convex hull for higher-order geometries). Like the previous function, this one also returns an adjacency list, since a point may lie on a facet, edge, or vertex shared by multiple cells
Finally, we want the code below to run in parallel when the mesh is distributed across multiple processors. In this case, not every point in points is guaranteed to be present on each processor. To handle this, we create a subset points_on_proc
that only contains the points found on the current processor
= []
cells = []
points_on_proc
# Find cells whose bounding-box collide with the the points
= geometry.compute_collisions_points(bb_tree, points.T)
cell_candidates
# Choose one of the cells that contains the point
= geometry.compute_colliding_cells(domain, cell_candidates, points.T)
colliding_cells for i, point in enumerate(points.T):
if len(colliding_cells.links(i)) > 0:
points_on_proc.append(point)0]) cells.append(colliding_cells.links(i)[
We now have a list of points assigned to the processor, along with the cell each point belongs to. Using this information, we can call uh.eval
and pressure.eval
to obtain the values of these functions at all points
= np.array(points_on_proc, dtype=np.float64)
points_on_proc
= uh.eval(points_on_proc, cells)
u_values = pressure.eval(points_on_proc, cells) p_values
Now that we have an array of coordinates together with two arrays of function values, we can use matplotlib
to plot the results
import matplotlib.pyplot as plt
= plt.figure()
fig
1], 50 *u_values,
plt.plot(points_on_proc[:, "k", linewidth=2, label="Deflection ($\\times 50$)")
1], p_values,
plt.plot(points_on_proc[:, "b--", linewidth=2, label="Load")
True)
plt.grid("y")
plt.xlabel( plt.legend()
Saving functions to file
To visualize the solution in ParaView
, we can save it to a file as follows:
from pathlib import Path
import dolfinx.io
= Path("fenicsx/fundamentals")
results_folder =True, parents=True)
results_folder.mkdir(exist_ok
= results_folder / "membrane"
filename
= "Load"
pressure.name = "Deflection"
uh.name
with dolfinx.io.VTXWriter(
MPI.COMM_WORLD, / "membrane_pressure.bp",
results_folder
[pressure], ="BP4") as vtx:
engine0.0)
vtx.write(
with dolfinx.io.VTXWriter(
MPI.COMM_WORLD, / "membrane_deflection.bp",
results_folder
[uh], ="BP4") as vtx:
engine0.0) vtx.write(
The aim of this chapter is to show how a variety of important PDEs from science and engineering can be solved with just a few lines of DOLFINx code. We begin with the heat equation, then move on to the nonlinear Poisson equation, the equations of linear elasticity, and the Navier–Stokes equations. Finally, we consider systems of nonlinear advection–diffusion–reaction equations. These examples demonstrate how to tackle time-dependent problems, nonlinear problems, vector-valued problems, and systems of PDEs. For each case, we derive the variational formulation and implement the problem in Python in a way that closely mirrors the underlying mathematics
Authors: Hans Petter Langtangen, Anders Logg, Jørgen S. Dokken
As a first extension of the Poisson problem introduced in the previous chapter, we now turn to the time-dependent heat equation (also known as the time-dependent diffusion equation). This equation can be viewed as the natural generalization of the Poisson equation, which describes the stationary distribution of heat in a body, to the case where the distribution evolves over time. By discretizing time into small intervals and applying standard time-stepping methods, we can solve the heat equation as a sequence of variational problems, in much the same way as we solved the Poisson equation
The PDE problem
The model problem for the time-dependent PDE is given by
\[ \begin{aligned} \frac{\partial u}{\partial t}&=\nabla^2 u + f && \text{in } \Omega \times (0, T] \\ u &= u_D && \text{on } \partial\Omega \times (0,T] \\ u &= u_0 && \text{at } t=0 \end{aligned}\]
Here, \(u\) depends on both space and time; for example, \(u = u(x,y,t)\) if the spatial domain \(\Omega\) is two-dimensional. The source term \(f\) and the boundary condition \(u_D\) may also vary with space and time, while the initial condition \(u_0\) is a function of space alone
The variational formulation
A simple approach to solving time-dependent PDEs with the finite element method is to first discretize the time derivative using a finite difference approximation. This reduces the problem to a sequence of stationary equations, each of which can then be written in variational form. We use the superscript \(n\) to denote a quantity at time \(t_n\), where \(n\) indexes the discrete time levels. For example, \(u^n\) represents the value of \(u\) at time level \(n\). The first step in a finite difference discretization of time is to evaluate the PDE at a chosen time level, such as \(t_{n+1}\)
\[ \begin{aligned} \left(\frac{\partial u }{\partial t}\right)^{n+1}= \nabla^2 u^{n+1}+ f^{n+1} \end{aligned}\]
The time derivative can be approximated by a difference quotient. For reasons of both simplicity and stability, we adopt the backward difference scheme:
\[\begin{aligned} \left(\frac{\partial u }{\partial t}\right)^{n+1}\approx \frac{u^{n+1}-u^n}{\Delta t} \end{aligned}\]
where \(\Delta t\) denotes the time-step size. Substituting this expression into the equation at time level \(n +1\) gives
\[\begin{aligned} \frac{u^{n+1}-u^n}{\Delta t}= \nabla^2 u^{n+1}+ f^{n+1}. \end{aligned}\]
This is the time-discrete form of the heat equation, known as the backward Euler or implicit Euler scheme
We rearrange the equation so that the left-hand side contains only the unknown \(u^{n+1}\), while the right-hand side contains terms that are already known. This yields a sequence of stationary problems for \(u^{n+1}\), given that \(u^{n}\) is available from the previous time step:
\[\begin{aligned} u^0 &= u_0\\ u^{n+1} - \Delta t \nabla^2 u^{n+1} &= u^{n} + \Delta t f^{n+1}, \quad n = 0,1,2,\dots \end{aligned}\]
Starting from the initial condition \(u_0\), we can then compute \(u^0\), \(u^1\), \(u^2\), and so forth
Next, we apply the finite element method. To do this, we first derive the weak formulation of the equation: we multiply by a test function \(v \in \hat{V}\) and integrate the second-order derivatives by parts. For simplicity, we now denote \(u^{n+1}\) by \(u\). The resulting weak formulation can be written as
\[\begin{aligned} a(u,v) &= L_{n+1}(v) \end{aligned}\]
where
\[\begin{aligned} a(u,v) &= \int_{\Omega} \big( u v + \Delta t \nabla u \cdot \nabla v \big),\mathrm{d}x \\ L_{n+1}(v) &= \int_{\Omega} \big(u^n + \Delta t f^{n+1}\big), v ,\mathrm{d}x \end{aligned}\]
Projection or interpolation of the initial condition
In addition to the variational problem solved at each time step, we also need to approximate the initial condition. This can likewise be expressed as a variational problem:
\[\begin{aligned} a_0(u,v) &= L_0(v) \end{aligned}\]
with
\[\begin{aligned} a_0(u,v) &= \int_{\Omega} u v ,\mathrm{d}x\\ L_0(v) &= \int_{\Omega} u_0 v ,\mathrm{d}x \end{aligned}\]
Solving this variational problem gives \(u^0\), which is the \(L^2\) projection of the prescribed initial condition \(u_0\) onto the finite element space
An alternative way to construct \(u^0\) is by directly interpolating the initial value \(u_0\). We discussed the use of interpolation in DOLFINx
in the membrane deflection
In DOLFINx
, the initial condition can be obtained either by projection or by interpolation. The most common approach is projection, which provides an approximation of \(u_0\). However, in applications where we want to verify the implementation against exact solutions, interpolation must be used. In this chapter, we will consider such a case
Author: Jørgen S. Dokken
Let us now consider a more interesting problem: the diffusion of a Gaussian hill. We take the initial condition to be
\[\begin{aligned} u_0(x,y) &= e^{-a (x^2 +y^2)} \end{aligned}\]
with \(a = 5\) on the domain \([-2,2]\times[-2,2]\). For this problem, we impose homogeneous Dirichlet boundary conditions (\(u_D = 0\))
The first difference from the previous problem is that the domain is no longer the unit square. We create the rectangular domain using dolfinx.mesh.create_rectangle
import numpy as np
import matplotlib as mpl
from petsc4py import PETSc
from mpi4py import MPI
from dolfinx import fem, mesh, io, plot
from dolfinx.fem.petsc import (assemble_vector, assemble_matrix,
create_vector, apply_lifting, set_bc)
import ufl
import pyvista
# Define temporal parameters
= 0.0 # Start time
t = 1.0 # Final time
T = 50
num_steps = T /num_steps # time step size
dt
# Define mesh
= 50, 50
nx, ny = mesh.create_rectangle(MPI.COMM_WORLD,
domain -2, -2]), np.array([2, 2])],
[np.array([
[nx, ny], mesh.CellType.triangle)
= fem.functionspace(domain, ("Lagrange", 1)) V
= domain.topology.dim
tdim = pyvista.UnstructuredGrid(*vtk_mesh(domain, tdim))
grid
= Path("fenicsx/heat")
results_folder =True, parents=True)
results_folder.mkdir(exist_ok
= pyvista.Plotter(off_screen=True)
plotter =True)
plotter.add_mesh(grid, show_edges
plotter.add_axes()
plotter.view_xy()
# if not pyvista.OFF_SCREEN:
# plotter.show()
# HTML 저장
"fenicsx/heat/heat_mesh.html") plotter.export_html(
Note that we have used a much higher resolution than before to better capture the features of the solution. We can also easily update the initial and boundary conditions. Instead of defining the initial condition using a class, we simply use a function
# Create initial condition
def initial_condition(x, a=5):
return np.exp(-a *(x[0]**2 +x[1]**2))
= fem.Function(V)
u_n = "u_n"
u_n.name
u_n.interpolate(initial_condition)
# Create boundary condition
= tdim -1
fdim = mesh.locate_entities_boundary(
boundary_facets
domain,
fdim, lambda x: np.full(x.shape[1], True, dtype=bool)
)
= fem.dirichletbc(
bc 0),
PETSc.ScalarType(
fem.locate_dofs_topological(V, fdim, boundary_facets),
V )
Time-dependent output
To visualize the solution in an external program such as Paraview
, we create an XDMFFile
, which can store multiple solutions. The main advantage of using an XDMFFile is that the mesh only needs to be stored once, and multiple solutions can be appended to the same grid, thereby reducing the storage requirements
The first argument to XDMFFile
specifies the communicator used to store the data. Since we want a single output file regardless of the number of processors, we use the COMM_WORLD
. The second argument is the name of the output file, and the third argument specifies the file mode, which can be read ("r"
), write ("w"
) or append ("a"
)
= io.XDMFFile(domain.comm, "fenicsx/heat/diffusion.xdmf", "w")
xdmf
xdmf.write_mesh(domain)
# Define solution variable,
# and interpolate initial solution for visualization in Paraview
= fem.Function(V)
uh = "uh"
uh.name
uh.interpolate(initial_condition) xdmf.write_function(uh, t)
Variational problem and solver
As in the previous example, we set up the necessary objects for the time-dependent problem so that we do not need to recreate the data structures
= ufl.TrialFunction(V), ufl.TestFunction(V)
u, v
= fem.Constant(domain, PETSc.ScalarType(0))
f
= u *v *ufl.dx +dt *ufl.dot(ufl.grad(u), ufl.grad(v)) *ufl.dx
a = (u_n +dt *f) *v *ufl.dx L
Preparing linear algebra structures for time dependent problems
Even though u_n
depends on time, we use the same function for both f
and u_n
at each time step. We then call dolfinx.fem.form
to create the assembly kernels for the matrix and vector
= fem.form(a)
bilinear_form = fem.form(L) linear_form
The left-hand side of the system, the matrix A
, does not change between time steps, so it only needs to be assembled once. The right-hand side, however, depends on the previous solution u_n
and must be assembled at every time step. For this reason, we create the vector b
from L
and reuse it at each step
= assemble_matrix(bilinear_form, bcs=[bc])
A
A.assemble()= create_vector(linear_form) b
Using petsc4py
to create a linear solver
Since we have already assembled a
into the matrix A
, we cannot use the dolfinx.fem.petsc.LinearProblem
class. Instead, we create a PETSc
linear solver, assign the matrix A
to it, and select a solution strategy
= PETSc.KSP().create(domain.comm)
solver
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY) solver.getPC().setType(PETSc.PC.Type.LU)
Visualization of time dependent problem using pyvista
We use the DOLFINx
plotting tools, based on PyVista
, to plot the solution at every 15-th time step. We also display a colorbar showing the minimum and maximum values of u
at each step. For this, we use the convenience function plot_function
:
= pyvista.Plotter()
plotter # conda install -c conda-forge imageio
"fenicsx/heat/heat_animation.gif", fps=10)
plotter.open_gif(
"uh"] = uh.x.array
grid.point_data[= grid.warp_by_scalar("uh", factor=3)
warped
= mpl.colormaps.get_cmap("viridis").resampled(25)
viridis = dict(
sargs =25,
title_font_size=20,
label_font_size="%.2e",
fmt="black",
color=0.1,
position_x=0.8,
position_y=0.8,
width=0.1
height
)
= plotter.add_mesh(
renderer
warped, =True,
show_edges=False,
lighting=viridis,
cmap=sargs,
scalar_bar_args=[0, max(uh.x.array)]
clim )
Updating the solution and right hand side per time step
To solve the variational problem at each time step, we must assemble the right-hand side and apply the boundary conditions before calling solver.solve(b, uh.x.petsc_vec)
. We begin by resetting the values in b
, since the vector is reused at every step. Next, we assemble the vector with dolfinx.fem.petsc.assemble_vector(b, L)
, which inserts the linear form L(v)
into b
Note that boundary conditions are not supplied during this assembly, unlike for the left-hand side. Instead, we apply them later using lifting, which preserves the symmetry of the matrix in the bilinear form without Dirichlet conditions. After applying the boundary conditions, we solve the linear system and update any values that may be shared across processors. Finally, before advancing to the next step, we update the previous solution with the current one
for i in range(num_steps):
+= dt
t
# Update the right hand side reusing the initial vector
with b.localForm() as loc_b:
set(0)
loc_b.
assemble_vector(b, linear_form)
# Apply Dirichlet boundary condition to the vector
apply_lifting(b, [bilinear_form], [[bc]])
b.ghostUpdate(=PETSc.InsertMode.ADD_VALUES,
addv=PETSc.ScatterMode.REVERSE
mode
)
set_bc(b, [bc])
# Solve linear problem
solver.solve(b, uh.x.petsc_vec)
uh.x.scatter_forward()
# Update solution at previous time step (u_n)
= uh.x.array
u_n.x.array[:]
# Write solution to file
xdmf.write_function(uh, t)
# Update plot
= grid.warp_by_scalar("uh", factor=3)
new_warped = new_warped.points
warped.points[:, :] "uh"][:] = uh.x.array
warped.point_data[
plotter.write_frame()
plotter.close() xdmf.close()
Author: Jørgen S. Dokken
Just as in the Poisson problem, we construct a test case that makes it straightforward to verify the correctness of the computations
Since our first-order time-stepping scheme is exact for linear functions in time, we design a problem with linear temporal variation combined with quadratic spatial variation. Accordingly, we choose the analytical solution
\[\begin{aligned} u = 1 + x^2+\alpha y^2 + \beta t \end{aligned}\]
which ensures that the computed values at the degrees of freedom are exact, regardless of the mesh size or time step \(\Delta t\), provided that the mesh is uniformly partitioned
Substituting this expression into the original PDE yields the right hand side \(f = \beta -2 -2\alpha\). The corresponding boundary condition is \(u_d(x,y,t)= 1 +x^2 +\alpha y^2 +\beta t\) and the initial condition is \(u_0(x,y)=1+x^2+\alpha y^2\)
We start by defining the temporal discretization parameters, along with the parameters for \(\alpha\) and \(\beta\)
import numpy as np
from petsc4py import PETSc
from mpi4py import MPI
from dolfinx import mesh, fem
from dolfinx.fem.petsc import (assemble_matrix, assemble_vector,
apply_lifting, create_vector, set_bc)
import ufl
= 0 # Start time
t = 2 # End time
T
= 20 # Number of time steps
num_steps = (T -t) /num_steps # Time step size
dt
= 3
alpha = 1.2 beta
As in the previous problem, we define the mesh and the appropriate function spaces
= 5, 5
nx, ny = mesh.create_unit_square(
domain
MPI.COMM_WORLD,
nx, ny,
mesh.CellType.triangle
)
= fem.functionspace(domain, ("Lagrange", 1)) V
Defining the exact solution
We implement a Python class to represent the exact solution
class exact_solution():
def __init__(self, alpha, beta, t):
self.alpha = alpha
self.beta = beta
self.t = t
def __call__(self, x):
return 1 +x[0]**2 +self.alpha *x[1]**2 +self.beta *self.t
= exact_solution(alpha, beta, t) u_exact
Defining the boundary condition
As in the previous sections, we define a Dirichlet boundary condition over the whole boundary
= fem.Function(V)
u_D
u_D.interpolate(u_exact)
= domain.topology.dim
tdim = tdim -1
fdim
domain.topology.create_connectivity(fdim, tdim)= mesh.exterior_facet_indices(domain.topology)
boundary_facets
= fem.dirichletbc(
bc
u_D,
fem.locate_dofs_topological(V, fdim, boundary_facets) )
Defining the variational formualation
Since we set \(t=0\) in u_exact
, we can reuse this variable to obtain \(u_n\) for the first time step
= fem.Function(V)
u_n u_n.interpolate(u_exact)
Because \(f\) is time-independent, we can treat it as a constant
= fem.Constant(domain, beta -2 -2 *alpha) f
We can now create our variational formulation, with the bilinear form a
and linear form L
= ufl.TrialFunction(V), ufl.TestFunction(V)
u, v = u *v *ufl.dx +dt *ufl.dot(ufl.grad(u), ufl.grad(v)) *ufl.dx -(u_n +dt *f) *v *ufl.dx
F
= fem.form(ufl.lhs(F))
a = fem.form(ufl.rhs(F)) L
Create the matrix and vector for the linear problem
To ensure that the variational problem is solved efficiently, we construct several structures that allow data reuse, such as matrix sparisty patterns. In particular, since the bilinear form a
is independent of time, the matrix only needs to be assembled once
= assemble_matrix(a, bcs=[bc])
A
A.assemble()= create_vector(L)
b
= fem.Function(V) uh
Define a linear variational solver
The resulting linear algebra problem is solved with PETSc, using the Python API petsc4py
to define a linear solver
= PETSc.KSP().create(domain.comm)
solver
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY) solver.getPC().setType(PETSc.PC.Type.LU)
Solving the time-dependent problem
With these structures in place, we construct our time-stepping loop. Within this loop, we first update the Dirichlet boundary condition by interpolating the updated expression u_exact
into V
. Next, we reassemble the vector b
using the current solution u_n
. The boundary condition is then applied to this vector via a lifting operation, which preserves the symmetry of the matrix. Finally, we solve the problem using PETSc
and update u_n
with the values from uh
for n in range(num_steps):
# Update Diriclet boundary condition
+= dt
u_exact.t
u_D.interpolate(u_exact)
# Update the right hand side reusing the initial vector
with b.localForm() as loc_b:
set(0)
loc_b.
assemble_vector(b, L)
# Apply Dirichlet boundary condition to the vector
apply_lifting(b, [a], [[bc]])
b.ghostUpdate(=PETSc.InsertMode.ADD_VALUES,
addv=PETSc.ScatterMode.REVERSE
mode
)
set_bc(b, [bc])
# Solve linear problem
solver.solve(b, uh.x.petsc_vec)
uh.x.scatter_forward()
# Update solution at previous time step (u_n)
= uh.x.array u_n.x.array[:]
Verifying the numerical solution
we compute the \(L^2\)-error and the maximum error at the mesh vertices for the final time step. This allows us to verify the correctness of our implementation
# Compute L2 error and error at nodes
= fem.functionspace(domain, ("Lagrange", 2))
V_ex = fem.Function(V_ex)
u_ex
u_ex.interpolate(u_exact)
= np.sqrt(
error_L2
domain.comm.allreduce(-u_ex)**2 *ufl.dx)),
fem.assemble_scalar(fem.form((uh =MPI.SUM)
op
)if domain.comm.rank == 0:
print(f"L2-error: {error_L2:.2e}")
# Compute values at mesh vertices
= domain.comm.allreduce(
error_max max(np.abs(uh.x.array -u_D.x.array)),
np.=MPI.MAX)
opif domain.comm.rank == 0:
print(f"Error_max: {error_max:.2e}")
L2-error: 2.83e-02
Error_max: 1.78e-15
Author: Jørgen S. Dokken
In this example, we solve the singular Poisson problem by incorporating information about the nullspace of the discretized system into the matrix formulation
The problem is defined as
\[\begin{aligned} -\Delta u &= f &&\text{in } \Omega\\ -\nabla u \cdot \mathbf{n} &= \mathbf{g} &&\text{on } \partial\Omega \end{aligned}\]
This problem possesses a nullspace: if \(\tilde u\) is a solution, then for any constant \(c\), \(u_c = \tilde u + c\) is also a solution
To investigate this problem, we consider a manufactured solution on the unit square, given by
\[\begin{aligned} u(x, y) &= \sin(2\pi x)\\ f(x, y) &= -4\pi^2\sin(2\pi x)\\ g(x, y) &= \begin{cases} -2\pi & \text{if } x=0,\\ \phantom{-}2\pi & \text{if } x=1,\\ \phantom{-}0 & \text{otherwise} \end{cases} \end{aligned}\]
As in previous tutorials on the Poisson problem, we define a simple wrapper function to set up the variational problem for a given manufactured solution
import typing
import numpy as np
from mpi4py import MPI
import dolfinx.fem.petsc
from dolfinx import fem, mesh
import ufl
def u_ex(mod, x):
return mod.sin(2 *mod.pi *x[0])
def setup_problem(N: int) -> typing.Tuple[fem.FunctionSpace, fem.Form, fem.Form]:
"""
Set up bilinear and linear form of the singular Poisson problem
Args:
N (int): Number of elements in each direction of the mesh
Returns:
The function space, the bilinear form and the linear form of the problem
"""
= dolfinx.mesh.create_unit_square(
domain
MPI.COMM_WORLD,
N, N, =mesh.CellType.quadrilateral
cell_type
)= fem.functionspace(domain, ("Lagrange", 1))
V
= ufl.TrialFunction(V)
u = ufl.TestFunction(V)
v
= ufl.SpatialCoordinate(domain)
x
= u_ex(ufl, x)
u_exact
= -ufl.div(ufl.grad(u_exact))
f = ufl.FacetNormal(domain)
n = -ufl.dot(ufl.grad(u_exact), n)
g
= ufl.dot(ufl.grad(u), ufl.grad(v)) *ufl.dx
F += ufl.inner(g, v) *ufl.ds
F -= f *v *ufl.dx
F
return V, *fem.form(ufl.system(F))
With the convenience function defined above, we can now address the nullspace. To this end, we use PETSc
, attaching additional information to the assembled matrices. PETSc
provides a built-in function for creating constant nullspaces, which we employ here
from petsc4py import PETSc
= PETSc.NullSpace().create(constant=True, comm=MPI.COMM_WORLD) nullspace
Direct solver
We begin by solving the singular problem using a direct solver (MUMPS). MUMPS provides additional options to handle singular matrices, which we utilize here
= {
petsc_options "ksp_error_if_not_converged": True,
"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type": "mumps",
"ksp_monitor": None,
}
Next, we configure the KSP solver
= PETSc.KSP().create(MPI.COMM_WORLD)
ksp "singular_direct")
ksp.setOptionsPrefix(= PETSc.Options()
opts
opts.prefixPush(ksp.getOptionsPrefix())
for key, value in petsc_options.items():
= value
opts[key]
ksp.setFromOptions()
for key, value in petsc_options.items():
del opts[key]
opts.prefixPop()
We then assemble the bilinear and linear forms and construct the matrix A
and the right-hand side vector b
= setup_problem(40)
V, a, L
= fem.petsc.assemble_matrix(a)
A
A.assemble()= fem.petsc.assemble_vector(L)
b =PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
b.ghostUpdate(addv ksp.setOperators(A)
We begin by verifying that this is indeed the nullspace of A
, after which we attach it to the matrix
assert nullspace.test(A)
A.setNullSpace(nullspace)
We can then solve the linear system of equations
= fem.Function(V)
uh
ksp.solve(b, uh.x.petsc_vec)
uh.x.scatter_forward()
ksp.destroy()
Residual norms for singular_direct solve.
0 KSP Residual norm 1.553142231547e+00
1 KSP Residual norm 1.542072413513e-14
<petsc4py.PETSc.KSP at 0x16d9abba0>
The \(L^2\)-error can now be evaluated against the analytical solution
def compute_L2_error(uh: fem.Function) -> float:
= uh.function_space.mesh
mesh = u_ex(ufl, ufl.SpatialCoordinate(mesh))
u_exact
= fem.form(ufl.inner(uh -u_exact, uh -u_exact) *ufl.dx)
error_L2 = fem.assemble_scalar(error_L2)
error_local return np.sqrt(mesh.comm.allreduce(error_local, op=MPI.SUM))
print(f"Direct solver L2 error: {compute_L2_error(uh):.5e}")
Direct solver L2 error: 1.59184e-03
We additionally confirm that the solution’s mean value coincides with that of the manufactured solution
= u_ex(ufl, ufl.SpatialCoordinate(V.mesh))
u_exact = V.mesh.comm.allreduce(
ex_mean *ufl.dx)), op=MPI.SUM
fem.assemble_scalar(fem.form(u_exact
)= V.mesh.comm.allreduce(
approx_mean *ufl.dx)), op=MPI.SUM
fem.assemble_scalar(fem.form(uh
)
print(f"Mean value of manufactured solution: {ex_mean:.5e}")
print(f"Mean value of computed solution (direct solver): {approx_mean:.5e}")
assert np.isclose(ex_mean, approx_mean), "Mean values do not match!"
Mean value of manufactured solution: -1.17019e-15
Mean value of computed solution (direct solver): -4.16960e-15
Iterative solver
We can also solve the problem using an iterative solver, such as GMRES
with AMG
preconditioning. To do this, we select a new set of PETSc
options and create a new KSP
solver
= PETSc.KSP().create(MPI.COMM_WORLD)
ksp_iterative "singular_iterative")
ksp_iterative.setOptionsPrefix(
= {
petsc_options_iterative "ksp_error_if_not_converged": True,
"ksp_monitor": None,
"ksp_type": "gmres",
"pc_type": "hypre",
"pc_hypre_type": "boomeramg",
"pc_hypre_boomeramg_max_iter": 1,
"pc_hypre_boomeramg_cycle_type": "v",
"ksp_rtol": 1.0e-13,
}
opts.prefixPush(ksp_iterative.getOptionsPrefix())for key, value in petsc_options_iterative.items():
= value
opts[key]
ksp_iterative.setFromOptions()
for key, value in petsc_options_iterative.items():
del opts[key]
opts.prefixPop()
Rather than defining the nullspace explicitly, we provide it as a near-nullspace to the multigrid preconditioner
= fem.petsc.assemble_matrix(a)
A_iterative
A_iterative.assemble()
A_iterative.setNearNullSpace(nullspace)
ksp_iterative.setOperators(A_iterative)
= fem.Function(V) uh_iterative
ksp_iterative.solve(b, uh_iterative.x.petsc_vec) uh_iterative.x.scatter_forward()
Residual norms for singular_iterative solve.
0 KSP Residual norm 2.661001756726e+01
1 KSP Residual norm 6.492588815947e-01
2 KSP Residual norm 1.847006602521e-02
3 KSP Residual norm 4.134969756499e-04
4 KSP Residual norm 9.983094877085e-06
5 KSP Residual norm 1.583296725038e-07
6 KSP Residual norm 4.552318917228e-09
7 KSP Residual norm 9.468343662156e-11
8 KSP Residual norm 2.090974426858e-12
When using the iterative solver, we correct the solution by subtracting its mean value and adding the mean value of the manufactured solution before evaluating the error
= V.mesh.comm.allreduce(
approx_mean *ufl.dx)), op=MPI.SUM
fem.assemble_scalar(fem.form(uh_iterative
)print(
"Mean value of computed solution (iterative solver):",
approx_mean
)
+= ex_mean -approx_mean
uh_iterative.x.array[:]
= V.mesh.comm.allreduce(
approx_mean *ufl.dx)), op=MPI.SUM
fem.assemble_scalar(fem.form(uh_iterative
)print(
"Mean value of computed solution (iterative solver) post normalization:",
approx_mean,
)print(f"Iterative solver L2 error: {compute_L2_error(uh_iterative):.5e}")
=1e-10, atol=1e-12) np.testing.assert_allclose(uh.x.array, uh_iterative.x.array, rtol
Mean value of computed solution (iterative solver): -0.1310056163277967
Mean value of computed solution (iterative solver) post normalization: -8.316194887131138e-16
Iterative solver L2 error: 1.59184e-03