|
| 1 | +import numpy as np |
| 2 | +import scipy.sparse as sp |
| 3 | +import scipy.sparse.linalg as spl |
| 4 | +import matplotlib.pyplot as plt |
| 5 | + |
| 6 | +import pySDC.helpers.transfer_helper as th |
| 7 | + |
| 8 | +np.random.seed(32) |
| 9 | + |
| 10 | +nX = 2 ** 3 - 1 |
| 11 | +nXc = 2 ** 2 - 1 |
| 12 | +dx = 1/(nX + 1) |
| 13 | +dxc = 1/(nXc + 1) |
| 14 | + |
| 15 | +stencil = [1, -2, 1] |
| 16 | +A = sp.diags(stencil, [-1, 0, 1], shape=(nX, nX), format='csc') |
| 17 | +A *= - 1.0 / (dx ** 2) |
| 18 | +b = np.zeros(nX) |
| 19 | +# %% |
| 20 | +PGS = spl.inv(sp.csc_matrix(np.tril(A.todense()))) |
| 21 | +PJ = sp.dia_matrix(np.diag(1/np.diag(A.todense()))) |
| 22 | + |
| 23 | +RGS = sp.eye(nX) - PGS @ A |
| 24 | +RJ = sp.eye(nX) - PJ @ A |
| 25 | +sRJ = sp.eye(nX) - 0.5 * PJ @ A |
| 26 | + |
| 27 | +Ac = sp.diags(stencil, [-1, 0, 1], shape=(nX//2, nX//2), format='csc') |
| 28 | +Ac *= - 1.0 / (dxc ** 2) |
| 29 | +Acinv = spl.inv(Ac) |
| 30 | +PJc = sp.dia_matrix(np.diag(1/np.diag(Ac.todense()))) |
| 31 | + |
| 32 | +fine_grid = np.linspace(dx, 1, nX, endpoint=False) |
| 33 | +coarse_grid = np.linspace(dxc, 1, nXc, endpoint=False) |
| 34 | +Pr = th.interpolation_matrix_1d(fine_grid, coarse_grid, k=2, periodic=False, equidist_nested=True) |
| 35 | +Prinj = th.interpolation_matrix_1d(fine_grid, coarse_grid, k=0, periodic=False, equidist_nested=True) |
| 36 | +Re = 0.5 * Pr.T |
| 37 | +# Re = Prinj.T |
| 38 | + |
| 39 | +TG = (sp.eye(nX) - Pr @ Acinv @ Re @ A) @ sRJ.power(3) |
| 40 | +TJ = (sp.eye(nX) - Pr @ PJc @ Re @ A) @ sRJ.power(3) |
| 41 | + |
| 42 | +# %% |
| 43 | + |
| 44 | +nIter = 100 |
| 45 | +dK = 10 |
| 46 | + |
| 47 | +def coarse(x, dK=dK): |
| 48 | + xK = x |
| 49 | + for k in range(dK): |
| 50 | + xK = RJ @ xK #+ RJ @ b |
| 51 | + return xK |
| 52 | + |
| 53 | + |
| 54 | +def fine(x, dK=dK): |
| 55 | + xK = x |
| 56 | + for k in range(dK): |
| 57 | + xK = TG @ xK #+ PGS @ b |
| 58 | + return xK |
| 59 | + |
| 60 | +def initBlocks(x): |
| 61 | + x0 = np.sin(np.linspace(dx, 1*2*np.pi, nX)) |
| 62 | + # x0 = np.random.rand(nX) |
| 63 | + for l in range(nB+1): |
| 64 | + x[l, :] = x0 |
| 65 | + |
| 66 | +nB = nIter//dK |
| 67 | +nK = nB+1 |
| 68 | + |
| 69 | + |
| 70 | +uSeq = np.zeros((nB+1, nX)) |
| 71 | +uNum = np.zeros((nK+1, nB+1, nX)) |
| 72 | + |
| 73 | +initBlocks(uNum[0]) |
| 74 | + |
| 75 | +uSeq[0] = uNum[0, 0] |
| 76 | +uNum[:, 0, :] = uNum[0, 0, :] |
| 77 | + |
| 78 | +for k in range(nK): |
| 79 | + for l in range(nB): |
| 80 | + |
| 81 | + uF = fine(uNum[k, l]) |
| 82 | + uGk = coarse(uNum[k, l], dK=1) |
| 83 | + uGkp1 = coarse(uNum[k+1, l], dK=1) |
| 84 | + |
| 85 | + uNum[k+1, l+1] = uF + 1*(uGkp1 - uGk) |
| 86 | + |
| 87 | + |
| 88 | +for l in range(nB): |
| 89 | + uSeq[l+1] = fine(uSeq[l]) |
| 90 | + |
| 91 | +# %% |
| 92 | + |
| 93 | +err = np.linalg.norm(uNum, axis=-1)[:, -1] |
| 94 | +errSeq = np.linalg.norm(uSeq, axis=-1) |
| 95 | + |
| 96 | +plt.figure() |
| 97 | +plt.semilogy(err, label='Parareal') |
| 98 | +plt.semilogy(errSeq, label='Sequential') |
| 99 | +plt.legend() |
| 100 | +plt.show() |
0 commit comments