Hi, I am recently looking at the mixed FEM in FEniCSx. I tried running the code in: Mixed formulation for the Poisson equation — DOLFINx 0.10.0.0 documentation, and compare the results to the code in: Poisson equation — DOLFINx 0.10.0.0 documentation. Here I changed the domain of the ‘demo_poisson’ to msh = mesh.create_unit_square(MPI.COMM_WORLD, 32, 32, mesh.CellType.quadrilateral) and facets to facets = mesh.locate_entities_boundary( msh, dim=(msh.topology.dim - 1), marker=lambda x: np.isclose(x[0], 0.0) | np.isclose(x[0], 1.0), ). Just to make sure they study the same problem. However, the mixed poisson code gives me completely different results, as you can see from the figures (left: mixed poisson, right: poisson). Does anyone have the same issue?
Are you sure you’re handling the boundary conditions appropriately?
The standard FE discretisation of the Poisson problem enforces the solution’s Dirichlet data in a strong sense, and the flux data weakly.
The H(\mathrm{Div}) conforming FE discretisation of the mixed Poisson problem enforces the solution’s Dirichlet data weakly, and the flux data in a strong sense.
Hi nate, thanks for the prompt response. I think so, because I didn’t change any other parts of the code.
Please provide the code you are running here with proper formatting, to ensure that we can reproduce it.
Here is the code I am using:
from mpi4py import MPI
import numpy as np
from basix.ufl import element, mixed_element
from dolfinx import default_real_type, fem, io, mesh
from dolfinx.fem.petsc import LinearProblem
from ufl import Measure, SpatialCoordinate, TestFunctions, TrialFunctions, div, exp, inner
msh = mesh.create_unit_square(MPI.COMM_WORLD, 32, 32, mesh.CellType.quadrilateral)
k = 1
Q_el = element("BDMCF", msh.basix_cell(), k, dtype=default_real_type)
P_el = element("DG", msh.basix_cell(), k - 1, dtype=default_real_type)
V_el = mixed_element([Q_el, P_el])
V = fem.functionspace(msh, V_el)
(sigma, u) = TrialFunctions(V)
(tau, v) = TestFunctions(V)
x = SpatialCoordinate(msh)
f = 10.0 * exp(-((x[0] - 0.5) * (x[0] - 0.5) + (x[1] - 0.5) * (x[1] - 0.5)) / 0.02)
dx = Measure("dx", msh)
a = inner(sigma, tau) * dx + inner(u, div(tau)) * dx + inner(div(sigma), v) * dx
L = -inner(f, v) * dx
# Get subspace of V
V0 = V.sub(0)
fdim = msh.topology.dim - 1
facets_top = mesh.locate_entities_boundary(msh, fdim, lambda x: np.isclose(x[1], 1.0))
Q, _ = V0.collapse()
dofs_top = fem.locate_dofs_topological((V0, Q), fdim, facets_top)
def f1(x):
values = np.zeros((2, x.shape[1]))
values[1, :] = np.sin(5 * x[0])
return values
f_h1 = fem.Function(Q)
f_h1.interpolate(f1)
bc_top = fem.dirichletbc(f_h1, dofs_top, V0)
facets_bottom = mesh.locate_entities_boundary(msh, fdim, lambda x: np.isclose(x[1], 0.0))
dofs_bottom = fem.locate_dofs_topological((V0, Q), fdim, facets_bottom)
def f2(x):
values = np.zeros((2, x.shape[1]))
values[1, :] = -np.sin(5 * x[0])
return values
f_h2 = fem.Function(Q)
f_h2.interpolate(f2)
bc_bottom = fem.dirichletbc(f_h2, dofs_bottom, V0)
bcs = [bc_top, bc_bottom]
problem = LinearProblem(
a,
L,
bcs=bcs,
petsc_options={
"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type": "superlu_dist",
},
)
try:
w_h = problem.solve()
except PETSc.Error as e: # type: ignore
if e.ierr == 92:
print("The required PETSc solver/preconditioner is not available. Exiting.")
print(e)
exit(0)
else:
raise e
sigma_h, u_h = w_h.split()
with io.XDMFFile(msh.comm, "out_mixed_poisson/u.xdmf", "w") as file:
file.write_mesh(msh)
file.write_function(u_h)
from mpi4py import MPI
import numpy as np
from dolfinx import default_scalar_type
import ufl
from dolfinx import fem, io, mesh, plot
from dolfinx.fem.petsc import LinearProblem
from ufl import ds, dx, grad, inner
msh = mesh.create_unit_square(MPI.COMM_WORLD, 32, 32, mesh.CellType.quadrilateral)
V = fem.functionspace(msh, ("Lagrange", 1))
facets = mesh.locate_entities_boundary(
msh,
dim=(msh.topology.dim - 1),
marker=lambda x: np.isclose(x[0], 0.0) | np.isclose(x[0], 1.0),
)
dofs = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facets)
bc = fem.dirichletbc(value=default_scalar_type(0), dofs=dofs, V=V)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
x = ufl.SpatialCoordinate(msh)
f = 10 * ufl.exp(-((x[0] - 0.5) ** 2 + (x[1] - 0.5) ** 2) / 0.02)
g = ufl.sin(5 * x[0])
a = inner(grad(u), grad(v)) * dx
L = inner(f, v) * dx + inner(g, v) * ds
problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
with io.XDMFFile(msh.comm, "out_poisson/poisson.xdmf", "w") as file:
file.write_mesh(msh)
file.write_function(uh)
Thanks! I checked my code and it seems like the problem with the solver. I reproduced the same results as you after changing the solver into mumps:
problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu",
"pc_factor_mat_solver_type": "mumps"})
However, I don’t know why “superlu_dist” is not working on my side. It might be because of the dolfinx version I am using is 0.8.0.
Likely it’s not detecting the zero block of the mixed Poisson saddle point system. I think MUMPS automatically checks for this.