Nest data structure and dolfinx.fem.petsc.LinearProblem

dolfinx
1

I have two issues:

Using dolfinx.fem.petsc.LinearProblem with kind=”nest” and explicitly specifying the u argument, like

problem = dolfinx.fem.petsc.LinearProblem(
    a_ufl,
    L_ufl,
    u = [u_h, p_h],
    kind="nest",
    bcs=bcs,
    P=a_p,
    petsc_options_prefix="demo_stokes__nested_iterative_solver_high_level_",
    petsc_options={
        "ksp_type": "minres",
        "ksp_rtol": 1e-9,
        "pc_type": "fieldsplit",
        "pc_fieldsplit_type": "additive",
        "fieldsplit_0_ksp_type": "preonly",
        "fieldsplit_0_pc_type": "gamg",
        "fieldsplit_1_ksp_type": "preonly",
        "fieldsplit_1_pc_type": "jacobi",
    },
)

does not work (nest is not used) if the Functions are defined with a name (e.g.

u_h, p_h = dolfinx.fem.Function(V, name='u'), dolfinx.fem.Function(Q, name='p')

), while it works if the optional name argument is not used (e.g.

u_h, p_h = dolfinx.fem.Function(V), Function(Q)

). The following MRE, adapter from the Stokes equation demo,

from mpi4py import MPI
from petsc4py import PETSc

import numpy as np

import ufl
import basix
import dolfinx
import dolfinx.fem.petsc

msh = dolfinx.mesh.create_rectangle(
    MPI.COMM_WORLD, [np.array([0, 0]), np.array([1, 1])], (32, 32), dolfinx.mesh.CellType.triangle
)


def noslip_boundary(x):
    return np.isclose(x[0], 0.0) | np.isclose(x[0], 1.0) | np.isclose(x[1], 0.0)


def lid(x):
    return np.isclose(x[1], 1.0)


def lid_velocity_expression(x):
    return np.stack((np.ones(x.shape[1]), np.zeros(x.shape[1])))


P2 = basix.ufl.element(
    "Lagrange", msh.basix_cell(), degree=2, shape=(msh.geometry.dim,), dtype=dolfinx.default_real_type
)
P1 = basix.ufl.element("Lagrange", msh.basix_cell(), degree=1, dtype=dolfinx.default_real_type)
V, Q = dolfinx.fem.functionspace(msh, P2), dolfinx.fem.functionspace(msh, P1)

noslip = np.zeros(msh.geometry.dim, dtype=PETSc.ScalarType)  # type: ignore
facets = dolfinx.mesh.locate_entities_boundary(msh, 1, noslip_boundary)
bc0 = dolfinx.fem.dirichletbc(noslip, dolfinx.fem.locate_dofs_topological(V, 1, facets), V)

lid_velocity = dolfinx.fem.Function(V)
lid_velocity.interpolate(lid_velocity_expression)
facets = dolfinx.mesh.locate_entities_boundary(msh, 1, lid)
bc1 = dolfinx.fem.dirichletbc(lid_velocity, dolfinx.fem.locate_dofs_topological(V, 1, facets))

bcs = [bc0, bc1]

(u, p) = ufl.TrialFunction(V), ufl.TrialFunction(Q)
(v, q) = ufl.TestFunction(V), ufl.TestFunction(Q)
f = dolfinx.fem.Constant(msh, (PETSc.ScalarType(0), PETSc.ScalarType(0)))  # type: ignore

a_ufl = [
    [ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx, ufl.inner(p, ufl.div(v)) * ufl.dx],
    [ufl.inner(ufl.div(u), q) * ufl.dx, None],
]


a = dolfinx.fem.form(a_ufl)
L_ufl = [ufl.inner(f, v) * ufl.dx, ufl.ZeroBaseForm((q,))]
L = dolfinx.fem.form(L_ufl)

a_p11 = dolfinx.fem.form(ufl.inner(p, q) * ufl.dx)
a_p = [[a[0][0], None], [None, a_p11]]

# prevent nesting
u_h, p_h = dolfinx.fem.Function(V, name='u'), dolfinx.fem.Function(Q, name='p')
# does not prevent nesting
# u_h, p_h = dolfinx.fem.Function(V), Function(Q)


problem = dolfinx.fem.petsc.LinearProblem(
    a_ufl,
    L_ufl,
    u = [u_h, p_h],
    kind="nest",
    bcs=bcs,
    P=a_p,
    petsc_options_prefix="demo_stokes__nested_iterative_solver_high_level_",
    petsc_options={
        "ksp_type": "minres",
        "ksp_rtol": 1e-9,
        "pc_type": "fieldsplit",
        "pc_fieldsplit_type": "additive",
        "fieldsplit_0_ksp_type": "preonly",
        "fieldsplit_0_pc_type": "gamg",
        "fieldsplit_1_ksp_type": "preonly",
        "fieldsplit_1_pc_type": "jacobi",
    },
)

