Variable Coefficient Poisson#

We want to solve an equation of the form $\nabla \cdot (α \nabla ϕ) = f$

We’ll do this with periodic boundary conditions. Consider the coefficient $α$ of the form:

α = 2 + \cos (2 π x) \cos (2 π x)

and the source, $f$ :

f = - 16 π^{2} (\cos (2 π x) \cos (2 π y) + 1) \sin (2 π x) \sin (2 π y)

The solution to this (with periodic BCs) is:

ϕ = \sin (2 π x) \sin (2 π y)

[1]:

import numpy as np
import matplotlib.pyplot as plt

Setting up the solver#

[2]:

import pyro.multigrid.variable_coeff_MG as MG

[3]:

def true(x, y):
    return np.sin(2.0*np.pi*x)*np.sin(2.0*np.pi*y)

[4]:

def alpha(x, y):
    return 2.0 + np.cos(2.0*np.pi*x)*np.cos(2.0*np.pi*y)

[5]:

def f(x, y):
    return -16.0*np.pi**2*(np.cos(2*np.pi*x)*np.cos(2*np.pi*y) + 1) * \
        np.sin(2*np.pi*x)*np.sin(2*np.pi*y)

Let’s create a patch to store the coefficient $α$

[6]:

from pyro.mesh import patch
import pyro.mesh.boundary as bnd

[7]:

N = 128

g = patch.Grid2d(N, N, ng=1)
d = patch.CellCenterData2d(g)
bc_alpha = bnd.BC(xlb="periodic", xrb="periodic",
                  ylb="periodic", yrb="periodic")
d.register_var("alpha", bc_alpha)
d.create()

Now we can fill the coefficient

[8]:

a = d.get_var("alpha")
a[:, :] = alpha(g.x2d, g.y2d)

With periodic BCs, solvability requires that $f$ sum to zero over the domain. Let’s check that.

[9]:

rhs = f(g.x2d, g.y2d)
print(f"rhs sum: {np.sum(rhs[g.ilo:g.ihi+1, g.jlo:g.jhi+1]):20.10g}")

rhs sum:       2.07633187e-12

Now we can create the multigrid object

[10]:

mg = MG.VarCoeffCCMG2d(N, N,
                       xl_BC_type="periodic", yl_BC_type="periodic",
                       xr_BC_type="periodic", yr_BC_type="periodic",
                       coeffs=a, coeffs_bc=bc_alpha,
                       verbose=1, vis=0, true_function=true)

cc data: nx = 2, ny = 2, ng = 1
         nvars = 4
         variables:
               v: min:    0.0000000000    max:    0.0000000000
                  BCs: -x: periodic     +x: periodic     -y: periodic     +y: periodic
               f: min:    0.0000000000    max:    0.0000000000
                  BCs: -x: periodic     +x: periodic     -y: periodic     +y: periodic
               r: min:    0.0000000000    max:    0.0000000000
                  BCs: -x: periodic     +x: periodic     -y: periodic     +y: periodic
          coeffs: min:    0.0000000000    max:    0.0000000000
                  BCs: -x: periodic     +x: periodic     -y: periodic     +y: periodic

cc data: nx = 4, ny = 4, ng = 1
         nvars = 4
         variables:
               v: min:    0.0000000000    max:    0.0000000000
                  BCs: -x: periodic     +x: periodic     -y: periodic     +y: periodic
               f: min:    0.0000000000    max:    0.0000000000
                  BCs: -x: periodic     +x: periodic     -y: periodic     +y: periodic
               r: min:    0.0000000000    max:    0.0000000000
                  BCs: -x: periodic     +x: periodic     -y: periodic     +y: periodic
          coeffs: min:    0.0000000000    max:    0.0000000000
                  BCs: -x: periodic     +x: periodic     -y: periodic     +y: periodic

cc data: nx = 8, ny = 8, ng = 1
         nvars = 4
         variables:
               v: min:    0.0000000000    max:    0.0000000000
                  BCs: -x: periodic     +x: periodic     -y: periodic     +y: periodic
               f: min:    0.0000000000    max:    0.0000000000
                  BCs: -x: periodic     +x: periodic     -y: periodic     +y: periodic
               r: min:    0.0000000000    max:    0.0000000000
                  BCs: -x: periodic     +x: periodic     -y: periodic     +y: periodic
          coeffs: min:    0.0000000000    max:    0.0000000000
                  BCs: -x: periodic     +x: periodic     -y: periodic     +y: periodic

