Diffusion operator

This is an example of various preconditioners for the 2D scalar diffusion problem.

Problem formulation

The scalar diffusion problem is given by

\[\nabla^2 u = f\]

with a weak form of

\[\int_\Omega \nabla u \nabla v = \int_\Omega v f, \forall v \in V\]

for an appropriate test space $V \subseteq H_0^1 \left( \Omega \right)$ on the domain. In this weak formulation, boundary terms have been omitted, as they are not present on the infinite grid for Local Fourier Analysis.

Plot for the symbol of the finite element operator for the 2D scalar diffusion problem with cubic basis.

LFAToolkit code

The diffusion operator is the classical test case for multigrid preconditioners and solvers.

Jacobi

This is an example of a Jacobi smoother.

# ------------------------------------------------------------------------------
# Jacobi smoother example
# ------------------------------------------------------------------------------

using LFAToolkit
using LinearAlgebra

# setup
mesh = Mesh2D(1.0, 1.0)
p = 3

# diffusion operator
diffusion = GalleryOperator("diffusion", p + 1, p + 1, mesh)

# Jacobi smoother
jacobi = Jacobi(diffusion)

# compute operator symbols
A = computesymbols(jacobi, [1.0], [π, π])
eigenvalues = real(eigvals(A))

# ------------------------------------------------------------------------------

Plot for the symbol of a Jacobi smoother for the 2D scalar diffusion problem with cubic basis.

Chebyshev

This is an example of a Chebyshev smoother.

# ------------------------------------------------------------------------------
# Chebyshev smoother example
# ------------------------------------------------------------------------------

using LFAToolkit
using LinearAlgebra

# setup
mesh = Mesh2D(1.0, 1.0)
p = 3

# diffusion operator
diffusion = GalleryOperator("diffusion", p + 1, p + 1, mesh)

# Chebyshev smoother
chebyshev = Chebyshev(diffusion)

# compute operator symbols
A = computesymbols(chebyshev, [3], [π, π])
eigenvalues = real(eigvals(A))

# ------------------------------------------------------------------------------

Plot for the symbol of a cubic Chebyshev smoother for the 2D scalar diffusion problem with cubic basis.

P-Multigrid

This is an example of a two level p-multigrid preconditioner.

# ------------------------------------------------------------------------------
# p-multigrid example
# ------------------------------------------------------------------------------

using LFAToolkit
using LinearAlgebra

# setup
mesh = Mesh2D(1.0, 1.0)
finep = 2
coarsep = 1
numbercomponents = 1
dimension = 2
ctofbasis = TensorH1LagrangeBasis(
    coarsep + 1,
    finep + 1,
    numbercomponents,
    dimension,
    collocatedquadrature = true,
)

# diffusion operators
finediffusion = GalleryOperator("diffusion", finep + 1, finep + 1, mesh)
coarsediffusion = GalleryOperator("diffusion", coarsep + 1, finep + 1, mesh)

# Chebyshev smoother
chebyshev = Chebyshev(finediffusion)

# p-multigrid preconditioner
multigrid = PMultigrid(finediffusion, coarsediffusion, chebyshev, [ctofbasis])

# compute operator symbols
A = computesymbols(multigrid, [3], [1, 1], [π, π])
eigenvalues = real(eigvals(A))

# ------------------------------------------------------------------------------

Example plot for the symbol of p-multigrid with a cubic Chebyshev smoother for the 2D scalar diffusion problem with cubic basis.

This is an example of a multilevel p-multigrid preconditioner.

# ------------------------------------------------------------------------------
# p-multigrid multilevel example
# ------------------------------------------------------------------------------

using LFAToolkit
using LinearAlgebra

# setup
mesh = Mesh2D(1.0, 1.0)
finep = 4
midp = 2
coarsep = 1
numbercomponents = 1
dimension = 2
ctombasis =
    TensorH1LagrangePProlongationBasis(coarsep + 1, midp + 1, numbercomponents, dimension)
mtofbasis =
    TensorH1LagrangePProlongationBasis(midp + 1, finep + 1, numbercomponents, dimension)

# diffusion operators
finediffusion = GalleryOperator("diffusion", finep + 1, finep + 1, mesh)
middiffusion = GalleryOperator("diffusion", midp + 1, finep + 1, mesh)
coarsediffusion = GalleryOperator("diffusion", coarsep + 1, finep + 1, mesh)

# Chebyshev smoothers
finechebyshev = Chebyshev(finediffusion)
midchebyshev = Chebyshev(middiffusion)

# p-multigrid preconditioner
midmultigrid = PMultigrid(middiffusion, coarsediffusion, midchebyshev, [ctombasis])
multigrid = PMultigrid(finediffusion, midmultigrid, finechebyshev, [mtofbasis])

