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

**URL:** <https://fenicsproject.discourse.group/t/how-to-solve-a-heat-equation-whose-analysis-solution-is-known-and-expressed-by-spatialcoordinate/15222>\
**Category:** General\
**Created:** [July 11, 2024, 12:37am UTC](https://fenicsproject.discourse.group/t/how-to-solve-a-heat-equation-whose-analysis-solution-is-known-and-expressed-by-spatialcoordinate/15222 "2024-07-11T00:37:33Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![OsbertWang](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/osbertwang/32/2138_2.png) [@OsbertWang](https://fenicsproject.discourse.group/u/OsbertWang)\
**Post date:** [July 11, 2024, 12:37am UTC](https://fenicsproject.discourse.group/t/how-to-solve-a-heat-equation-whose-analysis-solution-is-known-and-expressed-by-spatialcoordinate/15222/1 "2024-07-11T00:37:33Z")

</div>

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

\begin{equation} \kappa = \kappa (t). \end{equation}

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

\begin{equation} u\_{ex} = 1 + x^3 + 2y^3 + t^2. \end{equation}

The diffusion coefficient

\begin{equation} \kappa = \sqrt{t} \end{equation}

in this demo.

I tried this code:

```auto
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]])
    b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    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:

```auto
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]])
    b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    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!

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [July 11, 2024, 6:41am UTC](https://fenicsproject.discourse.group/t/how-to-solve-a-heat-equation-whose-analysis-solution-is-known-and-expressed-by-spatialcoordinate/15222/2 "2024-07-11T06:41:33Z")

</div>

> [@OsbertWang](#):
>
> ```auto
> 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)
> 
> ```

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:

```python
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]])
    b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    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

```bash
L2-error: 1.85e-02

```

---

<div class="post-metadata">

**Author:** ![OsbertWang](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/osbertwang/32/2138_2.png) [@OsbertWang](https://fenicsproject.discourse.group/u/OsbertWang)\
**Post date:** [July 12, 2024, 10:43am UTC](https://fenicsproject.discourse.group/t/how-to-solve-a-heat-equation-whose-analysis-solution-is-known-and-expressed-by-spatialcoordinate/15222/3 "2024-07-12T10:43:23Z")

</div>

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

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

```

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

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [July 12, 2024, 11:05am UTC](https://fenicsproject.discourse.group/t/how-to-solve-a-heat-equation-whose-analysis-solution-is-known-and-expressed-by-spatialcoordinate/15222/4 "2024-07-12T11:05:48Z")

</div>

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