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.
%matplotlib notebook
import numpy as np
from pymor.basic import *
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.
from interpolations import (interpolate_fv, interpolate_lagrange_p1)
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.
from visualizations import (visualize_fv, visualize_grid, visualize_lagrange_p1)
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
.omega = RectDomain(([0, 0], [1, 1]))
grid = TriaGrid(domain=omega.domain, num_intervals=(2, 2))
visualize_grid(grid)
CGVectorSpace
aus dem pymor.operators.cg
Modul und vergleichen Sie dessen Dimension mit der Anzahl der Knoten des Gitters.from pymor.operators.cg import CGVectorSpace
cg_space = CGVectorSpace(grid)
print('S_h^1 = {}'.format(cg_space))
print('number of grid nodes: {}'.format(grid.size(grid.dim)))
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?
from pymor.operators.cg import DiffusionOperatorP1
# => Hilfe anzeigen mit DiffusionOperatorP1?
# => Dieser Operator beinhaltet Randwertbehandlung:
# Für den Leplace Operator bezüglich H^1 benötigen wir also eine boundary_info,
# die keinen Dirichlet-Rand definiert.
h1_semi_product = DiffusionOperatorP1(grid, boundary_info=EmptyBoundaryInfo(grid))
print('h1_semi_product is a {}'.format(h1_semi_product))
print('the type of h1_semi_product is {}'.format(type(h1_semi_product)))
if h1_semi_product.linear:
print('h1_semi_product is linear')
cg_space
übereinstimmen.print(h1_semi_product.source)
assert h1_semi_product.source == cg_space
print(h1_semi_product.range)
assert h1_semi_product.range == cg_space
from interpolations import interpolate_lagrange_p1
from visualizations import visualize_lagrange_p1
f = ExpressionFunction('x[..., 0]', dim_domain=2, shape_range=())
f_h = interpolate_lagrange_p1(f, grid)
print(np.sqrt(h1_semi_product.apply2(f_h, f_h)))
visualize_lagrange_p1(f_h, grid, 'f(x,y) = x')
Es gilt $\nabla f (x, y) = (1, 0)^\top$, also
$$|f_h|_{H^1(\Omega)}^2 = \int_\Omega \nabla f_h \cdot \nabla f_h \text{d}x = \int_\Omega (1, 0)^\top \cdot (1, 0)^\top\text{d}x = \int_\Omega \text{d}x = 1$$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.
h1_semi_product = h1_semi_product.assemble()
print('h1_semi_product is a {}'.format(h1_semi_product))
print('the type of h1_semi_product is {}'.format(type(h1_semi_product)))
# NumpyMatrixOperator?
h1_semi_mat = h1_semi_product.matrix
print('the matrix representation of h1_semi_product is a {}'.format(type(h1_semi_mat)))
print('the matrix representation of h1_semi_product is \n{}'.format(h1_semi_mat))
# https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csc_matrix.html
h1_semi_mat = h1_semi_mat.todense()
print(type(h1_semi_mat))
print('the matrix representation of h1_semi_product is \n{}'.format(h1_semi_mat))
Berechnen Sie $b(f_h, f_h)$, indem Sie die Matrixdarstellung nutzen.
f_h
anzeigen lassen, Dokumentation nutzen).# - Anwendung des Produkts über die Matrixdarstellung/DoF Vektor
print(f_h.data.dot(h1_semi_mat.dot(f_h.data.T)))
# - moderne Syntax
print(f_h.data @ (h1_semi_mat @ f_h.data.T))
# - Anwedung des Produkts über die Bilinearform
print(h1_semi_product.apply2(f_h, f_h))
# Was sollte ein Produkt erfüllen?
# * Symmetrie
if np.all(h1_semi_mat == h1_semi_mat.T):
print('=> is symmetric')
else:
print('=> is NOT symmetric!')
# * Positive definitheit:
# - zB alle Eigenwerte > 0
min_ev = np.linalg.eigvals(h1_semi_mat).min()
if np.allclose(min_ev, 0.) or min_ev < 0.:
print('=> the minimum eigenvalue {} is 0 or negative!'.format(min_ev))
# - zB Determinante != 0
det = np.linalg.det(h1_semi_mat)
if np.allclose(det, 0.):
print('=> det is 0!')
# * Wenn b_h(v_h, v_h) != 0 so muss schon v_h = 0 gelten
one = h1_semi_product.source.from_data(np.ones(cg_space.dim))
print('(1, 1)_H^1 = {}'.format(h1_semi_product.apply2(one, one)))
# => h1_semi_product ist kein Produkt bezüglich H^1 (wie erwartet)
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.boundary_info = AllDirichletBoundaryInfo(grid)
print(boundary_info.boundary_types)
vertices = grid.centers(grid.dim)
print(len(vertices))
boundary_indices = np.where(boundary_info.mask('dirichlet', codim=2))[0]
print(boundary_indices)
boundary_vertices = vertices[boundary_indices]
print(len(boundary_vertices))
print(boundary_vertices)
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.
f_in_H1 = interpolate_lagrange_p1(f, grid)
visualize_lagrange_p1(f_in_H1, grid, 'f in H^1')
f_in_H10 = interpolate_lagrange_p1(f, grid, boundary_info)
visualize_lagrange_p1(f_in_H10, grid, 'f in H^1_0')
Nutzen Sie den DiffusionOperatorP1
, um das $H^1_0$ Produkt bezüglich $S_h^1$ zu assemblieren.
h10_product = DiffusionOperatorP1(grid, boundary_info=boundary_info, dirichlet_clear_columns=True)
h10_product = h10_product.assemble()
h10_mat = h10_product.matrix.todense()
print('the matrix representation of h10_product is \n{}'.format(h10_mat))
print(np.allclose(h10_mat[boundary_indices, :], np.eye(h10_mat.shape[0])[boundary_indices, :]))
print(np.allclose(h10_mat[:, boundary_indices], np.eye(h10_mat.shape[0])[:, boundary_indices]))
# Was sollte ein Produkt erfüllen?
# * Symmetrie
if np.all(h10_mat == h10_mat.T):
print('=> is symmetric')
else:
print('=> is NOT symmetric!')
# * Positive definitheit:
# - zB alle Eigenwerte > 0
min_ev = np.linalg.eigvals(h10_mat).min()
if np.allclose(min_ev, 0.) or min_ev < 0.:
print('=> the minimum eigenvalue {} is 0 or negative!'.format(min_ev))
else:
print('=> eigenvalues are positive')
# - zB Determinante != 0
det = np.linalg.det(h10_mat)
if np.allclose(det, 0.):
print('=> det is 0!')
else:
print('=> seems invertible')
# * Wenn b_h(v_h, v_h) != 0 so muss schon v_h = 0 gelten
one = interpolate_lagrange_p1(ConstantFunction(1, 2), grid, boundary_info)
print('(1, 1)_H^1 = {}'.format(h10_product.apply2(one, one)))
# => h1_semi_product ist kein Produkt bezüglich H^1 (wie erwartet)
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
omega = RectDomain(([-1, -1], [1, 1]))
grid = TriaGrid(domain=omega.domain, num_intervals=(8, 8))
boundary_info = AllDirichletBoundaryInfo(grid)
visualize_grid(grid)
Legen Sie Datenfunktionen für $A$ und $f$ an, visualisieren Sie $f$.
Hinweis: nutzen Sie nicht die ExpressionFunction
für $A$.
diffusion = ConstantFunction(1, dim_domain=omega.dim)
force = ExpressionFunction('0.5 * np.pi**2 * np.cos(0.5 * np.pi * x[..., 0]) * np.cos(0.5 * np.pi * x[..., 1])', 2, ())
visualize_lagrange_p1(interpolate_lagrange_p1(force, grid), grid, 'f')
b_h = DiffusionOperatorP1(grid, boundary_info=boundary_info, diffusion_function=diffusion)
b_h = b_h.assemble()
Nutzen Sie das L2ProductFunctionalP1
um die rechte Seite $l$ zu assemblieren, achten Sie auf korrekte Randwertbehandlung.
Von welchem Typ ist das assemblierte Funktional?
from pymor.operators.cg import L2ProductFunctionalP1
l_h = L2ProductFunctionalP1(grid, force, dirichlet_clear_dofs=True, boundary_info=boundary_info)
print('l_h is a {}'.format(l_h))
print('the type of l_h is {}'.format(type(l_h)))
l_h = l_h.assemble()
print('l_h is a {}'.format(l_h))
print('the type of l_h is {}'.format(type(l_h)))
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.
u_h = b_h.apply_inverse(l_h.as_vector())
visualize_lagrange_p1(u_h, grid, 'u_h')
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.
u = interpolate_lagrange_p1(ExpressionFunction('np.cos(0.5 * np.pi * x[..., 0]) * np.cos(0.5 * np.pi * x[..., 1])', 2, ()), grid)
visualize_lagrange_p1(u, grid, 'u')
h10_product = DiffusionOperatorP1(grid, boundary_info=boundary_info, dirichlet_clear_columns=True)
h10_product = h10_product.assemble()
print('approx. relative H^1_0 error: {}'.format(np.sqrt(h10_product.apply2(u - u_h, u - u_h)[0] / np.sqrt(h10_product.apply2(u, u)[0]))))
visualize_lagrange_p1(u - u_h, grid, 'u - u_h')
Es gilt
$$\|u - u_h\| \leq \|u - \Pi_{S_h^1}[u]\| + \|\Pi_{S_h^1}[u] - u_h\|$$und wir haben nur den letzten Term berechnet.