MG Cycle

.gitignore

.*.swp

fig1_x0.png

fig1_x1.png

fig1_x2.png

mg_cycle.py

import numpy as np
import scipy as sp
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.mlab import bivariate_normal
import matplotlib.pyplot as plt
import matplotlib.colors as col

import pyamg

plt.ion()

startcolor = '#AA0000'  # a dark olive 
midcolor = '#00AA00'    # a bright yellow
endcolor = '#0000AA'    # medium dark red
cmap2 = col.LinearSegmentedColormap.from_list('own2',[startcolor,midcolor,endcolor])
plt.cm.register_cmap(cmap=cmap2)

def plotit(X, Y, U, s):
    nx = X.shape[0]
    fig = plt.figure()
    ax = fig.gca(projection='3d')
    ax.plot_surface(X, Y, U.reshape((nx, nx)), rstride=1, cstride=1,
                    cmap=cmap2,
                    linewidth=1,
                    edgecolor='white',
                    )
    ax.set_xlim(-3, 3)
    ax.set_ylim(-3, 3)
    ax.set_zlim(-3, 3)
    ax.axis('off')
    f = 'fig1_%s.png' % s
    plt.savefig(f, bbox_inches='tight')
    #f = 'fig1_%s.pdf' % plt.gcf().number
    #plt.savefig(f, bbox_inches='tight')
    #import os
    #os.system('pdfcrop %s %s'%(f,f))
nx = 24
nxc = 12

stencil = [[-1, -1, -1], [-1, 8, -1], [-1, -1, -1.]]
A = pyamg.gallery.stencil_grid(stencil, (nx, nx), format='csr')
#A = pyamg.gallery.poisson((nx, nx), format='csr')
n = A.shape[0]

x = y = np.arange(-3.0, 3.0, 6.0 / nx)
X, Y = np.meshgrid(x, y)

Z1 = bivariate_normal(X, Y, 1.0, 1.0, 0.0, 0.0)
Z2 = bivariate_normal(X, Y, 1.5, 0.5, 1, 1)
Z = Z2 - Z1
X = X
Y = Y
Z = 10 * Z

u = Z.ravel()
pert = np.random.random(n)
u0 = u + pert
b = A * u

u1 = u0.copy()
print(np.linalg.norm(b - A * u0))
pyamg.relaxation.relaxation.jacobi(A, u1, b, iterations=2, omega=2.0/3.0)
print(np.linalg.norm(b - A * u1))

ml = pyamg.ruge_stuben_solver(A, strength=('classical', {'theta': 0.0}))
P = ml.levels[0].P
Ac = ml.levels[1].A

e1 = u - u1

xc = yc = np.arange(-3.0, 3.0, 6.0 / nxc)
Xc, Yc = np.meshgrid(xc, yc)
e2 = P.T * e1
e3 = 0. * e2.copy()
rc = P.T * (b - A * u1)

print(np.linalg.norm(rc - Ac * e3))
pyamg.relaxation.relaxation.jacobi(Ac, e3, rc, iterations=2, omega=2.0/3.0)
print(np.linalg.norm(rc - Ac * e3))

e4 = P * e3
u5 = u1 + e4
u6 = u5.copy()
print(np.linalg.norm(b - A * u6))
pyamg.relaxation.relaxation.jacobi(A, u6, b, iterations=2, omega=2.0/3.0)
print(np.linalg.norm(b - A * u6))

plotit(X, Y, u0, 'x0')
plotit(X, Y, u1, 'x1')
#plotit(X, Y, e1, 'e1')
#plotit(Xc, Yc, e2, 'e1c')
#plotit(Xc, Yc, e3, 'e1c')
plotit(X, Y, u5, 'x1correction')
plotit(X, Y, u6, 'x2')
plt.show()

This creates figures to make the following