Institut für Numerische und Angewandte Mathematik
FB Mathematik und Informatik der Universität Münster
Prof. Dr. Mario Ohlberger, Dr. Stephan Rave und Tobias Wedemeier
Damit die Visualisierungen der Lösungen in das Jupyter Notebook übernommen werden können, muss zu Beginn das folgende IPython Spezialkommando ("magic") ausgeführt werden:
%matplotlib notebook
Aufgabe 1 a)¶
Zunächst importieren wir die gängigsten pyMOR Klassen und Methoden mit der folgenden Anweisung:
from pymor.basic import *
Zur Spezifikation des zu lösenden Problems mit pyMORs Diskretisierungstoolkit müssen wir zunächst das Gebiet $\Omega$ definieren, auf dem die Gleichung zu lösen ist. Hierfür können verschiedene Klassen in pyMOR verwendet werden, die alle von pymor.domaindescriptions.interfaces.DomainDescriptionInterface ableiten.
In unserem Fall, können wir eine RectDomain verwenden:
domain = RectDomain([[0.,0.], [1.,1.]])
Datenfunktionen werden in pyMORs Diskretisierungstoolkit mit Hilfe von Klassen definiert, die von pymor.functions.interfaces.FunctionInterface ableiten. Die konstante Diffusivität geben wir mit Hilfe einer ConstantFunction an:
diffusion = ConstantFunction(1, 2)
Das erste Argument ist hierbei der konstante Wert, den die Funktion annehmen soll. Das zweite argument ist die räumliche Dimension des Gebiets, auf dem die Funktion definiert ist.
Zur Definition der rechten Seite verwenden wir eine ExpressionFunction, die einen beliebigen, die Zuweisungsvorschrift definierenden Python-Ausdruck als String erwartet. In diesem String sind die Koordinaten, an denen die Funktion ausgewertet werden soll, über die Variable x verfügbar. Die wichtigsten NumPy-Funktionen können in dem Ausruck direkt verwendet werden. Es kann aber auch das gesamte NumPy-Modul unter dem namen np angesprochen werden.
Eine Indikatorfunktion, die auf der um $(0.5, 0.5)$ zentrierten Kreisscheibe mit Radius $0.3$ den Wert $1.$ annimmt, ließe sich also wie folgt definieren:
(sqrt((x[0]-0.5)**2 + (x[1]-0.5)**2) <= 0.3) * 1.Allerdings ist zu beachten, dass alle in pyMOR definierten Funktionen vektorisert auswertbar sein müssen. D.h., dass x im Allgemeinen ein beliebig-dimensionales NumPy-Array von meheren Koordinaten ist, bei dem der Index der jeweiligen Koordinate über die letzte Array-Achse spezifiziert wird. Die Korrekte Definition der rechten Seite lautet daher:
rhs = ExpressionFunction('(np.sqrt( (x[...,0]-0.5)**2 + (x[...,1]-0.5)**2) <= 0.3) * 1.', 2, ())
Auch hier wird als zweites Argument die Dimension des Definitionsbereichs der Funktion angegeben. Da dies nicht automatisch von dem angegebenen String-Ausdruck abgeleitet werden kann, muss zusätzlich als drittes Argument die Form, der von der Funktion zurückgegebenen Werte, spezifiziert werden. Für skalare Funktionion ist dies das leere Tupel (). Für Funktionen, die einen dreidimensionalen Vektor zurückgeben, muss (3,) angegeben werden. Für Funktionen, die 2x2-Matrizen zurückgeben, müsste z.B (2,2) angegeben werden.
Schließlich werden alle Daten in einer Problembeschreibungs-Klasse (siehe pymor.analyticalproblems) zusammengefasst. Für stationäre Probleme ist dies StationaryProblem:
problem = StationaryProblem(
domain=domain,
diffusion=diffusion,
rhs=rhs,
)
Die analytischen Probleme dienen als Eingabedaten für 'discretizer' (siehe pymor.discretizers), die aus der gegebenen Problembeschreibung eine Diskrtisierungs-Klasse bauen (siehe pymor.discretizations.interfaces.DiscretizationInterface), welche als Operatoren die assemblierten Finite-Elemente- / Finite-Volumen-Matrizen enthalten.
Für Finitie-Elemente-Verfahren verwenden wir discretize_stationary_cg, welches den maximalen Element-Durchmesser als diameter-Argument erhält:
d, data = discretize_stationary_cg(problem, diameter=1./4.)
Die erhaltene Diskretisierung d können wir lösen und erhalten ein VectorArray (siehe pymor.vectorarrays.interfaces.VectorArrayInterface) zurück:
U = d.solve()
Schließlich visualisieren wir die erhaltene Lösung:
d.visualize(U)
Soll ein bestimmter Gitter-Typ verwendet werden (RectGrid oder TriaGrid), muss die entsprechende Klasse dem Discretizer als grid_type-Argument übergeben werden:
d, data = discretize_stationary_cg(problem, diameter=1./4., grid_type=RectGrid)
d.visualize(d.solve())
Eine Finite-Volumen-Diskretisierung erhalten wir mit discretize_stationary_fv:
d, data = discretize_stationary_fv(problem, diameter=1/32, grid_type=TriaGrid)
d.visualize(d.solve())
Aufgabe 1 b)¶
Um die ausführlichen Log-Meldungen von pyMOR im Folgenden zu unterdrücken verwenden wir die set_log_levels-Methode:
set_log_levels({'pymor': 'WARN'})
Wesentliche Neuerung bei dieser Aufgabe sind die gemischten Randbedingungen, die wir in dem Problem spezifizieren müssen. Hierfür geben wir zunächst bei der Definition von $\Omega$ an, welcher Teil des Randes mit Neumann-Randbedingungen versehen sein soll (Default ist Dirichlet-Randbedingung):
domain = RectDomain(bottom='neumann')
Soll die Funktion $u_N$ nicht Konstant 0 sein, müssen wir diese ebenfalls als pyMOR-Funktion angeben:
neumann_data = ConstantFunction(-1., 2)
Die Diffusivität definieren wir ähnlich wie oben:
diffusion = ExpressionFunction('1. - (np.sqrt( (x[...,0]-0.5)**2 + (x[...,1]-0.5)**2) <= 0.3) * 0.999' , 2, ())
Schließlich fassen wir alles wieder in einem StationaryProblem zusammen und diskretisieren:
problem = StationaryProblem(
domain=domain,
diffusion=diffusion,
neumann_data=neumann_data
)
d, data = discretize_stationary_cg(problem, diameter=1/32)
d.visualize(d.solve())
Aufgabe 1 c)¶
Hier ist die einzige Schwierigkeit die korrekte Definition der Diffusivität:
diffusion = ExpressionFunction(
'1. - (sqrt( (np.mod(x[...,0],1./K)-0.5/K)**2 + (np.mod(x[...,1],1./K)-0.5/K)**2) <= 0.3/K) * 0.999',
2, (),
values={'K': 10}
)
Hierbei haben wir den values-Parameter verwendet, um K als zuätzliche Konstante in dem Ausdruck verfügbar zu machen.
problem = StationaryProblem(
domain=domain,
diffusion=diffusion,
neumann_data=neumann_data
)
d, data = discretize_stationary_cg(problem, diameter=1/100)
d.visualize(d.solve())
Aufgabe 1 d)¶
Die BitmapFunction verwendet die Python Imaging Library (PIL) um ein Graustufenbild in üblichen Grafik-Dateiformaten einzulesen. Wir verwenden hier die folgende Grafik: 
Das resultierende zweidimensionale NumpyArray von Pixel-Werten definiert dann eine stückweise konstante Datenfunktion, wobei der Wertebereich (von schwarzem Pixel zu weißen Pixel) mit dem range-Parameter angegeben wird:
diffusion = BitmapFunction('RB.png', range=[0.001, 1])
problem = StationaryProblem(
domain=domain,
diffusion=diffusion,
neumann_data=neumann_data
)
d, data = discretize_stationary_cg(problem, diameter=1/100)
d.visualize(d.solve())
Die angezeigte Warnung rührt daher, dass die verwendete Grafik einen zusätzlichen Transparenz-Kanal (alpha-Kanal) besitzt, und kann igoriert werden.
Aufgabe 2 a)¶
In Aufgabe 2 wollen wir nun die Differentialgleichung um Parameter ergänzen. Ein Parameter in pyMOR ist immer ein Dictionary von NumPy-Arrays. Die einzelnen Einträge des Dictionarys werden auch als 'Parameter-Komponenten' bezeichnet. Zu jedem Parameter ist ein parameter_type assoziiert, der dem Namen jeder Parameter-Komponent das Shape des zugehörigen Numpy-Arrays zuordnet.
Soll eine parametrische ExpressionFunction definiert werden, wie in dieser Teilaufgabe bei den Neumann-Randwerten, muss dieser parameter_type zusätzlich spezifiziert werden. In diesem Fall haben wir einen Parameter mit nur einer, skalaren Komponente 'neum':
domain = RectDomain(bottom='neumann')
neumann_data = ExpressionFunction('-cos(pi*x[...,0])**2*neum', 2, (), parameter_type= {'neum': ()})
Die Parametrisierung der Datenfunktionen überträgt sich automatisch auf die resultierende Diskretisierung:
diffusion = ExpressionFunction(
'1. - (sqrt( (np.mod(x[...,0],1./K)-0.5/K)**2 + (np.mod(x[...,1],1./K)-0.5/K)**2) <= 0.3/K) * 0.999',
2, (),
values={'K': 10}
)
problem = StationaryProblem(
domain=domain,
diffusion=diffusion,
neumann_data=neumann_data
)
d, data = discretize_stationary_cg(problem, diameter=1/100)
d.parameter_type
Beim Lösen muss dann ein entsprechender Parameter angeben werden:
d.visualize(d.solve({'neum': 1.}))
Im Falle der solve-Methode kann der Parameter oft auch in einer vereinfachten Form angegeben werden:
d.visualize(d.solve(-100))
Aufgabe 2 b)¶
Analog wie oben definieren wir:
diffusion = ExpressionFunction(
'1. - (sqrt( (np.mod(x[...,0],1./K)-0.5/K)**2 + (np.mod(x[...,1],1./K)-0.5/K)**2) <= 0.3/K) * (1 - diffu)',
2, (),
values={'K': 10},
parameter_type= {'diffu': ()}
)
problem = StationaryProblem(
domain=domain,
diffusion=diffusion,
neumann_data=neumann_data
)
d, data = discretize_stationary_cg(problem, diameter=1/100)
d.parameter_type
pyMOR leitet hier automatisch ab, dass die Lösung der Diskretisierung von zwei Parameter-Komponenten abhängt. Diese können bei solve direkt in alphabetischer Reihenfolge angegeben werden:
d.visualize(d.solve({'diffu': 0.001, 'neum': 1}))
d.visualize(d.solve([1, -1]))
Aufgabe 2 c)¶
Hier benötigen wir einzelne Grafiken für die Buchstaben 'R' und 'B': 
diffusion_R = BitmapFunction('R.png', range=[1, 0])
diffusion_B = BitmapFunction('B.png', range=[1, 0])
Um diese Funktionen zu einer Gesamt-Funktion zusammenzusetzen verwenden wir eine LincombFunction die eine Liste von Funkionen und eine Liste von Linearkoeffizieten erhält. Parametrische Linearkoeffizienten werden mit Hilfe von ParameterFunktionalen (siehe pymor.parameters.functionals) realisiert.
diffusion = LincombFunction(
[ConstantFunction(1., 2), diffusion_R, diffusion_B],
[1., ExpressionParameterFunctional('-(1 - R)', {'R': ()}), ExpressionParameterFunctional('-(1 - B)', {'B': ()})]
)
diffusion.parameter_type
problem = StationaryProblem(
domain=domain,
diffusion=diffusion,
neumann_data=ConstantFunction(-1, 2)
)
d, data = discretize_stationary_cg(problem, diameter=1/100)
d.visualize(d.solve([1., 0.001]))
d.visualize(d.solve([0.001, 1]))
