Skip to content

Commit

Permalink
Update documentation
Browse files Browse the repository at this point in the history
  • Loading branch information
markuswess committed May 7, 2024
1 parent ac9dff5 commit f72c70a
Show file tree
Hide file tree
Showing 31 changed files with 766 additions and 40 deletions.
8 changes: 6 additions & 2 deletions _sources/explicit_methods/dg.md
Original file line number Diff line number Diff line change
Expand Up @@ -101,10 +101,10 @@ Following the reasoning for the 1d case we obtain the general DG formulation

````{card}
```{math}
:label: wave_eq_1o
:label: wave_eq_1o_dg
\begin{aligned}
\int_\Omega\partial_tv w+b(p,w) &= \int_\Omega f w,\\
\partial_t \int_\Omega pq-b(q,v) &=0,\\
\int_\Omega \partial_t pq-b(q,v) &=0,\\
p(0,\cdot)&=p_0,\\
v(0,\cdot)&=v_0,
\end{aligned}
Expand All @@ -114,3 +114,7 @@ for all $w\in\mathcal V_h,q\in\mathcal W_h$ with
b(p,w) = \sum_{T\in\mathcal T} \int_T \nabla p\cdot w + \int_{\partial T}(\{p\}-p)w
```
````

````{prf:theorem}
The DG formulation {eq}`wave_eq_1o` is consistent, i.e., the weak solution $p,v$ of {eq}`wave_eq_1o` also solves {eq}`wave_eq_1o_dg`
````
75 changes: 57 additions & 18 deletions _sources/explicit_methods/dg_in_ngsolve.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@
},
{
"cell_type": "code",
"execution_count": 29,
"execution_count": 7,
"id": "3bc53f84",
"metadata": {},
"outputs": [
Expand Down Expand Up @@ -47,6 +47,43 @@
"print(fes_p.ndof, fes_v.ndof)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "7279ff76",
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "2971b4273a0d4da996c8cb3978f868ea",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"WebGuiWidget(value={'gui_settings': {}, 'ngsolve_version': '6.2.2402-37-g45ed24bde', 'mesh_dim': 2, 'order2d':…"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"BaseWebGuiScene"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"gfp = GridFunction(fes_p)\n",
"gfp.vec[100]=1.\n",
"Draw(gfp)"
]
},
{
"cell_type": "markdown",
"id": "efdc411c",
Expand All @@ -59,16 +96,16 @@
},
{
"cell_type": "code",
"execution_count": 24,
"execution_count": 9,
"id": "5995a776",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Assembling and factorization: 3.9310574531555176\n",
"20 matrix applications: 3.2477149963378906\n"
"Assembling and factorization: 4.310660362243652\n",
"20 matrix applications: 3.2749738693237305\n"
]
}
],
Expand Down Expand Up @@ -102,26 +139,28 @@
},
{
"cell_type": "code",
"execution_count": 25,
"execution_count": 10,
"id": "70d2b014",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Assembling and factorization: 0.08060789108276367\n",
"20 matrix applications: 0.3024129867553711\n"
"Assembling and factorization: 0.09638404846191406\n",
"20 matrix applications: 0.30399298667907715\n"
]
}
],
"source": [
"from time import time\n",
"now = time()\n",
"\n",
"mass_p_inv = fes_p.Mass(1).Inverse()\n",
"mass_v_inv = fes_v.Mass(1).Inverse()\n",
"print(\"Assembling and factorization: \",time()-now)\n",
"\n",
"tmp_p = mass_p_inv.CreateVector()\n",
"tmp_v = mass_v_inv.CreateVector()\n",
"tmp_p.SetRandom()\n",
"tmp_v.SetRandom()\n",
"now = time()\n",
Expand All @@ -144,16 +183,16 @@
},
{
"cell_type": "code",
"execution_count": 26,
"execution_count": 11,
"id": "90687b31",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Assembling and factorization: 0.0005583763122558594\n",
"20 matrix applications: 3.0518743991851807\n"
"Assembling and factorization: 0.0006022453308105469\n",
"20 matrix applications: 3.425299644470215\n"
]
}
],
Expand Down Expand Up @@ -195,16 +234,16 @@
},
{
"cell_type": "code",
"execution_count": 27,
"execution_count": 12,
"id": "006c879d",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Assembling and factorization: 1.4245126247406006\n",
"20 matrix applications: 1.156559705734253\n"
"Assembling and factorization: 0.9520788192749023\n",
"20 matrix applications: 0.7076926231384277\n"
]
}
],
Expand All @@ -223,7 +262,7 @@
"Bel.Assemble()\n",
"\n",
"Btr = BilinearForm(trialspace=fes_tr, testspace=fes_v)\n",
"Btr += 0.5 * pT * (v_*n) * dx(element_boundary=True)\n",
"Btr += pT * (v_*n) * dx(element_boundary=True)\n",
"Btr.Assemble()\n",
"\n",
"emb_p = fes.embeddings[0]\n",
Expand All @@ -249,16 +288,16 @@
},
{
"cell_type": "code",
"execution_count": 28,
"execution_count": 13,
"id": "08d7d460",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Assembling and factorization: 0.04546666145324707\n",
"20 matrix applications: 0.25931668281555176\n"
"Assembling and factorization: 0.03208279609680176\n",
"20 matrix applications: 0.16531062126159668\n"
]
}
],
Expand Down
47 changes: 47 additions & 0 deletions _sources/explicit_methods/mass_lumping.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,47 @@
(mass_lumping)=
# Mass lumping for the wave equation

