Stokes eigeproblem MWE

stokes-eigenproblem-mwe.py

import numpy as np
import firedrake
from firedrake import inner, grad, div, dx, FiniteElement
from petsc4py import PETSc

nx, ny = 16, 16
mesh = firedrake.UnitSquareMesh(nx, ny, diagonal="crossed")

cg2 = FiniteElement("CG", "triangle", 2)
cg1 = FiniteElement("CG", "triangle", 1)
V = firedrake.VectorFunctionSpace(mesh, cg2)
Q = firedrake.FunctionSpace(mesh, cg1)
Z = V * Q

u, p = firedrake.TrialFunctions(Z)
v, q = firedrake.TestFunctions(Z)

A = -(inner(u, v) + inner(grad(u), grad(v)) - p * div(v) - q * div(u)) * dx
M = p * q * dx

problem = firedrake.LinearEigenproblem(A, M, restrict=False)

opts = PETSc.Options()
solver_type = opts.getString("solver")

if solver_type == "lu":
    solver_params = {
        "st_ksp_type": "gmres",
        "st_pc_type": "lu",
        "st_pc_factor_mat_solver_type": "mumps",
    }
elif solver_type == "fieldsplit":
    solver_params = {
        "st_ksp_type": "fgmres",
        "st_pc_type": "fieldsplit",
        "st_pc_fieldsplit_type": "schur",
        "st_fieldsplit_0": {
            "ksp_type": "preonly",
            "pc_type": "lu",
        },
        "st_fieldsplit_1": {
            "ksp_type": "cg",
            "ksp_rtol": 1e-8,
            "pc_type": "none",
        },
    }
else:
    raise ValueError("solver type must be either `lu` or `fieldsplit`!")

eps_params = {
    "eps_tol": 1e-8,
    "eps_gen_hermitian": None,
    "eps_target_real": None,
    "eps_smallest_magnitude": None,
    "st_type": "sinvert",
}

params = {
    "n_evals": 1, "solver_parameters": dict(**eps_params, **solver_params)
}
solver = firedrake.LinearEigensolver(problem, **params)
num_converged = solver.solve()
print(f"LBB constant: {np.sqrt(solver.eigenvalue(0).real):.2f}")