Mixed formulation of Poisson problem in dolfinx

dolfinx
7

The problem with the code above is the usage of the default petsc “lu” solver.
Changing it to mumps or umfpack works:

from mpi4py import MPI
from dolfinx import fem, mesh
import basix.ufl
import numpy as np
import ufl
from dolfinx.fem.petsc import LinearProblem


domain = mesh.create_unit_square(comm=MPI.COMM_WORLD,
                                 nx=32, ny=32,
                                 cell_type=mesh.CellType.triangle)
el0 = basix.ufl.element("BDM", domain.topology.cell_name(),1)
el1 = basix.ufl.element("DG", domain.topology.cell_name(), 0)
W = fem.functionspace(mesh=domain, element=basix.ufl.mixed_element([el0, el1]))

# test and trial functions
tau, v = ufl.TestFunctions(function_space=W)
sigma, u = ufl.TrialFunctions(function_space=W)


# function determining if a node is on the tray top
def on_left_and_right(x):
    return np.isclose(x[0], 0) | np.isclose(x[0], 1)


W0, _ = W.sub(0).collapse()
u_zero = fem.Function(W0)
u_zero.x.array[:] = 0.0  # ----------------------------------------------------   Changed here
boundary_facets = mesh.locate_entities_boundary(mesh=domain,
                                                dim=domain.topology.dim-1,
                                                marker=on_left_and_right)
boundary_dofs = fem.locate_dofs_topological(V=(W.sub(0), W0),
                                            entity_dim=domain.topology.dim-1,
                                            entities=boundary_facets)

# apply Dirichlet BC
bc = fem.dirichletbc(value=u_zero,
                     dofs=boundary_dofs,
                     V=W.sub(0))


# apply Neumann BC
x = ufl.SpatialCoordinate(domain=domain)
f = 10.0 * ufl.exp(-((x[0] - 0.5)**2 + (x[1] - 0.5)**2) / 0.02)
g = ufl.sin(5 * x[0])

# Variational form
a = (ufl.dot(sigma, tau) + ufl.div(tau) * u + ufl.div(sigma) * v) * ufl.dx
L = -ufl.inner(f, v) * ufl.dx + ufl.inner(g, v) * ufl.ds

# Solve
problem = LinearProblem(a, L, [bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu",
                                                    "pc_factor_mat_solver_type": "mumps"})
uh = problem.solve()
print(uh.x.array, problem.solver.getConvergedReason())

This code is written from dolfinx 0.7.0. The only difference is the usage of

el0 = basix.ufl.element("BDM", domain.topology.cell_name(),1)
el1 = basix.ufl.element("DG", domain.topology.cell_name(), 0)
W = fem.functionspace(mesh=domain, element=basix.ufl.mixed_element([el0, el1]))

instead of

I’ll aim to add a singular Poisson demo to the FEniCSx tutorial over the weekend.