Hello anyone,
I once again do a numerical experiment based on My previous question.
The problem is:
find u = u(\boldsymbol{x}, t) obeying
\begin{alignat}{2}
\frac{\partial u}{\partial t} - \nabla \cdot (\kappa \nabla u) &= f, && (\boldsymbol{x}, t) \in \Omega \times [0,T],\\
-\kappa \frac{\partial u}{\partial \boldsymbol{n}} &= g, &\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),\\
g &=
\begin{cases}
0, & y=0,\\
-6 \kappa, & y=1,\\
0, & x=0,\\
-3 \kappa. & x=1,
\end{cases}\\
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 follow this tutorial to code as:
from petsc4py import PETSc
from mpi4py import MPI
import ufl
from dolfinx import mesh, fem, default_scalar_type
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_matrix, create_vector, set_bc
import numpy as np
import time
start = time.time()
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))
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_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)
g = fem.Function(V)
boundaries = [
(1, lambda x: np.isclose(x[1], 0)),
(2, lambda x: np.isclose(x[1], 1)),
(3, lambda x: np.isclose(x[0], 0)),
(4, lambda x: np.isclose(x[0], 1))
]
facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for (marker, locator) in boundaries:
facets = mesh.locate_entities(domain, fdim, locator)
if locator == 1:
g.x.array[facets] = np.full_like(facets, 0.0, dtype = default_scalar_type)
elif locator == 2:
g.x.array[facets] = np.full_like(facets, -6.0 * kappa_expression, dtype = default_scalar_type)
elif locator == 3:
g.x.array[facets] = np.full_like(facets, 0.0, dtype = default_scalar_type)
elif locator == 4:
g.x.array[facets] = np.full_like(facets, -3.0 * kappa_expression, dtype = default_scalar_type)
else:
g.x.array[facets] = np.full_like(facets, 0.0, dtype = default_scalar_type)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = mesh.meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
ds = ufl.Measure("ds", domain = domain, subdomain_data = facet_tag)
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 + dt * g * v * ds(1) + dt * g * v * ds(2) + dt * g * v * ds(3) + dt * g * v * ds(4)
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
f.interpolate(f_expression)
A.zeroEntries()
assemble_matrix(A, a)
A.assemble()
with b.localForm() as loc_b:
loc_b.set(0)
assemble_vector(b, L)
apply_lifting(b, [a], [[]])
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, [])
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}")
end = time.time()
print(f"Execution time: {(end - start) * 10 ** 3:.6f} ms")
then the system tells me
Traceback (most recent call last):
File "/media/math/个人文件空间/WangRan/work-fenics/SIR/demo4/ex.py", line 69, in <module>
g.x.array[facets] = np.full_like(facets, 0.0, dtype = default_scalar_type)
IndexError: index 37 is out of bounds for axis 0 with size 36
I guess the error results from the definition of g since it is defined on the square, but how to define a function only on facet?
Moreover, is there any simple way to define the function g?
Thanks very much.