# How to solve a heat equation whose analysis solution is known and expressed by SpatialCoordinate

My problem is:
find u = u(\boldsymbol{x}, t) obeying

\begin{alignat}{2} \frac{\partial u}{\partial t} - \kappa \Delta u &= f, && (\boldsymbol{x}, t) \in \Omega \times [0,T],\\ u &= u_D, &\quad& (\boldsymbol{x}, t) \in \partial \Omega \times [0,T],\\ u(0) &= u_0, && \boldsymbol{x} \in \Omega, \end{alignat}

where \Omega is unit square, T =2, and the diffusion coefficient

$$\kappa = \kappa (t).$$

The sources

\begin{align} f &= 2t - \kappa (6x + 12y),\\ u_D &= 1 + x^3 + 2y^3 + t^2,\\ u_0 &= 1 + x^3 + 2y^3 \end{align}

are derived from exact solution

$$u_{ex} = 1 + x^3 + 2y^3 + t^2.$$

The diffusion coefficient

$$\kappa = \sqrt{t}$$

in this demo.

I tried this 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, set_bc
import numpy as np

t = 0.0
T = 2.0
num_steps = 20
dt = (T - t) / num_steps

nx, ny = 5, 5
domain = mesh.create_unit_square(MPI.COMM_WORLD, nx, ny, mesh.CellType.triangle)
V = fem.functionspace(domain, ("Lagrange", 1))

x = ufl.SpatialCoordinate(domain)
t_c = fem.Constant(domain, t)
u_exact = 1 + x[0]**3 + 2 * x[1]**3 + t_c ** 2

u_D = fem.Function(V)
u_D.interpolate(fem.Expression(u_exact, V.element.interpolation_points()))

tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets))

u_n = fem.Function(V)
u_n = u_exact

kappa = ufl.sqrt(t_c)

f = fem.Function(V)
f = 2 * t_c - kappa * (6 * x[0] + 12 * x[1])

u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
F = u * v * ufl.dx + dt * kappa * ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx - u_n * v * ufl.dx - dt * f * v * ufl.dx
a = fem.form(ufl.lhs(F))
L = fem.form(ufl.rhs(F))

A = create_matrix(a)
b = create_vector(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)

errors_L2 = []
for n in range(num_steps):
t += dt
t_c.value = t
u_D.interpolate(fem.Expression(u_exact, V.element.interpolation_points()))

A.zeroEntries()
assemble_matrix(A, a, bcs=[bc])
A.assemble()

with b.localForm() as loc_b:
loc_b.set(0)
assemble_vector(b, L)

apply_lifting(b, [a], [[bc]])
set_bc(b, [bc])

solver.solve(b, uh.vector)
uh.x.scatter_forward()

u_n = uh

V_ex = fem.functionspace(domain, ("Lagrange", 2))
u_ex = fem.Function(V_ex)
u_ex = u_exact
error_L2 = domain.comm.allreduce(fem.assemble_scalar(fem.form((uh - u_ex)**2 * ufl.dx)), op=MPI.SUM)
errors_L2.append(error_L2)

L2_error = np.sqrt(np.trapz(errors_L2, dx=dt))
print(f"L2-error: {L2_error:.2e}")


However, the error is about e^{-1}. Then I tried another 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, set_bc
import numpy as np

t = 0
T = 2
num_steps = 20
dt = float((T - t) / num_steps)

nx, ny = 5, 5
domain = mesh.create_unit_square(MPI.COMM_WORLD, nx, ny, mesh.CellType.triangle)
V = fem.functionspace(domain, ("Lagrange", 1))

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

def __call__(self, x):
return 1 + x[0] ** 3 + 2 * x[1] ** 3 + self.t ** 2

u_exact = exact_solution(t)

u_D = fem.Function(V)
u_D.interpolate(u_exact)
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets))

u_n = fem.Function(V)
u_n.interpolate(u_exact)

kappa_expression = np.sqrt(t)
kappa = fem.Constant(domain, kappa_expression)

class source_term():
def __init__(self, t, kappa):
self.t = t
self.kappa = kappa

def __call__(self, x):
return 2 * self.t - self.kappa * (6 * x[0] + 12 * x[1])

f_expression = source_term(t, kappa_expression)

f = fem.Function(V)
f.interpolate(f_expression)

