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 conditions.
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)\]
[1]:
import numpy as np
import matplotlib.pyplot as plt
Setting up the solver
[2]:
import pyro.mesh.boundary as bnd
import pyro.mesh.patch as patch
import pyro.multigrid.general_MG as gMG
For reference, we’ll define a function providing the analytic solution
[3]:
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
[4]:
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
[5]:
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
[6]:
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
[7]:
nx = 128
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
[8]:
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
[9]:
a.init_zeros()
a.init_RHS(f(a.x2d, a.y2d))
Source norm = 1.775181492337501
Solving the system
[10]:
a.solve(rtol=1.e-10)
source norm = 1.775181492337501
<<< beginning V-cycle (cycle 1) >>>
level = 6, nx = 128, residual change: 1.77518 → 188.933
level = 5, nx = 64, residual change: 129.938 → 56.2871
level = 4, nx = 32, residual change: 38.8869 → 18.7228
level = 3, nx = 16, residual change: 12.9261 → 6.74186
level = 2, nx = 8, residual change: 4.64648 → 2.06513
level = 1, nx = 4, residual change: 1.37453 → 0.0224452
bottom solve
level = 1, nx = 4, residual change: 0.0312525 → 0.0312525
level = 2, nx = 8, residual change: 2.80598 → 2.80598
level = 3, nx = 16, residual change: 8.7724 → 8.7724
level = 4, nx = 32, residual change: 19.591 → 19.591
level = 5, nx = 64, residual change: 50.4641 → 50.4641
level = 6, nx = 128, residual change: 160.213 → 160.213
cycle 1: relative err = 0.9999999999999981, residual err = 2.323786088373021
<<< beginning V-cycle (cycle 2) >>>
level = 6, nx = 128, residual change: 4.12514 → 2.42473
level = 5, nx = 64, residual change: 1.69154 → 1.04862
level = 4, nx = 32, residual change: 0.728342 → 0.455482
level = 3, nx = 16, residual change: 0.316533 → 0.221286
level = 2, nx = 8, residual change: 0.153325 → 0.0747197
level = 1, nx = 4, residual change: 0.0497494 → 0.000813357
bottom solve
level = 1, nx = 4, residual change: 0.00113252 → 0.00113252
level = 2, nx = 8, residual change: 0.101526 → 0.101526
level = 3, nx = 16, residual change: 0.298147 → 0.298147
level = 4, nx = 32, residual change: 0.521885 → 0.521885
level = 5, nx = 64, residual change: 0.991063 → 0.991063
level = 6, nx = 128, residual change: 2.04419 → 2.04419
cycle 2: relative err = 0.036315310129800826, residual err = 0.03283823443993396
<<< beginning V-cycle (cycle 3) >>>
level = 6, nx = 128, residual change: 0.0582938 → 0.0417201
level = 5, nx = 64, residual change: 0.0292467 → 0.0233563
level = 4, nx = 32, residual change: 0.0163063 → 0.0129066
level = 3, nx = 16, residual change: 0.00901111 → 0.00731526
level = 2, nx = 8, residual change: 0.0050815 → 0.00256253
level = 1, nx = 4, residual change: 0.00170641 → 2.79124e-05
bottom solve
level = 1, nx = 4, residual change: 3.88653e-05 → 3.88653e-05
level = 2, nx = 8, residual change: 0.00348191 → 0.00348191
level = 3, nx = 16, residual change: 0.010065 → 0.010065
level = 4, nx = 32, residual change: 0.0160323 → 0.0160323
level = 5, nx = 64, residual change: 0.0243037 → 0.0243037
level = 6, nx = 128, residual change: 0.0377753 → 0.0377753
cycle 3: relative err = 0.0012532978372415558, residual err = 0.0006216334987521017
<<< beginning V-cycle (cycle 4) >>>
level = 6, nx = 128, residual change: 0.00110351 → 0.000889832
level = 5, nx = 64, residual change: 0.00062574 → 0.00060774
level = 4, nx = 32, residual change: 0.000426042 → 0.000397674
level = 3, nx = 16, residual change: 0.000278462 → 0.000242683
level = 2, nx = 8, residual change: 0.000168818 → 8.63435e-05
level = 1, nx = 4, residual change: 5.75013e-05 → 9.40799e-07
bottom solve
level = 1, nx = 4, residual change: 1.30997e-06 → 1.30997e-06
level = 2, nx = 8, residual change: 0.000117324 → 0.000117324
level = 3, nx = 16, residual change: 0.000338509 → 0.000338509
level = 4, nx = 32, residual change: 0.000524953 → 0.000524953
level = 5, nx = 64, residual change: 0.000708087 → 0.000708087
level = 6, nx = 128, residual change: 0.000918517 → 0.000918517
cycle 4: relative err = 4.257466296364851e-05, residual err = 1.4967652930826935e-05
<<< beginning V-cycle (cycle 5) >>>
level = 6, nx = 128, residual change: 2.65703e-05 → 2.30982e-05
level = 5, nx = 64, residual change: 1.62749e-05 → 1.79061e-05
level = 4, nx = 32, residual change: 1.25859e-05 → 1.28807e-05
level = 3, nx = 16, residual change: 9.03506e-06 → 8.103e-06
level = 2, nx = 8, residual change: 5.6415e-06 → 2.90121e-06
level = 1, nx = 4, residual change: 1.93217e-06 → 3.16168e-08
bottom solve
level = 1, nx = 4, residual change: 4.40233e-08 → 4.40233e-08
level = 2, nx = 8, residual change: 3.94227e-06 → 3.94227e-06
level = 3, nx = 16, residual change: 1.14059e-05 → 1.14059e-05
level = 4, nx = 32, residual change: 1.7696e-05 → 1.7696e-05
level = 5, nx = 64, residual change: 2.28172e-05 → 2.28172e-05
level = 6, nx = 128, residual change: 2.72045e-05 → 2.72045e-05
cycle 5: relative err = 1.437223355768636e-06, residual err = 4.2910353907176844e-07
<<< beginning V-cycle (cycle 6) >>>
level = 6, nx = 128, residual change: 7.61737e-07 → 6.88796e-07
level = 5, nx = 64, residual change: 4.8583e-07 → 5.69884e-07
level = 4, nx = 32, residual change: 4.01145e-07 → 4.28873e-07
level = 3, nx = 16, residual change: 3.01132e-07 → 2.72291e-07
level = 2, nx = 8, residual change: 1.89676e-07 → 9.77049e-08
level = 1, nx = 4, residual change: 6.50717e-08 → 1.06486e-09
bottom solve
level = 1, nx = 4, residual change: 1.48271e-09 → 1.48271e-09
level = 2, nx = 8, residual change: 1.32767e-07 → 1.32767e-07
level = 3, nx = 16, residual change: 3.85631e-07 → 3.85631e-07
level = 4, nx = 32, residual change: 6.03884e-07 → 6.03884e-07
level = 5, nx = 64, residual change: 7.68242e-07 → 7.68242e-07
level = 6, nx = 128, residual change: 8.86509e-07 → 8.86509e-07
cycle 6: relative err = 4.849259894834445e-08, residual err = 1.3530556515124825e-08
<<< beginning V-cycle (cycle 7) >>>
level = 6, nx = 128, residual change: 2.40192e-08 → 2.21253e-08
level = 5, nx = 64, residual change: 1.56138e-08 → 1.88696e-08
level = 4, nx = 32, residual change: 1.32927e-08 → 1.44857e-08
level = 3, nx = 16, residual change: 1.01772e-08 → 9.19808e-09
level = 2, nx = 8, residual change: 6.40947e-09 → 3.30184e-09
level = 1, nx = 4, residual change: 2.19906e-09 → 3.59875e-11
bottom solve
level = 1, nx = 4, residual change: 5.01092e-11 → 5.01092e-11
level = 2, nx = 8, residual change: 4.48679e-09 → 4.48679e-09
level = 3, nx = 16, residual change: 1.30812e-08 → 1.30812e-08
level = 4, nx = 32, residual change: 2.0705e-08 → 2.0705e-08
level = 5, nx = 64, residual change: 2.62808e-08 → 2.62808e-08
level = 6, nx = 128, residual change: 2.99444e-08 → 2.99444e-08
cycle 7: relative err = 1.6392149576904378e-09, residual err = 4.458207725000789e-10
<<< beginning V-cycle (cycle 8) >>>
level = 6, nx = 128, residual change: 7.91413e-10 → 7.35586e-10
level = 5, nx = 64, residual change: 5.1922e-10 → 6.36466e-10
level = 4, nx = 32, residual change: 4.4855e-10 → 4.92822e-10
level = 3, nx = 16, residual change: 3.4637e-10 → 3.11941e-10
level = 2, nx = 8, residual change: 2.17418e-10 → 1.11945e-10
level = 1, nx = 4, residual change: 7.45572e-11 → 1.22015e-12
bottom solve
level = 1, nx = 4, residual change: 1.69894e-12 → 1.69894e-12
level = 2, nx = 8, residual change: 1.52121e-10 → 1.52121e-10
level = 3, nx = 16, residual change: 4.44914e-10 → 4.44914e-10
level = 4, nx = 32, residual change: 7.10977e-10 → 7.10977e-10
level = 5, nx = 64, residual change: 9.034e-10 → 9.034e-10
level = 6, nx = 128, residual change: 1.0238e-09 → 1.0238e-09
cycle 8: relative err = 5.555097426033948e-11, residual err = 1.5072807373286882e-11
Checking the result
We can compare to the true solution
[11]:
v = a.get_solution()
b = true(a.x2d, a.y2d)
e = v - b
The norm of the error is
[12]:
print(f"{e.norm():20.10g}")
1.671934405e-05