# compute operator symbols
A = computesymbols(multigrid, [3], [1, 1], [π, π])
eigenvalues = real(eigvals(A))

# ------------------------------------------------------------------------------

Example plot for the symbol of multilevel p-multigrid with a cubic Chebyshev smoother for the 2D scalar diffusion problem with cubic basis.

H-Multigrid

This is an example of a two level h-multigrid preconditioner.

# ------------------------------------------------------------------------------
# h-multigrid example
# ------------------------------------------------------------------------------

using LFAToolkit
using LinearAlgebra

# setup
mesh = Mesh2D(1.0, 1.0)
p = 1
numbercomponents = 1
numberfineelements1d = 2
dimension = 2
ctofbasis = TensorH1LagrangeHProlongationBasis(
    p + 1,
    numbercomponents,
    dimension,
    numberfineelements1d,
);

# operators
finediffusion =
    GalleryMacroElementOperator("diffusion", p + 1, p + 2, numberfineelements1d, mesh);
coarsediffusion = GalleryOperator("diffusion", p + 1, p + 2, mesh);

# Chebyshev smoother
chebyshev = Chebyshev(finediffusion)

# h-multigrid preconditioner
multigrid = HMultigrid(finediffusion, coarsediffusion, chebyshev, [ctofbasis])

# compute operator symbols
A = computesymbols(multigrid, [3], [1, 1], [π, π])
eigenvalues = real(eigvals(A))

# ------------------------------------------------------------------------------

Example plot for the symbol of h-multigrid with a cubic Chebyshev smoother for the 2D scalar diffusion problem with linear basis.

This is an example of a multilevel h-multigrid preconditioner.

# ------------------------------------------------------------------------------
# h-multigrid multilevel example
# ------------------------------------------------------------------------------

using LFAToolkit
using LinearAlgebra

# setup
mesh = Mesh2D(1.0, 1.0)
p = 1
numberfineelements1d = 4
numbermidelements1d = 2
numbercomponents = 1
dimension = 2
ctombasis = TensorH1LagrangeHProlongationBasis(
    p + 1,
    numbercomponents,
    dimension,
    numbermidelements1d,
);
mtofbasis = TensorH1LagrangeHProlongationMacroBasis(
    p + 1,
    numbercomponents,
    dimension,
    numbermidelements1d,
    numberfineelements1d,
);

# operators
finediffusion =
    GalleryMacroElementOperator("diffusion", p + 1, p + 2, numberfineelements1d, mesh);
middiffusion =
    GalleryMacroElementOperator("diffusion", p + 1, p + 2, numbermidelements1d, mesh);
coarsediffusion = GalleryOperator("diffusion", p + 1, p + 2, mesh);

# Chebyshev smoothers
finechebyshev = Chebyshev(finediffusion)
midchebyshev = Chebyshev(middiffusion)

# h-multigrid preconditioner
midmultigrid = HMultigrid(middiffusion, coarsediffusion, midchebyshev, [ctombasis])
multigrid = HMultigrid(finediffusion, midmultigrid, finechebyshev, [mtofbasis])

# compute operator symbols
A = computesymbols(multigrid, [3], [1, 1], [π, π])
eigenvalues = real(eigvals(A))

# ------------------------------------------------------------------------------

Example plot for the symbol of multilevel h-multigrid with a cubic Chebyshev smoother for the 2D scalar diffusion problem with linear basis.

Lumped BDDC

This is an example of a simple lumped BDDC preconditioner.

# ------------------------------------------------------------------------------
# lumped BDDC example
# ------------------------------------------------------------------------------

using LFAToolkit
using LinearAlgebra

# setup
mesh = Mesh2D(1.0, 1.0)
p = 3

# diffusion operator
diffusion = GalleryOperator("diffusion", p + 1, p + 1, mesh)

# lumped BDDC preconditioner
bddc = LumpedBDDC(diffusion)

# compute operator symbols
A = computesymbols(bddc, [0.2], [π, π])
eigenvalues = real(eigvals(A))

# ------------------------------------------------------------------------------

This is an example of a simple Dirichlet BDDC preconditioner.

# ------------------------------------------------------------------------------
# Dirichlet BDDC example
# ------------------------------------------------------------------------------

using LFAToolkit
using LinearAlgebra

# setup
mesh = Mesh2D(1.0, 1.0)
p = 3

# diffusion operator
diffusion = GalleryOperator("diffusion", p + 1, p + 1, mesh)

# Dirichlet BDDC preconditioner
bddc = DirichletBDDC(diffusion)

# compute operator symbols
A = computesymbols(bddc, [0.2], [π, π])
eigenvalues = real(eigvals(A))

# ------------------------------------------------------------------------------