# How to define a time-dependent Neumann boundary condition if I have known the exact solution

Hello,

I want to solve a heat equation like this:

\frac{\partial u}{\partial t} - \nabla \cdot (\kappa \nabla u) = f

with the Neumann boundary condition:

-\frac{\partial u}{\partial \boldsymbol{n}} = g

I want to compute the error between exact solution and FEM solution.

u_ex = lambda x: 1 + x[0]**2 + 2*x[1]**2
x = SpatialCoordinate(mesh)
n = FacetNormal(mesh)


Here is my question:
My exact solution is time-dependent, e.g. u_{ex} = te^{-x-y}. How to define g as the tutorial.

My code is

from mpi4py import MPI
import numpy as np
from petsc4py import PETSc
from dolfinx.mesh import create_unit_square
from dolfinx.fem import functionspace
from ufl import dot, FacetNormal, grad, SpatialCoordinate

t = 0.0

domain = create_unit_square(MPI.COMM_WORLD, 10, 10)
V = functionspace(domain, ("Lagrange", 1))

class u_ExactSolution():
def __init__(self, t):
self.t = t

def __call__(self, x):
value = np.zeros((2, x.shape[1]), PETSc.ScalrType)
value[0] = self.t * np.exp(-x[0] - x[1])
return value

u_ex = u_ExactSolution(t)



Of course, it fails since ‘SpatialCoordinate’ object has no attribute ‘shape’

If I were you, I’d rather do the math on paper (or with some symbolic algebra system) and compute myself -\partial u_{ex} / \partial \mathbf{n}. It will be a math expression depending on x and y, and once you have it, you can type it up using ufl.SpatialCoordinate for (x, y).

This may be known as “method of manufactured solutions” in some textbooks.

You can automate the solution with something along the lines of

t_c = dolfinx.fem.Constant(domain, 0.0)
x = ufl.SpatialCoordinate(mesh)
u_ex = t_c * ufl.exp(-x-y)


then if you want to update t_c, you can do that with
t_c.value += float(dt)

Thanks, @francesco-ballarin. I have already computed these terms on paper, and I think I should define different markers for each boundary.

Thank you very much!
Here is the MWE:

from mpi4py import MPI
from dolfinx.mesh import create_unit_square
from dolfinx.fem import functionspace
domain = create_unit_square(MPI.COMM_WORLD, 10, 10)
V = functionspace(domain, ("Lagrange", 1))

from dolfinx.fem import Constant
from ufl import SpatialCoordinate, exp
t = 0.0
t_c = Constant(domain, t)
x = SpatialCoordinate(domain)
u_exact = t_c * exp(-x[0] - x[1])

from ufl import dot, FacetNormal, grad
n = FacetNormal(domain)


Moreover, I tried another method whose MWE is

from mpi4py import MPI
import numpy as np
from petsc4py import PETSc
from dolfinx.mesh import create_unit_square
from dolfinx.fem import functionspace, Function
from ufl import dot, FacetNormal, grad

t = 0.0

domain = create_unit_square(MPI.COMM_WORLD, 10, 10)
V = functionspace(domain, ("Lagrange", 1))

class u_ExactSolution():
def __init__(self, t):
self.t = t

def __call__(self, x):
return self.t * np.exp(-x[0] - x[1])

u_exact = Function(V)
u_ex = u_ExactSolution(t)
u_exact.interpolate(u_ex)


I remember that the tutorial says two methods, one is SpatialCoordinate and one is interpolate.

Agreed, but that is shown in the tutorial you linked, search for meshtags. If I were, I would still try doing it (also) this way, because that will be the way you’ll set Neumann boundary conditions in any practical case (i.e., a case in which you do not have the exact solution, and can ask ufl to do grad(u) times n).