Blatt 07, 19.06.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 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.

Aufgabe 1: RB-Projektion¶

Schreiben Sie eine Funktion 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.operatorund d.rhs, um das reduzierte Problem zu assemblieren.

Schreiben Sie eine Funktion 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.

Diskretisieren Sie das parametrische Problem von Blatt 4 Aufgabe 1 und berechnen Sie mittels 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.

Bestimmen Sie empirisch den maximalen Modellreduktionsfehler aus 100 Parametern und plotten Sie diesen zusammen mit dem maximalen Fehler der orthogonalen Projektion auf den RB-Raum in Abhängigkeit von der Basisgröße.

Aufgabe 2: Offline-Online-Zerlegung¶

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.

Benutzen Sie 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öße

Hinweis: 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)

Wiederholen Sie das selbe Experiment für das parametrische Problem von Blatt 5 Aufabe 2. Was fällt Ihnen auf?

Passen Sie ihre Definition von Blatt 5 Aufgabe 2 so an, dass die Datenfunktionen nur aus 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?

Definieren Sie ihre Diskretisierung so, dass auch das 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.

Aufgabe 3: Reduktoren und Fehlerschätzer¶

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.

Vereinfachen Sie den Code aus Aufgabe 2, indem Sie 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

Verwenden Sie nun 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()

Versuchen Sie im Quellcode von 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?