I wonder how to solve a heat equation which parameter is time dependent

General
1

hello, anyone concerned,

I am a freshmen for FEniCSx.
Now I want to solve this problem: Find S= S(x,y,t) satisfying：

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

where \kappa(t) is time dependent, e.g. \kappa(t) = t.
I compute the source f = \exp(-x - y) - 2t^2 \exp(-x - y) - 2t^3 \exp(-x - y) from the exact solution S_{ex} = t \exp(-x - y). I first tried Dirichlet boundary condition in my code, and also compute the error between FEM solution and exact solution.
Here is my code:

import numpy as np

from mpi4py import MPI

from dolfinx import fem, mesh, io
import ufl

t = 0.0  
T = 1.0  
num_steps = 50
dt = T / num_steps  

nx, ny = 50, 50
domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([0, 0]), np.array([1, 1])], [nx, ny], mesh.CellType.triangle)
V = fem.FunctionSpace(domain, ("Lagrange", 1))

def initial_condition(x):
    return np.exp(-t * (x[0]**2 + x[1]**2))

S_n = fem.Function(V)
S_n.name = "S_n"
S_n.interpolate(initial_condition)

with io.XDMFFile(domain.comm, "diffusion.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    Sh = fem.Function(V)
    Sh.name = "Sh"
    Sh.interpolate(initial_condition)
    xdmf.write_function(Sh, t)

S = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

t_c = fem.Constant(domain, t)
x = ufl.SpatialCoordinate(domain)
S_exact = t_c * ufl.exp(-x[0] - x[1])

def boundary_D(x):
    return np.logical_or(np.isclose(x[0], 0), np.isclose(x[0], 1))
dofs_D = fem.locate_dofs_geometrical(V, boundary_D)

f = ufl.exp(-x[0] - x[1]) - 2 * (t_c**2) * ufl.exp(-x[0] - x[1]) - 2*(t_c**3) * ufl.exp(-x[0] - x[1])

a = S * v * ufl.dx + t_c * dt * ufl.dot(ufl.grad(S), ufl.grad(v)) * ufl.dx
L = (S_n + dt * f) * v * ufl.dx

from dolfinx.fem.petsc import LinearProblem

S_exact_func = fem.Function(V)
S_exact_func.name = "S_exact"

errors = []

for i in range(num_steps):
    t += dt
    t_c.value = t

    S_bc = fem.Function(V)
    S_bc.interpolate(lambda x: t * np.exp(-x[0] - x[1]))
    bc = fem.dirichletbc(S_bc, dofs_D)

    problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
    Sh = problem.solve()

    S_exact_func.interpolate(lambda x: t * np.exp(-x[0] - x[1]))

    error_L2 = fem.assemble_scalar(fem.form((Sh - S_exact_func)**2 * ufl.dx))**0.5
    errors.append(error_L2)

    with io.XDMFFile(domain.comm, "diffusion.xdmf", "a") as xdmf:
        xdmf.write_function(Sh, t)

print("error_L2:", errors)

However, the error seems too big. I trust there is something wrong, but I cannot discover it. Can anyone help me? Thanks very mush!

2

If S_ex is t e^{-x-y}
the initial condition for S(t=0) is 0, not

Secondly,

\frac{\partial S}{\partial t} = e^{-x-y}
\nabla \cdot (t\nabla S) = t \nabla \cdot S = t^2(2e^{-x-y})

thus your source is wrong.

3

Thank you very much.

Firstly, I noticed the error in the code and have corrected it.
Secondly, I wrote the wrong exact solution, mistakenly using the exact solution from another piece of code. Now, I have modified it and here is the new code.

Could you please take a look and see if there are any other issues?

import numpy as np
from mpi4py import MPI
from dolfinx import fem, mesh
import ufl
from petsc4py import PETSc

t = 0.0  
T = 1.0  
num_steps = 10
dt = float(T / num_steps)

nx = ny = 10
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 (t ** 2 + t + 1) * np.exp(-x[0] - x[1])

S_n = fem.Function(V)
S_ex = exact_solution(t)
S_n.interpolate(S_ex)

S = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

def boundary_D(x):
    return np.logical_or(np.isclose(x[0], 0), np.isclose(x[0], 1))
dofs_D = fem.locate_dofs_geometrical(V, boundary_D)

class source_term():
    def __init__(self, t):
        self.t = t
    
    def __call__(self, x):
        return 2 * t * (-t ** 2 - t - 1) * np.exp(-x[0] - x[1]) + (2 * t + 1) * np.exp(-x[0] - x[1])

f = fem.Function(V)
source = source_term(t)
f.interpolate(source)

kappa = fem.Constant(domain, t)
time_step = fem.Constant(domain, dt)

a = S * v * ufl.dx + kappa * time_step * ufl.dot(ufl.grad(S), ufl.grad(v)) * ufl.dx
L = (S_n + time_step * f) * v * ufl.dx

from dolfinx.fem.petsc import LinearProblem

errors = []

for i in range(num_steps):
    t += dt
    S_ex.t = t
    source.t = t
    kappa.value += float(t)

    S_bc = fem.Function(V)
    S_bc.interpolate(S_ex)
    bc = fem.dirichletbc(S_bc, dofs_D)

    problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
    Sh = problem.solve()

    S_n.x.array[:] = Sh.x.array

    error_each = fem.assemble_scalar(fem.form((Sh - S_bc) ** 2 * ufl.dx))
    errors.append(error_each)

print(errors)