Praktikum zu Vorlesung Modellreduktion parametrisierter Systeme
Mario Ohlberger, Felix Schindler
Blatt 03, 24.04.2019¶
Aktivieren Sie wie gewohnt ihre Arbeitsumgebung und starten Sie den Jupyter Notebook server, siehe zB Blatt 1, Aufgabe 0.
Erstellen Sie ein neues
Python 3Notebook oder laden Sie dieses von der Homepage herunter.Importieren Sie
numpyundpymor.basicund machen Siematplotlibfür das Notebook nutzbar.
%matplotlib notebook
import numpy as np
from pymor.basic import *
Aufgabe 1: vereinfachtes Diffusionsproblem mit pyMOR diskretisieren¶
Wir betrachten wieder das vereinfachte Diffusionsproblem aus Blatt 02, mit dem Gebiet und den Datenfunktionen aus Blatt 02, Aufgabe 4.
Modellieren Sie das vereinfachte Diffusionsproblem mit Hilfe eines analytischen Problems in pyMOR und nutzen Sie einen Diskretizer, um eine Diskretisierung des Problems zu erhalten.
- Machen Sie sich mit der Klasse
StationaryProblemvertraut und legen Sie einproblem_1an, welches die gewünschten Datenfunktionen enthält.
#StationaryProblem?
problem_1 = StationaryProblem(
domain=RectDomain([[-1., -1.], [1., 1.]]),
diffusion=ConstantFunction(1, 2),
rhs=ExpressionFunction('0.5 * np.pi**2 * np.cos(0.5 * np.pi * x[..., 0]) * np.cos(0.5 * np.pi * x[..., 1])', 2, ()),
)
- Visualisieren Sie die Diffusion und rechte Seite des Problems.
- Legen Sie dazu ein feines
TriaGridan. - Nutzen Sie die Funktionen aus
interpolations.pyundvisualizations.py.
- Legen Sie dazu ein feines
from interpolations import (interpolate_fv, interpolate_lagrange_p1)
from visualizations import (visualize_fv, visualize_grid, visualize_lagrange_p1)
visualization_grid = TriaGrid(domain=problem_1.domain.domain, num_intervals=(128, 128))
visualize_lagrange_p1(interpolate_lagrange_p1(problem_1.diffusion, visualization_grid), visualization_grid, 'diffusion')
visualize_lagrange_p1(interpolate_lagrange_p1(problem_1.rhs, visualization_grid), visualization_grid, 'force')
- Diskretisieren Sie das stationäre Problem mit einer CG Diskretisierung, indem Sie die entsprechende Funktion aus dem pymor.discretizers Paket verwenden.
- Rufen Sie diese Funktion nur mit dem analytischen Problem und $\tfrac{1}{2}$ als maximalem Durchmesser für die Gitterelemente auf.
- Was ist der Rückgabewert des Discretizers?
#discretize_stationary_cg?
d, data = discretize_stationary_cg(problem_1, diameter=1/2)
print(data.keys())
- Lassen Sie sich von der resultierenden Diskretisierung den assemblierten Operator $B$, die assemblierte rechte Seite $l$ und den Approximationsraum $S_h^1$ anzeigen.
print(d.operator)
print(d.rhs)
print(d.solution_space)
- Lassen Sie sich außerdem alle assemblierten Produkte anzeigen.
for id, prod in d.products.items():
print('{}: {}'.format(id, prod))
- Nutzen Sie die Diskretisierung, um das vereinfachte Diffusionsproblem zu lösen.
- Was erhalten Sie als Lösung?
- Stellen Sie sicher, dass die Lösung im korrekten Vektorraum enthalten ist!
U = d.solve()
print(type(U))
print(len(U))
print(U.dim)
assert U.space == d.solution_space
- Visualisieren Sie die Lösung mit Hilfe der Diskretisierung.
d.visualize(U)
- Nutzen Sie denselben Diskretizer, um das Problem auf einem Vierecksgitter zu diskretisieren und zu lösen, visualisieren Sie die Lösung.
d, _ = discretize_stationary_cg(problem_1, diameter=1/2, grid_type=RectGrid)
d.visualize(d.solve())
Diskretisieren Sie dasselbe Problem mit Hilfe einer Finite Volumen Diskretisierung für verschiedene Gittertypen und dieselbe Gitterweite $\tfrac{1}{2}$, visualisieren Sie die Lösungen.
Falls Sie die Meldungen von pyMOR einschränken wollen, können Sie
set_log_levels({'pymor': 'WARN'})
verwenden.
Vergleichen Sie die Lösungen für CG und FV.
set_log_levels({'pymor': 'WARN'})
d, _ = discretize_stationary_fv(problem_1, diameter=1/2, grid_type=TriaGrid)
d.visualize(d.solve())
d, _ = discretize_stationary_fv(problem_1, diameter=1/2, grid_type=RectGrid)
d.visualize(d.solve())
Aufgabe 2: Neumann Ränder¶
Wir betrachten parameter-unabhängige Diffusionsprobleme. Sei dazu
- $\Omega \subset \mathbb{R}^d$, $d = 1, 2, 3$ eine beschränktes zusammenhängendes Gebiet mit
Lipschitz-Rand $\partial \Omega$, der sich in einen
- nichtleeren Dirichlet Rand $\Gamma_\text{D} \subseteq \partial \Omega$ und einen
- Neumann Rand $\Gamma_\text{N} := \partial \Omega \backslash \Gamma_\text{D}$
mit $\Gamma_\text{D} \cap \Gamma_\text{N} = \emptyset$ und $\Gamma_\text{D} \cup \Gamma_\text{N} = \partial \Omega$ aufteilt.
Sei außerdem
- eine Diffusion $A \in L^\infty(\Omega)$,
- eine rechte Seite $f \in L^2(\Omega)$,
- Dirichlet Randwerte $g_\text{D} \in H^{-1/2}(\Gamma_\text{D})$ und
- Neumann Randwerte $g_\text{N} \in H^{-1/2}(\Gamma_\text{N})$,
gegeben. Gesucht ist eine schwache Lösung $u \in H^1(\Omega)$, sodass
$$\begin{align} -\nabla\cdot( A \nabla u ) &= f &&\text{in } \Omega\\ u &= g_\text{D} &&\text{auf } \Gamma_\text{D}\\ - (A \nabla u) \cdot n &= g_\text{N} &&\text{auf } \Gamma_\text{N}\\ \end{align}$$im schwachen Sinne, wobei $n \in \mathbb{R}^d$ die aüßere Normale an den Gebietsrand $\partial \Omega$ sei.
Approximieren Sie die Lösung dieses Problems auf dem Gebiet
- $\Omega = [0, 1]^2$
- $\Gamma_\text{N} = [0, 1] \times \{0\}$
- $\Gamma_\text{D} = \partial\Omega \backslash \Gamma_\text{N}$
für die Datenfunktionen
- $A(x) := \begin{cases}0.001, &&\big|x - (\tfrac{1}{2}, \tfrac{1}{2})\big| < 0.3\\1, &&\text{sonst}\end{cases}$
- $f = 0$
- $g_\text{D} = 0$
- $g_\text{N} = -1$
- Legen Sie ein entsprechendes analytisches Problem
problem_2an. Beachten Sie dass die Definition der Randwerte schon im Problem erfolgen muss.
problem_2 = StationaryProblem(
domain=RectDomain(bottom='neumann'),
diffusion=ExpressionFunction('1. - (np.sqrt( (x[...,0]-0.5)**2 + (x[...,1]-0.5)**2) <= 0.3) * 0.999' , 2, ()),
neumann_data=ConstantFunction(-1, 2)
)
- Diskretisieren Sie das analytische Problem mit dem CG Verfahren mit einer Gitterweite von höchstens $\tfrac{1}{32}$ und visualisieren Sie die Lösung.
d, data = discretize_stationary_cg(problem_2, diameter=1/32)
d.visualize(d.solve())
- Visualisieren Sie $A$ mit Hilfe der Diskretisierung.
d.visualize(interpolate_lagrange_p1(problem_2.diffusion, data['grid']))
- Lassen Sie sich die $L^\infty$-Norm des Vektors des assemblierten Rechte-Seite-Funktionals $l$ der Diskretisierung anzeigen.
print(d.rhs.as_vector().sup_norm())
Extrahieren Sie aus den bei der Diskretisierung anfallenden Daten das Rechte-Seite Funktional $l$ der nicht assemblierten Diskretisierung.
- Was ist $l$?
- Lassen Sie sich die Koeffizienten der Linearkombination von $l$ anzeigen und extrahieren Sie das zu $f$ und das zu $g_\text{N}$ assoziierte Funktional.
- Assemblieren Sie diese und lassen Sie sich die jeweilige $L^\infty$-Norm des assemblierten Vektors anzeigen.
ud = data['unassembled_d']
print(type(ud.rhs))
print(ud.rhs.coefficients)
print(ud.rhs.operators)
l_f = ud.rhs.operators[0].assemble()
l_gN = ud.rhs.operators[1].assemble()
print(l_f.as_vector().sup_norm())
print(l_gN.as_vector().sup_norm())
Aufgabe 3: parametrische Neumannwerte¶
Wir betrachten das Problem aus Aufgabe 2 mit parametrischen Neumann Randwerten $g_{\text{N}, \mu}$, d.h. für jeden Parameter $\mu \in \mathbb{R}$ ist $g_{\text{N}, \mu} \in H^{-1/2}(\Gamma_\text{N})$ und wir suchen zu einem Parameter $\mu \in \mathbb{R}$ eine schwache Lösung $u_\mu \in H^1(\Omega)$, sodass
$$\begin{align} -\nabla\cdot( A \nabla u_\mu ) &= f &&\text{in } \Omega\\ u_\mu &= g_\text{D} &&\text{auf } \Gamma_\text{D}\\ - (A \nabla u_\mu) \cdot n &= g_{\text{N}, \mu} &&\text{auf } \Gamma_\text{N}\\ \end{align}$$Approximieren Sie die Lösungen dieses Problems für verschiedene Parameter auf dem Gebiet wie oben und zu den Datenfunktionen wie oben, aber mit
- $g_{\text{N}, \mu}(x, y) := -\cos(\pi x_0)^2 \mu$
Legen Sie eine
ExpressionFunctionfür $g_{\text{N}, \mu}$ an- Machen Sie sich mit dem ersten Absatz der pyMOR Dokumentation zu parametrischen Objekten vertraut, was ist ein Parameter (zur Modellierung von $\mu$) in pyMOR, was ein Parametertyp?
- Lesen Sie die Dokumentation zur
ExpressionFunction, um die Neumann Randwerte als parametrische Funktion mit einem skalaren Parameter zu definieren (nennen Sie diesen z.B.'neum'). - Lassen Sie sich die Repräsentation der resultierende Funktion, ihren Paramtertypen sowie dessen Python-Typen ausgeben.
#ExpressionFunction?
g_N = ExpressionFunction('-cos(pi*x[...,0])**2 * neum', 2, (), parameter_type={'neum': ()})
print(g_N)
g_N.__repr__()
print(g_N.parameter_type)
print(type(g_N.parameter_type))
Werten Sie die Funktion für verschiedene Punkte auf dem Neumann Rand und verschiedene Parameter aus.
- Legen Sie dazu einen
Parameter-Objectmuan, welches $\mu = 1$ modelliert, lassen Sie sich den Paramter anzeigen und werten Sie die Funktion für dieses $\mu$ aus. - Finden Sie einfachere Wege, die Funktion für diesen und andere Parameterwerte auszuwerten.
- Legen Sie dazu einen
mu = Parameter({'neum': 1})
print(mu)
bottom_points = np.array([[0, 0], [0.5, 0], [1, 0]])
print(g_N.evaluate(bottom_points, mu))
print(g_N.evaluate(bottom_points, {'neum': 1}))
print(g_N.evaluate(bottom_points, {'neum': -2}))
print(g_N.evaluate(bottom_points, -2))
- Versuchen Sie mit dem parametrischen $g_{\text{N}, \mu}$ die Neumann Randwerte des
problem_2Objektes zu ersetzen, um ein neues analytisches Problem zu erhalten. Warum schlägt dies fehl? (Hinweis: nutzen Sie diemroFunktion aus deminspectModul, um sich alle Basisklassen des Typs vonproblem_2anzeigen zu lassen)
problem_3 = problem_2
#problem_3.neumann_data = g_N # <= funktioniert nicht
import inspect
inspect.getmro(type(problem_2))
- Machen Sie sich mit den entsprechenden Konzepten in pyMOR vertraut und legen Sie ein
problem_3als veränderte Kopie desproblem_2an, wobei nur die Neumann Randwerte ausgetauscht wurde.
problem_3 = problem_2.with_(neumann_data=g_N)
- Diskretisieren Sie
problem_3mit dem CG Verfahren mit einer Gitterweite von höchstens $\tfrac{1}{32}$. Ist die Diskretisierung parametrisch? Wenn ja, lassen Sie sich den entsprechenden Parametertypen anzeigen.
d, data = discretize_stationary_cg(problem_3, diameter=1/32)
print(d.parametric)
print(d.parameter_type)
- Lassen Sie sich den Linke-Seite-Operator und das Rechte-Seite-Funktional der Diskretisierung anzeigen. Sind diese parametrisch? Lassen Sie sich die entsprechenden Parametertypen anzeigen.
lhs = d.operator
print(lhs)
print(lhs.parametric)
print(lhs.parameter_type)
rhs = d.rhs
print(rhs)
print(rhs.parametric)
print(rhs.parameter_type)
- Ist das Rechte-Seite-Funktional assembliert? Lassen Sie sich gegebenenfalls den Paramtertypen aller Funktionale der Linearkombination anzeigen.
print(type(rhs)) # <= nicht assembliert
for func in rhs.operators:
print(func)
- Berechnen Sie die Lösung für verschiedene Paramter $\mu$ (analog zur Auswertung von $g_{\text{N}, \mu}$) und visualisieren Sie diese.
for mu in -10, 1, 2:
d.visualize(d.solve(mu))
- Lassen Sie sich die Summe der Lösungen $u_{\mu_1} + u_{\mu_2}$, für die Parameter $\mu_1 = 1$ und $\mu_2 = -1$ anzeigen und verifizieren Sie das Ergebnis.
U_1 = d.solve(mu={'neum': 1})
U_2 = d.solve(mu={'neum': -1})
d.visualize(U_1 + U_2)
Da $g_{\text{N}, \mu_1} = -g_{\text{N}, \mu_2}$ und $f = 0$ gilt $u_{\mu_1} = - u_{\mu_2}$, also $u_{\mu_1} + u_{\mu_2} = 0$.