MixedElement functionspace, DirichletBC with Function as value

dolfinx
1

Hi all,

I’m trying to solve a problem involving a mixed element and apply a DirichletBC with a function defined on the collapsed functionspace.

It works fine if instead I use a Constant as a value but when trying to use the fem.Function the solver diverges.

Any pointers?

Thanks in advance!

MWE:

import dolfinx
import numpy as np
import ufl
from dolfinx.fem.petsc import NonlinearProblem
from petsc4py import PETSc
from mpi4py import MPI
from dolfinx import mesh
from dolfinx import fem
import basix.ufl
from basix.ufl import element

domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)


cg_element = element("Lagrange", domain.basix_cell(), degree=1)

mixed_element = basix.ufl.mixed_element([cg_element, cg_element])

V = fem.functionspace(domain, mixed_element)

u = fem.Function(V)

cm, ct = ufl.split(u)

cm_post, ct_post = u.split()

tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)


def boundary_left(x):
    return np.isclose(x[0], 0)


def boundary_right(x):
    return np.isclose(x[0], 1)


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


V0, submap = V.sub(0).collapse()


dofs_right = fem.locate_dofs_geometrical((V.sub(0), V0), boundary_right)
dofs_left = fem.locate_dofs_geometrical((V.sub(0), V0), boundary_left)
dofs_top = fem.locate_dofs_geometrical((V.sub(0), V0), boundary_top)

# FIXME this makes the solver diverge (reason -5)
uD = fem.Function(V0)
uD.interpolate(lambda x: 1 + x[0] ** 2 + 2 * x[1] ** 2)
bc_left = fem.dirichletbc(uD, dofs_left[0])

# this works
# bc_left = fem.dirichletbc(fem.Constant(domain, 2.0), dofs_left[0], V.sub(0))
bc_right = fem.dirichletbc(fem.Constant(domain, 1.0), dofs_right[0], V.sub(0))

v_cm, v_ct = ufl.TestFunctions(V)

F1 = ufl.dot(ufl.grad(cm), ufl.grad(v_cm)) * ufl.dx
F2 = ufl.dot(ufl.grad(ct), ufl.grad(v_ct)) * ufl.dx

F = F1 + F2

# taken from https://github.com/FEniCS/dolfinx/blob/5fcb988c5b0f46b8f9183bc844d8f533a2130d6a/python/demo/demo_cahn-hilliard.py#L279C1-L286C28
use_superlu = PETSc.IntType == np.int64  # or PETSc.ScalarType == np.complex64
sys = PETSc.Sys()  # type: ignore
if sys.hasExternalPackage("mumps") and not use_superlu:
    linear_solver = "mumps"
elif sys.hasExternalPackage("superlu_dist"):
    linear_solver = "superlu_dist"
else:
    linear_solver = "petsc"

petsc_options = {
    "snes_type": "newtonls",
    "snes_linesearch_type": "none",
    "snes_stol": np.sqrt(np.finfo(dolfinx.default_real_type).eps) * 1e-2,
    "snes_atol": 1e-10,
    "snes_rtol": 1e-10,
    "snes_max_it": 100,
    "snes_divergence_tolerance": "PETSC_UNLIMITED",
    "ksp_type": "preonly",
    "pc_type": "lu",
    "pc_factor_mat_solver_type": linear_solver,
}

problem = NonlinearProblem(
    F,
    u,
    bcs=[bc_left, bc_right],
    petsc_options=petsc_options,
    petsc_options_prefix="Poisson",
)
problem.solve()
converged = problem.solver.getConvergedReason()
num_iter = problem.solver.getIterationNumber()
assert converged > 0, f"Solver did not converge, got {converged}."
print(
    f"Solver converged after {num_iter} iterations with converged reason {converged}."
)

Produces:

    assert converged > 0, f"Solver did not converge, got {converged}."
           ^^^^^^^^^^^^^
AssertionError: Solver did not converge, got -5.
2

The correct line is:

bc_left = fem.dirichletbc(uD, dofs_left,  V.sub(0))

The point here is is that your function uD is a member of the collapsed (stand-alone) space, while you’re aiming to set a dirichlet contraint on a function in the mixedspace. In order to do that, FEniCS needs to know how to connect which DoFs to which DoFs. That is why you need to pass both arrays in dofs_left (you’ll see dofs_left is a list with two arrays, one of dofs of V_sub(0) and one of dofs of V0).

1 Like