Skip to content

Commit

Permalink
Update documentation
Browse files Browse the repository at this point in the history
  • Loading branch information
markuswess committed May 10, 2024
1 parent f72c70a commit 53998f5
Show file tree
Hide file tree
Showing 39 changed files with 2,243 additions and 10 deletions.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
21 changes: 21 additions & 0 deletions _sources/explicit_methods/exercises.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,21 @@
(exercises_expl)=
# Exercises for explicit methods

## Exercise 1
Estimate the CFL condition of the Leap-Frog timestepping for the DG system. To this end estimate the largest eigenvalue of $\mathbf M^{-1}\mathbf K$ by $\mu_N$ for large enough $N$ given by a power iteration
```{math}
\mathbf x_{n+1}=\mathbf M^{-1}\mathbf K\mathbf x_n,\quad \mu_n:=\frac{\mathbf x_n^\top\mathbf M^{-1}\mathbf K\mathbf x_n}{\|\mathbf x_n\|}
```
with a random starting vector $\mathbf x_n$.


## Exercise 2
Implement the DG-Method for the acoustic wave equation on the `unit_square` using suitable initial and boundary conditions and the Leap-Frog time-stepping. Choose the time-step as large as possible to satisfy the CFL-condition.


## Exercise 3
Implement the first order mass-lumping method by supplying the according `IntegrationRule` to the differential symbol `dx`. Explore the sparsity pattern of the resulting mass matrix. Experiment with the argument `diagonal = True` for the mass `BilinearForm`.
Replace the `H1` space by `H1LumpingFESpace` (only for orders $1,2$). The `IntegrationRule`s can by obtained by `H1LumpingFESpace.GetIntegrationRules()`.

## Exercise 4
Use a suitable example to compare the efficiency of the explicit methods (Mass-Lumping, DG) to the efficiency of the conforming method (second order system, implicit time-stepping). Plot the respective errors against the computation times.
528 changes: 528 additions & 0 deletions _sources/explicit_methods/intrules.ipynb

Large diffs are not rendered by default.

