Multigrid examples¶
[1]:
%matplotlib inline
import matplotlib.pyplot as plt
[2]:
from __future__ import print_function
import numpy as np
import mesh.boundary as bnd
import mesh.patch as patch
import multigrid.MG as MG
Constant-coefficent Poisson equation¶
We want to solve
\[\phi_{xx} + \phi_{yy} = -2[(1-6x^2)y^2(1-y^2) + (1-6y^2)x^2(1-x^2)]\]
on
\[[0,1]\times [0,1]\]
with homogeneous Dirichlet boundary conditions (this example comes from “A Multigrid Tutorial”).
This has the analytic solution
\[u(x,y) = (x^2 - x^4)(y^4 - y^2)\]
We start by setting up a multigrid object–this needs to know the number of zones our problem is defined on
[3]:
nx = ny = 256
mg = MG.CellCenterMG2d(nx, ny,
xl_BC_type="dirichlet", xr_BC_type="dirichlet",
yl_BC_type="dirichlet", yr_BC_type="dirichlet", verbose=1)
cc data: nx = 2, ny = 2, ng = 1
nvars = 3
variables:
v: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
f: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
r: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
cc data: nx = 4, ny = 4, ng = 1
nvars = 3
variables:
v: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
f: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
r: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
cc data: nx = 8, ny = 8, ng = 1
nvars = 3
variables:
v: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
f: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
r: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
cc data: nx = 16, ny = 16, ng = 1
nvars = 3
variables:
v: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
f: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
r: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
cc data: nx = 32, ny = 32, ng = 1
nvars = 3
variables:
v: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
f: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
r: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
cc data: nx = 64, ny = 64, ng = 1
nvars = 3
variables:
v: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
f: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
r: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
cc data: nx = 128, ny = 128, ng = 1
nvars = 3
variables:
v: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
f: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
r: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
cc data: nx = 256, ny = 256, ng = 1
nvars = 3
variables:
v: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
f: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
r: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
Next, we initialize the RHS. To make life easier, the CellCenterMG2d object has the coordinates of the solution grid (including ghost cells) as mg.x2d and mg.y2d (these are two-dimensional arrays).
[4]:
def rhs(x, y):
return -2.0*((1.0-6.0*x**2)*y**2*(1.0-y**2) + (1.0-6.0*y**2)*x**2*(1.0-x**2))
mg.init_RHS(rhs(mg.x2d, mg.y2d))
Source norm = 1.097515813669473
The last setup step is to initialize the solution–this is the starting point for the solve. Usually we just want to start with all zeros, so we use the init_zeros() method
[5]:
mg.init_zeros()
we can now solve – there are actually two different techniques we can do here. We can just do pure smoothing on the solution grid using mg.smooth(mg.nlevels-1, N), where N is the number of smoothing iterations. To get the solution N will need to be large and this will take a long time.
Multigrid accelerates the smoothing. We can do a V-cycle multigrid solution using mg.solve()
[6]:
mg.solve()
source norm = 1.097515813669473
<<< beginning V-cycle (cycle 1) >>>
level: 7, grid: 256 x 256
before G-S, residual L2: 1.097515813669473
after G-S, residual L2: 1.502308451578657
level: 6, grid: 128 x 128
before G-S, residual L2: 1.0616243965458263
after G-S, residual L2: 1.4321452257629033
level: 5, grid: 64 x 64
before G-S, residual L2: 1.011366277976364
after G-S, residual L2: 1.281872470375375
level: 4, grid: 32 x 32
before G-S, residual L2: 0.903531158162907
after G-S, residual L2: 0.9607576999783505
level: 3, grid: 16 x 16
before G-S, residual L2: 0.6736112182020367
after G-S, residual L2: 0.4439774050299674
level: 2, grid: 8 x 8
before G-S, residual L2: 0.30721142286171554
after G-S, residual L2: 0.0727215591269748
level: 1, grid: 4 x 4
before G-S, residual L2: 0.04841813458618458
after G-S, residual L2: 3.9610700301811246e-05
bottom solve:
level: 0, grid: 2 x 2
level: 1, grid: 4 x 4
before G-S, residual L2: 3.925006722484123e-05
after G-S, residual L2: 1.0370099820862674e-09
level: 2, grid: 8 x 8
before G-S, residual L2: 0.07010129273961899
after G-S, residual L2: 0.0008815704830693547
level: 3, grid: 16 x 16
before G-S, residual L2: 0.4307377377402105
after G-S, residual L2: 0.007174899576794818
level: 4, grid: 32 x 32
before G-S, residual L2: 0.911086486792154
after G-S, residual L2: 0.01618756602227813
level: 5, grid: 64 x 64
before G-S, residual L2: 1.1945438349788615
after G-S, residual L2: 0.022021327892004925
level: 6, grid: 128 x 128
before G-S, residual L2: 1.313456560108626
after G-S, residual L2: 0.02518650395173617
level: 7, grid: 256 x 256
before G-S, residual L2: 1.3618314516335004
after G-S, residual L2: 0.026870007568672097
cycle 1: relative err = 0.999999999999964, residual err = 0.02448256984911586
<<< beginning V-cycle (cycle 2) >>>
level: 7, grid: 256 x 256
before G-S, residual L2: 0.026870007568672097
after G-S, residual L2: 0.025790216249923552
level: 6, grid: 128 x 128
before G-S, residual L2: 0.018218080204017304
after G-S, residual L2: 0.023654310121915274
level: 5, grid: 64 x 64
before G-S, residual L2: 0.01669077991582338
after G-S, residual L2: 0.01977335201785163
level: 4, grid: 32 x 32
before G-S, residual L2: 0.013922595404814862
after G-S, residual L2: 0.013577568890182053
level: 3, grid: 16 x 16
before G-S, residual L2: 0.009518306167970536
after G-S, residual L2: 0.006115159484497302
level: 2, grid: 8 x 8
before G-S, residual L2: 0.004244630812032651
after G-S, residual L2: 0.0010674120586864006
level: 1, grid: 4 x 4
before G-S, residual L2: 0.0007108144252738053
after G-S, residual L2: 5.818246254772703e-07
bottom solve:
level: 0, grid: 2 x 2
level: 1, grid: 4 x 4
before G-S, residual L2: 5.765281065294632e-07
after G-S, residual L2: 1.5231212503339452e-11
level: 2, grid: 8 x 8
before G-S, residual L2: 0.0010291471590693868
after G-S, residual L2: 1.2950948742201083e-05
level: 3, grid: 16 x 16
before G-S, residual L2: 0.006239446983842889
after G-S, residual L2: 0.00010483463130232172
level: 4, grid: 32 x 32
before G-S, residual L2: 0.014573363314854
after G-S, residual L2: 0.00026233988398787004
level: 5, grid: 64 x 64
before G-S, residual L2: 0.021564270263952755
after G-S, residual L2: 0.0003944827058086955
level: 6, grid: 128 x 128
before G-S, residual L2: 0.02579092712136628
after G-S, residual L2: 0.00048636495715121916
level: 7, grid: 256 x 256
before G-S, residual L2: 0.028051324215592862
after G-S, residual L2: 0.0005440874957950154
cycle 2: relative err = 13.739483825281054, residual err = 0.0004957445615074047
<<< beginning V-cycle (cycle 3) >>>
level: 7, grid: 256 x 256
before G-S, residual L2: 0.0005440874957950154
after G-S, residual L2: 0.0005095844930046698
level: 6, grid: 128 x 128
before G-S, residual L2: 0.0003597879816772893
after G-S, residual L2: 0.00044648485218937167
level: 5, grid: 64 x 64
before G-S, residual L2: 0.0003147892995472901
after G-S, residual L2: 0.0003492541721056348
level: 4, grid: 32 x 32
before G-S, residual L2: 0.0002457276904804801
after G-S, residual L2: 0.00022232862524233384
level: 3, grid: 16 x 16
before G-S, residual L2: 0.0001558932199490972
after G-S, residual L2: 9.511093023364078e-05
level: 2, grid: 8 x 8
before G-S, residual L2: 6.616899520585456e-05
after G-S, residual L2: 1.711006102346096e-05
level: 1, grid: 4 x 4
before G-S, residual L2: 1.139522901981679e-05
after G-S, residual L2: 9.33004809910226e-09
bottom solve:
level: 0, grid: 2 x 2
level: 1, grid: 4 x 4
before G-S, residual L2: 9.245125097272049e-09
after G-S, residual L2: 2.442311694447821e-13
level: 2, grid: 8 x 8
before G-S, residual L2: 1.64991725637487e-05
after G-S, residual L2: 2.0771258971860784e-07
level: 3, grid: 16 x 16
before G-S, residual L2: 0.00010097720436460624
after G-S, residual L2: 1.7241727900979902e-06
level: 4, grid: 32 x 32
before G-S, residual L2: 0.0002575410544503153
after G-S, residual L2: 4.766282851613449e-06
level: 5, grid: 64 x 64
before G-S, residual L2: 0.00041133882050328275
after G-S, residual L2: 7.600616845344458e-06
level: 6, grid: 128 x 128
before G-S, residual L2: 0.0005232809692242086
after G-S, residual L2: 9.860758095018993e-06
level: 7, grid: 256 x 256
before G-S, residual L2: 0.0005945070122423073
after G-S, residual L2: 1.1466134915427874e-05
cycle 3: relative err = 34.347638624909216, residual err = 1.0447352805871284e-05
<<< beginning V-cycle (cycle 4) >>>
level: 7, grid: 256 x 256
before G-S, residual L2: 1.1466134915427874e-05
after G-S, residual L2: 1.054466722279011e-05
level: 6, grid: 128 x 128
before G-S, residual L2: 7.442814693866286e-06
after G-S, residual L2: 8.955050475722099e-06
level: 5, grid: 64 x 64
before G-S, residual L2: 6.311313968968047e-06
after G-S, residual L2: 6.734553609148436e-06
level: 4, grid: 32 x 32
before G-S, residual L2: 4.737984987500691e-06
after G-S, residual L2: 4.091799997658277e-06
level: 3, grid: 16 x 16
before G-S, residual L2: 2.871028473858937e-06
after G-S, residual L2: 1.6319551993366253e-06
level: 2, grid: 8 x 8
before G-S, residual L2: 1.1372178077508109e-06
after G-S, residual L2: 2.961040430099916e-07
level: 1, grid: 4 x 4
before G-S, residual L2: 1.9721864323458624e-07
after G-S, residual L2: 1.61503943872384e-10
bottom solve:
level: 0, grid: 2 x 2
level: 1, grid: 4 x 4
before G-S, residual L2: 1.6003411195777404e-10
after G-S, residual L2: 4.2274326344473505e-15
level: 2, grid: 8 x 8
before G-S, residual L2: 2.855691101825338e-07
after G-S, residual L2: 3.5961118754371857e-09
level: 3, grid: 16 x 16
before G-S, residual L2: 1.7893831203170535e-06
after G-S, residual L2: 3.1136282101831173e-08
level: 4, grid: 32 x 32
before G-S, residual L2: 4.97129807196115e-06
after G-S, residual L2: 9.544819739422644e-08
level: 5, grid: 64 x 64
before G-S, residual L2: 8.281644276572538e-06
after G-S, residual L2: 1.56637783149839e-07
level: 6, grid: 128 x 128
before G-S, residual L2: 1.0888850082357996e-05
after G-S, residual L2: 2.0777271327080248e-07
level: 7, grid: 256 x 256
before G-S, residual L2: 1.2717522622400765e-05
after G-S, residual L2: 2.464531349025277e-07
cycle 4: relative err = 0.17409776671446628, residual err = 2.24555429482631e-07
<<< beginning V-cycle (cycle 5) >>>
level: 7, grid: 256 x 256
before G-S, residual L2: 2.464531349025277e-07
after G-S, residual L2: 2.2491138140311698e-07
level: 6, grid: 128 x 128
before G-S, residual L2: 1.5874562191875262e-07
after G-S, residual L2: 1.886249099391391e-07
level: 5, grid: 64 x 64
before G-S, residual L2: 1.3294481979637655e-07
after G-S, residual L2: 1.397710191717015e-07
level: 4, grid: 32 x 32
before G-S, residual L2: 9.836928907527788e-08
after G-S, residual L2: 8.269030961692836e-08
level: 3, grid: 16 x 16
before G-S, residual L2: 5.8062531341283565e-08
after G-S, residual L2: 3.034725896415429e-08
level: 2, grid: 8 x 8
before G-S, residual L2: 2.116912379336852e-08
after G-S, residual L2: 5.467519592468213e-09
level: 1, grid: 4 x 4
before G-S, residual L2: 3.6418116003284676e-09
after G-S, residual L2: 2.982625229812215e-12
bottom solve:
level: 0, grid: 2 x 2
level: 1, grid: 4 x 4
before G-S, residual L2: 2.955484162036181e-12
after G-S, residual L2: 7.806739482450516e-17
level: 2, grid: 8 x 8
before G-S, residual L2: 5.273610709946236e-09
after G-S, residual L2: 6.642323465658688e-11
level: 3, grid: 16 x 16
before G-S, residual L2: 3.4146989205844565e-08
after G-S, residual L2: 6.052228076583688e-10
level: 4, grid: 32 x 32
before G-S, residual L2: 1.031248597196911e-07
after G-S, residual L2: 2.0541497445308587e-09
level: 5, grid: 64 x 64
before G-S, residual L2: 1.7585349306604133e-07
after G-S, residual L2: 3.421022608879089e-09
level: 6, grid: 128 x 128
before G-S, residual L2: 2.3383756442516674e-07
after G-S, residual L2: 4.552170797983864e-09
level: 7, grid: 256 x 256
before G-S, residual L2: 2.7592842790687426e-07
after G-S, residual L2: 5.41488950707315e-09
cycle 5: relative err = 0.005391244339065405, residual err = 4.933769007818501e-09
<<< beginning V-cycle (cycle 6) >>>
level: 7, grid: 256 x 256
before G-S, residual L2: 5.41488950707315e-09
after G-S, residual L2: 4.948141362729419e-09
level: 6, grid: 128 x 128
before G-S, residual L2: 3.4929583962703016e-09
after G-S, residual L2: 4.154445183511443e-09
level: 5, grid: 64 x 64
before G-S, residual L2: 2.9288841397931397e-09
after G-S, residual L2: 3.074779198797186e-09
level: 4, grid: 32 x 32
before G-S, residual L2: 2.164991235492634e-09
after G-S, residual L2: 1.788028730183651e-09
level: 3, grid: 16 x 16
before G-S, residual L2: 1.2562223343389894e-09
after G-S, residual L2: 6.021983813990021e-10
level: 2, grid: 8 x 8
before G-S, residual L2: 4.2028073688787063e-10
after G-S, residual L2: 1.0655724637281067e-10
level: 1, grid: 4 x 4
before G-S, residual L2: 7.097871736854444e-11
after G-S, residual L2: 5.813506543301849e-14
bottom solve:
level: 0, grid: 2 x 2
level: 1, grid: 4 x 4
before G-S, residual L2: 5.760611936011378e-14
after G-S, residual L2: 1.521555112430923e-18
level: 2, grid: 8 x 8
before G-S, residual L2: 1.027891920456506e-10
after G-S, residual L2: 1.294879454701896e-12
level: 3, grid: 16 x 16
before G-S, residual L2: 6.914011940773812e-10
after G-S, residual L2: 1.2453691230551983e-11
level: 4, grid: 32 x 32
before G-S, residual L2: 2.2570491487662195e-09
after G-S, residual L2: 4.639035392364569e-11
level: 5, grid: 64 x 64
before G-S, residual L2: 3.908967396962745e-09
after G-S, residual L2: 7.803740782474827e-11
level: 6, grid: 128 x 128
before G-S, residual L2: 5.196394306272565e-09
after G-S, residual L2: 1.033274523100204e-10
level: 7, grid: 256 x 256
before G-S, residual L2: 6.117636729623554e-09
after G-S, residual L2: 1.2199402602477584e-10
cycle 6: relative err = 7.51413991329132e-05, residual err = 1.111546863428753e-10
<<< beginning V-cycle (cycle 7) >>>
level: 7, grid: 256 x 256
before G-S, residual L2: 1.2199402602477584e-10
after G-S, residual L2: 1.121992266879251e-10
level: 6, grid: 128 x 128
before G-S, residual L2: 7.921861122082639e-11
after G-S, residual L2: 9.493449600138316e-11
level: 5, grid: 64 x 64
before G-S, residual L2: 6.694993398453784e-11
after G-S, residual L2: 7.050995614737483e-11
level: 4, grid: 32 x 32
before G-S, residual L2: 4.9666563586565975e-11
after G-S, residual L2: 4.045094776680348e-11
level: 3, grid: 16 x 16
before G-S, residual L2: 2.843147343834713e-11
after G-S, residual L2: 1.2576313722677801e-11
level: 2, grid: 8 x 8
before G-S, residual L2: 8.777954081387978e-12
after G-S, residual L2: 2.170559196862902e-12
level: 1, grid: 4 x 4
before G-S, residual L2: 1.445876195415056e-12
after G-S, residual L2: 1.1842925278593641e-15
bottom solve:
level: 0, grid: 2 x 2
level: 1, grid: 4 x 4
before G-S, residual L2: 1.1735184729034125e-15
after G-S, residual L2: 3.0994757710835167e-20
level: 2, grid: 8 x 8
before G-S, residual L2: 2.094012660676073e-12
after G-S, residual L2: 2.6382579574150587e-14
level: 3, grid: 16 x 16
before G-S, residual L2: 1.466147487151147e-11
after G-S, residual L2: 2.6760553592700965e-13
level: 4, grid: 32 x 32
before G-S, residual L2: 5.130705216489902e-11
after G-S, residual L2: 1.0810419626613159e-12
level: 5, grid: 64 x 64
before G-S, residual L2: 9.001551103280705e-11
after G-S, residual L2: 1.8342879121275396e-12
level: 6, grid: 128 x 128
before G-S, residual L2: 1.1914921193827463e-10
after G-S, residual L2: 2.4124327865487605e-12
level: 7, grid: 256 x 256
before G-S, residual L2: 1.3907209384461257e-10
after G-S, residual L2: 2.8429898342353533e-12
cycle 7: relative err = 7.062255558417692e-07, residual err = 2.590386214782638e-12
We can access the solution on the finest grid using get_solution()
[7]:
phi = mg.get_solution()
[8]:
plt.imshow(np.transpose(phi.v()), origin="lower")
[8]:
<matplotlib.image.AxesImage at 0x7efe36035070>
we can also get the gradient of the solution
[9]:
gx, gy = mg.get_solution_gradient()
[10]:
plt.subplot(121)
plt.imshow(np.transpose(gx.v()), origin="lower")
plt.subplot(122)
plt.imshow(np.transpose(gy.v()), origin="lower")
[10]:
<matplotlib.image.AxesImage at 0x7efe2ddb6ca0>
General linear elliptic equation¶
The GeneralMG2d class implements support for a general elliptic equation of the form:
\[\alpha \phi + \nabla \cdot (\beta \nabla \phi) + \gamma \cdot \nabla \phi = f\]
with inhomogeneous boundary condtions.
It subclasses the CellCenterMG2d class, and the basic interface is the same
We will solve the above with
\begin{align}
\alpha &= 10 \\
\beta &= xy + 1 \\
\gamma &= \hat{x} + \hat{y}
\end{align}
and
\begin{align}
f = &-\frac{\pi}{2}(x + 1)\sin\left(\frac{\pi y}{2}\right) \cos\left(\frac{\pi x}{2}\right ) \\
&-\frac{\pi}{2}(y + 1)\sin\left(\frac{\pi x}{2}\right) \cos\left(\frac{\pi y}{2}\right ) \\
&+\left(\frac{-\pi^2 (xy+1)}{2} + 10\right) \cos\left(\frac{\pi x}{2}\right) \cos\left(\frac{\pi y}{2}\right)
\end{align}
on \([0, 1] \times [0,1]\) with boundary conditions:
\begin{align}
\phi(x=0) &= \cos(\pi y/2) \\
\phi(x=1) &= 0 \\
\phi(y=0) &= \cos(\pi x/2) \\
\phi(y=1) &= 0
\end{align}
This has the exact solution:
\[\phi = \cos(\pi x/2) \cos(\pi y/2)\]
[11]:
import multigrid.general_MG as gMG
For reference, we’ll define a function providing the analytic solution
[12]:
def true(x,y):
return np.cos(np.pi*x/2.0)*np.cos(np.pi*y/2.0)
Now the coefficents–note that since \(\gamma\) is a vector, we have a different function for each component
[13]:
def alpha(x,y):
return 10.0*np.ones_like(x)
def beta(x,y):
return x*y + 1.0
def gamma_x(x,y):
return np.ones_like(x)
def gamma_y(x,y):
return np.ones_like(x)
and the righthand side function
[14]:
def f(x,y):
return -0.5*np.pi*(x + 1.0)*np.sin(np.pi*y/2.0)*np.cos(np.pi*x/2.0) - \
0.5*np.pi*(y + 1.0)*np.sin(np.pi*x/2.0)*np.cos(np.pi*y/2.0) + \
(-np.pi**2*(x*y+1.0)/2.0 + 10.0)*np.cos(np.pi*x/2.0)*np.cos(np.pi*y/2.0)
Our inhomogeneous boundary conditions require a function that can be evaluated on the boundary to give the value
[15]:
def xl_func(y):
return np.cos(np.pi*y/2.0)
def yl_func(x):
return np.cos(np.pi*x/2.0)
Now we can setup our grid object and the coefficients, which are stored as a CellCenter2d object. Note, the coefficients do not need to have the same boundary conditions as \(\phi\) (and for real problems, they may not). The one that matters the most is \(\beta\), since that will need to be averaged to the edges of the domain, so the boundary conditions on the coefficients are important.
Here we use Neumann boundary conditions
[16]:
import mesh.patch as patch
nx = ny = 128
g = patch.Grid2d(nx, ny, ng=1)
d = patch.CellCenterData2d(g)
bc_c = bnd.BC(xlb="neumann", xrb="neumann",
ylb="neumann", yrb="neumann")
d.register_var("alpha", bc_c)
d.register_var("beta", bc_c)
d.register_var("gamma_x", bc_c)
d.register_var("gamma_y", bc_c)
d.create()
a = d.get_var("alpha")
a[:,:] = alpha(g.x2d, g.y2d)
b = d.get_var("beta")
b[:,:] = beta(g.x2d, g.y2d)
gx = d.get_var("gamma_x")
gx[:,:] = gamma_x(g.x2d, g.y2d)
gy = d.get_var("gamma_y")
gy[:,:] = gamma_y(g.x2d, g.y2d)
Now we can setup the multigrid object
[17]:
a = gMG.GeneralMG2d(nx, ny,
xl_BC_type="dirichlet", yl_BC_type="dirichlet",
xr_BC_type="dirichlet", yr_BC_type="dirichlet",
xl_BC=xl_func,
yl_BC=yl_func,
coeffs=d,
verbose=1, vis=0, true_function=true)
cc data: nx = 2, ny = 2, ng = 1
nvars = 7
variables:
v: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
f: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
r: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
alpha: min: 0.0000000000 max: 0.0000000000
BCs: -x: neumann +x: neumann -y: neumann +y: neumann
beta: min: 0.0000000000 max: 0.0000000000
BCs: -x: neumann +x: neumann -y: neumann +y: neumann
gamma_x: min: 0.0000000000 max: 0.0000000000
BCs: -x: neumann +x: neumann -y: neumann +y: neumann
gamma_y: min: 0.0000000000 max: 0.0000000000
BCs: -x: neumann +x: neumann -y: neumann +y: neumann
cc data: nx = 4, ny = 4, ng = 1
nvars = 7
variables:
v: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
f: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
r: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
alpha: min: 0.0000000000 max: 0.0000000000
BCs: -x: neumann +x: neumann -y: neumann +y: neumann
beta: min: 0.0000000000 max: 0.0000000000
BCs: -x: neumann +x: neumann -y: neumann +y: neumann
gamma_x: min: 0.0000000000 max: 0.0000000000
BCs: -x: neumann +x: neumann -y: neumann +y: neumann
gamma_y: min: 0.0000000000 max: 0.0000000000
BCs: -x: neumann +x: neumann -y: neumann +y: neumann
cc data: nx = 8, ny = 8, ng = 1
nvars = 7
variables:
v: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
f: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
r: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
alpha: min: 0.0000000000 max: 0.0000000000
BCs: -x: neumann +x: neumann -y: neumann +y: neumann
beta: min: 0.0000000000 max: 0.0000000000
BCs: -x: neumann +x: neumann -y: neumann +y: neumann
gamma_x: min: 0.0000000000 max: 0.0000000000
BCs: -x: neumann +x: neumann -y: neumann +y: neumann
gamma_y: min: 0.0000000000 max: 0.0000000000
BCs: -x: neumann +x: neumann -y: neumann +y: neumann
cc data: nx = 16, ny = 16, ng = 1
nvars = 7
variables:
v: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
f: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
r: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
alpha: min: 0.0000000000 max: 0.0000000000
BCs: -x: neumann +x: neumann -y: neumann +y: neumann
beta: min: 0.0000000000 max: 0.0000000000
BCs: -x: neumann +x: neumann -y: neumann +y: neumann
gamma_x: min: 0.0000000000 max: 0.0000000000
BCs: -x: neumann +x: neumann -y: neumann +y: neumann
gamma_y: min: 0.0000000000 max: 0.0000000000
BCs: -x: neumann +x: neumann -y: neumann +y: neumann
cc data: nx = 32, ny = 32, ng = 1
nvars = 7
variables:
v: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
f: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
r: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
alpha: min: 0.0000000000 max: 0.0000000000
BCs: -x: neumann +x: neumann -y: neumann +y: neumann
beta: min: 0.0000000000 max: 0.0000000000
BCs: -x: neumann +x: neumann -y: neumann +y: neumann
gamma_x: min: 0.0000000000 max: 0.0000000000
BCs: -x: neumann +x: neumann -y: neumann +y: neumann
gamma_y: min: 0.0000000000 max: 0.0000000000
BCs: -x: neumann +x: neumann -y: neumann +y: neumann
cc data: nx = 64, ny = 64, ng = 1
nvars = 7
variables:
v: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
f: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
r: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
alpha: min: 0.0000000000 max: 0.0000000000
BCs: -x: neumann +x: neumann -y: neumann +y: neumann
beta: min: 0.0000000000 max: 0.0000000000
BCs: -x: neumann +x: neumann -y: neumann +y: neumann
gamma_x: min: 0.0000000000 max: 0.0000000000
BCs: -x: neumann +x: neumann -y: neumann +y: neumann
gamma_y: min: 0.0000000000 max: 0.0000000000
BCs: -x: neumann +x: neumann -y: neumann +y: neumann
cc data: nx = 128, ny = 128, ng = 1
nvars = 7
variables:
v: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
f: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
r: min: 0.0000000000 max: 0.0000000000
BCs: -x: dirichlet +x: dirichlet -y: dirichlet +y: dirichlet
alpha: min: 0.0000000000 max: 0.0000000000
BCs: -x: neumann +x: neumann -y: neumann +y: neumann
beta: min: 0.0000000000 max: 0.0000000000
BCs: -x: neumann +x: neumann -y: neumann +y: neumann
gamma_x: min: 0.0000000000 max: 0.0000000000
BCs: -x: neumann +x: neumann -y: neumann +y: neumann
gamma_y: min: 0.0000000000 max: 0.0000000000
BCs: -x: neumann +x: neumann -y: neumann +y: neumann
just as before, we specify the righthand side and initialize the solution
[18]:
a.init_zeros()
a.init_RHS(f(a.x2d, a.y2d))
Source norm = 1.775181492337501
and we can solve it
[19]:
a.solve(rtol=1.e-10)
source norm = 1.775181492337501
<<< beginning V-cycle (cycle 1) >>>
level: 6, grid: 128 x 128
before G-S, residual L2: 1.775181492337501
after G-S, residual L2: 188.9332667507471
level: 5, grid: 64 x 64
before G-S, residual L2: 129.93801550392874
after G-S, residual L2: 56.28708770794368
level: 4, grid: 32 x 32
before G-S, residual L2: 38.88692621665778
after G-S, residual L2: 18.722754099081875
level: 3, grid: 16 x 16
before G-S, residual L2: 12.92606814051491
after G-S, residual L2: 6.7418584016115615
level: 2, grid: 8 x 8
before G-S, residual L2: 4.646478379380238
after G-S, residual L2: 2.0651261541465855
level: 1, grid: 4 x 4
before G-S, residual L2: 1.3745334259197386
after G-S, residual L2: 0.02244519721859255
bottom solve:
level: 0, grid: 2 x 2
level: 1, grid: 4 x 4
before G-S, residual L2: 0.03125252087247784
after G-S, residual L2: 8.232822131629766e-05
level: 2, grid: 8 x 8
before G-S, residual L2: 2.8059768631102897
after G-S, residual L2: 0.07481536016729919
level: 3, grid: 16 x 16
before G-S, residual L2: 8.77240243659538
after G-S, residual L2: 0.2436194269452686
level: 4, grid: 32 x 32
before G-S, residual L2: 19.59101132435104
after G-S, residual L2: 0.5448263647958954
level: 5, grid: 64 x 64
before G-S, residual L2: 50.4641088994847
after G-S, residual L2: 1.3597629173942345
level: 6, grid: 128 x 128
before G-S, residual L2: 160.2131163846867
after G-S, residual L2: 4.125142056231144
cycle 1: relative err = 0.9999999999999981, residual err = 2.323786088373021
<<< beginning V-cycle (cycle 2) >>>
level: 6, grid: 128 x 128
before G-S, residual L2: 4.125142056231144
after G-S, residual L2: 2.4247311846143984
level: 5, grid: 64 x 64
before G-S, residual L2: 1.6915411385849388
after G-S, residual L2: 1.048624109440286
level: 4, grid: 32 x 32
before G-S, residual L2: 0.728341635357186
after G-S, residual L2: 0.4554818109365305
level: 3, grid: 16 x 16
before G-S, residual L2: 0.3165327512850202
after G-S, residual L2: 0.22128563126748022
level: 2, grid: 8 x 8
before G-S, residual L2: 0.15332496186655523
after G-S, residual L2: 0.07471968817844267
level: 1, grid: 4 x 4
before G-S, residual L2: 0.04974939187294444
after G-S, residual L2: 0.000813357286041041
bottom solve:
level: 0, grid: 2 x 2
level: 1, grid: 4 x 4
before G-S, residual L2: 0.0011325179143730278
after G-S, residual L2: 2.9833778391774223e-06
level: 2, grid: 8 x 8
before G-S, residual L2: 0.10152627387884025
after G-S, residual L2: 0.00270070470024105
level: 3, grid: 16 x 16
before G-S, residual L2: 0.2981467241559525
after G-S, residual L2: 0.008199107952269268
level: 4, grid: 32 x 32
before G-S, residual L2: 0.5218848114624626
after G-S, residual L2: 0.014956130961989951
level: 5, grid: 64 x 64
before G-S, residual L2: 0.9910630869231989
after G-S, residual L2: 0.028422939317571477
level: 6, grid: 128 x 128
before G-S, residual L2: 2.0441877458177586
after G-S, residual L2: 0.05829382601881069
cycle 2: relative err = 0.036315310129800826, residual err = 0.03283823443993396
<<< beginning V-cycle (cycle 3) >>>
level: 6, grid: 128 x 128
before G-S, residual L2: 0.05829382601881069
after G-S, residual L2: 0.04172011870726864
level: 5, grid: 64 x 64
before G-S, residual L2: 0.029246699093099682
after G-S, residual L2: 0.02335632639759113
level: 4, grid: 32 x 32
before G-S, residual L2: 0.01630629679281779
after G-S, residual L2: 0.012906629461195187
level: 3, grid: 16 x 16
before G-S, residual L2: 0.009011110787953677
after G-S, residual L2: 0.00731526293890866
level: 2, grid: 8 x 8
before G-S, residual L2: 0.005081499522859446
after G-S, residual L2: 0.0025625265171556363
level: 1, grid: 4 x 4
before G-S, residual L2: 0.0017064130732666084
after G-S, residual L2: 2.7912387046731846e-05
bottom solve:
level: 0, grid: 2 x 2
level: 1, grid: 4 x 4
before G-S, residual L2: 3.88652692543315e-05
after G-S, residual L2: 1.0238217009469722e-07
level: 2, grid: 8 x 8
before G-S, residual L2: 0.0034819145217790757
after G-S, residual L2: 9.252096659805304e-05
level: 3, grid: 16 x 16
before G-S, residual L2: 0.010064990348703503
after G-S, residual L2: 0.000274405441825591
level: 4, grid: 32 x 32
before G-S, residual L2: 0.016032310448839227
after G-S, residual L2: 0.0004558226543272719
level: 5, grid: 64 x 64
before G-S, residual L2: 0.02430374388018733
after G-S, residual L2: 0.0007098551729200968
level: 6, grid: 128 x 128
before G-S, residual L2: 0.03777531891587048
after G-S, residual L2: 0.001103512282001738
cycle 3: relative err = 0.0012532978372415558, residual err = 0.0006216334987521017
<<< beginning V-cycle (cycle 4) >>>
level: 6, grid: 128 x 128
before G-S, residual L2: 0.001103512282001738
after G-S, residual L2: 0.0008898317346982837
level: 5, grid: 64 x 64
before G-S, residual L2: 0.0006257398720757915
after G-S, residual L2: 0.0006077401190832001
level: 4, grid: 32 x 32
before G-S, residual L2: 0.00042604165447805513
after G-S, residual L2: 0.0003976740182571413
level: 3, grid: 16 x 16
before G-S, residual L2: 0.0002784624522915077
after G-S, residual L2: 0.00024268300992448264
level: 2, grid: 8 x 8
before G-S, residual L2: 0.0001688184030128213
after G-S, residual L2: 8.63435240004183e-05
level: 1, grid: 4 x 4
before G-S, residual L2: 5.750132804421135e-05
after G-S, residual L2: 9.407985171394705e-07
bottom solve:
level: 0, grid: 2 x 2
level: 1, grid: 4 x 4
before G-S, residual L2: 1.309971480329214e-06
after G-S, residual L2: 3.4508339509814223e-09
level: 2, grid: 8 x 8
before G-S, residual L2: 0.0001173242104275028
after G-S, residual L2: 3.115753146780207e-06
level: 3, grid: 16 x 16
before G-S, residual L2: 0.0003385086711958763
after G-S, residual L2: 9.177601888021123e-06
level: 4, grid: 32 x 32
before G-S, residual L2: 0.0005249527904445416
after G-S, residual L2: 1.4651643231018942e-05
level: 5, grid: 64 x 64
before G-S, residual L2: 0.0007080871923403828
after G-S, residual L2: 2.0290645679962866e-05
level: 6, grid: 128 x 128
before G-S, residual L2: 0.0009185166830467631
after G-S, residual L2: 2.657030046653513e-05
cycle 4: relative err = 4.257466296364851e-05, residual err = 1.4967652930826935e-05
<<< beginning V-cycle (cycle 5) >>>
level: 6, grid: 128 x 128
before G-S, residual L2: 2.657030046653513e-05
after G-S, residual L2: 2.3098223935313785e-05
level: 5, grid: 64 x 64
before G-S, residual L2: 1.62748573956126e-05
after G-S, residual L2: 1.790614264207289e-05
level: 4, grid: 32 x 32
before G-S, residual L2: 1.258588239889598e-05
after G-S, residual L2: 1.2880701434507449e-05
level: 3, grid: 16 x 16
before G-S, residual L2: 9.035061893240921e-06
after G-S, residual L2: 8.103003189063974e-06
level: 2, grid: 8 x 8
before G-S, residual L2: 5.641504288200945e-06
after G-S, residual L2: 2.901212906819645e-06
level: 1, grid: 4 x 4
before G-S, residual L2: 1.9321695178564805e-06
after G-S, residual L2: 3.161675602297454e-08
bottom solve:
level: 0, grid: 2 x 2
level: 1, grid: 4 x 4
before G-S, residual L2: 4.402332099930838e-08
after G-S, residual L2: 1.1596974315664115e-10
level: 2, grid: 8 x 8
before G-S, residual L2: 3.9422658752906826e-06
after G-S, residual L2: 1.0466257646978307e-07
level: 3, grid: 16 x 16
before G-S, residual L2: 1.1405869021996882e-05
after G-S, residual L2: 3.081954658995012e-07
level: 4, grid: 32 x 32
before G-S, residual L2: 1.769602521372454e-05
after G-S, residual L2: 4.853326075346603e-07
level: 5, grid: 64 x 64
before G-S, residual L2: 2.2817221850081978e-05
after G-S, residual L2: 6.339093027063977e-07
level: 6, grid: 128 x 128
before G-S, residual L2: 2.7204506619363593e-05
after G-S, residual L2: 7.617366608567251e-07
cycle 5: relative err = 1.437223355768636e-06, residual err = 4.2910353907176844e-07
<<< beginning V-cycle (cycle 6) >>>
level: 6, grid: 128 x 128
before G-S, residual L2: 7.617366608567251e-07
after G-S, residual L2: 6.887955156426563e-07
level: 5, grid: 64 x 64
before G-S, residual L2: 4.858303576730231e-07
after G-S, residual L2: 5.698844687563798e-07
level: 4, grid: 32 x 32
before G-S, residual L2: 4.0114485957443323e-07
after G-S, residual L2: 4.288730517202583e-07
level: 3, grid: 16 x 16
before G-S, residual L2: 3.011320287772163e-07
after G-S, residual L2: 2.722913600393885e-07
level: 2, grid: 8 x 8
before G-S, residual L2: 1.8967555906447816e-07
after G-S, residual L2: 9.770491560795584e-08
level: 1, grid: 4 x 4
before G-S, residual L2: 6.507167362727163e-08
after G-S, residual L2: 1.0648579124114833e-09
bottom solve:
level: 0, grid: 2 x 2
level: 1, grid: 4 x 4
before G-S, residual L2: 1.4827137305196198e-09
after G-S, residual L2: 3.905880555246777e-12
level: 2, grid: 8 x 8
before G-S, residual L2: 1.327670548521078e-07
after G-S, residual L2: 3.5242457966015008e-09
level: 3, grid: 16 x 16
before G-S, residual L2: 3.856314492307069e-07
after G-S, residual L2: 1.0398885089175414e-08
level: 4, grid: 32 x 32
before G-S, residual L2: 6.038836851891428e-07
after G-S, residual L2: 1.6338312488394662e-08
level: 5, grid: 64 x 64
before G-S, residual L2: 7.682416354755789e-07
after G-S, residual L2: 2.077211623429353e-08
level: 6, grid: 128 x 128
before G-S, residual L2: 8.865085868685839e-07
after G-S, residual L2: 2.4019193506676187e-08
cycle 6: relative err = 4.849259894834445e-08, residual err = 1.3530556515124825e-08
<<< beginning V-cycle (cycle 7) >>>
level: 6, grid: 128 x 128
before G-S, residual L2: 2.4019193506676187e-08
after G-S, residual L2: 2.2125281372771198e-08
level: 5, grid: 64 x 64
before G-S, residual L2: 1.561381038807939e-08
after G-S, residual L2: 1.886960673167899e-08
level: 4, grid: 32 x 32
before G-S, residual L2: 1.329268816246029e-08
after G-S, residual L2: 1.4485741753521445e-08
level: 3, grid: 16 x 16
before G-S, residual L2: 1.0177211567137172e-08
after G-S, residual L2: 9.198083184898649e-09
level: 2, grid: 8 x 8
before G-S, residual L2: 6.409466908197649e-09
after G-S, residual L2: 3.3018376959918963e-09
level: 1, grid: 4 x 4
before G-S, residual L2: 2.199060578574703e-09
after G-S, residual L2: 3.598749904926513e-11
bottom solve:
level: 0, grid: 2 x 2
level: 1, grid: 4 x 4
before G-S, residual L2: 5.0109192227870255e-11
after G-S, residual L2: 1.3200150083511782e-13
level: 2, grid: 8 x 8
before G-S, residual L2: 4.48679191788079e-09
after G-S, residual L2: 1.1908944666026909e-10
level: 3, grid: 16 x 16
before G-S, residual L2: 1.308116183600633e-08
after G-S, residual L2: 3.5229822636126385e-10
level: 4, grid: 32 x 32
before G-S, residual L2: 2.0705036251434032e-08
after G-S, residual L2: 5.546643278736882e-10
level: 5, grid: 64 x 64
before G-S, residual L2: 2.628082249674873e-08
after G-S, residual L2: 6.964954416636588e-10
level: 6, grid: 128 x 128
before G-S, residual L2: 2.99444121990791e-08
after G-S, residual L2: 7.914127842417476e-10
cycle 7: relative err = 1.6392149576904378e-09, residual err = 4.458207725000789e-10
<<< beginning V-cycle (cycle 8) >>>
level: 6, grid: 128 x 128
before G-S, residual L2: 7.914127842417476e-10
after G-S, residual L2: 7.355859728039993e-10
level: 5, grid: 64 x 64
before G-S, residual L2: 5.192197657358953e-10
after G-S, residual L2: 6.364658784609479e-10
level: 4, grid: 32 x 32
before G-S, residual L2: 4.485501516615104e-10
after G-S, residual L2: 4.92822276848544e-10
level: 3, grid: 16 x 16
before G-S, residual L2: 3.463701923819941e-10
after G-S, residual L2: 3.1194086454251064e-10
level: 2, grid: 8 x 8
before G-S, residual L2: 2.1741762390347935e-10
after G-S, residual L2: 1.1194487678586295e-10
level: 1, grid: 4 x 4
before G-S, residual L2: 7.455716542482062e-11
after G-S, residual L2: 1.2201470093856338e-12
bottom solve:
level: 0, grid: 2 x 2
level: 1, grid: 4 x 4
before G-S, residual L2: 1.6989396365976684e-12
after G-S, residual L2: 4.475476506307932e-15
level: 2, grid: 8 x 8
before G-S, residual L2: 1.5212108616548206e-10
after G-S, residual L2: 4.037425067492746e-12
level: 3, grid: 16 x 16
before G-S, residual L2: 4.449137567321805e-10
after G-S, residual L2: 1.1972449146203713e-11
level: 4, grid: 32 x 32
before G-S, residual L2: 7.109768080523844e-10
after G-S, residual L2: 1.8912282007395565e-11
level: 5, grid: 64 x 64
before G-S, residual L2: 9.034003713243436e-10
after G-S, residual L2: 2.3606440301264504e-11
level: 6, grid: 128 x 128
before G-S, residual L2: 1.0238015900738368e-09
after G-S, residual L2: 2.6756968686627095e-11
cycle 8: relative err = 5.555097426033948e-11, residual err = 1.5072807373286882e-11
We can compare to the true solution
[20]:
v = a.get_solution()
b = true(a.x2d, a.y2d)
e = v - b
The norm of the error is
[22]:
print(f"{e.norm():20.10g}")
1.671934405e-05
[ ]: