Solving Stokes with in a domain with high aspect-ratio

[miguelcoolchips]

I want to solve the Navier-Stokes equations on a 3D domain with dimensions 1 x 0.1 x 0.01 and using an iterative solver. I started first with the Stokes equations to make my life easier.

import firedrake as fd
from firedrake import inner, grad, div, dx
from firedrake.petsc import PETSc


simulation_case = "elongated"
# simulation_case = "cube"

if simulation_case == "elongated":
    Lx = 1
    Ly = 0.01
    Lz = 0.1
    nx = 100
    ny = int(nx * Ly / Lx * 5)
    nz = int(nx * Lz / Lx * 5)
else:
    Lx = 1
    Ly = 1
    Lz = 1
    nx = 20
    ny = nx
    nz = nx
mesh = fd.BoxMesh(nx, ny, nz, Lx, Ly, Lz)

V = fd.VectorFunctionSpace(mesh, "CG", 2)
W = fd.FunctionSpace(mesh, "CG", 1)
Z = V * W
PETSc.Sys.Print(f"DOFs: {Z.dim()}")

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


def a(u, v, p, q):
    return (inner(grad(u), grad(v)) - inner(p, div(v)) + inner(div(u), q)) * dx


L = inner(fd.Constant((0, 0, 0)), v) * dx

max_vel = 16
x, y, z = fd.SpatialCoordinate(mesh)
u_inlet = fd.as_vector(
    [
        max_vel * (1 - y / Ly) * y / Ly * z / Lz * (1 - z / Lz),
        0,
        0,
    ]
)


bcs = [
    fd.DirichletBC(Z.sub(0), u_inlet, (1,)),
    fd.DirichletBC(Z.sub(0), fd.Constant((0, 0, 0)), (3, 4, 5, 6)),
]

up = fd.Function(Z)
a_form = a(u, v, p, q)

direct_solver_settings = {
    "ksp_type": "preonly",
    "ksp_monitor_true_residual": None,
    "ksp_monitor": None,
    "pc_type": "lu",
    "pc_factor_mat_solver_type": "mumps",
    "mat_mumps_icntl_14": 200,
    "mat_mumps_icntl_24": 1,
}


class MassInvPC(fd.AssembledPC):

    _prefix = "Mp_"

    def form(self, pc, test, trial):
        _, bcs = super(MassInvPC, self).form(pc)

        appctx = self.get_appctx(pc)
        mu = appctx.get("mu", 1.0)
        a = inner((1 / mu) * trial, test) * dx
        return a, bcs

    def set_nullspaces(self, pc):
        # the mass matrix does not have a nullspace
        pass


assembled_hypre_settings = {
    "ksp_type": "cg",
    # "ksp_converged_reason": None,
    "ksp_max_it": 100,
    "ksp_atol": 1e-15,
    "ksp_rtol": 1e-15,
    "pc_type": "python",
    "pc_python_type": "firedrake.AssembledPC",
    "assembled_pc_type": "hypre",
    "assembled_pc_hypre_type": "boomeramg",
    "assembled_pc_hypre_boomeramg_coarsen_type": "PMIS",
    "assembled_pc_hypre_boomeramg_interp_type": "ext+i",
    "assembled_pc_hypre_boomeramg_strong_threshold": 0.1,
    "assembled_pc_hypre_boomeramg_agg_nl": 1,
    "assembled_pc_hypre_boomeramg_relax_type_all": "l1-Gauss-Seidel",
    "assembled_pc_hypre_boomeramg_smooth_type": "ILU",
    "assembled_pc_hypre_boomeramg_truncfactor": 0.3,
    "assembled_pc_hypre_boomeramg_max_levels": 25,
    "assembled_pc_hypre_boomeramg_max_iter": 1,
}

mass_inv_settings = {
    "ksp_type": "gmres",
    "ksp_converged_reason": None,
    "ksp_monitor_true_residual": None,
    "ksp_max_it": 100,
    "ksp_atol": 1e-12,
    "ksp_rtol": 1e-12,
    "pc_type": "python",
    "pc_python_type": "__main__.MassInvPC",
    "Mp_pc_type": "lu",
}
parameters_mass_inv = {
    "ksp_type": "fgmres",
    "ksp_atol": 1e-7,
    "ksp_rtol": 1e-6,
    "ksp_monitor_true_residual": None,
    "pc_type": "fieldsplit",
    "pc_fieldsplit_type": "schur",
    "pc_fieldsplit_schur_fact_type": "full",
    "pc_fieldsplit_schur_precondition": "selfp",
    "fieldsplit_0": assembled_hypre_settings,
    "fieldsplit_1": mass_inv_settings,
}

