Praktikum zu Vorlesung Modellreduktion parametrisierter Systeme
Mario Ohlberger, Felix Schindler
Aktivieren Sie wie gewohnt ihre Arbeitsumgebung und starten Sie den Jupyter Notebook server, siehe zB Blatt 1, Aufgabe 0.
Erstellen Sie ein neues Python 3 notebook oder laden Sie dieses von der Homepage herunter.
Importieren Sie numpy und pymor.basic und machen Sie matplotlib für das Notebook nutzbar.
Erstellen Sie eine neue Textdatei und benennen Sie diese in interpolations.py um.
Kopieren Sie die Funktionen interpolate_lagrange_p1 und interpolate_fv aus den letzten Blättern in diese neue Datei, sodass diese Funktionen in Zukunft wiederverwendet werden können.
Kopieren Sie insbesondere alle dazu nötigen imports.
Ändern Sie die interpolate_lagrange_p1, sodass anstelle des bisherigen NumpyVectorSpace der CGVectorSpace aus dem pymor.operators.cg Modul benutzt wird; entsprechend in interpolate_fv der FVVectorSpace aus dem pymor.operators.fv Modul.
Führen Sie zum Test from interpolations import (interpolate_fv, interpolate_lagrange_p1) in Ihrem Notebook aus.
Erstellen Sie eine neue Textdatei und benennen Sie diese in visualizations.py um.
Kopieren Sie analog die Funktionen visualize_fv, visualize_grid, visualize_lagrange_p1 sowie alle dazu nötigen Imports und testen Sie den Import dieser Fnuktionen in Ihrem Notebook.
Wir betrachten ein vereinfachtes parameter-unabhängiges Diffusionsproblem. Sei dazu
gegeben. Gesucht ist eine schwache Lösung $u \in H^1_0(\Omega)$, sodass
$$\begin{align} -\nabla\cdot( A \nabla u ) &= f &&\text{in } \Omega\\ u &= 0 &&\text{auf } \partial \Omega, \end{align}$$im schwachen Sinne erfüllt ist, also äquivalent: wir suchen $u \in H^1_0(\Omega)$ als Lösung des Variationsproblems
$$\begin{align} b(u, v) &= l(v) &&\text{für alle } v \in H^1_0(\Omega) \end{align}$$mit der Bilinearform $b: H^1(\Omega) \times H^1(\Omega) \to \mathbb{R}$ definiert durch
$$b(u, v) := \int_\Omega (A \nabla u) \cdot \nabla v \text{d}x$$und dem linearen Funktional $l: H^1(\Omega) \to \mathbb{R}$ definiert durch
$$l(v) := \int_\Omega f v \text{d}x.$$Beachte, dass sowohl $b$ als auch $l$ bezüglich $H^1(\Omega)$ definiert werden können, die Lösung des Variationsproblems aber nur in $H^1_0(\Omega)$ gesucht wird.
In der Praxis werden wir die Integrale in der Definition von $b$ und $l$ durch numerische Quadraturen ersetzen, die für allgemeine (nicht-polinomielle Datenfunktionen $A$ und $f$) nicht exakt integrieren. Wir erhalten dadurch Approximationen $b_h$ von $b$ und $l_h$ von $l$. Wir werden aber im Folgenden nicht zwischen $b$ und $l$ und ihren Approximationen unterscheiden!
Seien $V$, $W$ reelle Hilberträume. Dann kann jede Bilinearform $b: W \times V \to \mathbb{R}$ eindeutig mit einem linearen Operator $L: V \to W'$ identifiziert werden, wobei $W'$ der Dualraum zu $W$ ist, d.h.
$$B: V \to W', \quad v\mapsto B[v] \in W'$$wobei das Funktional $B[v]$ definiert ist durch
$$B[v]: W \to \mathbb{R},\quad w \mapsto B[v](w) := b(w, v).$$Äquivalent kann jeder lineare Operator $B: V \to W'$ eindeutig mit einer Bilinearform
$$b: W \times V \to \mathbb{R},\quad (w, v) \mapsto b(w, v) := B[v](w)$$identifiziert werden.
Zum Beispiel entsprict $b$ wie oben definiert für konstante Diffusion $A = 1$ dem (schwachen) Laplace Operator $\Delta$ bezüglich der Hilberträume $V = W = H^1(\Omega)$.
Für konstante Diffusion $A = 1$ entspricht die Bilinearform $b$ dem $H^1$-semi Produkt und induziert die $H^1$ Halbnorm
$$|v|_{H^1} := \sqrt{\int_\Omega \nabla v \cdot \nabla v \;\text{d}x} = \sqrt{b(v, v)}\quad\quad\text{für } v \in H^1(\Omega).$$Eingeschränkt auf $H_0^1(\Omega) \subset H^1(\Omega)$ hingegen entspricht $b$ (für $A = 1$) dem $H_0^1$ Produkt und induziert damit die $H^1_0$ Norm (die aufgrund der Poincare Ungleichung eine Norm ist), also
$$\|v\|_{H^1_0} := |v|_{H^1}\quad\quad\text{für } v \in H^1_0(\Omega).$$Für konstante Diffusion $A = 1$ entspricht die Bilinearform $b$ dem $H^1$-semi Produkt (bzw. der associierte Operator $B$ dem schwachen Laplace Operator $\Delta$) und induziert die $H^1$ Halbnorm
$$|v|_{H^1} := \sqrt{\int_\Omega \nabla v \cdot \nabla v \;\text{d}x} = \sqrt{b(v, v)}\quad\quad\text{für } v \in H^1(\Omega).$$Eingeschränkt auf $H_0^1(\Omega) \subset H^1(\Omega)$ hingegen entspricht $b$ (für $A = 1$) dem $H_0^1$ Produkt und induziert damit die $H^1_0$ Norm (die aufgrund der Poincare Ungleichung eine Norm ist), also
$$\|v\|_{H^1_0} := |v|_{H^1}\quad\quad\text{für } v \in H^1_0(\Omega).$$visualizations.py.CGVectorSpace aus dem pymor.operators.cg Modul und vergleichen Sie dessen Dimension mit der Anzahl der Knoten des Gitters.Generell werden in pyMOR Operatoren, Bilinearformen und Produkte durch Operatoren dargestellt, die von pymor.operators.interfaces.OperatorInterface ableiten.
Auf einem TriaGrid wird die Bilinearform $b$ von oben für beliebige Diffusion $A$ durch den DiffusionOperatorP1 aus dem pymor.operators.cg Modul repräsentiert.
Nutzen Sie den DiffusionOperatorP1, um den Laplace Operator $\Delta$ (und damit das $H^1$-semi Produkt) bezüglich $S_h^1 \subset H^1(\Omega)$ zu definieren (nicht $H^1_0(\Omega)$!). Nutzen Sie dazu in der pyMOR Dokumentation in der API Dokumentation im operators Paket im cg Modul die entsprechende Hilfe oder zeigen Sie sich die Dokumentation im Notebook an.
Achtung: der DiffusionOperatorP1 beinhaltet eine automatische Randwertbehandlung (mehr dazu unten), stellen Sie sicher, dass Sie den Operator bezüglch $H^1(\Omega)$ und nicht $H^1_0(\Omega)$ anlegen.
Suchen Sie dazu in der pyMOR Dokumentation in der API Dokumentation im grids Paket im boundaryinfos Modul nach einer geeigneten Möglichkeit.
Lassen Sie sich den erstellten Operator sowie dessen Typ anzeigen. Ist dieser Operator linear?
cg_space übereinstimmen.Jeder lineare Operator $B: V \to W$ zwischen endlichdimensionalen Hilberträumen
$V$ mit Dimension $J := \dim(V)$ und Basis $\big\{ \varphi_1, \dots, \varphi_J \big\}$ (der Ansatzraum)
$W$ mit Dimension $I := \dim(W)$ und Basis $\big\{ \psi_i, \dots, \psi_I \big\}$ (der Testraum)
(bzw. jede Bilineaform $b: W \times V \to \mathbb{R}$) lässt sich eindeutig mit einer Matrix $\underline{b} \in \mathbb{R}^{I \times J}$ identifizieren (der Basisdarstellung von $b$):
$$\underline{b}_{i, j} := b(\psi_i, \varphi_j)$$Analog lässt sich jedes Funktional $l: W \to \mathbb{R}$ eindeutig mit einem Vektor $\underline{l} \in \mathbb{R}^I$, gegeben durch $\underline{l}_i := l(\psi_i)$, identifizieren.
Gegeben diskrete Funktionen $v_h \in V$ und $w_h \in W$ mit DoF-Vektoren $\underline{v_h} \in \mathbb{R}^J$ und $\underline{w_h} \in \mathbb{R}^I$ lässt sich außerdem die Anwendung von $B$ bzw. $b$ durch Matrix/Vektor Multiplikation darstellen:
$$B[v_h](w_h) = b(w_h, v_h) = \underline{w_h}^\top\; \underline{b}\; \underline{v_h}$$Wir betrachten als Approximaiton des Hilbertraums $V = W = H^1(\Omega)$ den endlichdimensionalen Hilbertraum $S_h^1 \subset H^1(\Omega)$ von Blatt 1, Aufgabe 2 zu einem gegebenen Gitter $\mathcal{T}_h$ von $\Omega$, mit Dimension $I := |\mathcal{T}_h^d|$ (Anzahl der Knotenpunkte des Gitters) und entsprechender Lagrange Basis.
Viele Operator in pyMOR bieten die Möglichkeit, ihre Matrixdarstellung zu assemblieren (wenn Sie nicht von einem Parameter abhängen). Suchen Sie in der Dokumentation von OperatorInterface (alle Operatoren in pyMOR leiten davon ab) nach einer entsprechenden möglichkeit.
Berechnen Sie $b(f_h, f_h)$, indem Sie die Matrixdarstellung nutzen.
f_h anzeigen lassen, Dokumentation nutzen).Gegeben $S_h^1$ als Approximation von $H^1(\Omega)$ zu einem Gitter $\mathcal{T}_h$, wie erhalten wir eine Approximation von $H^1_0(\Omega)$?
Sei dazu $\mathcal{T}_h^d$ die Menge der Knoten des Gitters und sei $\mathcal{T}_h^{\partial\Omega} \subset \mathcal{T}_h^d$ die Teilmenge derjenigen Knoten, die auf dem Gebietsrand $\partial \Omega$ liegen.
Haben wir zu $v_h \in S_h^1$ den DoF Vektor $\underline{v_h} \in \mathbb{R}^I$ gegeben, so gilt
$$v_h \in S_h^1 \cap H^1_0(\Omega) \quad\Leftrightarrow\quad \underline{v_h}_i = 0\quad\text{für alle } \nu_i \in \mathcal{T}_h^{\partial \Omega}$$
Haben wir die Matrix $\underline{b_h}$ zu einer Bilinearform $b_h$ bezüglich der Basis von $S_h^1 \subset H^1(\Omega)$ gegeben, so sind die Zeilen der Matrix mit dem zweiten Argument von $b_h$ assoziiert (Funktionen aus dem Ansatzraum) und die Spalten der Matrix mit dem ersten Argument von $b_h$ assoziiert (Funktionen aus dem Testraum), siehe oben.
=> Wir erhalten also die Matrix der Einschränkung von $b_h$ bezüglich des ersten (bzw. zweiten) Argumentes auf $S_h^1 \cap H^1_0(\Omega)$, indem wir all jene Spalten (bzw. Zeilen) von $\underline{b_h}$ durch Einheitsspalten (bzw. Einheitszeilen) ersetzen, die zu einem Randknoten des Gitters gehören.
Zur Modellierung des Dirichlet-Randes bietet pyMOR sogenannten boundary infos an, die von BoundaryInfoInterface ableiten.
Diese erlauben es, zu einem gegeben Gitter diejenigen Randknoten zu erhalten, die mit dem entsprechenden Teil des Gebietsrandes assoziiert zind.
pymor.grids.boundaryinfos Moduls nach einer Möglichkeit, den gesamten Rand von \Omega als Dirichlet-Rand zu definieren und legen Sie ein entsprechendes Objekt boundary_info an.boundary_info beschreibt.boundary_info, um die Indices aller Randknoten des Gitters zu erhalten, und vergleichen Sie diese mit der Dokumentation des Gitters.Implementieren Sie die Interpolation $\Pi_{S_h^1\cap H^1_0(\Omega)}: C^0(\Omega) \to S_h^1\cap H^1_0(\Omega)$, die zusätzlich zur Lagrange-Interpolation noch die korrekte Randwertbehandlung vornimmt.
interpolate_lagrange_p1 Funktion, sodass eine optionale boundary_info übergeben werden kann, die den Dirichlet Rand beschreibt.Bemerkung: durch die Änderung der Funktion interpolate_lagrange_p1 in der interpolations.py Datei stimmt das bisher importierte Modul nicht mehr mit dieser Datei überein! Starten Sie zB den Kernel des Notebooks neu, um das Modul neu zu importieren.
Nutzen Sie den DiffusionOperatorP1, um das $H^1_0$ Produkt bezüglich $S_h^1$ zu assemblieren.
Für gegebenes $b$ und $l$ und deren Matrix und Vektordarstellung $\underline{b_h}$ und $\underline{l_h}$ ist folgendes äquivalent:
$u_h \in S_h^1 \cap H^1_0(\Omega)$ ist die Lösung des Variationsproblems
$$b(u_h, v_h) = l(v_h)\quad\quad\text{für alle } v_h \in S_h^1$$
$\underline{u_h} \in \mathbb{R}^I$ ist die Lösung des linearen Gleichungssystems
$$\underline{b_h}\; \underline{u_h} = \underline{l_h},$$
wobei
Was folgt dadurch für die Einträge von $\underline{u_h}$, die mit einem Dirichlet-Knoten assoziiert sind?
Wir betrachten das parameter-unabhängige Diffusionsproblem von oben mit
Legen Sie Datenfunktionen für $A$ und $f$ an, visualisieren Sie $f$.
Hinweis: nutzen Sie nicht die ExpressionFunction für $A$.
Nutzen Sie das L2ProductFunctionalP1 um die rechte Seite $l$ zu assemblieren, achten Sie auf korrekte Randwertbehandlung.
Von welchem Typ ist das assemblierte Funktional?
Lösen Sie das lineare Gleichungssystem, um $u_h$ zu erhalten und visualisieren Sie die Lösung.
Suchen Sie dazu in der Dokumentation nach einer Möglichkeit, um die Vektordarstellung von $l$ zu erhalten, und um den zu $b$ zugehörigen Operator zu invertieren.
Die exakte Lösung des Problems ist durch
$$u(x, y) = \frac{1}{2}\pi^2\; \cos(\frac{1}{2}\pi x)\; \cos(\frac{1}{2}\pi y)$$gegeben.