Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update of fundamentals.md #169

Merged
merged 5 commits into from
Jan 16, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
59 changes: 30 additions & 29 deletions chapter1/fundamentals.md
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@ Authors: Hans Petter Langtangen, Anders Logg

Adapted to FEniCSx by Jørgen S. Dokken

The goal of this tutorial is to solve one of the most basic PDEs, the Poisson equations, with a few lines of code in FEniCSx. We start by introducing the most fundamental FEniCSx objects, such as `Function`, `FunctionSpace`, `TrialFunction` and `TestFunction`, and learn how to write a basic PDE solver.
The goal of this tutorial 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
Expand All @@ -17,11 +17,12 @@ The Poisson equation is the following boundary-value problem
u(\mathbf{x}) &= u_D(\mathbf{x})&& \mathbf{x} \in \partial\Omega
\end{align}

Here, $u=u(\mathbf{x})$ is the unknown function, $f=f(\mathbf{x})$ is a prescribed function, $\nabla^2$ the Laplace operator, often written as $\Delta$, $\Omega$ the spatial domain, and $\partial\Omega$ is 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 state before it makes sense to start solving it numerically with FEniCSx.
Here, $u=u(\mathbf{x})$ is the unknown function, $f=f(\mathbf{x})$ a prescribed function, $\nabla^2$ 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.

In the two dimensional space with coordinates $x$ and $y$, we can expand the Poisson equation as
In the two-dimensional space with coordinates $x$ and $y$, we can expand the Poisson equation as

$ -\frac{\partial^2 u}{\partial x^2} - \frac{\partial^2 u}{\partial y^2} = f(x,y)$
$-\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$.

Expand All @@ -35,57 +36,57 @@ Solving a boundary value problem in FEniCSx consists of the following steps:

1. Identify the computational domain $\Omega$, the PDE, and its corresponding boundary conditions and source terms $f$.
2. Reformulate the PDE as a finite element variational problem.
3. Write a Python program defining the computational domain, the boundary conditions, the variational problem and the source terms, using FEniCSx.
3. Write a Python program defining the computational domain, the boundary conditions, the variational problem, and the source terms, using FEniCSx.
4. Run the Python program to solve the boundary-value problem. Optionally, you can extend the program to derive quantities such as fluxes and averages,
and visualize the results.

As we have already covered step 1, we shall now cover step 2-4.
As we have already covered step 1, we shall now cover steps 2-4.

## Finite element variational formulation

FEniCSx is based on the finite element method, which is a general and
efficient mathematical machinery for the numerical solution of
PDEs. The starting point for the finite element methods is a PDE
expressed in _variational form_. For readers not familiar with variational problems, it is suggested to reading a proper book on the finite element method in addition, as this tutorial is meant as a brief introduction to the subject. See the original tutorial {cite}`FenicsTutorial` (Chapter 1.6.2).
efficient mathematical technique for the numerical solution of
PDEs. The starting point for finite element methods is a PDE
expressed in _variational form_. For readers not familiar with variational problems, we suggest reading a proper treatment on the finite element method, as this tutorial is meant as a brief introduction to the subject. See the original tutorial {cite}`FenicsTutorial` (Chapter 1.6.2).

The basic recipe for turning a PDE into a variational problem is:
- Multiply the PDE by a function $v$
- Integrate the resulting equation over the domain $\Omega$
- Perform integration by parts of terms with second order derivatives.
- Perform integration by parts of those terms with second order derivatives

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 equation by a test function $v$ and integrate over $\Omega$:
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.$
$\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 the boundary of $\Omega$.
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 is that when we derive variational formulations is that one tries to keep the order of derivatives of $u$ and $v$ as small as possible.
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](https://en.wikipedia.org/wiki/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_\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 n$ is the derivative of $u$ in the outward normal direction $n$ on the boundary.
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. See for instance {cite}`Langtangen_Mardal_FEM_2019`.

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 is zero, and we have that
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.$
$\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:
$\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}.$
$\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 test and trial spaces $V$ and $\hat{V}$ is defined as
For the present problem, the trial and test spaces $V$ and $\hat{V}$ are defined as
\begin{align}
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 \}.
Expand All @@ -99,28 +100,28 @@ 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
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{align}
\int_\Omega \nabla u_h \cdot \nabla v ~\mathrm{d} x &= \int_\Omega fv~ \mathrm{d} x && \forall v \in \hat{V}_h.
\int_\Omega \nabla u_h \cdot \nabla v~\mathrm{d} x &= \int_\Omega fv~\mathrm{d} x && \forall v \in \hat{V}_h.
\end{align}
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 looks the same.
but means that the finite element variational problem and the continuous variational problem look the same.

