It seems like you are using a slightly out of date dolfinx.
I’ve made some changes to your code, including using locate_dofs_topological for the boundary conditions, as you cannot use locate_dofs_geometrical on higher order Serendipity spaces, as the dofs are defined through intergral moments.
Here is my altered version, that runs for me:
import numpy as np
import ufl
from dolfinx import log, fem
from dolfinx import __version__ as version
from dolfinx.common import git_commit_hash
from dolfinx.mesh import CellType, create_unit_square, locate_entities_boundary
from mpi4py import MPI
from petsc4py import PETSc
# Save all logging to file
log.set_output_file("log.txt")
#
# Next, various model parameters are defined::
print(f"Running DOLFINx version {version}, git commit hash {git_commit_hash}")
def solve_dh(h, dt=1.0e-2, T=8, decoupled=False, k=2):
domain = create_unit_square(MPI.COMM_WORLD, h, h, CellType.quadrilateral)
P2 = ufl.FiniteElement("S", domain.ufl_cell(), k)
H = fem.FunctionSpace(domain, P2)
x = ufl.SpatialCoordinate(domain)
t = fem.Constant(domain, 0.0)
def phi_e(x):
return np.sin(np.pi * x[0]) * np.cos(np.pi * x[1]) ** 2
def f_3():
return 3 * ufl.pi ** 2 * ufl.cos(ufl.pi * x[1]) ** 2 * ufl.sin(ufl.pi * x[0]) \
- 2 * ufl.pi ** 2 * \
ufl.sin(ufl.pi * x[0]) * ufl.sin(ufl.pi * x[1]) ** 2
# geometric boundary
def boundary_D(x):
return np.logical_or(np.logical_or(np.isclose(x[0], 0), np.isclose(x[0], 1)),
np.logical_or(np.isclose(x[1], 0), np.isclose(x[1], 1)))
facets_D = locate_entities_boundary(
domain, domain.topology.dim-1, boundary_D)
dofs_D = fem.locate_dofs_topological(H, domain.topology.dim-1, facets_D)
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
phi_e_bc = fem.Function(H)
phi_e_bc.interpolate(phi_e)
# Weak statement of the equations
# terms for Eq. 13
phi_e = ufl.TrialFunction(H)
psi = ufl.TestFunction(H)
# Eq. 13
F3 = ufl.inner(ufl.grad(phi_e), ufl.grad(psi)) * ufl.dx \
- f_3() * psi * ufl.dx
a3 = fem.form(ufl.lhs(F3))
L3 = fem.form(ufl.rhs(F3))
A3 = fem.petsc.assemble_matrix(a3, bcs=[fem.dirichletbc(phi_e_bc, dofs_D)])
A3.assemble()
b3 = fem.petsc.assemble_vector(L3)
phi_e_m1, phi_e_m2 = fem.Function(H), fem.Function(H)
phi_e = fem.Function(H)
# numbered by step, not equation
solver1 = PETSc.KSP().create(domain.comm)
solver1.setOperators(A3)
solver1.setType(PETSc.KSP.Type.PREONLY)
pc1 = solver1.getPC()
pc1.setType(PETSc.PC.Type.LU)
# interpolate exact solution at initial time
t.value = 0.0
phi_e_m2.interpolate(phi_e)
t.value += dt
phi_e_m1.interpolate(phi_e)
# time stepping
print(f"starting timestepping...")
while t.value < T:
t.value += dt
phi_e_bc.interpolate(phi_e)
A3.zeroEntries()
fem.petsc.assemble_matrix(
A3, a3, bcs=[fem.dirichletbc(phi_e_bc, dofs_D)])
A3.assemble()
solver1.setOperators(A3)
with b3.localForm() as loc_b:
loc_b.set(0)
fem.petsc.assemble_vector(b3, L3)
# Apply Dirichlet boundary condition to the vector
fem.apply_lifting(b3, [a3], [[fem.dirichletbc(phi_e_bc, dofs_D)]])
b3.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES,
mode=PETSc.ScatterMode.REVERSE)
fem.set_bc(b3, [fem.dirichletbc(phi_e_bc, dofs_D)])
r = solver1.solve(b3, phi_e.vector)
phi_e.x.scatter_forward()
# Compute L2 error and error at nodes
error_e_L2 = np.sqrt(domain.comm.allreduce(fem.assemble_scalar(
fem.form((phi_e - phi_e_bc) ** 2 * ufl.dx)), op=MPI.SUM))
print(f"completed solve with h = {h}")
return None
ns = [2, 4, 8, 16, 32]
errors = [solve_dh(h=N, T=2, dt=1.0/(10*N), decoupled=True, k=2) for N in ns]
and returns
Running DOLFINx version 0.3.1.0, git commit hash 2f63364f479899b1500490d9f27869a83d6a5d8e
starting timestepping...
completed solve with h = 2
starting timestepping...
completed solve with h = 4
starting timestepping...
completed solve with h = 8
starting timestepping...
completed solve with h = 16
starting timestepping...
completed solve with h = 32
starting timestepping...