Praktikum zur Vorlesung Numerik partieller Differentialgleichungen I

Mario Ohlberger, Felix Schindler

Blatt 01, 23.10.2019

Aufgabe 0: jupyter Notebook Server starten (Erinnerung)

Aktivieren Sie die virtuelle Umgebung zum Praktikum und starten Sie den Notebook server. Zur Erinnerung:

  • Terminal starten (super-Taste/Windows-Taste, terminal eingeben, Enter)

  • ins Projektverzeichnis wechseln

    cd ~/NPDGL1

  • die virtuelle Umgebung aktivieren

    source python_umgebung_NPDGL1/bin/activate

  • Notebook server starten

    jupyter-notebook --notebook-dir=notebooks
  • Laden Sie dieses Blatt von der Homepage herunter und speichern Sie es unter ~/NPDGL1/notebooks/blatt_01.ipynb.
  • Importieren Sie numpy und pymor.basic und machen Sie matplotlib für das Notebook nutzbar.
  • unterdrücken Sie die automatische Ausgabe des Diskretizers mit Hilfe von set_log_levels({'pymor': 'WARN'})
In [1]:
%matplotlib notebook
import numpy as np
from pymor.basic import *
set_log_levels({'pymor': 'WARN'})

Definition: stationäres Diffusionsproblem mit verallgemeinerten Randwerten

Wir betrachten ein stationäres Diffusionsproblem. Sei dazu

  • ein beschränktes polygonales Gebiet $\Omega \subset \mathbb{R}^d$ für $d = 1, 2, 3$ mit äußerer Normale $n \in \mathbb{R}^d$ und Lipschitz-Rand $\partial \Omega$, aufgeteilt in
  • einen nicht-leeren Dirichlet-Rand $\Gamma_\text{D} \subset \partial \Omega$, $\Gamma_\text{D} \neq \emptyset$ und
  • einen Neumann-Rand $\Gamma_\text{N} := \partial \Omega \backslash \Gamma_\text{D}$,
  • eine Diffusion $A: \Omega \to \mathbb{R}$,
  • eine rechte Seite (Quellen und Senken) $f: \Omega \to \mathbb{R}$,
  • Dirichlet-Randwerte $g_\text{D}: \Gamma_\text{D} \to \mathbb{R}$ und
  • Neumann-Randwerte $g_\text{N}: \Gamma_\text{N} \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 &= 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, siehe Kapitel 1) erfüllt ist.

Aufgabe 1: homogene Dirichlet-Randwerte (Wiederholung)

Wir betrachten obiges Problem mit $d = 2$, $\Omega = [0, 1]^2$, $\Gamma_\text{D} = \partial\Omega$, $\Gamma_\text{N} = \emptyset$, $A = f = 1$ und $g_\text{D} = 0$. Modellieren, diskretisieren und lösen Sie dieses Problem mit Hilfe von pyMOR und visualisieren Sie die Lösung. Nehmen Sie dafür die Musterlösung zu Blatt 00 und die API Dokumentation zu Hilfe. Zur Erinnerung:

  • modellieren Sie $\Omega$ als RectDomain
  • legen Sie ein StationaryProblem mit den entsprechenden Datenfunktionen an
  • diskretisieren Sie mit discretize_stationary_cg und speichern Sie die zurückgegebene Diskretisierung als disc ab
  • lösen Sie mit Hilfe der Diskretisierung
  • visualisieren Sie mit Hilfe der Diskretisierung
In [2]:
omega = RectDomain(([0, 0], [1, 1]))
problem = StationaryProblem(
    domain=omega,
    diffusion=ConstantFunction(1, 2),
    rhs=ConstantFunction(1, 2),
)
disc, _ = discretize_stationary_cg(problem)
u_h = disc.solve()
disc.visualize(u_h)

Aufgabe 2: Quellen und senken

(2.i) eine Quelle

Wir betrachten obiges Problem mit $d = 2$, $\Omega = [0, 1]^2$, $\Gamma_\text{D} = \partial\Omega$, $\Gamma_\text{N} = \emptyset$, $A = 1$, $g_\text{D} = 0$ und einer Quelle, modelliert durch

$$f(x) = \begin{cases}1, &x \in [0.2, 0.4]^2\\0, &\text{sonst}\end{cases}.$$
  • Modellieren, diskretisieren und lösen Sie dieses Problem mit Hilfe von pyMOR und visualisieren Sie die Lösung.
  • Vergleichen Sie das Maximum von $u$ mit dem aus Aufgabe 1.
In [3]:
problem = StationaryProblem(
    domain=omega,
    diffusion=ConstantFunction(1, 2),
    rhs=ExpressionFunction(
        '(0.2 <= x[..., 0]) * (x[..., 0] <= 0.4) * (0.2 <= x[..., 1]) * (x[..., 1] <= 0.4) * 1.', dim_domain=2, shape_range=()),
)
disc, _ = discretize_stationary_cg(problem)
u_h = disc.solve()
disc.visualize(u_h)

