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:
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:
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!