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

variational formulation
1

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.
I noticed https://jsdokken.com/dolfinx-tutorial/chapter3/robin_neumann_dirichlet.html which reads

u_ex = lambda x: 1 + x[0]**2 + 2*x[1]**2
x = SpatialCoordinate(mesh)
n = FacetNormal(mesh)
g = -dot(n, grad(u_ex(x)))

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)

g = -dot(FacetNormal(domain), grad(u_ex(SpatialCoordinate(domain))))

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

2

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.

3

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)

4

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

5

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)
g = -dot(n, grad(u_exact))

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.

6

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).