null_vec = dolfinx.fem.petsc.create_vector(dolfinx.fem.extract_function_spaces(L), "nest")

null_vecs = null_vec.getNestSubVecs()
null_vecs[0].set(0.0), null_vecs[1].set(1.0)

null_vec.normalize()
nsp = PETSc.NullSpace().create(vectors=[null_vec])
problem.A.setNullSpace(nsp)

A00 = problem.A.getNestSubMatrix(0, 0)
A00.setOption(PETSc.Mat.Option.SPD, True)

P00, P11 = problem.P_mat.getNestSubMatrix(0, 0), problem.P_mat.getNestSubMatrix(1, 1)
P00.setOption(PETSc.Mat.Option.SPD, True)
P11.setOption(PETSc.Mat.Option.SPD, True)

problem.solve()

leads to

WARNING! There are options you set that were not used!
WARNING! could be spelling mistake, etc!
There are 4 unused database options. They are:
Option left: name:-demo_stokes__nested_iterative_solver_high_level_fieldsplit_0_ksp_type value: preonly source: code
Option left: name:-demo_stokes__nested_iterative_solver_high_level_fieldsplit_0_pc_type value: gamg source: code
Option left: name:-demo_stokes__nested_iterative_solver_high_level_fieldsplit_1_ksp_type value: preonly source: code
Option left: name:-demo_stokes__nested_iterative_solver_high_level_fieldsplit_1_pc_type value: jacobi source: code

No warning is issued if the functions are defined without a name, as in the commented line of the MRE.

LinearProblem takes a list of bcs, also for the case of kind=”nest”. I fear that this could lead to issues (for example, if two blocks are built over the same vector space and more than one bc needs to be applied to one of them). Am I wrong and/or misunderstanding something? Same issue with bcs_by_block?

Thank you in advance,

Marco

2

Some more data: the first issue seems to be because of these lines starting from python/dolfinx/fem/petsc.py:871

        fieldsplit_IS = tuple(
            [
                (f"{u.name + '_' if u.name != 'f' else ''}{i}", IS)
                for i, (u, IS) in enumerate(zip(self.u, nest_IS[0]))
            ]
        )

because of these, if the Functions have a name, we end up with

(('u_0', <petsc4py.PETSc.IS object at 0x7fe9f780abb0>), ('p_1', <petsc4py.PETSc.IS object at 0x7fe9f780ac00>))

for fieldsplit_IS rather than

(('0', <petsc4py.PETSc.IS object at 0x7f05ee7eaa70>), ('1', <petsc4py.PETSc.IS object at 0x7f05ee7eaac0>))

and from there this leads to

Option left: name:-demo_stokes__nested_iterative_solver_high_level_fieldsplit_0_ksp_type value: preonly source: code
Option left: name:-demo_stokes__nested_iterative_solver_high_level_fieldsplit_0_pc_type value: gamg source: code
Option left: name:-demo_stokes__nested_iterative_solver_high_level_fieldsplit_1_ksp_type value: preonly source: code
Option left: name:-demo_stokes__nested_iterative_solver_high_level_fieldsplit_1_pc_type value: jacobi source: code

thus I should have written

problem = LinearProblem.LinearProblem(
    a_ufl,
    L_ufl,
    u = [u_h, p_h],
    kind="nest",
    bcs=bcs,
    P=a_p,
    petsc_options_prefix="demo_stokes__nested_iterative_solver_high_level_",
    petsc_options={
        "ksp_type": "minres",
        "ksp_rtol": 1e-9,
        "pc_type": "fieldsplit",
        "pc_fieldsplit_type": "additive",
        "fieldsplit_u_0_ksp_type": "preonly",
        "fieldsplit_u_0_pc_type": "gamg",
        "fieldsplit_p_1_ksp_type": "preonly",
        "fieldsplit_p_1_pc_type": "jacobi",
    },
)

although I find this somewhat counter-intuitive (at this point, why `u_0` and not u since it has an explicit name?).

As for the second issue, I’m still convinced that the interface is wrong, i.e. that LinearProblem should explicitly get a list of list of bcs, and that bcs_by_block should perhaps be nuked and not used within LinearProblem. Am I wrong?

3

Change the options to:

        f"fieldsplit_{u_h.name}_ksp_type": "preonly",
        f"fieldsplit_{u_h.name}_pc_type": "gamg",
        f"fieldsplit_{p_h.name}_ksp_type": "preonly",
        f"fieldsplit_{p_h.name}_pc_type": "jacobi",

If you want to use the “same” function space for multiple blocks, you need to use dolfinx.fem.FunctionSpace.clone (dolfinx.fem — DOLFINx 0.10.0.post5 documentation)
i.e.

V0 = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
V1 = V0.clone()

otherwise there is no way for the block assembler to distinguish between what is the 0th, 1st, 2nd, 3rd block etc.