Difference in performance between FEniCS and FEniCSX for higher order

General
1

Hi community, I am comparing some libraries, and FEniCSX seems to have issues for higher order. The test is solving

-\Delta u + u = 1, \quad \in (0,1)^3

with homogeneous Neumann BCs. The number of iterations using CG+Jacobi for P1 and P2 are equal between FEniCS and FEniCSX, but for P3 they do 227 and 247 resp. To confim, I have checked with FreeFEM and Firedrake, and they all do 227, so it is probably a problem with FEniCSX. I have checked, and varying the quadrature degree from 4 to 12 gives the same result (less than that results in divergence).

If you have another idea please let me know, if not I will report a bug.

Best!

1 Like
2

Without the code corresponding to the results you are seeing, it is hard to give you any pointers.

1 Like
3

Sorry, dumb mistake. Here they are:

FEniCS:

from petsc4py.PETSc import Options
from dolfin import *
N = 20
P = 3
mesh = UnitCubeMesh(N, N, N)
V = FunctionSpace(mesh, "Lagrange", P)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v))*dx + u*v*dx
F = v*dx
u_sol = Function(V)
A = assemble(a)
FF = assemble(F)
solver_ops = Options()
solver_ops.setValue('ksp_type', 'cg')
solver_ops.setValue('pc_type', 'jacobi')
solver_ops.setValue('ksp_monitor', None)
solver_ops.setValue('ksp_atol', 1e-14)
solver_ops.setValue('ksp_rtol', 1e-6)
solver = PETScKrylovSolver()
solver.set_from_options()
solver.solve(A, u_sol.vector(), FF)

FEniCSX

from dolfinx import fem, io, mesh, plot
import ufl
from ufl import ds, dx, grad, inner
from mpi4py import MPI
N = 20
P = 3
msh = mesh.create_unit_cube(comm=MPI.COMM_WORLD, nx=N, ny=N, nz=N)
V = fem.FunctionSpace(msh, ("Lagrange", P))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = inner(grad(u), grad(v)) + u*v
a = a * dx
L = v * dx
problem = fem.petsc.LinearProblem(a, L, petsc_options={"ksp_type": "cg", "pc_type": "jacobi", "ksp_monitor": None, "ksp_atol": 1e-14, "ksp_rtol": 1e-6})
uh = problem.solve()
4

This is because DOLFINx uses GLL dof placement of higher order Lagrange dofs, see:
https://docs.fenicsproject.org/dolfinx/v0.5.1/python/demos/demo_lagrange_variants.html#equispaced-points-vs-gll-points

Thus, you would get different matrices, and different results, as your underlying polynomials and definitions of degrees of freedom would differ.

MWE. to produce the same iteration count as with DOLFIN

from dolfinx import fem, io, mesh, plot
import ufl
from ufl import ds, dx, grad, inner
from mpi4py import MPI
from basix import ufl_wrapper
import basix

N = 20
P = 3
msh = mesh.create_unit_cube(comm=MPI.COMM_WORLD, nx=N, ny=N, nz=N)
element = ufl_wrapper.create_element(
    "Lagrange", msh.topology.cell_type.name, P,
    lagrange_variant=basix.LagrangeVariant.equispaced)
V = fem.FunctionSpace(msh, element)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = inner(grad(u), grad(v)) + u*v
a = a * dx
L = v * dx
options = {"ksp_type": "cg",
           "pc_type": "jacobi",
           "ksp_monitor": None,
           "ksp_atol": 1e-14,
           "ksp_rtol": 1e-6}
problem = fem.petsc.LinearProblem(a, L, petsc_options=options)
uh = problem.solve()
1 Like
5

Amazing, thank you very much!