Multigrid examples
[1]:
%matplotlib inline
import matplotlib.pyplot as plt
[2]:
from __future__ import print_function
import numpy as np
import pyro.mesh.boundary as bnd
import pyro.mesh.patch as patch
import pyro.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 0x7ff7a8485ea0>
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 0x7ff7a81f26e0>
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{equation} 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{equation}
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 pyro.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 pyro.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.93801550392877
after G-S, residual L2: 56.28708770794367
level: 4, grid: 32 x 32
before G-S, residual L2: 38.88692621665777
after G-S, residual L2: 18.722754099081875
level: 3, grid: 16 x 16
before G-S, residual L2: 12.926068140514912
after G-S, residual L2: 6.7418584016115615
level: 2, grid: 8 x 8
before G-S, residual L2: 4.646478379380239
after G-S, residual L2: 2.0651261541465873
level: 1, grid: 4 x 4
before G-S, residual L2: 1.3745334259197388
after G-S, residual L2: 0.022445197218592287
bottom solve:
level: 0, grid: 2 x 2
level: 1, grid: 4 x 4
before G-S, residual L2: 0.03125252087247734
after G-S, residual L2: 8.232822131646219e-05
level: 2, grid: 8 x 8
before G-S, residual L2: 2.805976863110291
after G-S, residual L2: 0.07481536016730257
level: 3, grid: 16 x 16
before G-S, residual L2: 8.772402436595389
after G-S, residual L2: 0.24361942694526753
level: 4, grid: 32 x 32
before G-S, residual L2: 19.59101132435104
after G-S, residual L2: 0.5448263647958932
level: 5, grid: 64 x 64
before G-S, residual L2: 50.46410889948471
after G-S, residual L2: 1.3597629173942607
level: 6, grid: 128 x 128
before G-S, residual L2: 160.21311638468657
after G-S, residual L2: 4.125142056231161
cycle 1: relative err = 0.9999999999999981, residual err = 2.323786088373031
<<< beginning V-cycle (cycle 2) >>>
level: 6, grid: 128 x 128
before G-S, residual L2: 4.125142056231161
after G-S, residual L2: 2.424731184614396
level: 5, grid: 64 x 64
before G-S, residual L2: 1.6915411385849446
after G-S, residual L2: 1.048624109440289
level: 4, grid: 32 x 32
before G-S, residual L2: 0.7283416353571882
after G-S, residual L2: 0.4554818109365319
level: 3, grid: 16 x 16
before G-S, residual L2: 0.3165327512850212
after G-S, residual L2: 0.22128563126748188
level: 2, grid: 8 x 8
before G-S, residual L2: 0.15332496186655636
after G-S, residual L2: 0.07471968817844332
level: 1, grid: 4 x 4
before G-S, residual L2: 0.049749391872944894
after G-S, residual L2: 0.000813357286041042
bottom solve:
level: 0, grid: 2 x 2
level: 1, grid: 4 x 4
before G-S, residual L2: 0.0011325179143730421
after G-S, residual L2: 2.983377839180104e-06
level: 2, grid: 8 x 8
before G-S, residual L2: 0.10152627387884115
after G-S, residual L2: 0.0027007047002410205
level: 3, grid: 16 x 16
before G-S, residual L2: 0.2981467241559551
after G-S, residual L2: 0.008199107952269191
level: 4, grid: 32 x 32
before G-S, residual L2: 0.521884811462466
after G-S, residual L2: 0.014956130961989571
level: 5, grid: 64 x 64
before G-S, residual L2: 0.9910630869232032
after G-S, residual L2: 0.028422939317571935
level: 6, grid: 128 x 128
before G-S, residual L2: 2.044187745817753
after G-S, residual L2: 0.05829382601881336
cycle 2: relative err = 0.036315310129801166, residual err = 0.032838234439935464
<<< beginning V-cycle (cycle 3) >>>
level: 6, grid: 128 x 128
before G-S, residual L2: 0.05829382601881336
after G-S, residual L2: 0.04172011870730236
level: 5, grid: 64 x 64
before G-S, residual L2: 0.029246699093089436
after G-S, residual L2: 0.023356326397591245
level: 4, grid: 32 x 32
before G-S, residual L2: 0.016306296792817917
after G-S, residual L2: 0.012906629461195626
level: 3, grid: 16 x 16
before G-S, residual L2: 0.009011110787954003
after G-S, residual L2: 0.0073152629389098165
level: 2, grid: 8 x 8
before G-S, residual L2: 0.005081499522860255
after G-S, residual L2: 0.00256252651715607
level: 1, grid: 4 x 4
before G-S, residual L2: 0.0017064130732668981
after G-S, residual L2: 2.79123870467358e-05
bottom solve:
level: 0, grid: 2 x 2
level: 1, grid: 4 x 4
before G-S, residual L2: 3.886526925433716e-05
after G-S, residual L2: 1.0238217009498558e-07
level: 2, grid: 8 x 8
before G-S, residual L2: 0.0034819145217796638
after G-S, residual L2: 9.252096659806975e-05
level: 3, grid: 16 x 16
before G-S, residual L2: 0.010064990348705135
after G-S, residual L2: 0.00027440544182565064
level: 4, grid: 32 x 32
before G-S, residual L2: 0.016032310448840032
after G-S, residual L2: 0.00045582265432728727
level: 5, grid: 64 x 64
before G-S, residual L2: 0.024303743880187103
after G-S, residual L2: 0.0007098551729201522
level: 6, grid: 128 x 128
before G-S, residual L2: 0.037775318915905735
after G-S, residual L2: 0.0011035122820145216
cycle 3: relative err = 0.0012532978372416998, residual err = 0.000621633498759303
<<< beginning V-cycle (cycle 4) >>>
level: 6, grid: 128 x 128
before G-S, residual L2: 0.0011035122820145216
after G-S, residual L2: 0.0008898317346927453
level: 5, grid: 64 x 64
before G-S, residual L2: 0.000625739872070729
after G-S, residual L2: 0.000607740119080587
level: 4, grid: 32 x 32
before G-S, residual L2: 0.00042604165447622126
after G-S, residual L2: 0.000397674018255364
level: 3, grid: 16 x 16
before G-S, residual L2: 0.0002784624522902657
after G-S, residual L2: 0.00024268300992371118
level: 2, grid: 8 x 8
before G-S, residual L2: 0.00016881840301228077
after G-S, residual L2: 8.634352400013625e-05
level: 1, grid: 4 x 4
before G-S, residual L2: 5.7501328044023475e-05
after G-S, residual L2: 9.407985171363897e-07
bottom solve:
level: 0, grid: 2 x 2
level: 1, grid: 4 x 4
before G-S, residual L2: 1.3099714803249488e-06
after G-S, residual L2: 3.4508339509542035e-09
level: 2, grid: 8 x 8
before G-S, residual L2: 0.00011732421042711928
after G-S, residual L2: 3.1157531467700274e-06
level: 3, grid: 16 x 16
before G-S, residual L2: 0.00033850867119458294
after G-S, residual L2: 9.177601887988707e-06
level: 4, grid: 32 x 32
before G-S, residual L2: 0.0005249527904411394
after G-S, residual L2: 1.4651643230949486e-05
level: 5, grid: 64 x 64
before G-S, residual L2: 0.0007080871923349976
after G-S, residual L2: 2.0290645679874954e-05
level: 6, grid: 128 x 128
before G-S, residual L2: 0.0009185166831419814
after G-S, residual L2: 2.657030047801386e-05
cycle 4: relative err = 4.257466296345195e-05, residual err = 1.4967652937293164e-05
<<< beginning V-cycle (cycle 5) >>>
level: 6, grid: 128 x 128
before G-S, residual L2: 2.657030047801386e-05
after G-S, residual L2: 2.3098223934902124e-05
level: 5, grid: 64 x 64
before G-S, residual L2: 1.627485739129596e-05
after G-S, residual L2: 1.7906142640602113e-05
level: 4, grid: 32 x 32
before G-S, residual L2: 1.2585882397882331e-05
after G-S, residual L2: 1.2880701433409095e-05
level: 3, grid: 16 x 16
before G-S, residual L2: 9.035061892447504e-06
after G-S, residual L2: 8.103003187752153e-06
level: 2, grid: 8 x 8
before G-S, residual L2: 5.641504287281656e-06
after G-S, residual L2: 2.90121290636155e-06
level: 1, grid: 4 x 4
before G-S, residual L2: 1.9321695175514817e-06
after G-S, residual L2: 3.1616756017989616e-08
bottom solve:
level: 0, grid: 2 x 2
level: 1, grid: 4 x 4
before G-S, residual L2: 4.402332099236765e-08
after G-S, residual L2: 1.159697431375499e-10
level: 2, grid: 8 x 8
before G-S, residual L2: 3.942265874668158e-06
after G-S, residual L2: 1.0466257645320907e-07
level: 3, grid: 16 x 16
before G-S, residual L2: 1.1405869020197195e-05
after G-S, residual L2: 3.081954658477386e-07
level: 4, grid: 32 x 32
before G-S, residual L2: 1.7696025211153306e-05
after G-S, residual L2: 4.853326074732553e-07
level: 5, grid: 64 x 64
before G-S, residual L2: 2.281722184539092e-05
after G-S, residual L2: 6.339093026169473e-07
level: 6, grid: 128 x 128
before G-S, residual L2: 2.720450664755927e-05
after G-S, residual L2: 7.617366442287773e-07
cycle 5: relative err = 1.4372233556002303e-06, residual err = 4.291035297048683e-07
<<< beginning V-cycle (cycle 6) >>>
level: 6, grid: 128 x 128
before G-S, residual L2: 7.617366442287773e-07
after G-S, residual L2: 6.887955204981268e-07
level: 5, grid: 64 x 64
before G-S, residual L2: 4.858303555357937e-07
after G-S, residual L2: 5.69884466371041e-07
level: 4, grid: 32 x 32
before G-S, residual L2: 4.0114485792406196e-07
after G-S, residual L2: 4.288730508746916e-07
level: 3, grid: 16 x 16
before G-S, residual L2: 3.011320281769793e-07
after G-S, residual L2: 2.722913591694268e-07
level: 2, grid: 8 x 8
before G-S, residual L2: 1.8967555845982626e-07
after G-S, residual L2: 9.770491535790604e-08
level: 1, grid: 4 x 4
before G-S, residual L2: 6.507167346101259e-08
after G-S, residual L2: 1.064857909707231e-09
bottom solve:
level: 0, grid: 2 x 2
level: 1, grid: 4 x 4
before G-S, residual L2: 1.4827137267543486e-09
after G-S, residual L2: 3.905880545305021e-12
level: 2, grid: 8 x 8
before G-S, residual L2: 1.327670545122848e-07
after G-S, residual L2: 3.5242457874456e-09
level: 3, grid: 16 x 16
before G-S, residual L2: 3.856314482142684e-07
after G-S, residual L2: 1.0398885055660428e-08
level: 4, grid: 32 x 32
before G-S, residual L2: 6.038836839330361e-07
after G-S, residual L2: 1.6338312450885264e-08
level: 5, grid: 64 x 64
before G-S, residual L2: 7.682416327633833e-07
after G-S, residual L2: 2.077211614487108e-08
level: 6, grid: 128 x 128
before G-S, residual L2: 8.865086249123311e-07
after G-S, residual L2: 2.401917734271176e-08
cycle 6: relative err = 4.849259893743e-08, residual err = 1.3530547409597026e-08
<<< beginning V-cycle (cycle 7) >>>
level: 6, grid: 128 x 128
before G-S, residual L2: 2.401917734271176e-08
after G-S, residual L2: 2.212526516368908e-08
level: 5, grid: 64 x 64
before G-S, residual L2: 1.561380930193641e-08
after G-S, residual L2: 1.8869605253841384e-08
level: 4, grid: 32 x 32
before G-S, residual L2: 1.329268715576614e-08
after G-S, residual L2: 1.448574192361666e-08
level: 3, grid: 16 x 16
before G-S, residual L2: 1.0177211702855934e-08
after G-S, residual L2: 9.198083215964472e-09
level: 2, grid: 8 x 8
before G-S, residual L2: 6.409466931052197e-09
after G-S, residual L2: 3.301837694167828e-09
level: 1, grid: 4 x 4
before G-S, residual L2: 2.1990605773774708e-09
after G-S, residual L2: 3.5987499029820695e-11
bottom solve:
level: 0, grid: 2 x 2
level: 1, grid: 4 x 4
before G-S, residual L2: 5.010919220080045e-11
after G-S, residual L2: 1.3200150076369783e-13
level: 2, grid: 8 x 8
before G-S, residual L2: 4.486791915539427e-09
after G-S, residual L2: 1.1908944659498306e-10
level: 3, grid: 16 x 16
before G-S, residual L2: 1.3081161953516478e-08
after G-S, residual L2: 3.5229822904429924e-10
level: 4, grid: 32 x 32
before G-S, residual L2: 2.0705036864241283e-08
after G-S, residual L2: 5.54664336501412e-10
level: 5, grid: 64 x 64
before G-S, residual L2: 2.6280822545639532e-08
after G-S, residual L2: 6.964954084454355e-10
level: 6, grid: 128 x 128
before G-S, residual L2: 2.994436334203582e-08
after G-S, residual L2: 7.914114964975534e-10
cycle 7: relative err = 1.6392149521630824e-09, residual err = 4.4582004708456516e-10
<<< beginning V-cycle (cycle 8) >>>
level: 6, grid: 128 x 128
before G-S, residual L2: 7.914114964975534e-10
after G-S, residual L2: 7.355563798981241e-10
level: 5, grid: 64 x 64
before G-S, residual L2: 5.192187012540891e-10
after G-S, residual L2: 6.364663794896949e-10
level: 4, grid: 32 x 32
before G-S, residual L2: 4.4855051982056546e-10
after G-S, residual L2: 4.92823302521467e-10
level: 3, grid: 16 x 16
before G-S, residual L2: 3.463709288611066e-10
after G-S, residual L2: 3.1194153963395727e-10
level: 2, grid: 8 x 8
before G-S, residual L2: 2.1741809566511453e-10
after G-S, residual L2: 1.1194508137188934e-10
level: 1, grid: 4 x 4
before G-S, residual L2: 7.455730160141668e-11
after G-S, residual L2: 1.2201492319919709e-12
bottom solve:
level: 0, grid: 2 x 2
level: 1, grid: 4 x 4
before G-S, residual L2: 1.6989427313305885e-12
after G-S, residual L2: 4.47548465864942e-15
level: 2, grid: 8 x 8
before G-S, residual L2: 1.5212136435309373e-10
after G-S, residual L2: 4.037432498150632e-12
level: 3, grid: 16 x 16
before G-S, residual L2: 4.4491471941548957e-10
after G-S, residual L2: 1.1972476764962434e-11
level: 4, grid: 32 x 32
before G-S, residual L2: 7.109786940662245e-10
after G-S, residual L2: 1.891232351798468e-11
level: 5, grid: 64 x 64
before G-S, residual L2: 9.034023723354117e-10
after G-S, residual L2: 2.3606468222333935e-11
level: 6, grid: 128 x 128
before G-S, residual L2: 1.023747974109598e-09
after G-S, residual L2: 2.6771270744711738e-11
cycle 8: relative err = 5.555103471809045e-11, residual err = 1.5080864047010877e-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
[21]:
print(f"{e.norm():20.10g}")
1.671934405e-05
[ ]: