Hello,
For an exercise class in numerics I wanted to setup a simple example for the curl-curl operator (magneto-statics) but I cannot seem to get the correct convergence rates for the “N1curl” space, so I stripped to problem down to a simple L^2 projection into the space for some function with homogeneous zero tangential trace. It turns out that when I set the the boundary dofs (which in the case of N1curl are exactly the line integrals of the given vector field) to zero explicitly in the space, I get reduced h^{1/2} convergence rate in the L^2 norm and no convergence at all in the H^curl semi-norm. If I increase the degree of the discrete space in general I get h^{k-1/2} for the L^2 norm and h^{k-3/2} for the semi-norm. If I perform the projection in the full space (uncommenting line 88 in the attached code) the optimal rate is restored. All the documentation I found here in the forum indicates that I am setting the boundary conditions the correct way. What is the issue? I am trying to debug the solution in Paraview, but since I have to interpolate it in a DG space, it’s a bit useless for debugging what is happening in the conforming space. I cannot actually really visualise what I am computing…
Thanks in advance for any help.
import numpy as np
import math
from dolfinx import fem, mesh
from ufl import dx, SpatialCoordinate, \
sin, cos, \
curl, inner, \
TrialFunction, TestFunction, as_vector, \
FiniteElement
import matplotlib.pyplot as plt
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.io import VTXWriter
refinements = 8 # number of mesh refiments
maxh = np.zeros((refinements,1)) # maximum mesh size array
error_u_l2 = np.zeros((refinements,1))
error_u_hcurl = np.zeros((refinements,1))
meshsize = 0.5 # initial meshsize
for k in range(0,refinements) :
msh = mesh.create_rectangle(MPI.COMM_WORLD,
[np.array([-1, -1]), np.array([1, 1])],
[2**(k+1), 2**(k+1)],
mesh.CellType.triangle)
msh.topology.create_connectivity(msh.topology.dim-1, msh.topology.dim)
x = SpatialCoordinate(msh)
# Exact solutions and source term
u_exact = as_vector(( cos (math.pi*x[0])*sin (math.pi*x[1]),
-sin (math.pi*x[0])*cos (math.pi*x[1])))
curlu_exact = -2.0 * math.pi * cos (math.pi*x[0])*cos (math.pi*x[1])
# curlu_exact = fem.Constant(msh,0.0)
# Define boundary conditions
def gamma1(x):
return np.isclose(x[1], -1) # only the bottom boundary
def gamma2(x):
return np.isclose(x[0], 1) # only the right boundary
def gamma3(x):
return np.isclose(x[1], 1) # only the upper boundary
def gamma4(x):
return np.isclose(x[0], -1) # only the left boundary
V = fem.FunctionSpace(msh, FiniteElement("N1curl", msh.ufl_cell(), 1))
boundary_facets_gamma1 = mesh.locate_entities_boundary(msh, 1, gamma1)
boundary_dofs_gamma1 = fem.locate_dofs_topological(V, 1, boundary_facets_gamma1)
boundary_facets_gamma2 = mesh.locate_entities_boundary(msh, 1, gamma2)
boundary_dofs_gamma2 = fem.locate_dofs_topological(V, 1, boundary_facets_gamma2)
boundary_facets_gamma3 = mesh.locate_entities_boundary(msh, 1, gamma3)
boundary_dofs_gamma3 = fem.locate_dofs_topological(V, 1, boundary_facets_gamma3)
boundary_facets_gamma4 = mesh.locate_entities_boundary(msh, 1, gamma4)
boundary_dofs_gamma4 = fem.locate_dofs_topological(V, 1, boundary_facets_gamma4)
ubc = fem.Function(V)
class exact_expression():
def __call__(self, x):
values = np.zeros((2, x.shape[1]),dtype=PETSc.ScalarType)
values[0] = np.cos (math.pi*x[0])*np.sin (math.pi*x[1])
values[1] = np.sin (math.pi*x[0])*np.cos (math.pi*x[1])
return values
ubc.interpolate(exact_expression())
# ubc.vector.set(0.0) #should also work...
bcu_gamma1 = fem.dirichletbc(ubc, boundary_dofs_gamma1)
bcu_gamma2 = fem.dirichletbc(ubc, boundary_dofs_gamma2)
bcu_gamma3 = fem.dirichletbc(ubc, boundary_dofs_gamma3)
bcu_gamma4 = fem.dirichletbc(ubc, boundary_dofs_gamma4)
bc = [bcu_gamma1, bcu_gamma2, bcu_gamma3, bcu_gamma4]
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
# Assemble LHS matrix and RHS vector
a = fem.form( (inner(u,v)) * dx )
L = fem.form( inner(u_exact, v) * dx )
A = fem.petsc.assemble_matrix(a, bc) # this fails
# A = fem.petsc.assemble_matrix(a) # this gives the expected rate
A.assemble()
b = fem.petsc.assemble_vector(L)
# Create and configure solver
ksp = PETSc.KSP().create(msh.comm)
ksp.setOperators(A)
ksp.setType("preonly")
ksp.getPC().setType("lu")
ksp.getPC().setFactorSolverType("superlu_dist")
# Compute the solution
u_h = fem.Function(V)
ksp.solve(b, u_h.vector)
# Interpolate into the DG1 space and export to Paraview
DG1_vec = fem.VectorFunctionSpace(msh, ("DG", 1), dim=2)
u_dg1 = fem.Function(DG1_vec, name="A")
u_dg1.interpolate(u_h)
with VTXWriter(MPI.COMM_WORLD, "curlcurl_solution_paraview.bp", [u_dg1]) as f:
f.write(0.0)
# Compute errors
maxh[k] = meshsize
err_l2 = fem.assemble_scalar(fem.form( (inner(u_h - u_exact,u_h - u_exact)) * dx ) )
err_seminorm = fem.assemble_scalar(fem.form( (curl(u_h) - curlu_exact)**2 * dx ) )
error_u_l2[k] = np.sqrt(msh.comm.allreduce(err_l2, op=MPI.SUM))
error_u_hcurl[k] = np.sqrt(msh.comm.allreduce(err_l2+err_seminorm, op=MPI.SUM))
if k > 0:
print('Convergence rate in L2 norm at refinement step ',k,
' is ',
(np.log(error_u_l2[k]) - np.log(error_u_l2[k-1]) )/ \
(np.log(maxh[k]) - np.log(maxh[k-1]) ) )
print('Convergence rate in Hcurl norm at refinement step ',k,
' is ',
(np.log(error_u_hcurl[k]) - np.log(error_u_hcurl[k-1]) )/ \
(np.log(maxh[k]) - np.log(maxh[k-1]) ) )
meshsize = 0.5*meshsize
plt.plot(maxh[:], error_u_l2[:], label='$||u - u_h||_{L^2}$', marker='*', markersize='6.0', linestyle='-')
plt.plot(maxh[:], error_u_hcurl[:], label='$||u - u_h||_{Hcurl}$', marker='o', markersize='6.0', linestyle='-')
plt.loglog()
plt.xlabel('$h$')
plt.ylabel('Error')
plt.grid(True)
plt.legend()
plt.show()
plt.savefig('rates.png')