Hello !
A few weeks ago, I made a post about applying a gauge to a electromagnetic problem in potential formulation in 3D. I was recommended a mixed formulation to add another equation. As I explored this possibility I discovered some research papers warning about issues of spurious modes and interface conditions with conductors that could be solved with the introduction of Nedelec elements. This supposedly also solves the gauge setting issue that I faced.
However, despite a residual error similar to a Lagrange element formulation, the shape of the field (vector potential and magnetic) is completely wrong for PAIR orders of Nedelec elements even when applying the Method of Manufactured Solution. I now wonder if the issue comes from a misunderstanding in how interpolation works in FEniCSx.
Is this example a correct application of Nedelec elements ? Why are pair orders so wrong ?
Here is a MWE in a vaccum. I had to include 2 object in gmsh, otherwise I get a gmshio error, which is probably unrelated (?), L2-errors are comparable for both elements but the field is completely different.
import gmsh
import dolfinx
from dolfinx.fem.petsc import LinearProblem
import ufl
import numpy as np
from mpi4py import MPI
from ufl import dx, inner
bkg_tag = 0
sphere_tag = 1
# Mesh definition
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 1)
gmsh.model.add("MWE")
domain = gmsh.model.occ.addSphere(0, 0, 0, 1, tag=bkg_tag)
perturbator = gmsh.model.occ.addSphere(0, 0, 0, 1e-1, tag=sphere_tag)
whole_domain = gmsh.model.occ.fragment([(3, domain)], [(3, perturbator)])
gmsh.model.occ.synchronize()
gmsh.model.addPhysicalGroup(3, [1], tag=bkg_tag)
gmsh.model.addPhysicalGroup(3, [2], tag=sphere_tag)
gmsh.option.setNumber("Mesh.Algorithm", 2)
gmsh.model.mesh.generate(3)
gmsh.model.mesh.optimize("Netgen")
# gmsh.fltk.run()
mesh, ct, ft = dolfinx.io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=3)
gmsh.finalize()
# SOURCE
def background_potential(x):
return np.array([-x[1] / 2, x[0] / 2, x[0] * 0])
# -> uniform field
# PROBLEM DEFINITION
V = dolfinx.fem.FunctionSpace(mesh, ("N1curl", 1))
# V = dolfinx.fem.FunctionSpace(mesh, ("N1curl", 2))
# V = dolfinx.fem.VectorFunctionSpace(mesh, ("Lagrange", 1))
A_exact = dolfinx.fem.Function(V)
A_exact.interpolate(background_potential)
A_exact.x.scatter_forward()
f_source = dolfinx.fem.Function(V)
f_source.interpolate(
dolfinx.fem.Expression(
ufl.curl(ufl.curl(A_exact)), V.element.interpolation_points()))
f_source.x.scatter_forward()
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
a = inner(ufl.curl(u), ufl.curl(v)) * dx
L = inner(f_source, v) * dx
# SOLVING
A_h = dolfinx.fem.Function(V)
problem = LinearProblem(
a,
L,
u=A_h,
bcs=[],
petsc_options={
"ksp_type": "gmres",
"pc_type": "eisenstat",
"ksp_atol": 1e-20,
"ksp_max_it": 2000,
},
)
problem.solve()
V_E = dolfinx.fem.VectorFunctionSpace(mesh, ("Lagrange", 5))
A_h_VE = dolfinx.fem.Function(V_E)
A_exact_VE = dolfinx.fem.Function(V_E)
A_h_VE.interpolate(A_h)
A_exact_VE.interpolate(A_exact)
A_error = dolfinx.fem.Function(V_E)
A_error.x.array[:] = A_h_VE.x.array - A_exact_VE.x.array
comm = A_h.function_space.mesh.comm
error = dolfinx.fem.form(inner(A_error, A_error) * dx)
E = np.sqrt(comm.allreduce(dolfinx.fem.assemble_scalar(error), MPI.SUM))
#Residual errors are similar, but fields are different
#E_lagrange = 0.6279271187085357+0j
#E_nedelec = 0.6279271187085358+0j
print(E)
# POST-PROCESS
W = dolfinx.fem.VectorFunctionSpace(mesh, ("Lagrange", 1))
A_solution = dolfinx.fem.Function(W)
A_solution.interpolate(A_h)
B_solution = dolfinx.fem.Function(W)
B_solution.interpolate(
dolfinx.fem.Expression(
ufl.as_vector(
(
A_solution.sub(2).dx(1) - A_solution.sub(1).dx(2), # Dz/Dy - Dy/Dz
A_solution.sub(0).dx(2) - A_solution.sub(2).dx(0), # Dx/Dz - Dz/Dx
A_solution.sub(1).dx(0) - A_solution.sub(0).dx(1), # Dy/Dx - Dx/Dy
)
),
W.element.interpolation_points(),
)
)
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "potential.xdmf", "w") as file:
file.write_mesh(mesh)
file.write_function(A_solution)
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "field.xdmf", "w") as file:
file.write_mesh(mesh)
file.write_function(B_solution)
Any help is appreciated, thank you !