Praktikum zur Vorlesung Numerik partieller Differentialgleichungen I
Mario Ohlberger, Felix Schindler
Blatt 00, 16.10.2019¶
Aufgabe 0¶
(i) einmalige Einrichtung¶
Dazu müssen auf Ihrem Rechner die Befehle python3 und virtualenv verfügbar sein. Das können Sie überprüfen, indem Sie ein Terminal starten (unter Ubuntu: Super/Windows Taste, terminal eingeben, Return) und die folgenden Befehle eingeben:
python3 --version
virtualenv --version
Beide Befehle sollten Ihnen die jeweilige Version des Programms anzeigen (und keine Fehlermeldung). Zur Einrichtung gehen Sie folgendermaßen vor:
- Ins
HOMEVerzeichnis wechseln und einen Projekt-Ordner anlegen
cd # ins HOME Verzeichnis wechseln
pwd # Ort des aktuellen Verzeichnises anzeigen
ls -lh # Inhalt des aktuellen Verzeichnisses anzeigen
mkdir NPDGL1 # ein neues Verzeichnis erstellen
cd NPDGL1 # in dieses Verzeichnis wechseln
pwd
ls -lh
- eine virtuelle Umgebung für die Programmierung mit Python schaffen
virtualenv --python=python3 python_umgebung_NPDGL1
ls -lh
- die virtuelle Umgebung aktivieren
source python_umgebung_NPDGL1/bin/activate
- benötigte Software in der virtuellen Umgebung installieren
pip install --upgrade pip
pip install pymor[full]
pip install notebook
du -sh python_umgebung_NPDGL1
Dabei installiert pip install pymor[full] das neueste Release von pyMOR (momentan 0.5.2). Mehr dazu finden Sie auf https://pymor.org (gehen Sie dort auf View on Github und wählen anstelle von branch: master links oben Tags: 0.5.2 aus, um die entsprechende README.md anzuzeigen).
- Arbeitsverzeichnis erstellen
mkdir notebooks
(ii) ab jetzt jedes mal¶
- ins Projektverzeichnis wechseln
cd ~/NPDGL1
- die virtuelle Umgebung aktivieren
source python_umgebung_NPDGL1/bin/activate
- Notebook server starten
jupyter-notebook --notebook-dir=notebooks
(iii) Blatt 00 nutzbar machen¶
Laden Sie sich von der Homepage der Vorlesung blatt_00.ipynb herunter und speichern Sie dieses im Verzeichnis
NPDGL1/notebooks/alsblatt_00.ipynb.Starten Sie den Notebook server (siehe ii) und führen Sie
blatt_00.ipynbaus.
Ab jetzt können Sie direkt im jeweiligen Übungsblatt arbeiten, fügen Sie dazu einfach direkt nach der Aufgabenstellung mit + (oben links) eine neue Code Zelle hinzu. Eine Code Zelle können Sie mit STRG + Return ausführen. Hilfe zu einem Typ Foo oder einer Objekt foo können Sie mit Foo? oder foo? erhalten.
# führen Sie in einem Notebook immer die folgenden Befehle als erstes aus
%matplotlib notebook
import numpy as np
from pymor.basic import *
Ab jetzt ist die numpy Bibliothek unter np verfügbar und alles aus pymor.basic direkt verfügbar.
Definition: vereinfachtes Diffusionsproblem¶
Wir betrachten ein vereinfachtes stationäres Diffusionsproblem. Sei dazu
- ein beschränktes polygonales Gebiet $\Omega \subset \mathbb{R}^2$ mit Lipschitz-Rand $\partial \Omega$,
- eine Diffusion $A: \Omega \to \mathbb{R}$ und
- eine rechte Seite $f: \Omega \to \mathbb{R}$
gegeben. Gesucht ist die Lösung $u: \Omega \to \mathbb{R}$, sodass
$$\begin{align} -\nabla\cdot( A \nabla u ) &= f &&\text{in } \Omega\\ u &= 0 &&\text{auf } \partial \Omega, \end{align}$$(im schwachen Sinne, siehe Kapitel 1) erfüllt ist.
Aufgabe 1: Gebiete und Gitter¶
- Suchen Sie in der pyMOR Dokumentation auf
docs.pymor.org/nach einer Klasse, die das Gebiet $$\Omega = [0, 1]^2$$ mit Dirichlet-Randwerten modellieren kann und legen Sie ein entsprechendes Objekt mit Namenomegaan. Nutzen Sie dazu auch die Möglichkeit des jupyter Notebooks, sich zu einem beliebigen Objekt mit Hilfe des?-postfix die Dokumentation anzeigen zu lassen (z.B.np.array?).
# => API Dokumentation http://docs.pymor.org/en/0.5.2/index.html#api-documentation
# => Modul "domaindescriptions": http://docs.pymor.org/en/0.5.2/generated/pymor.domaindescriptions.html
# => Klasse RectDomain, diese ist Dank des * imports direkt verfügbar
# => Hilfe zu RectDomain anzeigen lassen: RectDomain?
# => enthält schon Dirichlet Randwerttyp überall
omega = RectDomain(([0, 0], [1, 1]))
Suchen Sie nach einer Möglichkeit ein konformes Dreiecksgitter $\mathcal{T}_h$ als Partitionierung des Gebiets $\Omega$ mit 16 Dreiecken zu erstellen, und geben Sie folgende Informationen aus:
- Dimension
ddes Gitters - Anzahl der Elemente (Kodimension 0 Entitäten, in diesem Fall Dreiecke) im Gitter
- Anzahl der Kanten (Kodimension 1 Entitäten, in diesem Fall Intervalle) im Gitter
- Anzahl der Knoten (Kodimension
dEntitäten, Punkte) im Gitter - die Position eines jeden Knoten des Gitters mit Hilfe einer
for-Schleife
Vergleichen Sie die Positionen der Knoten mit der Dokumentation.
- Dimension
# => API Dokumentation, Modul grid, Klasse TriaGrid
# => Hilfe anzeigen lassen: TriaGrid?
grid = TriaGrid(num_intervals=(2, 2), domain=omega.domain)
print(grid)
print()
print('grid lives in a {}-dimensional space'.format(grid.dim))
print('grid has {} simplices'.format(grid.size(0)))
print('grid has {} edges'.format(grid.size(1)))
print('grid has {} vertices'.format(grid.size(2)))
vertices = grid.centers(grid.dim)
print()
print(type(vertices))
print(vertices.shape)
print()
print('vertices:')
counter = 0
for vertex in vertices:
print('vertex {} has position {}'.format(counter, vertex))
counter += 1
Aufgabe 2: analytische Funktionen vektorisiert implementieren¶
Sei $f: \Omega \to \mathbb{R}$ gegeben durch $f(x_0, x_1) := x_0 \cdot x_1$.
- Schreiben Sie eine Funktion
f, die die Auswertung von $f$ and einem Punkt $x$ modelliert und werten Sie die Funktionffür jeden Knotenpunkt des Gitters einzeln mit Hilfe einerfor-Schleife aus.
def f(x):
return x[0]*x[1]
for x in vertices:
print('f({}) = {}'.format(x, f(x)))
- Erstellen Sie ein feineres Gitter mit $1024^2$ Intervallen und messen Sie die benötigte Zeit, um
fan jedem Knotenpunkt des feineren Gitters auszuwerten.
fine_grid = TriaGrid(num_intervals=(1024, 1024), domain=omega.domain)
print(fine_grid)
import time
tic = time.time()
values = f(fine_grid.centers(grid.dim))
toc = time.time() - tic
print()
print('evaluating non-vectorized function took {}s'.format(toc))
- Versuchen Sie,
fan allen Knotenpunkten des groben Gitters gleichzeitig auszuwerten. Warum funktioniert das nicht?
# => an allen Punkten gleichzeitig auswerten
values = f(vertices)
print(values)
print(len(vertices))
print(len(values))
# Warum kommt das raus?
# Beim Aufruf f(vertices) passiert
vertices[0]*vertices[1]
# also eine elementweise Multiplikation des 0ten und 1ten Knotens
- Schreiben Sie eine vektorisierte Variante von
f, welche die elementweisen Operationen desnumpy.ndarrayausnutzt und werten Sie diese Variante für alle Knotenpunkte des feineren Gitters gleichzeitig aus und messen Sie die benötigte Zeit.
# Gegeben eine Menge von Punkten als Matrix
print(vertices)
print(vertices.shape)
# extrahiere 0te und 1te Komponente für alle Punkte gleichzeitig
x_0s = vertices[..., 0] # 0te Spalte
x_1s = vertices[..., 1] # 1te Spalte
print()
print(type(x_0s))
print(x_0s)
print(x_1s)
# Das ganze kompakt in einer Funktion
def f(x):
return x[..., 0]*x[..., 1]
# Kann an einem Punkt ausgewertet werden ...
print(f(vertices[0]))
# ... oder an allen Punkten gleichzeitig
print(f(vertices))
tic = time.time()
values = f(fine_grid.centers(grid.dim))
toc = time.time() - tic
print('took {}s'.format(toc))
- Erstellen Sie eine
ExpressionFunction, um $f$ zu modellieren, werten Sie diese an allen Knotenpunkten des feinen Gitters aus und messen Sie die benötigte Zeit.
f = ExpressionFunction('x[..., 0]*x[..., 1]', dim_domain=2, shape_range=())
tic = time.time()
values = f.evaluate(fine_grid.centers(grid.dim))
toc = time.time() - tic
print('took {}s'.format(toc))
Aufgabe 3: Diffusionsproblem in pyMOR modellieren¶
Modellieren Sie das vereinfachte Diffusionsproblem mit
- $A = 1$
- $f = 1$
mit Hilfe eines analytischen Problems in pyMOR.
- Suchen Sie nach einer geeigneten Klasse, um $A$ und $f$ als Funktionen in pyMOR zu modellieren.
#ConstantFunction?
A = ConstantFunction(1, 2)
f = ConstantFunction(value=1, dim_domain=2)
- Machen Sie sich mit der Klasse
StationaryProblemvertraut und legen Sie einprobleman, welches die gewünschten Datenfunktionen enthält.
#StationaryProblem?
problem = StationaryProblem(
domain=omega,
diffusion=ConstantFunction(1, 2),
rhs=ConstantFunction(1, 2),
)
Aufgabe 4: Diffusionsproblem in pyMOR diskretisieren¶
- Diskretisieren Sie das stationäre Problem mit einer CG Diskretisierung (CG = continuous Galerkin = stetige Finite Elemente, siehe Kapitel 2), indem Sie die entsprechende Funktion aus dem pymor.discretizers Paket verwenden.
- Rufen Sie diese Funktion mit dem analytischen Problem auf.
- Was ist der Rückgabewert?
result = discretize_stationary_cg(problem)
type(result)
print()
for rr in result:
print(rr)
print()
disc, _ = discretize_stationary_cg(problem)
- Berechnen Sie eine approximative Lösung
u_hdes vereinfachten Diffusionsproblems mit Hilfe der Diskretisierung. Machen Sie sich dazu mit der Dokumentation vonStationaryDiscretizationvertraut oder finden Sie die entsprechende Funktion interaktiv (um die Methoden eines Objektesfooanzuzeigen, schreiben Siefoo.und drücken SieTab).
u_h = disc.solve()
- Visualisieren Sie die approximative Lösung mit Hilfe der Diskretisierung
disc.visualize(u_h)
Aufgabe 5: Approximationsgüte¶
Berechnen Sie Approximationen u_h für verschieden Gitterweiten h = 1/2, 1/4, 1/16, 1/64, indem Sie das diameter Argument der Funktion discretize_stationary_cg nutzen.
- Berechnen und visualisieren Sie diese Approximationen und vergleichen Sie sie.
- Unterdrücken Sie die automatische Ausgabe des Diskretizers mit Hilfe von
set_log_levels({'pymor': 'WARN'}) - Geben Sie den jeweiligen diskreten Lösungsraum aus.
set_log_levels({'pymor': 'WARN'})
def discretize_and_visualize(h):
disc, _ = discretize_stationary_cg(problem, diameter=h)
print(disc.solution_space)
u_h = disc.solve()
disc.visualize(u_h)
for h in (1/2, 1/4, 1/16, 1/64):
print('h = {}'.format(h))
discretize_and_visualize(h)
print()