(2.ii) eine Quelle und eine Senke

Wir betrachten obiges Problem mit $d = 2$, $\Omega = [0, 1]^2$, $\Gamma_\text{D} = \partial\Omega$, $\Gamma_\text{N} = \emptyset$, $A = 1$, $g_\text{D} = 0$ und je einer Quelle und Senke, modelliert durch

$$f(x) = \begin{cases}1, &x \in [0.2, 0.4]^2\\-1, &x \in [0.6, 0.8]^2\\0, &\text{sonst}\end{cases}.$$
  • Modellieren, diskretisieren und lösen Sie dieses Problem mit Hilfe von pyMOR und visualisieren Sie die Lösung.
  • Betrachten Sie den Fluss $- A \nabla u$ entlang der Diagonalen $x_0 = x_1$.
In [4]:
problem = StationaryProblem(
    domain=omega,
    diffusion=ConstantFunction(1, 2),
    rhs=(ExpressionFunction('((0.2 <= x[..., 0]) * (x[..., 0] <= 0.4) * (0.2 <= x[..., 1]) * (x[..., 1] <= 0.4))*1.', dim_domain=2, shape_range=())
         - ExpressionFunction('((0.6 <= x[..., 0]) * (x[..., 0] <= 0.8) * (0.6 <= x[..., 1]) * (x[..., 1] <= 0.8))*1.', dim_domain=2, shape_range=())),
)
disc, _ = discretize_stationary_cg(problem)
u_h = disc.solve()
disc.visualize(u_h)

Aufgabe 3: Dirichlet-Randwerte

(3.i)

Wir betrachten obiges Problem mit $d = 2$, $\Omega = [0, 1]^2$, $\Gamma_\text{D} = \partial\Omega$, $\Gamma_\text{N} = \emptyset$, $A = f = 1$ und $g_\text{D} = 1$.

  • Modellieren, diskretisieren und lösen Sie dieses Problem mit Hilfe von pyMOR und visualisieren Sie die Lösung.
  • Vergleichen Sie $u$ mit der Lösung aus Aufgabe 1.
In [5]:
problem = StationaryProblem(
    domain=omega,
    diffusion=ConstantFunction(1, 2),
    rhs=ConstantFunction(1, 2),
    dirichlet_data=ConstantFunction(1, 2)
)
disc, _ = discretize_stationary_cg(problem)
u_h = disc.solve()
disc.visualize(u_h)

(3.ii)

Wir betrachten obiges Problem mit $d = 2$, $\Omega = [0, 1]^2$, $\Gamma_\text{D} = \partial\Omega$, $\Gamma_\text{N} = \emptyset$, $A = f = 1$ und

$$g_\text{D}(x) = \begin{cases}1, &x_0 = 0\\-1, &x_0 = 1\\0, &\text{sonst}\end{cases}.$$

Modellieren, diskretisieren und lösen Sie dieses Problem mit Hilfe von pyMOR und visualisieren Sie die Lösung.

In [6]:
problem = StationaryProblem(
    domain=omega,
    diffusion=ConstantFunction(1, 2),
    rhs=ConstantFunction(1, 2),
    dirichlet_data=ExpressionFunction('(x[..., 0] < 1e-7)*1. + (x[..., 0] > (1 - 1e-7))*-1.', dim_domain=2, shape_range=())
)
disc, _ = discretize_stationary_cg(problem)
u_h = disc.solve()
disc.visualize(u_h)

Aufgabe 4: Neumann-Randwerte

(4.i)

Wir betrachten obiges Problem mit $d = 2$, $\Omega = [0, 1]^2$, $\Gamma_\text{N} = \{x \in \overline{\Omega} \;|\; x_0 = 1\}$, $\Gamma_\text{D} = \partial\Omega \backslash \Gamma_{N}$, $A = f = 1$ und $g_N = 1$. Modellieren, diskretisieren und lösen Sie dieses Problem mit Hilfe von pyMOR und visualisieren Sie die Lösung.

Beachten Sie dabei, dass die Definition der Art der Randwerte schon beim Anlegen des Gebietes erfolgen muss. Nehmen Sie im Zweifel die Dokumentation von RectDomain und StationaryProblem zu Hilfe.

In [7]:
problem = StationaryProblem(
    domain=RectDomain(([0, 0], [1, 1]), right='neumann'),
    diffusion=ConstantFunction(1, 2),
    rhs=ConstantFunction(1, 2),
    neumann_data=ConstantFunction(1, 2)
)
disc, _ = discretize_stationary_cg(problem)
u_h = disc.solve()
disc.visualize(u_h)

(4.ii)

