General Linear Elliptic Equation#

The GeneralMG2d class implements support for a general elliptic equation of the form:

α ϕ + \nabla \cdot (β \nabla ϕ) + γ \cdot \nabla ϕ = f

with inhomogeneous boundary conditions.

It subclasses the CellCenterMG2d class, and the basic interface is the same

We will solve the above with

$\begin{aligned} α & = 10 \\ β & = x y + 1 \\ γ & = \hat{x} + \hat{y} \end{aligned}$

and

$\begin{aligned} f = & - \frac{π}{2} (x + 1) \sin (\frac{π y}{2}) \cos (\frac{π x}{2}) \\ - \frac{π}{2} (y + 1) \sin (\frac{π x}{2}) \cos (\frac{π y}{2}) \\ + (\frac{- π^{2} (x y + 1)}{2} + 10) \cos (\frac{π x}{2}) \cos (\frac{π y}{2}) \end{aligned}$

on $[0, 1] \times [0, 1]$ with boundary conditions:

$\begin{aligned} ϕ (x = 0) & = \cos (π y / 2) \\ ϕ (x = 1) & = 0 \\ ϕ (y = 0) & = \cos (π x / 2) \\ ϕ (y = 1) & = 0 \end{aligned}$

This has the exact solution:

ϕ = \cos (π x / 2) \cos (π 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 $γ$ 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 $ϕ$ (and for real problems, they may not). The one that matters the most is $β$ , 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

General Linear Elliptic Equation

Contents

General Linear Elliptic Equation#

Setting up the solver#

Solving the system#

Checking the result#