125 changes: 124 additions & 1 deletion _sources/explicit_methods/mass_lumping.md
Original file line number Diff line number Diff line change
Expand Up @@ -44,4 +44,127 @@ which leads to an approximated mass matrix $\tilde{\mathbf M}$ with entries
0,&\mathrm{else.}\end{cases}
\end{aligned}
```
% ### Mass lumping in higher dimensions
### Mass lumping in higher dimensions

For higher dimensions the first order finite element spaces (for triangular or tetrahedral elements $T\in\mathcal T$) are again given by piecewise linear functions. A convenient basis is given by the according functions which vanish on all but one vertex points, i.e.,

```{math}
V_h = \{f\in H^1(\Omega): f|_T\in \mathcal P^1,\forall T\in\mathcal T\}
```
where the polynomial space $\mathcal P^k$ is defined by
```{math}
\mathcal P^k:=\{(x,y,z)\mapsto x^ny^mz^l:n+m+l\leq k\},
```
(and in two dimensions the exponent $l$ of the $z$ coordinate is always $0$).
If we denote the set of vertices by $V_0,\ldots,V_N$ the nodal basis is given by the condition
```{math}
b_j(V_i)=\delta_{j,i}.
```
Following the same reasoning as in one dimension we may construct integration rules locally on each element $T$ by
```{math}
I_T(f):=\frac{|T|}{d+1}\sum_{i=0}^{d+1} f(V_{T,i})
```
where $V_{T_i}$ are the vertices of the element $T$ and $d$ is the spacial dimension.
Again an approximation to the correct mass integrals is given by

```{math}
\tilde m_{i,j}:=\sum_{T\in\mathcal T} I_T(b_ib_j)=\delta_{i,j}\frac{1}{d+1}\sum_{k=1}^{K_i} |T_{i,k}|
```
where $T_{i,k},k=1,\ldots,K_i$ are the elements sharing the vertex $V_i$ and $K_i$ is the number of these elements.

## Higher order mass lumping

Usually finite element basis functions are defined on a reference element first. Then they are appropriately transformed to the physical elements and basis functions on neighbouring elements are glued together to obtain conforming basis functions. I.e., if $T$ is a physical element (interval, triangle, quadrilateral, tetrahedron, hexahedron,...) we have a reference element $\hat T$ and a diffeomorphism $\phi:\hat T\to T$.
If a local basis function $\hat b:\hat T\to\mathbb R$ is given then the basis function on the phyiscal element is given by $b := \hat b \circ \phi ^{-1}$.
We use the following reference elements:
```{math}
\begin{aligned}
\hat T_{\mathrm{segm}} &= \mathrm{conv}\{0,1\}=[0,1],\\
\hat T_{\mathrm{trig}} &= \mathrm{conv}\{(0,0),(1,0),(0,1)\},\\
\hat T_{\mathrm{tet}} &= \mathrm{conv}\{(0,0,0),(1,0,0),(0,1,0),(0,0,1)\} ,\\
\hat T_{\mathrm{quad}} &= \mathrm{conv}\{(0,0),(1,0),(0,1),(1,1)\}=[0,1]^2,\\
\hat T_{\mathrm{hex}} &= \mathrm{conv}\{(0,0,0),(1,0,0),(0,1,0),(1,1,0),(0,0,1),(1,0,1),(0,1,1),(1,1,1)\}=[0,1]^3,
\end{aligned}
```
where $\mathrm{conv}$ denotes the convex hull of a set.

For a given integration rule consisting of integration points $\hat x_j\in \hat T$ we use basis functions which vanish in all but one integration points, i.e.,
```{math}
\hat b_i(\hat x_j)=\delta_{i,j}.
```
The task is now to find suitable integration rules and spaces on the reference elements.

### Higher order mass lumping in one dimension

To be able to couple the mapped basis functions on two adjoining elements (intervals) we need to fix the endpoints $0,1$ of the reference interval $T_{\mathrm{segm}}$ as integration points. Then for a given total number of points $N+1$ the integration rule with highest order is the Gauss-Lobatto rule which is exact for polynomials of up to order $2N-1$.

As local spaces we may choose the spaces of polynomials of up to order $N$ with the Lagrangian basis
```{math}
\hat b_j(x):= \prod_{k=0,k\neq j}^N\frac{x-\hat x_k}{\hat x_j-\hat x_k},
```
where $\hat x_j$ denote the integration points on the reference element. The basis functions are then mapped to the physical elements as described above and the first basis function of an element is glued to the last basis function of the previous element to obtain globally continuous basis functions.

### Higher order mass lumping for quadrilaterals and hexahedra

In the case of quadrilaterals and hexahedra integration rules may be obtained by tensorizing the one dimensional Gauss-Lobatto rules. I.e., if a Gauss-Lobatto rule is given by the points $\hat x_0,\ldots,\hat x_N$ and weights $\hat w_0,\ldots, \hat w_N$ we define
```{math}
\hat {\mathbf x}_{i+(N+1)j}:=(\hat x_i,\hat x_j)^\top,\quad \hat {\mathbf w}_{i+(N+1)j}:=\hat w_i\hat w_j,
```
on the unit square.
The according Lagrangian polynomial basis is given by
```{math}
\hat{\mathbf b}_{i+(N+1)j}(x,y) = \hat b_i(x)\hat b_j(y).
```
Note that here we use the polynomial space
```{math}
\mathcal Q^N:=\{x^ky^l:k,l\leq N\}
```
which differs from $\mathcal P^N$.
For three dimensions the construction works similarly.

### Higher order mass lumping on simplexes

A naive ansatz to construct a second order mass-lumping method for triangles would be to follow the ideas from the first-order case: There we used the given basis and constructed a suitable integration rule. Since the dimension of $\mathcal P^2$ is $6$ and we need three symmetric degrees of freedom on each boundary to ensure continuity of the correctly coupled basis functions. The integration points on the reference triangle are fixed by
```{math}
\begin{aligned}
\hat{\mathbf x}_1&=(0,0),&
\hat{\mathbf x}_2&=(1/2,0),&
\hat{\mathbf x}_3&=(1,0),\\
\hat{\mathbf x}_4&=(1/2,1/2),&
\hat{\mathbf x}_5&=(0,1),&
\hat{\mathbf x}_6&=(0,1/2)\}.
\end{aligned}
```
It remains to specify the integration weights to obtain an optimal accuracy. Due to symmetry we only have to specify two weights $w_2=w_4=w_6=w_E$ for the mid-points of the edges and $w_1=w_3=w_5=w_V$ for the vertices, i.e., our integration rule is given by
```{math}
I_{\hat T}(f)=w_V\left(f(\hat{\mathbf x}_1)+f(\hat{\mathbf x}_3)+f(\hat{\mathbf x}_5)\right)+w_E\left(f(\hat{\mathbf x}_2)+f(\hat{\mathbf x}_4)+f(\hat{\mathbf x}_6)\right).
```
The condition that constants are integrated exactly yields
```{math}
\frac{1}{2}=3w_V+3w_E.
```
Imposing that also the function $(x,y)\mapsto x$ is integrated exactly yields
```{math}
\int_{\hat T}xdx=\frac{1}{6}=w_V\left(0+1+0\right)+w_E\left(\frac{1}{2}+\frac{1}{2}+0\right)=w_V+w_E,
```
which is the same condition as before. Imposing that also the function $(x,y)=x^2$ is integrated exactly yields
```{math}
\int_{\hat T} x^2 dx = \int_0^1\int_0^{1-x} x^2 dydx = \int_0^1x^2-x^3 dx = \frac{1}{12}
=w_V\left(0+1+0\right)+w_E\left(\frac{1}{4}+\frac{1}{4}+0\right)=w_V+\frac{1}{2}w_E.
```
Inserting the first condition $w_E=1/6-w_V$ yields
```{math}
w_V = 0,\quad w_E = \frac{1}{6}.
```
Using this quadrature for approximating our mass matrix is not an option since we would end up with a singular approximated mass matrix. We may settle for an integration rule which only integrates linear functions exactly, however it turns out that this will also reduce the overall convergence order of the resulting method.

The remedy to this problem is to add an additional integration point to the integration rule. From symmetry considerations this can only be the center $\hat{\mathbf x}_7=(1/3,1/3)$ of the reference triangle.

We also have to enlargen the local space of polynomials by another basis function.

Since this has to be a basis function corresponding to the integration point in the center we have to choose a bubble function namely $\hat b_7(x,y):= xy(1-x-y)\}.$
The resulting space
```{math}
\hat V_h:=\mathcal P^2\oplus \{\hat b_7\}.
```
now lies in between the polynomial spaces $\mathcal P^2$ and $\mathcal P^3$. To obtain a Lagrangian basis with respect to the integration points one may simply use the Lagrangian basis of $\mathcal P^2$ with respect to the original integration points and subtract the bubble function.
2 changes: 1 addition & 1 deletion _sources/time_integration/exercises.md
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
(exercises_th)=
(exercises_ti)=
# Exercises for time-stepping schemes

## Exercise 1
Expand Down
2 changes: 1 addition & 1 deletion _sources/time_integration/second_order.md
Original file line number Diff line number Diff line change
Expand Up @@ -256,7 +256,7 @@ which can be rewritten by
Moreover, let $\phi_n,\lambda_n$ be the normalized eigenpairs of
```{math}
:label: disc_evp
\mathbf K u=\lambda \mathbf M
\mathbf K \mathbf u=\lambda \mathbf M\mathbf u
```
Then we have
```{math}
Expand Down
2 changes: 2 additions & 0 deletions explicit_methods.html
Original file line number Diff line number Diff line change
Expand Up @@ -208,6 +208,8 @@
<li class="toctree-l2"><a class="reference internal" href="explicit_methods/dg.html">5.1. Discontinuous Galerkin methods for the wave equation</a></li>
<li class="toctree-l2"><a class="reference internal" href="explicit_methods/dg_in_ngsolve.html">5.2. DG Operators in NGSolve</a></li>
<li class="toctree-l2"><a class="reference internal" href="explicit_methods/mass_lumping.html">5.3. Mass lumping for the wave equation</a></li>
<li class="toctree-l2"><a class="reference internal" href="explicit_methods/intrules.html">5.4. Integration rules in NGSolve</a></li>
<li class="toctree-l2"><a class="reference internal" href="explicit_methods/exercises.html">5.5. Exercises for explicit methods</a></li>
</ul>
</li>
</ul>
Expand Down
2 changes: 2 additions & 0 deletions explicit_methods/dg.html
Original file line number Diff line number Diff line change
Expand Up @@ -210,6 +210,8 @@
<li class="toctree-l2 current active"><a class="current reference internal" href="#">5.1. Discontinuous Galerkin methods for the wave equation</a></li>
<li class="toctree-l2"><a class="reference internal" href="dg_in_ngsolve.html">5.2. DG Operators in NGSolve</a></li>
<li class="toctree-l2"><a class="reference internal" href="mass_lumping.html">5.3. Mass lumping for the wave equation</a></li>
<li class="toctree-l2"><a class="reference internal" href="intrules.html">5.4. Integration rules in NGSolve</a></li>
<li class="toctree-l2"><a class="reference internal" href="exercises.html">5.5. Exercises for explicit methods</a></li>
</ul>
</li>
</ul>
Expand Down
6 changes: 4 additions & 2 deletions explicit_methods/dg_in_ngsolve.html

Large diffs are not rendered by default.

Loading

0 comments on commit 53998f5

Please sign in to comment.