Multigrid Benchmarks

In multigrid methods, there are tradeoffs between the levels of coarsening and the number of cycles performed. This docs page provides a simple demonstrations of the relative error of different methods, and the time of execution for these methods. This can help you choose the parameters to be used by your multigrid solves.

This code demonstrates a simple pipeline of creating a mesh, executing a multigrid solve, and plotting statistics about those executions.

using CombinatorialSpaces
using CairoMakie
using GeometryBasics: Point3d
using LinearAlgebra: norm
const Point3D = Point3{Float64}

s = triangulated_grid(1, 1, 1/4, sqrt(3)/2*1/4, Point3D, false)

function plot_residuals(s, title; cycles=1:50, timeit=false)
  i = 0
  function residuals(s, scheme::AbstractSubdivisionScheme, level, sizes)
    series = PrimalGeometricMapSeries(s, scheme, level);
    md = MGData(series, sd -> ∇²(0, sd), 3, scheme)
    sd = finest_mesh(series)
    push!(sizes, nv(sd))
    L = first(md.operators)
    # Observe: putting into range of the Laplacian for solvability.
    b = L*rand(nv(sd))
    u0 = zeros(nv(sd))
    ress = map(cycles) do cyc
      if timeit
        i += 1
        @elapsed u = multigrid_vcycles(u0,b,md,cyc)
      else
        u = multigrid_vcycles(u0,b,md,cyc)
        norm(L*u-b)/norm(b)
      end
    end
  end
  bin_sizes = []
  bin_ress = map(lev -> residuals(s, BinarySubdivision(), lev, bin_sizes), 2:4)
  cub_sizes = []
  cub_ress = map(lev -> residuals(s, CubicSubdivision(), lev, cub_sizes), 2:4)
  f = Figure()
  ax = CairoMakie.Axis(f[1,1];
    title = "Multigrid V-cycles, $title",
    yscale = timeit ? identity : log10,
    ylabel = timeit ? "execution time [s]" : "log₁₀(relative error)",
    xlabel = "# cycles")
  colors = [:blue, :green, :orange]
  for lev in 1:3
    lines!(ax, bin_ress[lev],
      label="binary, $(lev+1) levels, $(bin_sizes[lev]) vertices",
      linestyle=:solid, color=colors[lev])
    lines!(ax, cub_ress[lev],
      label="cubic, $(lev+1) levels, $(cub_sizes[lev]) vertices",
      linestyle=:dash, color=colors[lev])
  end
  f[1,2] = Legend(f,ax,"Scheme")
  f
end

plot_residuals (generic function with 1 method)

plot_residuals(s, "Residuals")

We see that the log of the observed relative error is roughly in-line with the resolution of the base mesh used.

plot_residuals(s, "Times", timeit=true)

The spikes in these times are due to the fact that we are only executing the cycles a single time, for runtimes' sake. We observe that the timings are linear in the number of cycles. We observe that the slopes are greater for those meshes with more vertices, with the cubic-4-level method being the finest.