cc data: nx = 16, ny = 16, ng = 1
         nvars = 4
         variables:
               v: min:    0.0000000000    max:    0.0000000000
                  BCs: -x: periodic     +x: periodic     -y: periodic     +y: periodic
               f: min:    0.0000000000    max:    0.0000000000
                  BCs: -x: periodic     +x: periodic     -y: periodic     +y: periodic
               r: min:    0.0000000000    max:    0.0000000000
                  BCs: -x: periodic     +x: periodic     -y: periodic     +y: periodic
          coeffs: min:    0.0000000000    max:    0.0000000000
                  BCs: -x: periodic     +x: periodic     -y: periodic     +y: periodic

cc data: nx = 32, ny = 32, ng = 1
         nvars = 4
         variables:
               v: min:    0.0000000000    max:    0.0000000000
                  BCs: -x: periodic     +x: periodic     -y: periodic     +y: periodic
               f: min:    0.0000000000    max:    0.0000000000
                  BCs: -x: periodic     +x: periodic     -y: periodic     +y: periodic
               r: min:    0.0000000000    max:    0.0000000000
                  BCs: -x: periodic     +x: periodic     -y: periodic     +y: periodic
          coeffs: min:    0.0000000000    max:    0.0000000000
                  BCs: -x: periodic     +x: periodic     -y: periodic     +y: periodic

cc data: nx = 64, ny = 64, ng = 1
         nvars = 4
         variables:
               v: min:    0.0000000000    max:    0.0000000000
                  BCs: -x: periodic     +x: periodic     -y: periodic     +y: periodic
               f: min:    0.0000000000    max:    0.0000000000
                  BCs: -x: periodic     +x: periodic     -y: periodic     +y: periodic
               r: min:    0.0000000000    max:    0.0000000000
                  BCs: -x: periodic     +x: periodic     -y: periodic     +y: periodic
          coeffs: min:    0.0000000000    max:    0.0000000000
                  BCs: -x: periodic     +x: periodic     -y: periodic     +y: periodic

cc data: nx = 128, ny = 128, ng = 1
         nvars = 4
         variables:
               v: min:    0.0000000000    max:    0.0000000000
                  BCs: -x: periodic     +x: periodic     -y: periodic     +y: periodic
               f: min:    0.0000000000    max:    0.0000000000
                  BCs: -x: periodic     +x: periodic     -y: periodic     +y: periodic
               r: min:    0.0000000000    max:    0.0000000000
                  BCs: -x: periodic     +x: periodic     -y: periodic     +y: periodic
          coeffs: min:    0.0000000000    max:    0.0000000000
                  BCs: -x: periodic     +x: periodic     -y: periodic     +y: periodic

Initialize the solution to 0

[11]:

mg.init_zeros()

Now initialize the RHS

[12]:

rhs = f(mg.x2d, mg.y2d)
mg.init_RHS(rhs)

Source norm =  81.3868428575047

Solving the system#

[13]:

mg.solve(rtol=1.e-11)

source norm =  81.3868428575047
<<< beginning V-cycle (cycle 1) >>>

  level =  6, nx =  128, residual change:     81.3868 →     112.091
  level =  5, nx =   64, residual change:     79.2048 →      101.02
  level =  4, nx =   32, residual change:     71.2411 →     68.0705
  level =  3, nx =   16, residual change:     47.6577 →     14.4578
  level =  2, nx =    8, residual change:      9.8448 →   0.0171409
  level =  1, nx =    4, residual change:   0.0106141 → 1.69329e-16
  bottom solve
  level =  1, nx =    4, residual change: 3.07629e-16 → 3.07629e-16
  level =  2, nx =    8, residual change:   0.0168243 →   0.0168243
  level =  3, nx =   16, residual change:     14.4901 →     14.4901
  level =  4, nx =   32, residual change:      69.371 →      69.371
  level =  5, nx =   64, residual change:     103.884 →     103.884
  level =  6, nx =  128, residual change:     115.511 →     115.511