## Abstract finite element variational formulation

We will introduce the following notations for variational problems:
We will introduce the following notation for variational problems:
Find $u\in V$ such that
\begin{align}
a(u,v)&=L(v)&& \forall v \in \hat{V}.
\end{align}
For the Poisson equation, we have:
\begin{align}
a(u,v) &= \int_{\Omega} \nabla u \cdot \nabla v ~\mathrm{d} x,\\
L(v) &= \int_{\Omega} fv ~\mathrm{d} x.
a(u,v) &= \int_{\Omega} \nabla u \cdot \nabla v~\mathrm{d} x,\\
L(v) &= \int_{\Omega} fv~\mathrm{d} x.
\end{align}
From literature $a(u,v)$ is known as the _bilinear form_ and $L(v)$ as a _linear form_.
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:
Expand Down
12 changes: 6 additions & 6 deletions chapter1/membrane.md
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
# Deflection of a membrane
Authors: Hans Petter Langtangen and Anders Logg.

Modified to DOLFINx by Jørgen S. Dokken
Modified for DOLFINx by Jørgen S. Dokken

In the first FEniCSx program, we solved a simple problem which we could easily use to verify the implementation.
In this section, we will turn our attentition to a physically more relevant problem with solutions of a somewhat more exciting shape.
Expand All @@ -12,13 +12,13 @@ We would like to compute the deflection $D(x,y)$ of a two-dimensional, circular
\end{align}
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{align}
p(x,y)&=\frac{A}{2\pi\sigma}e^{-\frac{1}{2}\left(\frac{x-x_0}{\sigma}\right)^2-\frac{1}{2}\left(\frac{y-y_0}{\sigma}\right)^2}
p(x,y)&=\frac{A}{2\pi\sigma}e^{-\frac{1}{2}\left(\frac{x-x_0}{\sigma}\right)^2-\frac{1}{2}\left(\frac{y-y_0}{\sigma}\right)^2}.
\end{align}
The parameter $A$ is the amplitude of the pressure, $(x_0, y_0)$ the localization 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$.
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{align}
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)}
+\left(\frac{y-R_0}{\sigma}\right)^2\right)}.
\end{align}
## Scaling the equation

Expand All @@ -36,7 +36,7 @@ With $D_e=\frac{AR^2}{8\pi\sigma T}$ and dropping the bars we obtain the scaled
\begin{align}
-\nabla^2 w = 4e^{-\beta^2(x^2+(y-R_0)^2)}
\end{align}
to be solved over the unit disc with $w=0$ on the boundary. Now there are only two parameters to vary the dimensionless extent of the pressure, $\beta$, and the localization of the pressure peak, $R_0\in[0,1]$. As $\beta\to 0$, the solution will approach the special case $1-x^2-y^2$. Given a computed scaed solution $w$, the physical deflection can be computed by
to be solved over the unit disc with $w=0$ on the boundary. Now there are only two parameters which vary the dimensionless extent of the pressure, $\beta$, and the location of the pressure peak, $R_0\in[0,1]$. As $\beta\to 0$, the solution will approach the special case $w=1-x^2-y^2$. Given a computed scaed solution $w$, the physical deflection can be computed by
\begin{align}
D=\frac{AR^2}{8\pi\sigma T}w
D=\frac{AR^2}{8\pi\sigma T}w.
\end{align}
2 changes: 1 addition & 1 deletion chapter2/intro.md
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
# A Gallery of finite element solvers

