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’

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.

dokken
June 6, 2024, 11:19am
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)`

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

.

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`

).