cycle 1: relative err = 1.0000000000000007, residual err = 0.025573219961900512

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

  level =  6, nx =  128, residual change:     2.08132 →     2.02861
  level =  5, nx =   64, residual change:     1.43441 →     1.83684
  level =  4, nx =   32, residual change:     1.29871 →     1.24144
  level =  3, nx =   16, residual change:    0.877452 →    0.256041
  level =  2, nx =    8, residual change:    0.180588 → 0.000272285
  level =  1, nx =    4, residual change: 0.000187447 → 6.55867e-17
  bottom solve
  level =  1, nx =    4, residual change: 1.19217e-16 → 1.19217e-16
  level =  2, nx =    8, residual change: 0.000298554 → 0.000298554
  level =  3, nx =   16, residual change:     0.33718 →     0.33718
  level =  4, nx =   32, residual change:     1.72809 →     1.72809
  level =  5, nx =   64, residual change:     2.60971 →     2.60971
  level =  6, nx =  128, residual change:     2.89721 →     2.89721
cycle 2: relative err = 2.1803634390217064, residual err = 0.0006486396426301177

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

  level =  6, nx =  128, residual change:   0.0527907 →   0.0515129
  level =  5, nx =   64, residual change:   0.0364241 →   0.0468113
  level =  4, nx =   32, residual change:    0.033097 →   0.0318323
  level =  3, nx =   16, residual change:   0.0224975 →  0.00656631
  level =  2, nx =    8, residual change:  0.00463131 → 6.95479e-06
  level =  1, nx =    4, residual change:  4.7921e-06 → 3.47155e-16
  bottom solve
  level =  1, nx =    4, residual change: 6.30821e-16 → 6.30821e-16
  level =  2, nx =    8, residual change: 7.63309e-06 → 7.63309e-06
  level =  3, nx =   16, residual change:  0.00864876 →  0.00864876
  level =  4, nx =   32, residual change:   0.0442789 →   0.0442789
  level =  5, nx =   64, residual change:   0.0665472 →   0.0665472
  level =  6, nx =  128, residual change:   0.0736819 →   0.0736819
cycle 3: relative err = 0.04844393523115633, residual err = 1.659245815001406e-05

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

  level =  6, nx =  128, residual change:  0.00135041 →  0.00131762
  level =  5, nx =   64, residual change: 0.000931668 →  0.00119765
  level =  4, nx =   32, residual change: 0.000846751 → 0.000816239
  level =  3, nx =   16, residual change: 0.000576837 → 0.000168502
  level =  2, nx =    8, residual change:  0.00011884 → 1.78399e-07
  level =  1, nx =    4, residual change: 1.22925e-07 →  2.1827e-16
  bottom solve
  level =  1, nx =    4, residual change: 3.96622e-16 → 3.96622e-16
  level =  2, nx =    8, residual change: 1.95801e-07 → 1.95801e-07
  level =  3, nx =   16, residual change: 0.000221902 → 0.000221902
  level =  4, nx =   32, residual change:   0.0011347 →   0.0011347
  level =  5, nx =   64, residual change:  0.00170278 →  0.00170278
  level =  6, nx =  128, residual change:  0.00188597 →  0.00188597
cycle 4: relative err = 0.0012759605329324085, residual err = 4.2728946362388976e-07

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

  level =  6, nx =  128, residual change: 3.47757e-05 → 3.38943e-05
  level =  5, nx =   64, residual change: 2.39659e-05 →   3.073e-05
  level =  4, nx =   32, residual change:  2.1726e-05 → 2.09217e-05
  level =  3, nx =   16, residual change: 1.47845e-05 → 4.32098e-06
  level =  2, nx =    8, residual change: 3.04737e-06 → 4.57173e-09
  level =  1, nx =    4, residual change: 3.15043e-09 → 4.90471e-17
  bottom solve
  level =  1, nx =    4, residual change: 8.91242e-17 → 8.91242e-17
  level =  2, nx =    8, residual change: 5.01821e-09 → 5.01821e-09
  level =  3, nx =   16, residual change: 5.68972e-06 → 5.68972e-06
  level =  4, nx =   32, residual change: 2.90707e-05 → 2.90707e-05
  level =  5, nx =   64, residual change: 4.36992e-05 → 4.36992e-05
  level =  6, nx =  128, residual change: 4.85557e-05 → 4.85557e-05
