Praktikum zu Vorlesung Modellreduktion parametrisierter Systeme
Mario Ohlberger, Felix Schindler, Tim Keil
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.
solve_reduced(d, basis, mu)
, welche zu einer gegebenen Diskretisierung d
und einer reduzierten Basis basis
die RB-Approximation der Lösung zum Parameter mu
bestimmt. Verwenden Sie dabei die apply2
- und apply_adjoint
- Methoden von d.operator
und d.rhs
, um das reduzierte Problem zu assemblieren.build_basis(d, N)
, um zu einer Diskretisierung d
eine reduzierte Basis der Größe N
aus 2*N
zufälligen Lösungssnapshots zu berechnen. Verwenden Sie die Methode greedy_approximation
von Blatt 6.solve_reduced
reduzierte Lösungen für verschiedene Basisgrößen (mit oder ohne Greedy-Verfahren), und visualisieren Sie die reduzierte Lösung zusammen mit der hochdimensionalen Lösung sowie der Differenz der beiden Lösungen.pyMOR's project
-Methode ermöglicht die Reduzierte-Basis-Projektion von beliebigen pyMOR-Operatoren. Hierfür muss project
lediglich der zu projizierende Operator sowie reduzierte Basen für den range
- und source
-VectorSpace
des Operators übergeben werden. Als Ergebnis erhält man reduzierte pyMOR-Operatoren, welche die entsprechenden reduzierten Systemmatrizen enthalten.
Affine Zerlegungen, die im analytischen Problem mittels einer LincombFunction
dargestellt wurden, werden von discretize_stationary_cg
/ discretize_stationary_fv
automatisch in affine Zerlegungen der zugehörigen diskreten Operatoren überführt. Diese werden analog zur LincombFunction
durch einen LincombOperator
dargestellt:
print(type(d.operator))
print(d.operator.operators)
print(d.operator.coefficients)
Die project-Methode berücksichtigt automatisch solche affinen Operator-Zerlegungen und assembliert einen reduzierten LincombOperator, indem (rekursiv) die einzelnen Operatoren in d.operator.operators projiziert werden.
pymor.algorithms.projection.project
, um für das parametrische Problem von Blatt 4 Aufgabe 1 einen RB-projizierten, affin-zerlegten Systemoperator und ein RB-projiziertes Rechte-Seite-Funktional zu erhalten. Wiederholen Sie hiermit die Experimente von Aufgabe 1. Messen Sie dabei zusätzlich die durchschnittlich benötigte Zeit zur Lösung des reduzierten Problems (ohne hochdimensionale Rekonstruktion), und plotten Sie diese gegen die BasisgrößeHinweis: Nutzen Sie die projected_rhs.as_vector()
um eine VectorArray
- Darstellung der rechnten Seite zu erhalten.
def compute_errors(d, basis, test_mus):
#Your code
return errors, proj_errors, times
errors, proj_errors, times = compute_errors(d, basis, test_mus)
Wir plotten direkt den Speedup, den das Lösen des reduzierten Problems gegenüber dem Lösen des hochdimensionalen Problems bringt:
def plot_errors(proj_errors, errors, times, d):
tic = time.time()
d.solve(d.parameter_space.sample_randomly(1, seed=9876)[0])
fom_time = time.time() - tic
plt.figure()
plt.subplot(1, 2, 1)
plt.semilogy(np.max(proj_errors, axis=1), label='orth. proj.')
plt.semilogy(np.max(errors, axis=1), label='RB')
plt.legend()
plt.subplot(1, 2, 2)
plt.plot(np.median(fom_time / times, axis=1), label='speedup')
plt.legend()
plot_errors(proj_errors, errors, times, d)
Lincomb
Operatoren besteht, speichern Sie die Funktion als problem_B5_A2_parametric_decomb
und wiederholen Sie erneut das Experiment von oben. Was ist hier noch falsch? BoundaryL2ProductFunctionalP1
Offline-Online zerlegt werden kann und Wiederholen Sie das Experiment.print(d.rhs.coefficients)
print(d.rhs.operators)
print(d.rhs.operators[1])
Bemerkung: Im Falle von einem Problem welches nur eine Parameterabhängigkeit für die rechte Seite aufweist beobachten wir einen deutlich geringeren Speedup als für die übrigen Probleme. Dies liegt daran, dass hier die Systemmatrix nicht von mu
abhängt und der verwendete Gleichungssystem-Löser die LR-Zerlegung dieser Matrix im Speicher halten kann. Damit können, nach einmaliger Zerlegung, neue Lösungen für verschiedene rechte Seiten sehr effizient berechnet werden.
Bisher haben wir die RB-Projektionen von d.operator
und d.rhs
manuell durchgeführt. Die in pyMOR definierten Reduktoren
führen RB-Projektionen automatisch für alle in der gegebenen Diskretisierung enthaltenen Operatoren (einschließlich Skalarprodukte und etwaige Ausgabefunktionale) durch und liefern eine diese projizierten Operatoren enthaltende reduzierte Diskretisierung zurück.
Hierfür muss der Reduktor mit der zu reduzierenden Diskretisierungen und einer reduzierten Basis initialisiert werden. Die reduzierte Diskretisierung wird dann von der reduce
-Methode des Reduktors berechnet. Die hochdimensionale Rekonstruktion eines RB-Lösungsvektors ist mittels der reconstruct
-Methode möglich.
GenericRBReductor
verwenden, um eine RB-projizierte reduzierte Diskretisierung zu erhalten. Führen Sie das Experiment aus Aufgabe 2 erneut durch. Hinweis: Die reduce
Methode kann zusätzlich eine gewünschte reduzierte Dimension übergeben werden (die kleiner als die tatsächliche Dimension der reduzierten Basis sein muss). In diesem Fall wird aus den reduzierten Matrizen zur vollen reduzierten Basis ein entsprechender Teilblock ausgeschnitten, sodass nach erstmaliger Berechnung der reduzierten Diskretisierung, sehr schnell reduzierte Modelle für Anfangsstücke der reduzierten Basis berechnet werden können.
def compute_errors(d, basis, test_mus):
#Your code
return errors, proj_errors, times
CoerciveRBReductor
, um zusätzlich einen Fehlerschätzer für den Modellreduktionsfehler zu assemblieren. Leiten Sie hierfür zunächst eine geeignete untere Schranke für die Koerzivitätskonstante des jeweils betrachteten Problem her. Wiederholen Sie erneut das Experimant von oben und visualisieren Sie im selben Plot den Fehlerschätzer. Plotten Sie außerden empirisch bestimmte minimale und maximale Effizienz des Fehlerschätzers gegen die Basisgröße.def compute_errors(d, basis, test_mus):
#Your code
return errors, estimates, proj_errors, times
errors, estimates, proj_errors, times = compute_errors(d, basis, test_mus)
plt.figure()
plt.subplot(1, 2, 1)
plt.semilogy(np.max(proj_errors, axis=1), label='orth. proj.')
plt.semilogy(np.max(errors, axis=1), label='RB')
plt.semilogy(np.max(estimates, axis=1), label='est.')
plt.legend()
plt.subplot(1, 2, 2)
plt.plot(np.min(errors / estimates, axis=1), label='min. estimator effectivity')
plt.legend()
pyMOR
die Funktionsweise von CoerciveRBReductor.reduce()
nachzuvollziehen. Wie unterscheidet sich die Offline-Online-Zerlegung des Fehlerschätzers in CoerciveRBReductor
von der in der Vorlesung behandelten Zerlegung?