The goal of this chapter is to demonstrate how a range of important PDEs from science and engineering can be quickly solved with a few lines of DOLFINx code.
We will start with the heat equation, then continue with the non-linear Poisson equation, the equations for linear elasticity, the Navier-Stokes equations, and finally look at how to solve systems of nonlinear advection-diffusion-reaction equations. These problems illustrate how to solve time-dependent problems, non-linear problems, vector-values problems and systems of PDEs. For each problem, we derive the variational formulation and express the problem in Python in a way that closely resembels the mathematics.
We will start with the heat equation, then continue with the nonlinear Poisson equation, the equations for linear elasticity, the Navier-Stokes equations, and finally look at how to solve systems of nonlinear advection-diffusion-reaction equations. These problems illustrate how to solve time-dependent problems, nonlinear problems, vector-valued problems and systems of PDEs. For each problem, we derive the variational formulation and express the problem in Python in a way that closely resembles the mathematics.
20 changes: 10 additions & 10 deletions chapter2/nonlinpoisson.md
Original file line number Diff line number Diff line change
@@ -1,18 +1,18 @@
# A nonlinear Poisson equation
Authors: Anders Logg and Hans Petter Langtangen

We shall now address how to solve non-linear PDEs. We will see that non-linear problems introduce some subtle differences on how we define the variational form.
We shall now address how to solve nonlinear PDEs. We will see that nonlinear problems introduce some subtle differences in how we define the variational form.

## The PDE problem
As a model for the solution of non-linear PDEs, we take the following non-linear Poisson equation
As a model for the solution of nonlinear PDEs, we take the following nonlinear Poisson equation
\begin{align}
- \nabla \cdot (q(u) \nabla u)&=f && \text{in } \Omega,\\
u&=u_D && \text{on } \partial \Omega,
u&=u_D && \text{on } \partial \Omega.
\end{align}
and the coefficients $q(u)$ makes the problem non-linear (unless q(u) is constant in $u$).
The coefficients $q(u)$ make the problem nonlinear (unless $q(u)$ is constant in $u$).

## Variational formulation
As usual, we multiply the PDE by a test function $v\in \hat{V}$, integrate over the domain, and integrate second-order derivatives by parts. The boundary integrals arising from integration by parts vanishes wherever we employ Dirichlet conditions. The resulting variational formulation of our model problem becomes:
As usual, we multiply the PDE by a test function $v\in \hat{V}$, integrate over the domain, and integrate second-order derivatives by parts. The boundary integrals arising from integration by parts vanish wherever we employ Dirichlet conditions. The resulting variational formulation of our model problem becomes:

Find $u\in V$ such that
\begin{align}
Expand All @@ -24,14 +24,14 @@ where
\end{align}
and
\begin{align}
V&=\left\{v\in H^1(\Omega)\vert v=u_D \text{on } \partial \Omega \right\}\\
\hat{V}&=\left\{v\in H^1(\Omega)\vert v=0 \text{on } \partial \Omega \right\}
V&=\left\{v\in H^1(\Omega)\vert v=u_D \text{ on } \partial \Omega \right\}\\
\hat{V}&=\left\{v\in H^1(\Omega)\vert v=0 \text{ on } \partial \Omega \right\}
\end{align}

The discrete problem arises as usual by restricting $V$ and $\hat{V}$ to a pair of discrete spaces. The discrete non-linear problem can therefore be written as:
The discrete problem arises as usual by restricting $V$ and $\hat{V}$ to a pair of discrete spaces. The discrete nonlinear problem can therefore be written as:

Find $u_h \in V_h$ such that
\begin{align}
F(u_h, v) &=0 \quad \forall v \in \hat{V}_h,
F(u_h; v) &=0 \quad \forall v \in \hat{V}_h,
\end{align}
with $u_h=\sum_{j=1}^N U_j\phi_j$. Since $F$ is non-linear in $u$, the variational statement gives rise to a system of non-linear algebraic equation in the unknowns $U_1,\dots,U_N$.
with $u_h=\sum_{j=1}^N U_j\phi_j$. Since $F$ is nonlinear in $u$, the variational statement gives rise to a system of nonlinear algebraic equations in the unknowns $U_1,\dots,U_N$.