Multiphysics solvers (direct, AMG, GMG, matrix-free)

Overview

Every solve(PhysicsType(), mesh, material; …) accepts a method that selects the solver. This tutorial shows the options and the matrix-free geometric multigrid that applies the weak form on the fly.

using Rasmah

# 1. the geometry, 2. exact BRep + face selection, 3. mesh
part = box(1.0, 1.0, 1.0)
b = brep(part)
left  = faces_in_box(b, [-0.6, -1.0, -1.0], [-0.4, 1.0, 1.0])   # clamp
right = faces_in_box(b, [ 0.4, -1.0, -1.0], [ 0.6, 1.0, 1.0])   # pull
mesh = tetrahedralize(b, 6)

TetMesh{Matrix{Float64}, Matrix{Int64}, Rasmah.BRepFaceLabels}([-0.5 -0.33333333333333337 … 0.33333333333333337 0.5; -0.5 -0.5 … 0.5 0.5; -0.5 -0.5 … 0.5 0.5], [1 1 … 286 286; 2 9 … 335 336; 9 8 … 336 287; 58 58 … 343 343], Rasmah.BRepFaceLabels(FaceSet[FaceSet([19 17 … 23 8; 26 24 … 24 9; 27 25 … 31 16]), FaceSet([295 331 … 300 307; 302 332 … 301 314; 303 339 … 308 315]), FaceSet([1 52 … 246 3; 2 101 … 247 52; 51 102 … 296 53]), FaceSet([289 141 … 96 244; 290 190 … 97 245; 339 191 … 146 294]), FaceSet([127 232 … 15 253; 176 239 … 22 260; 183 288 … 71 309]), FaceSet([238 7 … 91 140; 287 14 … 140 147; 294 63 … 147 196])]))

Linear solvers

u_direct = solve(Elasticity(), mesh, steel;
                 fixed = left => [0.0, 0.0, 0.0], loads = right => [1e6, 0.0, 0.0])   # sparse Cholesky
u_lu     = solve(Elasticity(), mesh, steel;
                 fixed = left => [0.0, 0.0, 0.0], loads = right => [1e6, 0.0, 0.0], method = :lu) # sparse LU (non-SPD)
u_amg    = solve(Elasticity(), mesh, steel;
                 fixed = left => [0.0, 0.0, 0.0], loads = right => [1e6, 0.0, 0.0], method = :amg) # algebraic multigrid + CG
u_gmg    = solve(Elasticity(), mesh, steel;
                 fixed = left => [0.0, 0.0, 0.0], loads = right => [1e6, 0.0, 0.0], method = :gmg, levels = 3)       # geometric MG
u_mf     = solve(Elasticity(), mesh, steel;
                 fixed = left => [0.0, 0.0, 0.0], loads = right => [1e6, 0.0, 0.0], method = :gmg_mf, levels = 3)   # matrix-free GMG

println(isapprox(u_direct, u_amg; atol = 1e-6), " ",
        isapprox(u_direct, u_gmg; atol = 1e-6), " ",
        isapprox(u_direct, u_mf; atol = 1e-6))

true true true

:direct / :lu are sparse direct factorizations.
:amg preconditioned CG; :gmg box-coarsened geometric multigrid (the axis-aligned box BRep meshes to a uniform tetrahedralize_box lattice, which :gmg needs); :gmg_mf applies every level’s operator from the weak form without assembling it.

Matrix-free operators

The weak form itself supports an on-the-fly matrix–vector product and diagonal — the building blocks of matrix-free GMG:

U = FESpace(mesh, ReferenceFE(3, 1))
a = BilinearForm(U, U, heat_integrand(thermal_conductivity(steel)))
v = rand(ndofs(U))
@assert matvec(a, v) ≈ assemble_matrix(a) * v
@assert diagonal(a) ≈ [assemble_matrix(a)[i, i] for i in 1:ndofs(U)]
println("matvec and diagonal match the assembled operator")

matvec and diagonal match the assembled operator

Nonlinear solvers

Nonlinear physics are solved with Newton by default (or Jacobian-free Newton–Krylov / Anderson). We solve on a small mesh with a compressible elastomer (the compressible Neo-Hookean model needs a material with Poisson’s ratio well below 0.5):

elastomer = Material(name="Compressible elastomer",
                     youngs_modulus=1.0e6, poisson_ratio=0.3)
mh = tetrahedralize(b, 2)
u_newton = solve(HyperElasticity(), mh, elastomer;
                 fixed = [left => [0.0, 0.0, 0.0], right => [0.1, 0.0, 0.0]])
u_jfnk   = solve(HyperElasticity(), mh, elastomer;
                 fixed = [left => [0.0, 0.0, 0.0], right => [0.1, 0.0, 0.0]],
                 method = :newton_krylov)
println(isapprox(u_newton, u_jfnk; atol = 1e-4))

true

The solver choices are documented on each solve method; see also geometric-multigrid and sparse-linear-algebra.