up.assign(0)
# NOT CONVERGING
fd.solve(a_form == L, up, bcs=bcs, solver_parameters=parameters_mass_inv)

I started with the options in https://www.firedrakeproject.org/demos/stokes.py.html, but these were not working. I ended up switching to using hypre for both fields in the fieldsplit preconditioner.

hypre_settings = {
    "ksp_type": "cg",
    "ksp_max_it": 100,
    "ksp_atol": 1e-15,
    "ksp_rtol": 1e-15,
    "pc_type": "hypre",
    "pc_hypre_type": "boomeramg",
    "pc_hypre_boomeramg_coarsen_type": "PMIS",
    "pc_hypre_boomeramg_interp_type": "ext+i",
    "pc_hypre_boomeramg_strong_threshold": 0.1,
    "pc_hypre_boomeramg_agg_nl": 1,
    "pc_hypre_boomeramg_relax_type_all": "l1-Gauss-Seidel",
    "pc_hypre_boomeramg_smooth_type": "ILU",
    "pc_hypre_boomeramg_truncfactor": 0.3,
    "pc_hypre_boomeramg_max_levels": 25,
    "pc_hypre_boomeramg_max_iter": 1,
}


parameters_working = {
    "ksp_type": "fgmres",
    "ksp_atol": 1e-7,
    "ksp_rtol": 1e-6,
    "ksp_monitor_true_residual": None,
    "pc_type": "fieldsplit",
    "pc_fieldsplit_type": "schur",
    "pc_fieldsplit_schur_fact_type": "full",
    "pc_fieldsplit_schur_precondition": "selfp",
    "fieldsplit_0": assembled_hypre_settings,
    "fieldsplit_1": hypre_settings,
}

up.assign(0)
fd.solve(a_form == L, up, bcs=bcs, solver_parameters=parameters_working)

This was working, although the preconditioned norm of the KSP is suspiciously high. Once I had a Stokes solver working, I moved to solving the Navier-Stokes equation using the PCD preconditioner (we have modified the PCD solver to account for inflow boundary conditions). However, passing "pc_fieldsplit_schur_precondition": "selfp" is incompatible with having "mat_type": "matfree" which is required by PCDPC. So then, I went back to the Stokes problem and tried to build my own selfp preconditioner:

class SchurApprox(fd.AuxiliaryOperatorPC):
    def form(self, pc, test, trial):
        appctx = self.get_appctx(pc)
        up = appctx["state"]
        V = up.function_space()
        Vrestric = fd.RestrictedFunctionSpace(V.sub(0), boundary_set=[1, 3, 4, 5, 6])
        V = Vrestric * V.sub(1)

        up = fd.TrialFunction(V)
        vq = fd.TestFunction(V)
        u, p = fd.split(up)
        v, q = fd.split(vq)
        a = (inner(grad(u), grad(v)) - inner(p, div(v)) + inner(div(u), q)) * dx
        A = fd.Tensor(a)
        A_blocks = A.blocks

        A00 = A_blocks[0, 0]
        A10 = A_blocks[1, 0]
        A01 = A_blocks[0, 1]

        D00 = fd.DiagonalTensor(A00)
        Sp = -A10 * D00.inv * A01
        return (Sp, None)


parameters = {
    "ksp_type": "fgmres",
    "ksp_atol": 1e-7,
    "ksp_rtol": 1e-6,
    "ksp_monitor_true_residual": None,
    "pc_type": "fieldsplit",
    "pc_fieldsplit_type": "schur",
    "pc_fieldsplit_schur_fact_type": "full",
    "pc_fieldsplit_schur_precondition": "user",
    "fieldsplit_0": direct_solver_settings,
    "fieldsplit_1": {
        "pc_type": "python",
        "pc_python_type": "__main__.SchurApprox",
        "aux": direct_solver_settings,
    },
}

up.assign(0)
# NOT CONVERGING
fd.solve(a_form == L, up, bcs=bcs, solver_parameters=parameters)

but this is not converging either. These are the things that I tried for this long thin domain. Are there better ideas that I can try?