Institut für Numerische und Angewandte Mathematik
FB Mathematik und Informatik der Universität Münster
Prof. Dr. Mario Ohlberger und Dr. Stephan Rave
%matplotlib notebook
from pymor.basic import *
from problems import *
from algorithms import *
from matplotlib import pyplot as plt
import time
set_log_levels({'pymor': 'WARN',
'pymor.functions.bitmap': 'CRITICAL'})
Aufgabe 1 a)¶
Um die reduzierte Systemmatrix und die reduzierte rechte Seite zu assemblieren, müssen wir d.operator sowie d.rhs auf den Vektoren der gegebenen reduzierten Basis auswerten. d.operator fassen wir dabei als Bilinearform auf, die wir wie gewohnt mittels apply2 auswerten. Für das Rechte-Seite-Funktional verwenden wir die apply-Methode, welche ein neues VectorArray mit den Operatorauswertungen zurückliefert. Dabei haben Funktionale in pyMOR grundsätzlich als range den Raum NumpyVectorSpace(1). Mit d.rhs.apply(basis).data erhalten wir also den reduzierten Rechte-Seite-Vektor als ein len(basis) x 1 NumPy-Array.
Außerdem ist zu beachten, dass d.operator und d.rhs parametrisch sein können, weshalb wir apply bzw. apply2 zusätzlich den Parameter mu übergeben müssen:
def solve_reduced(d, basis, mu):
lhs = d.operator.apply2(basis, basis, mu=mu)
rhs = d.rhs.apply(basis, mu=mu).data
u = np.linalg.solve(lhs, rhs).ravel()
U = basis.lincomb(u)
return U
Aufgabe 1 b)¶
Wir verwenden hier die Methode greedy_approximation von Blatt 3, um eine reduzierte Basis der Größe N aus 2*N zufälligen Lösungssnapshots zu berechnen:
def build_basis(d, N):
snapshots = d.solution_space.empty()
for mu in d.parameter_space.sample_randomly(2*N):
snapshots.append(d.solve(mu))
basis, _ = greedy_approximation(snapshots, N, product=d.h1_0_semi_product, orthogonalize=True)
return basis
d, _ = discretize_stationary_cg(problem_1_2_c())
basis = build_basis(d, 30)
Wir berechnen gleichzeitig für zwei verschiedene Parameter die reduzierten Lösungen für verschiedene Größen der reduzierten Basis:
Us = [d.solution_space.empty(), d.solution_space.empty()]
U_rbs = [d.solution_space.empty(), d.solution_space.empty()]
for i_mu, mu in enumerate(d.parameter_space.sample_randomly(2, seed=5678)):
U = d.solve(mu)
for N in range(len(basis) + 1):
Us[i_mu].append(U)
U_rbs[i_mu].append(solve_reduced(d, basis[:N], mu))
d.visualize hat neben separate_colorbars zusätzlich den Parameter rescale_colorbars, der eine Reskalierung der Colorbar mit jedem Zeitschritt erzwingt:
d.visualize((Us[0], U_rbs[0], Us[0] - U_rbs[0],
Us[1], U_rbs[1], Us[1] - U_rbs[1]),
separate_colorbars=True, rescale_colorbars=True, columns=3)
Aufgabe 1 c)¶
Da die diskreten Lösungsvektoren des hochdimensionalen Problems für die gewählte Auflösung noch vergleichsweise klein sind, können wir es uns erlauben, diese für alle Test-Parameter im Speicher vorzuhalten:
test_mus = d.parameter_space.sample_randomly(100, seed=5678)
Us = d.solution_space.empty()
for mu in test_mus:
print('.', end='', flush=True)
Us.append(d.solve(mu))
Dies beschleunigt dann die Berechnung der Modellreduktionsfehler:
errors = np.empty((len(basis) + 1, len(test_mus)))
proj_errors = np.empty((len(basis) + 1, len(test_mus)))
for N in range(len(basis) + 1):
print('.', end='', flush=True)
for i_mu, mu in enumerate(test_mus):
U_rb = solve_reduced(d, basis[:N], mu)
errors[N, i_mu] = d.h1_0_semi_norm(Us[i_mu] - U_rb)
U_proj = orthogonal_projection(Us[i_mu], basis[:N], product=d.h1_0_semi_product)
proj_errors[N, i_mu] = d.h1_0_semi_norm(Us[i_mu] - U_proj)
Neben den maximalen Fehlern für orthogonale Projektion und reduzierte Lösung plotten wir zusätzlich noch den minimalen Quotienten zwischen beiden Fehlern als Maß für den Verlust an Genauigkeit, wenn die Bestapproximation im reduzierten Raum (mittels orthogonaler Projektion) durch die Galerkin-Projektion ersetzt wird:
plt.figure()
plt.subplot(1, 2, 1)
plt.semilogy(np.arange(len(basis) + 1), np.max(proj_errors, axis=1), label='orth. proj.')
plt.semilogy(np.arange(len(basis) + 1), np.max(errors, axis=1), label='RB')
plt.legend()
plt.subplot(1, 2, 2)
plt.plot(np.min(proj_errors / errors, axis=1), label='min. efficiency')
plt.legend()
Aufgabe 2 a)¶
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.
Um das reduzierte Problem zu lösen, kann der gleiche Code wie in StationaryDiscretization._solve verwendet werden: Zunächst muss das reduzierte Rechte-Seite-Funktional mittels as_source_array für einen gegebenen Parameter in ein den assemblierten Rechte-Seite-Vektor enthaltendes VectorArray überführt werden. Mittels der apply_inverse-Methode des reduzierten Systemoperators kann dann das reduzierte problem mit diesem VectorArray als rechter Seite gelöst werden.
Ansonsten verläuft die Berechnung der Modellreduktionsfehler analog zu Aufgabe 1 c):
def compute_errors(d, basis, test_mus):
Us = d.solution_space.empty()
for mu in test_mus:
print('.', end='', flush=True)
Us.append(d.solve(mu))
print()
errors = np.empty((len(basis) + 1, len(test_mus)))
proj_errors = np.empty((len(basis) + 1, len(test_mus)))
times = np.empty((len(basis) + 1, len(test_mus)))
for N in range(len(basis) + 1):
print('.', end='', flush=True)
op = project(d.operator, basis[:N], basis[:N])
rhs = project(d.rhs, None, basis[:N])
for i_mu, mu in enumerate(test_mus):
tic = time.time()
u_rb = op.apply_inverse(rhs.as_source_array(mu), mu=mu)
times[N, i_mu] = time.time() - tic
U_rb = basis[:N].lincomb(u_rb.data)
errors[N, i_mu] = d.h1_0_semi_norm(Us[i_mu] - U_rb)
U_proj = orthogonal_projection(Us[i_mu], basis[:N], product=d.h1_0_semi_product)
proj_errors[N, i_mu] = d.h1_0_semi_norm(Us[i_mu] - U_proj)
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)
Aufgabe 2 b)¶
In problems.py wurden bereits die Anpassungen von Aufgabe 2 c) durchgeführt. Die ursprünglichen Implementierungen sind weiterhin als problem_1_2_a_no_decomp und problem_1_2_b_no_decomp verfügbar:
d, _ = discretize_stationary_cg(problem_1_2_a_no_decomp())
basis = build_basis(d, 30)
test_mus = d.parameter_space.sample_randomly(100, seed=5678)
errors, proj_errors, times = compute_errors(d, basis, test_mus)
plot_errors(proj_errors, errors, times, d)
d, _ = discretize_stationary_cg(problem_1_2_b_no_decomp())
basis = build_basis(d, 30)
test_mus = d.parameter_space.sample_randomly(100, seed=5678)
errors, proj_errors, times = compute_errors(d, basis, test_mus)
plot_errors(proj_errors, errors, times, d)
In beiden Fällen ergibt sich aufgrund der fehlenden Offline-Online-Zerlegung nur ein sehr geringer Speedup. Die von pyMOR ausgegebene Warnung weist darauf hin, dass keine offline-online-zerlegte Projektion des jeweiligen Operators möglich ist und daher ein generischer ProjectedOperator angelegt wird, der zwar ein mathematisch korrektes Ergebnis liefert, in apply jedoch stets eine Auswertung des ursprünglichen hochdimensionalen Operators benötigt.
Aufgabe 2 c)¶
In problems.py wurden die Definitionen von problem_1_2_a und problem_1_2_b so angepasst, dass alle Datenfunktionen durch LincombFunctions gegeben sind, deren einzelnen Komponenten keine Parameterabhängigkeit mehr aufweisen:
%cat problems.py
d, _ = discretize_stationary_cg(problem_1_2_a())
basis = build_basis(d, 30)
test_mus = d.parameter_space.sample_randomly(100, seed=5678)
errors, proj_errors, times = compute_errors(d, basis, test_mus)
plot_errors(proj_errors, errors, times, d)
Im Falle von problem_1_2_a 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.
d, _ = discretize_stationary_cg(problem_1_2_b())
basis = build_basis(d, 30)
test_mus = d.parameter_space.sample_randomly(100, seed=5678)
errors, proj_errors, times = compute_errors(d, basis, test_mus)
plot_errors(proj_errors, errors, times, d)
Aufgabe 3 a)¶
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. reduce 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.
Damit lässt sich der Modellreduktionsfehler effizient, ohne Zwischenspeicherung der hochdimensionalen Lösungen wie folgt berechnen:
def compute_errors(d, basis, test_mus):
reductor = GenericRBReductor(d, basis)
errors = np.empty((len(basis) + 1, len(test_mus)))
proj_errors = np.empty((len(basis) + 1, len(test_mus)))
times = np.empty((len(basis) + 1, len(test_mus)))
for i_mu, mu in enumerate(test_mus):
print('.', end='', flush=True)
for N in range(len(basis) + 1):
rd = reductor.reduce(N)
tic = time.time()
u_rb = rd.solve(mu)
times[N, i_mu] = time.time() - tic
U_rb = reductor.reconstruct(u_rb)
errors[N, i_mu] = d.h1_0_semi_norm(Us[i_mu] - U_rb)
U_proj = orthogonal_projection(Us[i_mu], basis[:N], product=d.h1_0_semi_product)
proj_errors[N, i_mu] = d.h1_0_semi_norm(Us[i_mu] - U_proj)
return errors, proj_errors, times
d, _ = discretize_stationary_cg(problem_1_2_c())
basis = build_basis(d, 30)
test_mus = d.parameter_space.sample_randomly(100, seed=5678)
errors, proj_errors, times = compute_errors(d, basis, test_mus)
plot_errors(proj_errors, errors, times, d)
Aufgabe 3 b)¶
CoerciveRBReductor ist identisch zu GenericRBReductor, assembliert jedoch zusätzlich noch für koerzive Probleme einen Fehlerschätzer für den Modellreduktionsfehler. Dabei müssen wir das korrekte Skalarprodukt zur Berechnung der Riesz-Repräsentanten, sowie eine untere Abschätzung der Koerzivitätskonstante als ParameterFunctional angeben:
def compute_errors(d, basis, test_mus):
reductor = CoerciveRBReductor(
d, basis, product=d.h1_0_semi_product,
coercivity_estimator=ExpressionParameterFunctional('min([1, R, B])', d.parameter_type)
)
errors = np.empty((len(basis) + 1, len(test_mus)))
estimates = np.empty((len(basis) + 1, len(test_mus)))
proj_errors = np.empty((len(basis) + 1, len(test_mus)))
times = np.empty((len(basis) + 1, len(test_mus)))
for i_mu, mu in enumerate(test_mus):
print('.', end='', flush=True)
U = d.solve(mu)
for N in range(len(basis) + 1):
rd = reductor.reduce(N)
tic = time.time()
u_rb = rd.solve(mu)
estimates[N, i_mu] = rd.estimate(u_rb, mu)
times[N, i_mu] = time.time() - tic
U_rb = reductor.reconstruct(u_rb)
errors[N, i_mu] = d.h1_0_semi_norm(U - U_rb)
U_proj = orthogonal_projection(U, basis[:N], product=d.h1_0_semi_product)
proj_errors[N, i_mu] = d.h1_0_semi_norm(U - U_proj)
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()
