Manufactured solution for poisson type problem

dolfinx
1

Hello FEniCS experts!

I have this script as an MMS test for a Poisson-type problem \nabla(D \cdot \nabla u) + \Gamma = 0 where \Gamma is a source term for an exact solution u = 1 + \sin(2 \pi x). It has been adapted from the Error control tutorial on the FEniCSx tutorial and a Previous post.

from mpi4py import MPI
from dolfinx import fem
import ufl
import numpy as np
from dolfinx.mesh import create_unit_square
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver


def u_ex(mod):
    return lambda x: 1 + mod.sin(2 * mod.pi * x[0])


u_ufl = u_ex(ufl)
u_numpy = u_ex(np)

# mesh
my_mesh = create_unit_square(MPI.COMM_WORLD, 50, 50)
x = ufl.SpatialCoordinate(my_mesh)

elements = ufl.FiniteElement("P", my_mesh.ufl_cell(), 1)
V = fem.FunctionSpace(my_mesh, elements)
u = fem.Function(V)
v = ufl.TestFunction(V)

# define bcs
u_bc = fem.Function(V)
u_bc.interpolate(u_numpy)
dofs_L = fem.locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 0))
dofs_R = fem.locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 1))
bc_left = fem.dirichletbc(u_bc, dofs_L)
bc_right = fem.dirichletbc(u_bc, dofs_R)
bcs = [bc_left, bc_right]

D = 1.0

# source term
f = ufl.grad(D * ufl.grad(u_ufl(x)))[0, 0]

F = 0
F += ufl.dot(D * ufl.grad(u), ufl.grad(v)) * ufl.dx
F += f * v * ufl.dx

problem = NonlinearProblem(F, u, bcs=bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.solve(u)

I was just wondering if someone would be able to explain the significance of the index [0, 0] on the source term f. The script only seems to work with it there.

Plus I was wondering if there was any way to have this source term more generic i.e is it possible to interpolate it into a function for instance?

Thanks in advance,
James

2

This should be ufl.div(D*….) not gradient if it is a Poisson problem.

You can interpolate it into a function space, sure, But then you pay the price of approximating the source term in a space it does not lie in (as you are basing your analytical solution on sinus). Is there a specific motivation for wanting to do this?

3

Yeah you’re right, should’ve been ufl.div(D*...)! Thanks

And regarding the interpolation, it’s just for ease of handling the value of the source term within our application, just easier to allow passing of a value of type fem.Function rather than allowing a bunch of different ufl objects (ufl.differentiation.Div or ufl.algebra.Product for instance). But again you’re right it does seem interpolating it onto a function space does seem to introduce approximations, increasing the L2 error.

In some testing, I interpolated the value using an fem.Expression:

f_value = -ufl.div(D * ufl.grad(u_ufl(x)))
f_expr = fem.Expression(f_value, V.element.interpolation_points())
f = fem.Function(V)
f.interpolate(f_expr)

And in doing this it increased the L2 error compared to the analytical solution from 2.5e-08 to 4.7e-08 in this case. Not sure how much this may change with other examples or how acceptable this is…

4

As long as you are aware of any errors introduced by the interpolation, I guess you are fine.
For analytical problems you can compute the error by computing the L2 error between the interpolated value ad the value evaluated at quadrature points.

If you are concerned about type hints, the general type for such “expressions” is ufl.core.expr.Expr which can take in any ufl expression.

5

Ok noted, cheers @dokken !