u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
F = u * v * ufl.dx + dt * kappa * ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx - u_n * v * ufl.dx - dt * f * v * ufl.dx
a = fem.form(ufl.lhs(F))
L = fem.form(ufl.rhs(F))

A = create_matrix(a)
b = create_vector(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)

errors_L2 = []
for n in range(num_steps):

t += dt
kappa_expression = np.sqrt(t)
kappa.value = kappa_expression
u_exact.t = t
f_expression.t = t
f_expression.kappa = kappa_expression
u_D.interpolate(u_exact)
f.interpolate(f_expression)

A.zeroEntries()
assemble_matrix(A, a, bcs = [bc])
A.assemble()

with b.localForm() as loc_b:
loc_b.set(0)
assemble_vector(b, L)

apply_lifting(b, [a], [[bc]])
set_bc(b, [bc])

solver.solve(b, uh.vector)
uh.x.scatter_forward()

u_n.x.array[:] = uh.x.array

V_ex = fem.functionspace(domain, ("Lagrange", 2))
u_ex = fem.Function(V_ex)
u_ex.interpolate(u_exact)
error_L2 = domain.comm.allreduce(fem.assemble_scalar(fem.form((uh - u_ex)**2 * ufl.dx)), op=MPI.SUM)
errors_L2.append(error_L2)

L2_error = np.sqrt(np.trapz(errors_L2, dx = dt))
print(f"L2-error: {L2_error:.2e}")


and the error is about 5e^{-2}.
I think I misunderstood the usage of SpatialCoordinate. Could anyone tell me how to correct the first code ?
Thanks very much!

The issue is that you are reassigning variables, which are not reflected in updating your variational form.

Consider the following MWE that leverages DOLFINx form compilation:

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, set_bc
import numpy as np

t = 0.0
T = 2.0
num_steps = 10
dt = (T - t) / num_steps

nx, ny = 12, 12
domain = mesh.create_unit_square(MPI.COMM_WORLD, nx, ny, mesh.CellType.triangle)
V = fem.functionspace(domain, ("Lagrange", 1))

x = ufl.SpatialCoordinate(domain)
t_c = fem.Constant(domain, t)
u_exact = 1 + x[0]**3 + 2 * x[1]**3 + t_c ** 2
exact_expr = fem.Expression(u_exact, V.element.interpolation_points())
u_D = fem.Function(V)
u_D.interpolate(exact_expr)

tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets))

u_n = fem.Function(V)
u_n.interpolate(exact_expr)

kappa = ufl.sqrt(t_c)

f = fem.Function(V)
f = 2 * t_c - kappa * (6 * x[0] + 12 * x[1])

u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
F = u * v * ufl.dx + dt * kappa * ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx - u_n * v * ufl.dx - dt * f * v * ufl.dx
a = fem.form(ufl.lhs(F))
L = fem.form(ufl.rhs(F))

A = create_matrix(a)
b = create_vector(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)
error_form = fem.form((uh - u_exact)**2 * ufl.dx)
errors_L2 = []
for n in range(num_steps):
t += dt
t_c.value = t
u_D.interpolate(fem.Expression(u_exact, V.element.interpolation_points()))

A.zeroEntries()
assemble_matrix(A, a, bcs=[bc])
A.assemble()

with b.localForm() as loc_b:
loc_b.set(0)
assemble_vector(b, L)

apply_lifting(b, [a], [[bc]])
set_bc(b, [bc])

solver.solve(b, uh.vector)
uh.x.scatter_forward()

u_n.x.array[:] = uh.x.array

error_L2 = domain.comm.allreduce(fem.assemble_scalar(error_form), op=MPI.SUM)
errors_L2.append(error_L2)

L2_error = np.sqrt(np.trapz(errors_L2, dx=dt))
print(f"L2-error: {L2_error:.2e}")


yielding

L2-error: 1.85e-02


Thanks.Moreover, I just wonder that I can defint the error_form like

error_form = fem.form((uh - u_D) ** 2 * ufl.dx)


I tried this form and got a smaller error.
Why the error is smaller?

If you use spatial coordinate directly the error is computed between the approximate solution (uh) and the exact solution at the quadrature points. This means that your represent u_exact as accuratly as possible.

If you use u_D, u_exact is first interpolated into the approximation space and you are looking at the difference between the interpolated approximation of u_exact and uh