Hi,
I’m having some problems solving an elliptic PDE equation using the finite element method, and I’m finding that the error I’m getting is much larger compared to the expected error value. Below is the implementation of my code, including mesh generation, boundary condition setting and solving steps.
from petsc4py import PETSc
from mpi4py import MPI
import ufl
from dolfinx import mesh, fem
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_matrix, create_vector
import numpy as np
import time
start = time.time()
nx, ny = 10, 10
domain = mesh.create_unit_square(MPI.COMM_WORLD, nx, ny, mesh.CellType.triangle)
V = fem.FunctionSpace(domain, ("Lagrange", 1))
class exact_solution():
def __call__(self, x):
return 1 + x[0]**2 + x[1]**2
u_exact_expr = exact_solution()
def source_term(x):
return 4.0 + (1 + x[0]**2 + x[1]**2)
f = fem.Function(V)
f.interpolate(source_term)
u_D = fem.Function(V)
u_D.interpolate(u_exact_expr)
n = ufl.FacetNormal(domain)
g = ufl.dot(n, ufl.grad(u_D))
def left_boundary(x):
return np.isclose(x[0], 0.0)
bc = fem.dirichletbc(u_D, fem.locate_dofs_geometrical(V, left_boundary))
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
F = u * v * ufl.dx + ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx - f * v * ufl.dx + g * v * ufl.ds
a = ufl.lhs(F)
L = ufl.rhs(F)
A = fem.petsc.assemble_matrix(fem.form(a), bcs=[bc])
A.assemble()
b = fem.petsc.create_vector(fem.form(L))
uh = fem.Function(V)
solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)
fem.petsc.assemble_vector(b, fem.form(L))
apply_lifting(b, [fem.form(a)], bcs=[[bc]]) # bc needs to be wrapped in a list of lists
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
fem.set_bc(b, [bc])
solver.solve(b, uh.vector)
uh.x.scatter_forward()
u_exact_fem = fem.Function(V)
u_exact_fem.interpolate(u_exact_expr)
error_L2 = np.sqrt(domain.comm.allreduce(fem.assemble_scalar(fem.form((uh - u_exact_fem)**2 * ufl.dx)), op=MPI.SUM))
print(f"L2-error: {error_L2:.2e}")
end = time.time()
print(f"Execution time: {(end - start) * 10 ** 3:.6f} ms")
Can you please tell me what modifications I should make to reduce the error value of the result. I would be appreciated if you could help point out my problem and suggest specific modifications.
change your FE space to P2. That space is able to represent exactly the solution, so you must get an error which is zero (up to numerical precision). If you don’t get zero, then there is something wrong in your formulation and/or your code. If you get zero, continue to step 2.
go back to P1 as FE space. Carry out a convergence analysis by decreasing the mesh size, i.e. try nx=ny=10 and store in a spreadsheet the error you get, then try nx=ny=20 and store again the error you get, and so on. Compare the convergence rate to the one you expect from the FE error analysis.
As a side note:
this seems like an exercise for an introductory FE class. If that is the case, I’d rather not use the definition of g you have in your code. Just take some time to compute the gradient of the exact solution yourself, and the normal of each of the Neumann boundaries.
writing g * v * ufl.ds will integrate over the whole boundary, including left_boundary. Integrating, or not, over that part of the boundary won’t make any difference (because when you write the formulation on paper you assume v = 0 on the left boundary), but at least be aware of this.
there is no need to interpolate f in the FE space, just use a ufl expression for that. The way you are doing now makes comparing the rates you get to the expected ones more difficult, because your FE forcing term is a P1 approximation of the exact forcing rate (which is a P2 function). In layman term, you are solving a slightly different problem from the one you thought you were solving. Similarly for u_exact_fem when computing the error.
I changed my code as you suggested but the error result changed from 7.16e-01 to 9.93e-01.
Is this because my code or the equation itself is wrong? But I can’t see what the problem is, can you help me to modify my code? I will be very grateful!
Okay. Firstly, I made a mistake in my original f calculation and modified f after recalculating.
Secondly, I adjusted the plus and minus signs of the Newman boundary conditions and changed my FE space to P2.
Here’s the revised code:
from petsc4py import PETSc
from mpi4py import MPI
import ufl
from dolfinx import mesh, fem
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_matrix, create_vector
import numpy as np
import time
start = time.time()
nx, ny = 10, 10
domain = mesh.create_unit_square(MPI.COMM_WORLD, nx, ny, mesh.CellType.triangle)
V = fem.FunctionSpace(domain, ("Lagrange", 2))
def exact_solution(x):
return 1 + x[0]**2 + x[1]**2
u_exact_expr = exact_solution
def source_term(x):
return (1 + x[0]**2 + x[1]**2) - 4.0
f = fem.Function(V)
f.interpolate(source_term)
u_D = fem.Function(V)
u_D.interpolate(u_exact_expr)
n = ufl.FacetNormal(domain)
g = ufl.dot(n, ufl.grad(u_D))
def left_boundary(x):
return np.isclose(x[0], 0.0)
bc = fem.dirichletbc(u_D, fem.locate_dofs_geometrical(V, left_boundary))
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
F = u * v * ufl.dx + ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx - f * v * ufl.dx - g * v * ufl.ds
a = ufl.lhs(F)
L = ufl.rhs(F)
A = fem.petsc.assemble_matrix(fem.form(a), bcs=[bc])
A.assemble()
b = fem.petsc.create_vector(fem.form(L))
uh = fem.Function(V)
solver = PETSc.KSP().create(domain.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)
fem.petsc.assemble_vector(b, fem.form(L))
apply_lifting(b, [fem.form(a)], bcs=[[bc]])
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
fem.set_bc(b, [bc])
solver.solve(b, uh.vector)
uh.x.scatter_forward()
u_exact_fem = fem.Function(V)
u_exact_fem.interpolate(u_exact_expr)
error_L2 = np.sqrt(domain.comm.allreduce(fem.assemble_scalar(fem.form((uh - u_exact_fem)**2 * ufl.dx)), op=MPI.SUM))
print(f"L2-error: {error_L2:.2e}")
end = time.time()
print(f"Execution time: {(end - start) * 10 ** 3:.6f} ms")
The error results obtained from this code run are significantly smaller: L2-error: 8.38e-14