Hi everyone,
I am implementing a FEM framework in C++ and am using dolfinx to have a prototype that I can check against.
I am currently trying to implement periodic boundary conditions for the Poisson equation, with order 1 Lagrange elements in 1D and a 5th order quadrature (satisfying > 2p+1).
To study the convergence, I am setting the right hand side \pi^2 sin(\pi x) and solving in the domain [-1, 1]. This gives correct 2nd order convergence.
However, when I change the right hand side to \pi^2 cos(\pi x), it converges with a lower rate (order 1).
Does this mean there is something wrong in my implementation of the boundary conditions? I have printed out and checked a lot of values and it seems correct, and it does converge, just with a lower order. Or am I missing some mathematical nuance that prescribes this lower order of convergence for this problem?
I have also tried shifting the domain to [0,2] instead of [-1,1] with the right hand side as \pi^2 sin(\pi x) and I also get a lower order of convergence.
Here is the code I am using for the sin case (and I just change the exact solution and the right hand side to cos when checking that):
from dolfinx import fem, mesh, default_scalar_type, io
from mpi4py import MPI
import ufl
import numpy as np
from dolfinx_mpc import LinearProblem, MultiPointConstraint
from petsc4py import PETSc
comm = MPI.COMM_WORLD
rank = comm.Get_rank()
tol = 250 * np.finfo(default_scalar_type).resolution
def run(num_els=10,
space_order=1,
interpolate_rhs=False,
quadrature_degree=5,
write_vtk=False):
exact_sol = lambda x : np.sin(np.pi * (x[0]))
domain = mesh.create_interval(comm, num_els, [-1.0, 1.0])
dim = domain.topology.dim #dim=1
fdim = dim - 1# facet dimension
V = fem.functionspace(domain, ("Lagrange", space_order),)
u_exact = fem.Function(V, name="u_exact")
u_exact.interpolate( exact_sol )
def periodic_boundary(x):
return np.isclose(x[0], 1, atol=tol)
def periodic_relation(x):
out_x = np.zeros_like(x)
out_x[0] = -x[0]
out_x[1] = +x[1]
out_x[2] = +x[2]
return out_x
bcs = []
mpc = MultiPointConstraint(V)
mpc.create_periodic_constraint_geometrical(V, periodic_boundary, periodic_relation, bcs=None)
mpc.finalize()
(u, v) = (ufl.TrialFunction(V), ufl.TestFunction(V))
dx = ufl.dx(metadata={"quadrature_degree":quadrature_degree})
a_ufl = ufl.dot( ufl.grad(u), ufl.grad(v) ) * dx
x = ufl.SpatialCoordinate(domain)
if interpolate_rhs:
f = fem.Function(V, name="f_interpolated")
f.interpolate(lambda x : (np.pi**2)*exact_sol(x))
else:
f = ufl.pi**2 * ufl.sin(ufl.pi * (x[0]))
l_ufl = v * f * dx
# Solve
problem = LinearProblem(a_ufl, l_ufl, mpc, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
solver = problem.solver
solver_prefix = "dolfinx_mpc_solve_{}".format(id(solver))
solver.setOptionsPrefix(solver_prefix)
petsc_options = {
"ksp_type": "cg",
"ksp_rtol": 1e-6,
"pc_type": "hypre",
"pc_hypre_type": "boomeramg",
"pc_hypre_boomeramg_max_iter": 1,
"pc_hypre_boomeramg_cycle_type": "v", # ,
# "pc_hypre_boomeramg_print_statistics": 1
}
# Set PETSc options
opts = PETSc.Options() # type: ignore
opts.prefixPush(solver_prefix)
if petsc_options is not None:
for k, v in petsc_options.items():
opts[k] = v
opts.prefixPop()
solver.setFromOptions()
nullspace = PETSc.NullSpace().create(constant=True,comm=MPI.COMM_WORLD)
problem.A.setNullSpace(nullspace)
uh = problem.solve()
uh.name = "uh"
# Compute error
uex = ufl.sin(ufl.pi * (x[0])) # exact sol
L2_error = fem.form(ufl.inner(uh - uex, uh - uex) * dx)
error_local = fem.assemble_scalar(L2_error) # local to rank
error_L2 = np.sqrt(domain.comm.allreduce(error_local, op=MPI.SUM))
error_max = np.max(np.abs(u_exact.x.array-uh.x.array))
# Only print the error on one process
if domain.comm.rank == 0:
print(f"Ran with {num_els} elements, interpolate_rhs = {interpolate_rhs}, quadrature_degree = {quadrature_degree}, order {space_order} elements.")
print(f"L2 error : {error_L2:.2e}")
print(f"Max nodal error : {error_max:.2e}")
print("===========================")
if write_vtk:
with io.VTKFile(domain.comm, f"post/solution{num_els}_els.pvd", "w") as pf:
pf.write_function([uh,u_exact])
return error_L2, error_max
def get_convergence_rate():
space_order = 1
interpolate_rhs = True
l2_errors = []
quadrature_degree = 5
total_elements = np.array([4, 8, 16, 32, 64])
l2_errors = []
for num_els in total_elements:
l2_error, _ = run(num_els=num_els,
space_order=space_order,
interpolate_rhs=interpolate_rhs,
quadrature_degree=quadrature_degree,
write_vtk=True)
l2_errors.append( l2_error )
domain_size = 2.0
element_sizes = domain_size / total_elements
line = np.polyfit(np.log(element_sizes), np.log(l2_errors), 1) # fit degree 1 polynomial to the logs
convergence_rate = line[0]
if rank==0:
print(f"Ran with interpolate_rhs = {interpolate_rhs}, quadrature_degree = {quadrature_degree}, order {space_order} elements.")
print(f"Convergence rate = {convergence_rate}")
if __name__=="__main__":
get_convergence_rate()
Thanks in advance for any help!