A different approach to DG methods from {numref}`dg` to obtain easily invertible mass matrices is the so called mass lumping. We consider the one-dimensional case first

(mass_lumping_1o)=
## Mass lumping for linear finite elements

We focus on the simple case of linear finite elements, first in one and subsequently in higher dimensions.

### Mass lumping for linear 1d elements

Suppose our spacial domain is given by the interval $\Omega=(a,b)$ decomposed into $N\in\mathbb N$ subintervals $T_j:=(a_{j-1},a_j):=(a_{j-1},a_{j-1}+d_j$ with $a=a_0<a_j<\ldots<a_N=b$ and $d_j=a_{j}-a_{j-1}$ for $j=1,\ldots, N$.

Then the classical hat function basis of a first order finite element space is given by
```{math}
V_h = \mathrm{span}\{b_j, j = 0,\ldots, N\},\quad b_j =\begin{cases}\frac{x-a_{j-1}}{d_j},&x\in T_j\\\frac{a_{j+1}-x}{d_{j+1}},&x\in T_{j+1}\\0,&\mathrm{else}\end{cases}
```
with the appropriate restrictions for $b_0,b_N$.

The entries of the mass matrix can be easily computed by hand via
```{math}
m_{i,j}=\int_a^b b_j b_i = \begin{cases}
\int_{a_{j-1}}^{a_{j+1}}b_j^2=(d_j+d_{j+1})\int_0^1 x^2 dx=\frac{d_j+d_{j+1}}{3},&i=j,\\
\int_{a_j}^{a_{j+1}}b_j b_{j+1}=d_{j+1}\int_0^1 x(1-x)dx = \frac{d_{j+1}}{6},&i=j+1,\\
\int_{a_{j-1}}^{a_{j}}b_j b_{j-1}=d_{j}\int_0^1 x(1-x)dx = \frac{d_{j}}{6},&i=j-1,\\
0,&\mathrm{else,}\end{cases}
```
where we set $d_0=d_{N+1}=0$.
Instead of computing the integrals for the entries of the mass matrix exactly we may approximate them using the trapezoidal rule
```{math}
\int_a^b f \approx I_{a,b}f:= \frac{b-a}{2}(f(a)+f(b)).
```
which leads to an approximated mass matrix $\tilde{\mathbf M}$ with entries
```{math}
\begin{aligned}
\tilde m_{i,j}&=I_{a,b} (b_j b_i) \\
&=\sum_{k=1}^N I_{a_{k-1},a_k}(b_j,b_i)\\
&=\sum_{k=1}^N \frac{d_k}{2}(b_j(a_{k-1})b_i(a_{k-1})+b_j(a_{k})b_i(a_{k}))\\
&=\frac{d_1}{2}b_i(a_0)b_j(a_0)+\sum_{k=1}^{N-1}\frac{d_k+d_{k+1}}{2}b_j(a_k)b_i(a_k)+\frac{d_N}{2}b_j(a_N)b_i(a_N)\\
&=
\begin{cases}
\frac{d_i+d_{i+1}}{2},& i=j,\\
0,&\mathrm{else.}\end{cases}
\end{aligned}
```
% ### Mass lumping in higher dimensions
1 change: 1 addition & 0 deletions explicit_methods.html
Original file line number Diff line number Diff line change
Expand Up @@ -207,6 +207,7 @@
<li class="toctree-l1 current active has-children"><a class="current reference internal" href="#">5. Explicit methods for time-domain waves</a><input checked="" class="toctree-checkbox" id="toctree-checkbox-5" name="toctree-checkbox-5" type="checkbox"/><label class="toctree-toggle" for="toctree-checkbox-5"><i class="fa-solid fa-chevron-down"></i></label><ul>
<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>
</ul>
</li>
</ul>
Expand Down
14 changes: 10 additions & 4 deletions explicit_methods/dg.html
Original file line number Diff line number Diff line change
Expand Up @@ -209,6 +209,7 @@
<li class="toctree-l1 current active has-children"><a class="reference internal" href="../explicit_methods.html">5. Explicit methods for time-domain waves</a><input checked="" class="toctree-checkbox" id="toctree-checkbox-5" name="toctree-checkbox-5" type="checkbox"/><label class="toctree-toggle" for="toctree-checkbox-5"><i class="fa-solid fa-chevron-down"></i></label><ul class="current">
<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>
</ul>
</li>
</ul>
Expand Down Expand Up @@ -426,7 +427,7 @@ <h2> Contents </h2>
<div class="sd-card sd-sphinx-override sd-mb-3 sd-shadow-sm docutils">
<div class="sd-card-body docutils">
<div class="math notranslate nohighlight" id="equation-wave-eq-1o">
<span class="eqno">(5.7)<a class="headerlink" href="#equation-wave-eq-1o" title="Link to this equation">#</a></span>\[\begin{split}\begin{aligned}
<span class="eqno">(5.1)<a class="headerlink" href="#equation-wave-eq-1o" title="Link to this equation">#</a></span>\[\begin{split}\begin{aligned}
\partial_tv+\nabla p &amp;= f,&amp;\text{in }&amp;[0,T_0]\times \Omega,\\
\partial_t p+\nabla\cdot v &amp;=0,&amp;\text{in }&amp;[0,T_0]\times \Omega,\\
p(0,\cdot)&amp;=p_0,&amp;\text{in }&amp;\Omega,\\
Expand Down Expand Up @@ -501,10 +502,10 @@ <h2>DG in higher dimensions<a class="headerlink" href="#dg-in-higher-dimensions"
<p>Following the reasoning for the 1d case we obtain the general DG formulation</p>
<div class="sd-card sd-sphinx-override sd-mb-3 sd-shadow-sm docutils">
<div class="sd-card-body docutils">
<div class="math notranslate nohighlight" id="equation-wave-eq-1o">
<span class="eqno">(5.7)<a class="headerlink" href="#equation-wave-eq-1o" title="Link to this equation">#</a></span>\[\begin{split}\begin{aligned}
<div class="math notranslate nohighlight" id="equation-wave-eq-1o-dg">
<span class="eqno">(5.7)<a class="headerlink" href="#equation-wave-eq-1o-dg" title="Link to this equation">#</a></span>\[\begin{split}\begin{aligned}
\int_\Omega\partial_tv w+b(p,w) &amp;= \int_\Omega f w,\\
\partial_t \int_\Omega pq-b(q,v) &amp;=0,\\
\int_\Omega \partial_t pq-b(q,v) &amp;=0,\\
p(0,\cdot)&amp;=p_0,\\
v(0,\cdot)&amp;=v_0,
\end{aligned}\end{split}\]</div>
Expand All @@ -513,7 +514,12 @@ <h2>DG in higher dimensions<a class="headerlink" href="#dg-in-higher-dimensions"
\[b(p,w) = \sum_{T\in\mathcal T} \int_T \nabla p\cdot w + \int_{\partial T}(\{p\}-p)w\]</div>
</div>
</div>
<div class="proof theorem admonition" id="theorem-1">
<p class="admonition-title"><span class="caption-number">Theorem 5.1 </span></p>
<section class="theorem-content" id="proof-content">
<p>The DG formulation <a class="reference internal" href="../time_integration/mixed_methods.html#equation-wave-eq-1o">(4.10)</a> is consistent, i.e., the weak solution <span class="math notranslate nohighlight">\(p,v\)</span> of <a class="reference internal" href="../time_integration/mixed_methods.html#equation-wave-eq-1o">(4.10)</a> also solves <a class="reference internal" href="#equation-wave-eq-1o-dg">(5.7)</a></p>
</section>
</div></section>
</section>

<script type="text/x-thebe-config">
Expand Down
56 changes: 43 additions & 13 deletions explicit_methods/dg_in_ngsolve.html

Large diffs are not rendered by default.

Loading

0 comments on commit f72c70a

Please sign in to comment.