# Multigrid examples



In [1]:

%matplotlib inline
import matplotlib.pyplot as plt




In [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



In [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).



In [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.09751581367



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



In [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()



In [6]:

mg.solve()




source norm =  1.09751581367
<<< 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()



In [7]:

phi = mg.get_solution()




In [8]:

plt.imshow(np.transpose(phi.v()), origin="lower")




Out[8]:

<matplotlib.image.AxesImage at 0x7f7c47a0cc50>



we can also get the gradient of the solution



In [9]:




In [10]:

plt.subplot(121)
plt.imshow(np.transpose(gx.v()), origin="lower")
plt.subplot(122)
plt.imshow(np.transpose(gy.v()), origin="lower")




Out[10]:

<matplotlib.image.AxesImage at 0x7f7c47a345f8>



## 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)$$



In [11]:

import multigrid.general_MG as gMG



For reference, we'll define a function providing the analytic solution



In [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



In [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



In [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



In [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



In [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



In [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



In [18]:

a.init_zeros()
a.init_RHS(f(a.x2d, a.y2d))




Source norm =  1.77518149234



and we can solve it



In [19]:

a.solve(rtol=1.e-10)




source norm =  1.77518149234
<<< 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.741858401611561

level: 2, grid: 8 x 8
before G-S, residual L2: 4.646478379380238
after G-S, residual L2: 2.065126154146587

level: 1, grid: 4 x 4
before G-S, residual L2: 1.3745334259197384
after G-S, residual L2: 0.02244519721859215

bottom solve:
level: 0, grid: 2 x 2

level: 1, grid: 4 x 4
before G-S, residual L2: 0.031252520872477096
after G-S, residual L2: 8.232822131685586e-05

level: 2, grid: 8 x 8
before G-S, residual L2: 2.8059768631102893
after G-S, residual L2: 0.07481536016730024

level: 3, grid: 16 x 16
before G-S, residual L2: 8.772402436595382
after G-S, residual L2: 0.24361942694526875

level: 4, grid: 32 x 32
before G-S, residual L2: 19.591011324351037
after G-S, residual L2: 0.5448263647958976

level: 5, grid: 64 x 64
before G-S, residual L2: 50.4641088994847
after G-S, residual L2: 1.3597629173942398

level: 6, grid: 128 x 128
before G-S, residual L2: 160.2131163846867
after G-S, residual L2: 4.125142056231141

cycle 1: relative err = 0.9999999999999981, residual err = 2.3237860883730193

<<< beginning V-cycle (cycle 2) >>>

level: 6, grid: 128 x 128
before G-S, residual L2: 4.125142056231141
after G-S, residual L2: 2.4247311846143957

level: 5, grid: 64 x 64
before G-S, residual L2: 1.6915411385849393
after G-S, residual L2: 1.0486241094402862

level: 4, grid: 32 x 32
before G-S, residual L2: 0.7283416353571861
after G-S, residual L2: 0.45548181093652995

level: 3, grid: 16 x 16
before G-S, residual L2: 0.3165327512850198
after G-S, residual L2: 0.22128563126748008

level: 2, grid: 8 x 8
before G-S, residual L2: 0.15332496186655512
after G-S, residual L2: 0.0747196881784426

level: 1, grid: 4 x 4
before G-S, residual L2: 0.04974939187294444
after G-S, residual L2: 0.0008133572860410457

bottom solve:
level: 0, grid: 2 x 2

level: 1, grid: 4 x 4
before G-S, residual L2: 0.0011325179143730458
after G-S, residual L2: 2.98337783917788e-06

level: 2, grid: 8 x 8
before G-S, residual L2: 0.10152627387884022
after G-S, residual L2: 0.0027007047002410374

level: 3, grid: 16 x 16
before G-S, residual L2: 0.29814672415595245
after G-S, residual L2: 0.00819910795226899

level: 4, grid: 32 x 32
before G-S, residual L2: 0.5218848114624619
after G-S, residual L2: 0.014956130961989498

level: 5, grid: 64 x 64
before G-S, residual L2: 0.9910630869231989
after G-S, residual L2: 0.028422939317571984

level: 6, grid: 128 x 128
before G-S, residual L2: 2.044187745817752
after G-S, residual L2: 0.058293826018797935

cycle 2: relative err = 0.036315310129800826, residual err = 0.032838234439926776

<<< beginning V-cycle (cycle 3) >>>

level: 6, grid: 128 x 128
before G-S, residual L2: 0.058293826018797935
after G-S, residual L2: 0.0417201187072595

level: 5, grid: 64 x 64
before G-S, residual L2: 0.029246699093099564
after G-S, residual L2: 0.023356326397591495

level: 4, grid: 32 x 32
before G-S, residual L2: 0.016306296792818056
after G-S, residual L2: 0.012906629461195234

level: 3, grid: 16 x 16
before G-S, residual L2: 0.009011110787953703
after G-S, residual L2: 0.007315262938908486

level: 2, grid: 8 x 8
before G-S, residual L2: 0.005081499522859323
after G-S, residual L2: 0.002562526517155576

level: 1, grid: 4 x 4
before G-S, residual L2: 0.0017064130732665692
after G-S, residual L2: 2.7912387046731474e-05

bottom solve:
level: 0, grid: 2 x 2

level: 1, grid: 4 x 4
before G-S, residual L2: 3.886526925433118e-05
after G-S, residual L2: 1.0238217009484441e-07

level: 2, grid: 8 x 8
before G-S, residual L2: 0.0034819145217789937
after G-S, residual L2: 9.252096659805176e-05

level: 3, grid: 16 x 16
before G-S, residual L2: 0.01006499034870321
after G-S, residual L2: 0.0002744054418255884

level: 4, grid: 32 x 32
before G-S, residual L2: 0.016032310448838724
after G-S, residual L2: 0.0004558226543272663

level: 5, grid: 64 x 64
before G-S, residual L2: 0.024303743880186898
after G-S, residual L2: 0.0007098551729201239

level: 6, grid: 128 x 128
before G-S, residual L2: 0.037775318915862
after G-S, residual L2: 0.0011035122819927912

cycle 3: relative err = 0.0012532978372415335, residual err = 0.0006216334987470617

<<< beginning V-cycle (cycle 4) >>>

level: 6, grid: 128 x 128
before G-S, residual L2: 0.0011035122819927912
after G-S, residual L2: 0.0008898317346917108

level: 5, grid: 64 x 64
before G-S, residual L2: 0.0006257398720776081
after G-S, residual L2: 0.000607740119084607

level: 4, grid: 32 x 32
before G-S, residual L2: 0.00042604165447901086
after G-S, residual L2: 0.00039767401825608673

level: 3, grid: 16 x 16
before G-S, residual L2: 0.0002784624522907369
after G-S, residual L2: 0.00024268300992319052

level: 2, grid: 8 x 8
before G-S, residual L2: 0.0001688184030119159
after G-S, residual L2: 8.63435239999583e-05

level: 1, grid: 4 x 4
before G-S, residual L2: 5.750132804390505e-05
after G-S, residual L2: 9.407985171344554e-07

bottom solve:
level: 0, grid: 2 x 2

level: 1, grid: 4 x 4
before G-S, residual L2: 1.3099714803222558e-06
after G-S, residual L2: 3.450833950914012e-09

level: 2, grid: 8 x 8
before G-S, residual L2: 0.00011732421042687768
after G-S, residual L2: 3.1157531467636086e-06

level: 3, grid: 16 x 16
before G-S, residual L2: 0.00033850867119400885
after G-S, residual L2: 9.17760188796962e-06

level: 4, grid: 32 x 32
before G-S, residual L2: 0.0005249527904418192
after G-S, residual L2: 1.4651643230958405e-05

level: 5, grid: 64 x 64
before G-S, residual L2: 0.0007080871923387015
after G-S, residual L2: 2.0290645679943462e-05

level: 6, grid: 128 x 128
before G-S, residual L2: 0.0009185166830535544
after G-S, residual L2: 2.6570300453995103e-05

cycle 4: relative err = 4.2574662963457396e-05, residual err = 1.4967652923762853e-05

<<< beginning V-cycle (cycle 5) >>>

level: 6, grid: 128 x 128
before G-S, residual L2: 2.6570300453995103e-05
after G-S, residual L2: 2.3098223923757352e-05

level: 5, grid: 64 x 64
before G-S, residual L2: 1.6274857395354832e-05
after G-S, residual L2: 1.7906142642175535e-05

level: 4, grid: 32 x 32
before G-S, residual L2: 1.258588239896169e-05
after G-S, residual L2: 1.2880701433730278e-05

level: 3, grid: 16 x 16
before G-S, residual L2: 9.035061892671461e-06
after G-S, residual L2: 8.10300318788889e-06

level: 2, grid: 8 x 8
before G-S, residual L2: 5.641504287378599e-06
after G-S, residual L2: 2.9012129063955126e-06

level: 1, grid: 4 x 4
before G-S, residual L2: 1.932169517574082e-06
after G-S, residual L2: 3.161675601835735e-08

bottom solve:
level: 0, grid: 2 x 2

level: 1, grid: 4 x 4
before G-S, residual L2: 4.4023320992879136e-08
after G-S, residual L2: 1.1596974313938014e-10

level: 2, grid: 8 x 8
before G-S, residual L2: 3.9422658747144435e-06
after G-S, residual L2: 1.0466257645445924e-07

level: 3, grid: 16 x 16
before G-S, residual L2: 1.1405869020431955e-05
after G-S, residual L2: 3.0819546585464564e-07

level: 4, grid: 32 x 32
before G-S, residual L2: 1.7696025211842327e-05
after G-S, residual L2: 4.853326074858634e-07

level: 5, grid: 64 x 64
before G-S, residual L2: 2.281722184794443e-05
after G-S, residual L2: 6.339093026629609e-07

level: 6, grid: 128 x 128
before G-S, residual L2: 2.7204506586512792e-05
after G-S, residual L2: 7.61736677407384e-07

cycle 5: relative err = 1.4372233555992132e-06, residual err = 4.2910354839513e-07

<<< beginning V-cycle (cycle 6) >>>

level: 6, grid: 128 x 128
before G-S, residual L2: 7.61736677407384e-07
after G-S, residual L2: 6.887955287148536e-07

level: 5, grid: 64 x 64
before G-S, residual L2: 4.858303580829294e-07
after G-S, residual L2: 5.698844682533653e-07

level: 4, grid: 32 x 32
before G-S, residual L2: 4.011448592273346e-07
after G-S, residual L2: 4.2887305175998083e-07

level: 3, grid: 16 x 16
before G-S, residual L2: 3.011320287970724e-07
after G-S, residual L2: 2.7229135972437344e-07

level: 2, grid: 8 x 8
before G-S, residual L2: 1.8967555884605451e-07
after G-S, residual L2: 9.770491553515245e-08

level: 1, grid: 4 x 4
before G-S, residual L2: 6.507167357899105e-08
after G-S, residual L2: 1.0648579116334552e-09

bottom solve:
level: 0, grid: 2 x 2

level: 1, grid: 4 x 4
before G-S, residual L2: 1.4827137294363792e-09
after G-S, residual L2: 3.9058805523605475e-12

level: 2, grid: 8 x 8
before G-S, residual L2: 1.3276705475319977e-07
after G-S, residual L2: 3.524245793876337e-09

level: 3, grid: 16 x 16
before G-S, residual L2: 3.8563144896921417e-07
after G-S, residual L2: 1.0398885077513769e-08

level: 4, grid: 32 x 32
before G-S, residual L2: 6.038836850187365e-07
after G-S, residual L2: 1.6338312481157817e-08

level: 5, grid: 64 x 64
before G-S, residual L2: 7.682416346530921e-07
after G-S, residual L2: 2.0772116210685317e-08

level: 6, grid: 128 x 128
before G-S, residual L2: 8.865086230602598e-07
after G-S, residual L2: 2.401923227919822e-08

cycle 6: relative err = 4.8492598977484135e-08, residual err = 1.353057835656594e-08

<<< beginning V-cycle (cycle 7) >>>

level: 6, grid: 128 x 128
before G-S, residual L2: 2.401923227919822e-08
after G-S, residual L2: 2.2125290070425652e-08

level: 5, grid: 64 x 64
before G-S, residual L2: 1.5613809613835955e-08
after G-S, residual L2: 1.8869606239963252e-08

level: 4, grid: 32 x 32
before G-S, residual L2: 1.3292687837677291e-08
after G-S, residual L2: 1.4485742520315527e-08

level: 3, grid: 16 x 16
before G-S, residual L2: 1.0177212111802273e-08
after G-S, residual L2: 9.198083791538658e-09

level: 2, grid: 8 x 8
before G-S, residual L2: 6.409467335640698e-09
after G-S, residual L2: 3.3018379633629456e-09

level: 1, grid: 4 x 4
before G-S, residual L2: 2.1990607567876347e-09
after G-S, residual L2: 3.598750197454369e-11

bottom solve:
level: 0, grid: 2 x 2

level: 1, grid: 4 x 4
before G-S, residual L2: 5.010919630110133e-11
after G-S, residual L2: 1.3200151156453123e-13

level: 2, grid: 8 x 8
before G-S, residual L2: 4.48679228107323e-09
after G-S, residual L2: 1.1908945622999935e-10

level: 3, grid: 16 x 16
before G-S, residual L2: 1.3081162779667808e-08
after G-S, residual L2: 3.522982496836639e-10

level: 4, grid: 32 x 32
before G-S, residual L2: 2.0705037621548675e-08
after G-S, residual L2: 5.546643639307605e-10

level: 5, grid: 64 x 64
before G-S, residual L2: 2.6280822057541362e-08
after G-S, residual L2: 6.964954384251476e-10

level: 6, grid: 128 x 128
before G-S, residual L2: 2.994449911367404e-08
after G-S, residual L2: 7.914383325620475e-10

cycle 7: relative err = 1.6392150533299687e-09, residual err = 4.4583516444840087e-10

<<< beginning V-cycle (cycle 8) >>>

level: 6, grid: 128 x 128
before G-S, residual L2: 7.914383325620475e-10
after G-S, residual L2: 7.355629304356289e-10

level: 5, grid: 64 x 64
before G-S, residual L2: 5.19218220597571e-10
after G-S, residual L2: 6.364663261794707e-10

level: 4, grid: 32 x 32
before G-S, residual L2: 4.485504928535875e-10
after G-S, residual L2: 4.928237246176745e-10

level: 3, grid: 16 x 16
before G-S, residual L2: 3.4637122000064977e-10
after G-S, residual L2: 3.1194162913950586e-10

level: 2, grid: 8 x 8
before G-S, residual L2: 2.174181615639314e-10
after G-S, residual L2: 1.1194514367241423e-10

level: 1, grid: 4 x 4
before G-S, residual L2: 7.455734323986808e-11
after G-S, residual L2: 1.2201499216239134e-12

bottom solve:
level: 0, grid: 2 x 2

level: 1, grid: 4 x 4
before G-S, residual L2: 1.6989436916357301e-12
after G-S, residual L2: 4.475487188247986e-15

level: 2, grid: 8 x 8
before G-S, residual L2: 1.521214490284944e-10
after G-S, residual L2: 4.037434677870943e-12

level: 3, grid: 16 x 16
before G-S, residual L2: 4.4491498629640967e-10
after G-S, residual L2: 1.197248120085576e-11

level: 4, grid: 32 x 32
before G-S, residual L2: 7.109792371905777e-10
after G-S, residual L2: 1.8912335700376235e-11

level: 5, grid: 64 x 64
before G-S, residual L2: 9.034017109357381e-10
after G-S, residual L2: 2.3606466325271617e-11

level: 6, grid: 128 x 128
before G-S, residual L2: 1.0238349148814258e-09
after G-S, residual L2: 2.678477889744364e-11

cycle 8: relative err = 5.555107077431201e-11, residual err = 1.5088473495842003e-11



We can compare to the true solution



In [20]:

v = a.get_solution()
b = true(a.x2d, a.y2d)
e = v - b



The norm of the error is



In [21]:

e.norm()




Out[21]:

1.6719344048744095e-05




In [ ]: