In [ ]:
import sys
# Diese Zeile muss angepasst werden!
sys.path.append("/home/moser/MG_2016/pyMG-2016/")
# sys.path.append("/home/dima/lectures/MG/pyMG/")
import scipy as sp
import scipy.sparse.linalg as sprsla
import numpy as np
import matplotlib as mplt
import matplotlib.pyplot as plt
%matplotlib inline
In [ ]:
import pyamg
from project.gauss_seidel import GaussSeidel
from project.weighted_jacobi import WeightedJacobi
from project.pfasst.plot_tools import eigvalue_plot_list, matrix_plot, matrix_row_plot,heat_map
from project.pfasst.transfer_tools import to_dense
from project.pfasst.matrix_method_tools import matrix_power
from pymg.collocation_classes import CollGaussRadau_Right, CollSplineRight
from project.mymultigrid import MyMultigrid
from project.poisson1d import Poisson1D
from project.linear_transfer import LinearTransfer
from project.weighted_jacobi import WeightedJacobi
In [ ]:
def plot_mat_and_eigs(A):
A = to_dense(A)
fig1, (ax2, ax) = plt.subplots(ncols=2,figsize=(15,5))
plt_mat = ax.imshow(A, cmap=plt.cm.jet, interpolation='nearest')
ax.spines['left'].set_position(('outward', 10))
ax.spines['bottom'].set_position(('outward', 10))
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.yaxis.set_ticks_position('left')
ax.xaxis.set_ticks_position('bottom')
plt.colorbar(plt_mat)
eig=sp.linalg.eigvals(to_dense(A))
real_part = np.real(eig)
img_part = np.imag(eig)
ax2.plot(real_part,img_part,'bo')
ax2.set_xlabel("real part")
ax2.set_ylabel("img part")
ax2.set_title('eigenvalues')
fig1.tight_layout()
plt.show()
def plot_solution(sol,N_t):
heat_map(sol.reshape(N_t,-1))
def lineplot(vertices, indices, linewidths=1):
"""Plot 2D line segments"""
vertices = np.asarray(vertices)
indices = np.asarray(indices)
#3d tensor [segment index][vertex index][x/y value]
lines = vertices[np.ravel(indices),:].reshape((indices.shape[0],2,2))
col = mplt.collections.LineCollection(lines)
col.set_color('k')
col.set_linewidth(linewidths)
sub = mplt.pylab.gca()
sub.add_collection(col,autolim=True)
sub.autoscale_view()
In [ ]:
def system_matrix_poisson1d(n):
hh1d = Poisson1D(n)
return hh1d.A
def spec_rad(A):
return np.max(np.abs(sp.linalg.eigvals(to_dense(A))))
Wir beginnnen mit einem Vergleich zwischen geometrischem und algebraischen Multigrid. Zunächst geometrisches Multigrid.
In [ ]:
ndofs = 15
niter_list = []
nlevels = int(np.log2(ndofs+1))
prob = Poisson1D(ndofs=ndofs)
mymg = MyMultigrid(ndofs=ndofs, nlevels=nlevels)
mymg.attach_transfer(LinearTransfer)
mymg.attach_smoother(WeightedJacobi,prob.A,omega=2.0/3.0)
# mymg.attach_smoother(GaussSeidel,prob.A)
k = 6
xvalues = np.array([(i+1) * prob.dx for i in range(prob.ndofs)])
prob.rhs = (np.pi*k)**2 * np.sin(np.pi*k*xvalues)
uex = sprsla.spsolve(prob.A, prob.rhs)
res = 1
niter = 0
err = []
u = np.zeros(uex.size)
res_list_mg = []
print nlevels
while res > 1E-10 and niter < 10:
niter += 1
u = mymg.do_v_cycle(u, prob.rhs, 2, 2, 0)
res = np.linalg.norm(prob.A.dot(u)-prob.rhs, np.inf)
res_list_mg.append(res)
err.append(np.linalg.norm(u-uex, np.inf))
print(niter,res,err[-1])
In [ ]:
plot_mat_and_eigs(prob.A)
In [ ]:
mls = pyamg.ruge_stuben_solver(prob.A,max_levels = 10, max_coarse = 1,CF='RS',keep=True)
print mls
res_list_amg = []
sol = mls.solve(b=prob.rhs, x0=np.zeros(ndofs),tol=1e-10,residuals=res_list_amg)
In [ ]:
plt.semilogy(res_list_amg)
plt.semilogy(res_list_mg)
In [ ]:
# The CF splitting, 1 == C-node and 0 == F-node
splitting = mls.levels[0].splitting
C_nodes = splitting == 1
F_nodes = splitting == 0
plt.figure(figsize=(10,2))
V = np.zeros((ndofs,2))
V[:,0] = prob.domain
plt.scatter(V[:,0][C_nodes], V[:,1][C_nodes], c='r', s=200.0) #plot C-nodes in red
plt.scatter(V[:,0][F_nodes], V[:,1][F_nodes], c='b', s=200.0) #plot F-nodes in blue
plt.show()
Frage
AMG führt hier zum selben Ergebnis wie ein einfaches geometrischen Multigrid. Warum sollte man also AMG verwenden?
Anstatt
$$ \Delta u(x) = 0,\quad u(0)=u(1)=0 $$betrachten wir
$$ u_t(t,x)=\Delta u(t,x) ,\quad u(t,0)=u(t,1)=0, u(0,x) = u_0. $$Wir verwenden die Rothe Methode um in jedem Schritt wieder eine elliptische PDGL zu lösen.
$$ u_{i+1}(x)-u_{i}(x) = \Delta t \Delta u_{i+1}(x) $$Wir diskretisieren weiter ...
$$ u_{i+1}-u_{i} = \Delta t \mathbf{A} u_{i+1} \\ (\mathbf{I} - \Delta t \mathbf{A})u_{i+1} = u_i $$Mehrere Schritte dieser Art können als ein System geschrieben werden.
In [ ]:
def impl_euler_system(n_steps,dt,A):
I_L = sp.sparse.eye(n_steps,format='csc')
I_N = sp.sparse.eye(A.shape[0],format='csc')
if n_steps == 1:
return I_N - dt*A
else:
E_mat = sp.sparse.diags([1], [-1], shape=(n_steps, n_steps)).tocsc()
As = [[I_L,I_N],[-I_L*dt,-A],[-E_mat,I_N]]
return reduce(lambda a, b: a + b, map(lambda a_row: reduce(lambda a, b: sp.sparse.kron(a, b), a_row), As))
def impl_euler_rhs(n_steps,domain,k):
rhs = np.zeros(n_steps*domain.size)
init = np.sin(domain*k*np.pi)
rhs[:domain.size] = init
return rhs
def impl_euler_x0(n_steps,domain,k):
init = np.sin(domain*k*np.pi)
return np.kron(np.ones(n_steps), init)
In [ ]:
plot_mat_and_eigs(impl_euler_system(2,0.1,prob.A))
In [ ]:
n_steps = 20
dt = 1e-2
k = 2
M = impl_euler_system(n_steps,dt,prob.A)
rhs = impl_euler_rhs(n_steps,prob.domain,k)
x0 = impl_euler_x0(n_steps,prob.domain,k)
sol=sprsla.spsolve(M,rhs)
plot_solution(sol,n_steps)
In [ ]:
mls = pyamg.ruge_stuben_solver(M,max_levels = 10, max_coarse = 1,CF='RS',keep=True)
print mls
res_amg_list = []
sol = mls.solve(b=rhs, x0=x0,tol=1e-10,residuals=res_amg_list)
plt.semilogy(res_amg_list)
In [ ]:
# The CF splitting, 1 == C-node and 0 == F-node
V = np.zeros((ndofs*n_steps,2))
t_space = np.array([(i) * dt for i in range(n_steps)])
x_space = prob.domain
for j in range(n_steps):
for k in range(ndofs):
V[ndofs*j+k,0]=x_space[k]
V[ndofs*j+k,1]=t_space[j]
E = np.vstack((M.tocoo().row,M.tocoo().col)).T
splitting = mls.levels[0].splitting
C_nodes = splitting == 1
F_nodes = splitting == 0
plt.figure(figsize=(10,10))
lineplot(V,E,0.1)
plt.scatter(V[:,0][C_nodes], V[:,1][C_nodes], c='r', s=200.0) #plot C-nodes in red
plt.scatter(V[:,0][F_nodes], V[:,1][F_nodes], c='b', s=200.0) #plot F-nodes in blue
plt.show()
Aufgabe
Probieren Sie verschiedene $dt$ aus, wie verändert sich die Vergröberungsstrategien von AMG? Welche Vergröberungsstrategien erkennt man?
Wichtiger Bestandteil für unseren Formalismus ist die Quadratur-Matrix $Q$. Sei $v$ der Vektor der die Funktionsauswertungen der Funktion $f$ an den Quadratur-Knoten beinhaltet, so entspricht $Qv$ einer numerischen Approximation von $F(t)-F(t_0)$ ausgewertet an
In [ ]:
def Q_mat_and_nodes(N, which='radau_r'):
if which is 'radau_r':
cgr = CollGaussRadau_Right(N,0.0,1.0)
return cgr.Qmat[1:,1:], cgr.nodes
elif which is 'spline':
spl = CollSplineRight(N,0,1,4)
return spl.Qmat[1:,1:], spl.nodes
elif which is 'radau_delta':
cgr = CollGaussRadau_Right(N,0.0,1.0)
return cgr.QDmat, cgr.nodes
elif which is 'spline_delta':
spl = CollSplineRight(N,0,1,1)
return spl.QDmat, spl.nodes
Aufgabe:
Plotten Sie mithilfe von plot_mat_and_eigs die Q-Matrix für verschiedene $N$ mit und ohne splines an. Wann geht das ganze schief?
In [ ]:
Aufgabe
Testen sie nun für verschiedene Konfigurationen, wie gut Monome bis zur Ordnung 20 integriert werden.
In [ ]:
# Als Beispiel der lineare Fall
N = 3
Q,nodes = Q_mat_and_nodes(N)
v = nodes.copy()
e = Q.dot(v)-v**2/2
print np.linalg.norm(e)
In [ ]:
def test_monom_quadrature(Q,nodes,order):
v = nodes.copy()**order
e = # Zeile hier einfügen
return np.linalg.norm(e)
In [ ]:
Q_5,nodes_5 = Q_mat_and_nodes(5)
Q_10,nodes_10 = Q_mat_and_nodes(10)
Q_50,nodes_50 = Q_mat_and_nodes(50)
Q_s_5,nodes_s_5 = Q_mat_and_nodes(5,'spline')
Q_s_10,nodes_s_10 = Q_mat_and_nodes(10,'spline')
Q_s_50,nodes_s_50 = Q_mat_and_nodes(50,'spline')
In [ ]:
fig1, (ax1,ax2) = plt.subplots(ncols=2,figsize=(15,5),sharey=True)
ax1.semilogy(range(1,20),map(lambda order:test_monom_quadrature(Q_5,nodes_5,order),range(1,20)))
ax1.semilogy(range(1,20),map(lambda order:test_monom_quadrature(Q_10,nodes_10,order),range(1,20)))
ax1.semilogy(range(1,20),map(lambda order:test_monom_quadrature(Q_50,nodes_50,order),range(1,20)))
ax1.set_title("Gauss Radau")
ax2.semilogy(range(1,20),map(lambda order:test_monom_quadrature(Q_s_5,nodes_s_5,order),range(1,20)))
ax2.semilogy(range(1,20),map(lambda order:test_monom_quadrature(Q_s_10,nodes_s_10,order),range(1,20)))
ax2.semilogy(range(1,20),map(lambda order:test_monom_quadrature(Q_s_50,nodes_s_50,order),range(1,20)))
ax2.set_title("Splines")
plt.show()
Frage
Wo ist der Unterschied zwischen GaussRadau und Splines?
Wie wir in der Vorlesung gesehen haben, kann die Q-Matrix dazu genutzt werden um das Kollokationsproblem
$$ \left( \mathbf{I} - \Delta t \cdot \mathbf{Q}\otimes\mathbf{A} \right)U = U_0 $$aufzustellen. Wobei A hierbei die Systemmatrix für das periodische Problem ist.
Aufgabe
Schreiben Sie eine Funktion colloc_mat(Nt,Nx,sig,dt,typ='radau_r'), welche Kollokationsmatrix abhängig von $\sigma,dt,N_x,$ und $N_t$ aufstellt.
In [ ]:
def colloc_mat(Nt,Nx,dt,typ='radau_r'):
I = # Zeile hier einfügen
Q,nod = Q_mat_and_nodes(Nt,typ)
A = -system_matrix_poisson1d(Nx)
return # Zeile hier einfügen
Buntebilderaufgabe
Plotten Sie mithilfe von plot_mat_and_eigs einige Konfigurationen. Welche Strukturen sind zu erkennen?
In [ ]:
plot_mat_and_eigs(to_dense(colloc_mat(3,10,0.1)))
Nun brauchen wir noch eine rechte Seite um das lineare Gleichungssystem zu vervollständigen. Die rechte Seite entsteht durch einen "spread" auf die einzelnen Quadraturknoten. Wir werden nun Funktion $s_k(x) = \sin(k\pi x)$, ausgewertet auf den Raumknoten als Anfangswert verwenden.
In [ ]:
def colloc_rhs(Nt,Nx,k):
dx = 1.0/(Nx+1)
domain = np.array([(i+1) * dx for i in range(Nx)])
init = np.sin(domain*k*np.pi)
return np.kron(np.ones(Nt),init)
In [ ]:
plot_solution(colloc_rhs(10,16,2),10)
Aufgabe
Finden sie mithilfe von scipy die Lößung des Kollokationsproblems und plotten Sie diese für $N_t=10, N_x=32, dt=0.05$. Probieren sie verschiedene Konfigurationen aus. Was passiert für unterschiedliche $k$?
In [ ]:
Nx = 24
Nt = 10
k=2
M = # Zeile hier einfügen
rhs = # Zeile hier einfügen
sol = sprsla.spsolve(M,rhs)
plot_solution(sol,Nt)
Desweiteren wollen wir nun mehrere dieser Kollokationsprobleme hintereinander aufreihen, wir stellen dementsprechend die Matrix $$ \mathbf{I}_L\otimes \mathbf{I}_M \otimes \mathbf{I}_N -\mathbf{I}_L\otimes\mathbf{Q}\otimes\mathbf{A}-\mathbf{E}\otimes\mathbf{N}\otimes\mathbf{I}_N $$ auf, wobei $L$ die Anzahl der Subintervale, $N$ die Anzahl der Freiheitsgraden im Raum und $M$ die Anzahl der Quadraturknoten ist.
In [ ]:
def compound_colloc_mat(N,M,L,dt,typ='radau_r'):
if L == 1:
return colloc_mat(M,N,dt,typ)
else:
I_N = sp.sparse.eye(N,format='csc')
I_M = sp.sparse.eye(M,format='csc')
I_L = sp.sparse.eye(L,format='csc')
E_mat = sp.sparse.diags([1], [-1], shape=(L, L)).tocsc()
N_mat = np.zeros((M,M))
N_mat[:,-1] = 1.0
N_mat = N_mat
Q,nod = Q_mat_and_nodes(M,typ)
A = -system_matrix_poisson1d(N)
As = [[I_L,I_M,I_N],[-I_L*dt,Q,A],[-E_mat,N_mat,I_N]]
return reduce(lambda a, b: a + b, map(lambda a_row: reduce(lambda a, b: sp.sparse.kron(a, b), a_row), As))
In [ ]:
def compound_colloc_rhs(N,M,L,k):
rhs = np.zeros(N*M*L)
rhs[0:N*M] = colloc_rhs(M,N,k)
return rhs
In [ ]:
def compound_x0(N,M,L,rhs):
x0 = rhs.copy()
for k in range(L-1):
x0[(k+1)*N*M:(k+2)*N*M] = rhs[:N*M]
return x0
In [ ]:
plot_mat_and_eigs(to_dense(compound_colloc_mat(8,3,4,0.1)))
PyAMG bietet diverse Möglichkeiten zum Lösen des Systems (https://github.com/pyamg/pyamg/wiki/Examples).
In [ ]:
# Baue Kollokationsmatrix
N = 255
M = 5
L = 10
dt = 5e-5
A = compound_colloc_mat(N,M,L,dt,"radau_r")
rhs = compound_colloc_rhs(N,M,L,16)
x0 = compound_x0(N,M,L,rhs)
sol = sprsla.spsolve(A,rhs)
plot_solution(sol,M*L)
In [ ]:
bb_solution = pyamg.solve(A,rhs,verb=False, tol=1e-11,maxiter=10)
In [ ]:
plot_solution(bb_solution-sol,M*L)
In [ ]:
mls_rs = pyamg.ruge_stuben_solver(A, max_levels=10, CF='RS',keep=True)
print mls_rs
rs_residuals = []
rs_solution = mls_rs.solve(b=rhs,x0=x0, tol=1e-10, residuals=rs_residuals )
plt.semilogy(rs_residuals)
plot_solution(rs_solution,M*L)
Aufgabe
Finden Sie experimentell heraus wo die Grenzen der AMG Löser sind, probieren Sie insbesondere verschiedene $\Delta t$ aus. Finden Sie ein Setup welches den Kollokationsansatz und den Ansatz mit dem impl. Euler vergleichbar macht. Was fällt beim Vergleich auf?
In [ ]:
# Baue Kollokationsmatrix
N = 255
M = 5
L = 10
dt = 1e-2
k=2
A = compound_colloc_mat(N,M,L,dt,"radau_r")
rhs = compound_colloc_rhs(N,M,L,k)
x0 = compound_x0(N,M,L,rhs)
A_impl_euler = # Zeile hier einfügen
dx = 1.0/(N+1)
domain = np.asarray([(i + 1) * dx for i in range(N)])
rhs_impl_euler = impl_euler_rhs(L*M,domain,k)
sol = sprsla.spsolve(A,rhs)
sol_impl_euler = sprsla.spsolve(A_impl_euler,rhs_impl_euler)
plot_solution(sol_impl_euler,M*L)
plot_solution(sol,M*L)
print np.linalg.norm(sol_impl_euler[-N:] - sol[-N:])
In [ ]:
mls_rs = pyamg.ruge_stuben_solver(A, max_levels=10, CF='RS',keep=True)
print mls_rs
rs_residuals = []
rs_solution = mls_rs.solve(b=rhs,x0=x0, tol=1e-10, residuals=rs_residuals )
plt.semilogy(rs_residuals)
plot_solution(rs_solution,M*L)
In [ ]:
mls_rs = pyamg.ruge_stuben_solver(A_impl_euler, max_levels=10, CF='RS',keep=True)
print mls_rs
rs_residuals = []
rs_solution = mls_rs.solve(b=rhs_impl_euler,x0=x0, tol=1e-10, residuals=rs_residuals )
plt.semilogy(rs_residuals)
plot_solution(rs_solution,M*L)
Nun wollen wir uns etwas genauer anschauen was AMG macht.
In [ ]:
def compound_collocation_vertices(n_space,n_comp,q_nodes,dt=1e-2,plot_it=True):
n_time = q_nodes.shape[0]
V = np.zeros((n_space*n_time*n_comp,2))
dx = 1.0/(n_space+1)
x_space = np.array([(i+1) * dx for i in range(n_space)])
t_space = q_nodes.copy() * dt
for i in range(n_comp):
for j in range(n_time):
for k in range(n_space):
V[n_time*n_space*i+n_space*j+k,0]=x_space[k]
V[n_time*n_space*i+n_space*j+k,1]=t_space[j]+i*dt
if plot_it:
plt.figure(figsize=(10,10))
plt.scatter(V[:,0], V[:,1], c='r', s=100.0)
# if not with_initial_val:
# pylab.scatter(x_space,np.zeros(n_space),c='b',s=100.0)
plt.show()
return V
def CF_plotter_classical(N,M,L,dt,typ='radau_r',symmetric=False):
Q,nodes = Q_mat_and_nodes(M,typ)
A = compound_colloc_mat(N,M,L,dt,typ)
V = compound_collocation_vertices(N,L,nodes,dt,plot_it=False)
E = np.vstack((A.tocoo().row,A.tocoo().col)).T
if symmetric:
A = A.T.dot(A)
# Use Ruge-Stuben Splitting Algorithm (use 'keep' in order to retain the splitting)
mls = pyamg.ruge_stuben_solver(A, max_levels=2, max_coarse=1, CF='RS',keep=True)
# print mls
# The CF splitting, 1 == C-node and 0 == F-node
splitting = mls.levels[0].splitting
C_nodes = splitting == 1
F_nodes = splitting == 0
plt.figure(figsize=(20,10))
lineplot(V, E,linewidths=.05)
plt.scatter(V[:,0][C_nodes], V[:,1][C_nodes], c='r', s=200.0) #plot C-nodes in red
plt.scatter(V[:,0][F_nodes], V[:,1][F_nodes], c='b', s=200.0) #plot F-nodes in blue
plt.show()
return mls
Aufgabe
Probieren Sie die Funktion CF_plotter_classical mit verschiedenen Parametern aus.
Was fällt im Vergleich zum impliziten Euler Vorgehen auf?
In [ ]:
mls=CF_plotter_classical(7,5,2,1e1,typ='radau_r')
print mls
In der Vorlesung haben wir gelernt, dass die Theorie zu AMG auf sog. $\mathbf{M}$ Matrizen aufgebaut ist. Das Kollokationsproblem erfüllt die meisten Eigenschaften davon nicht. Zudem sehen wir die schlechten Konvergenzeigenschaften von AMG für das kombinierte Kollokationsproblem.
Aufgabe
Symmetrisieren Sie das Problem und wenden Sie SmoothedAggregation und Classical AMG an.
In [ ]:
# Baue Kollokationsmatrix
N = 255
M = 5
L = 4
dt = 1e-4
A = compound_colloc_mat(N,M,L,dt,"radau_r")
rhs = compound_colloc_rhs(N,M,L,16)
x0 = compound_x0(N,M,L,rhs)
A_sym = # Zeile hier einfügen
rhs_sym = # Zeile hier einfügen
x0_sym = # Zeile hier einfügen
sol = sprsla.spsolve(A_sym,rhs_sym)
plot_solution(sol,M*L)
In [ ]:
# Count of the non-zero entries
print A_sym.nnz, A.nnz, (N*M*L)**2
In [ ]:
mls_sa = pyamg.ruge_stuben_solver(A_sym, max_levels=5, max_coarse=1, CF='RS',keep=True)
print mls_sa
sa_residuals = []
sa_solution = mls_sa.solve(b=rhs_sym,x0=x0_sym, tol=1e-10, residuals=sa_residuals )
plt.semilogy(sa_residuals)
plot_solution(sa_solution, M*L)
In [ ]:
mls_sa = pyamg.ruge_stuben_solver(A, max_levels=5, max_coarse=1, CF='RS',keep=True)
print mls_sa
sa_residuals = []
sa_solution = mls_sa.solve(b=rhs,x0=x0, tol=1e-10, residuals=sa_residuals )
plt.semilogy(sa_residuals)
plot_solution(sa_solution, M*L)
Aufgabe
Vergleichen Sie die Wahl der Grobgitterpunkte zwischen dem Kollokationionsproblem und seiner symmetrisierten Version.
In [ ]:
N = 7
M = 5
L = 2
dt =1e-3
In [ ]:
mls=CF_plotter_classical(N,M,L,dt)
print mls
In [ ]:
mls=CF_plotter_classical(N,M,L,dt,symmetric=True)
print mls