cycle 5: relative err = 3.301203447716335e-05, residual err = 1.1068868945958364e-08

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

  level =  6, nx =  128, residual change:  9.0086e-07 → 8.76274e-07
  level =  5, nx =   64, residual change: 6.19593e-07 → 7.90474e-07
  level =  4, nx =   32, residual change: 5.58852e-07 → 5.36009e-07
  level =  3, nx =   16, residual change: 3.78756e-07 → 1.10732e-07
  level =  2, nx =    8, residual change: 7.80911e-08 → 1.17095e-10
  level =  1, nx =    4, residual change: 8.06977e-11 → 5.69861e-16
  bottom solve
  level =  1, nx =    4, residual change:  1.0355e-15 →  1.0355e-15
  level =  2, nx =    8, residual change: 1.28541e-10 → 1.28541e-10
  level =  3, nx =   16, residual change: 1.45795e-07 → 1.45795e-07
  level =  4, nx =   32, residual change:  7.4452e-07 →  7.4452e-07
  level =  5, nx =   64, residual change: 1.12439e-06 → 1.12439e-06
  level =  6, nx =  128, residual change: 1.25658e-06 → 1.25658e-06
cycle 6: relative err = 8.544249588823554e-07, residual err = 2.88200772432267e-10

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

  level =  6, nx =  128, residual change: 2.34558e-08 → 2.27531e-08
  level =  5, nx =   64, residual change:  1.6088e-08 → 2.03781e-08
  level =  4, nx =   32, residual change: 1.44068e-08 → 1.37252e-08
  level =  3, nx =   16, residual change: 9.69812e-09 → 2.83563e-09
  level =  2, nx =    8, residual change: 1.99971e-09 → 2.99732e-12
  level =  1, nx =    4, residual change: 2.06576e-12 → 1.74907e-17
  bottom solve
  level =  1, nx =    4, residual change: 3.17826e-17 → 3.17826e-17
  level =  2, nx =    8, residual change: 3.29051e-12 → 3.29051e-12
  level =  3, nx =   16, residual change: 3.73325e-09 → 3.73325e-09
  level =  4, nx =   32, residual change: 1.90594e-08 → 1.90594e-08
  level =  5, nx =   64, residual change: 2.89959e-08 → 2.89959e-08
  level =  6, nx =  128, residual change: 3.26638e-08 → 3.26638e-08
cycle 7: relative err = 2.210681933627904e-08, residual err = 7.534885150074738e-12

Visualizing the solution#

[14]:

v = mg.get_solution()

[15]:

fig, ax = plt.subplots()

im = ax.imshow(np.transpose(v.v()),
              interpolation="nearest", origin="lower",
              extent=[mg.xmin, mg.xmax, mg.ymin, mg.ymax])
fig.colorbar(im, ax=ax)

[15]:

<matplotlib.colorbar.Colorbar at 0x7fa45d0a1090>

_images/multigrid-variable-coeff_25_1.png

Comparing to the exact solution#

[16]:

phi = true(mg.x2d, mg.y2d)

With periodic BCs all around, there is nothing to normalize the solution, so we subtract off the average of $ϕ$ from the MG solution to ensure it is normalized (we’ll do the same with the true solution, just to be sure)

[17]:

e = v - np.sum(v.v()) / N**2 - (phi - np.sum(phi[mg.ilo:mg.ihi+1, mg.jlo:mg.jhi+1]) / N**2)

Now we can look at the norm of the error:

[18]:

error_norm = e.norm()
print(f"error = {error_norm:20.10g}")

error =      9.754984685e-05

and we can plot the error

[19]:

fig, ax = plt.subplots()

im = ax.imshow(np.transpose(e.v()),
               interpolation="nearest", origin="lower",
               extent=[mg.xmin, mg.xmax, mg.ymin, mg.ymax])
fig.colorbar(im, ax=ax)

[19]:

<matplotlib.colorbar.Colorbar at 0x7fa45d096f90>

_images/multigrid-variable-coeff_33_1.png

Variable Coefficient Poisson

Contents

Variable Coefficient Poisson#

Setting up the solver#

Solving the system#

Visualizing the solution#

Comparing to the exact solution#