Wir betrachten obiges Problem mit $d = 2$, $\Omega = [0, 1]^2$, $\Gamma_\text{N} = \{x \in \overline{\Omega} \;|\; x_0 = 0 \vee x_0 = 1\}$, $\Gamma_\text{D} = \partial\Omega \backslash \Gamma_{N}$, $A = f = 1$ und

$$g_\text{N}(x) = \begin{cases}1, &x_0 = 0\\-1, &x_0 = 1\\0, &\text{sonst}\end{cases}.$$
  • Modellieren, diskretisieren und lösen Sie dieses Problem mit Hilfe von pyMOR und visualisieren Sie die Lösung.
  • Vergleichen Sie die Lösung $u$ mit der aus Aufgabe 3.ii.
In [8]:
problem = StationaryProblem(
    domain=RectDomain(([0, 0], [1, 1]), left='neumann', right='neumann'),
    diffusion=ConstantFunction(1, 2),
    rhs=ConstantFunction(1, 2),
    neumann_data=ExpressionFunction('(x[..., 0] < 1e-7)*-1. + (x[..., 0] > (1 - 1e-7))*1.', dim_domain=2, shape_range=())
)
disc, _ = discretize_stationary_cg(problem)
u_h = disc.solve()
disc.visualize(u_h)

Aufgabe 5: Diffusion (Wärmeleitkoeffizient)

(i): isolierende Wand

Wir betrachten obiges Problem mit $d = 2$, $\Omega = [0, 1]^2$, $\Gamma_\text{D} = \partial\Omega$, $\Gamma_\text{N} = \emptyset$, $g_D = 0$ und

  • $f(x) = \begin{cases}1, &x \in [0.2, 0.4]^2\\-1, &x \in [0.6, 0.8]^2\\0, &\text{sonst}\end{cases},$
  • $A(x) = \begin{cases}0.01, &x_0 \in [0.45, 0.55]\\1, &\text{sonst}\end{cases}.$

Modellieren, diskretisieren und lösen Sie dieses Problem mit Hilfe von pyMOR und visualisieren Sie die Lösung.

In [9]:
problem = StationaryProblem(
    domain=omega,
    diffusion= ConstantFunction(1, 2) - 0.99*ExpressionFunction('((0.45 <= x[..., 0]) * (x[..., 0] <= 0.55))*1.', dim_domain=2, shape_range=()),
    rhs=ExpressionFunction('((0.2 <= x[..., 0]) * (x[..., 0] <= 0.4) * (0.2 <= x[..., 1]) * (x[..., 1] <= 0.4))*1. + ((0.6 <= x[..., 0]) * (x[..., 0] <= 0.8) * (0.6 <= x[..., 1]) * (x[..., 1] <= 0.8))*-1.', dim_domain=2, shape_range=()),
)
disc, _ = discretize_stationary_cg(problem)
u_h = disc.solve()
disc.visualize(u_h)

(5.2): Wärmeleitkanal

Wir betrachten obiges Problem mit $d = 2$, $\Omega = [0, 1]^2$, $\Gamma_\text{D} = \partial\Omega$, $\Gamma_\text{N} = \emptyset$, $g_D = 0$ und

  • $f(x) = \begin{cases}1, &x \in [0.2, 0.4]^2\\-1, &x \in [0.6, 0.8]^2\\0, &\text{sonst}\end{cases},$
  • $A(x) = \begin{cases}101, &x \in [0.25, 0.75] \times [0.25, 0.35] \cup [0.65, 0.75] \times [0.35, 0.75]\\1, &\text{sonst}\end{cases}.$

Modellieren, diskretisieren und lösen Sie dieses Problem mit Hilfe von pyMOR und visualisieren Sie die Lösung.

In [10]:
problem = StationaryProblem(
    domain=omega,
    diffusion=(
        ConstantFunction(1, 2)
        + 100*(ExpressionFunction('((0.25 <= x[..., 0]) * (x[..., 0] <= 0.75) * (0.25 <= x[..., 1]) * (x[..., 1] <= 0.35))*1.', dim_domain=2, shape_range=())
               + ExpressionFunction('((0.65 <= x[..., 0]) * (x[..., 0] <= 0.75) * (0.35 <= x[..., 1]) * (x[..., 1] <= 0.75))*1.', dim_domain=2, shape_range=()))),
    rhs=ExpressionFunction('((0.2 <= x[..., 0]) * (x[..., 0] <= 0.4) * (0.2 <= x[..., 1]) * (x[..., 1] <= 0.4))*1. + ((0.6 <= x[..., 0]) * (x[..., 0] <= 0.8) * (0.6 <= x[..., 1]) * (x[..., 1] <= 0.8))*-1.', dim_domain=2, shape_range=()),
)
disc, _ = discretize_stationary_cg(problem)
u_h = disc.solve()
disc.visualize(u_h)