Issue with DOLFINx solution when solving axisymmetric Maxwell's equations

[mikics]

Hi,

I am trying to solve the axisymmetric version of Maxwell's equations in DOLFINx for a simple problem: an electromagnetic wave scattered by a metallic sphere. I have already solved this problem in legacy DOLFIN, and I have tested the results against analytical ones, getting a very good agreement.

For this problem, I discretize the $\rho$ and $z$ components of the electric field with Nedelec elements, while the $\phi$ component is discretized with Lagrange elements.

The issue in DOLFINx is the following one: whenever I use degree = 2 for the elements, I get quite good results, as confirmed by a comparison between numerical and analytical efficiency. However, whenever I try to use degree = 3 or higher, the solution in DOLFINx is affected by artifacts, and the error between the numerical and the analytical efficiency gets really large. This issue does not affect legacy DOLFIN, where the solution looks clean and agrees with analytical calculation even up to degree = 6. The image here below shows a comparison between:

• $\operatorname{Re}(E_\phi)$ in legacy DOLFIN, for degree = 3
• $\operatorname{Re}(E_\phi)$ in DOLFINx, for degree = 2
• $\operatorname{Re}(E_\phi)$ in DOLFINx, for degree = 3, which has the aforementiond artifacts

In this repository I have uploaded two minimal working examples in legacy DOLFIN and DOLFINx, showing this issue.

[adeebkor]

Hi @mikics

I noticed you are using xdmf for the degree 3 dolfinx output. Have you tried to use vtk or adios.

[mikics]

Hi @mikics

I noticed you are using xdmf for the degree 3 dolfinx output. Have you tried to use vtk or adios.

Hi @adeebkor ! Thank you for your reply.

I am using XDMF for the legacy DOLFIN output, as you can see here, while for the DOLFINx output I am using VTXWriter, as you can see here, so I think I am already using ADIOS2. If it is not, could you be more specific? How should I change my code?

Also: if the problem is only with the visualization, I do not understand why the result for the efficiency changes so dramatically when passing from degree=2 to degree=3, since I calculate it directly from the DOLFINx output.

[jhale]

Can you try passing "pc_factor_mat_solver_type": "mumps" as part of petsc_options?

[mikics]

Can you try passing "pc_factor_mat_solver_type": "mumps" as part of petsc_options?

I've set the solver to mumps as shown here, but unfortunately the problem persists.

[chrisrichardson]

It looks like the artefacts are on the boundary. Do you have BCs set there for the Nedelec space? I wonder if there is a problem with higher order Nedelec and BCs?

[jhale]

I also thought about this, but I think the problem only has natural boundary conditions.

[mikics]

It looks like the artefacts are on the boundary. Do you have BCs set there for the Nedelec space? I wonder if there is a problem with higher order Nedelec and BCs?

I do not have BCs set on that boundary (corresponding to $\rho = 0$), I only have scattering boundary conditions (implemented in the weak form) on the other part of the boundary (a half-circle not visible in the image I posted above).

[mikics]

As a side note: the $\nabla\times$ operator in cylindrical coordinates has two components with a $\frac{1}{\rho}$ term, and since the axis corresponds to $\rho = 0$, I was wondering if this issue is arising due to a division-by-zero (or something close to that). However, as I already said, in legacy DOLFIN I implemented the same formula and I had no problems of this kind.

[chrisrichardson]

Does the interpolation work OK? Can you visualise Eb_m and/or the RHS before solve?

[mikics]

I can interpolate and visualise the different components of Eb_m. To do this, I have added the saving of the different components of Eb_m in this way, and I have compared them for degree = 2 and degree = 3 with paraview, as shown below, but it seems there are no relevant differences or artifacts in the background field. I can try to calculate the difference between the components for degree = 2 and degree = 3, just to be sure there are no subtle differences not visible at a quick visual inspection.

Real part:

Imaginary part:

[chrisrichardson]

Can you assemble the matrix and RHS separately, and check for artefacts (i.e. visualise the RHS vector before solve)?

[mikics]

As a side note: the ∇× operator in cylindrical coordinates has two components with a 1ρ term, and since the axis corresponds to ρ=0, I was wondering if this issue is arising due to a division-by-zero (or something close to that). However, as I already said, in legacy DOLFIN I implemented the same formula and I had no problems of this kind.

As a side note to this comment: I have also noticed that for m = 0, the simulation in DOLFINx works fine for whatever degree, therefore I suspect there might be something wrong with the terms in curl_axis "activated" by m = 1, see here. In general, the value of m only influences background_field and curl_axis, but since background_field has no weird issues as shown above, there might be something wrong with the curl_axis.

[mikics]

Can you assemble the matrix and RHS separately, and check for artefacts (i.e. visualise the RHS vector before solve)?

I will do that!

[jorgensd]

The issue is due to something with mixed spaces. A smaller example can be found here:

import sys

import gmsh
import numpy as np
from dolfinx import fem
from dolfinx.io import VTXWriter
from dolfinx.io.gmshio import model_to_mesh
from mpi4py import MPI
from ufl import FiniteElement, MixedElement, SpatialCoordinate, cos, dot, dx

degree = 3

um = 1
nm = um * 10**-3

radius_sph = 0.025 * um
radius_dom = 1 * um

in_sph_size = 2 * nm
on_sph_size = 0.8 * nm
scatt_size = 40 * nm

au_tag = 1
bkg_tag = 2
scatt_tag = 3

gmsh.initialize(sys.argv)

gmsh.model.add("geometry")

gmsh.model.occ.addCircle(0, 0, 0, radius_sph * 0.5,
                         angle1=-np.pi / 2, angle2=np.pi / 2, tag=1)
gmsh.model.occ.addCircle(0, 0, 0, radius_sph,
                         angle1=-np.pi / 2, angle2=np.pi / 2, tag=2)
gmsh.model.occ.addCircle(0, 0, 0, radius_dom,
                         angle1=-np.pi / 2, angle2=np.pi / 2, tag=3)

gmsh.model.occ.addLine(6, 4, tag=4)
gmsh.model.occ.addLine(4, 2, tag=5)
gmsh.model.occ.addLine(2, 1, tag=6)
gmsh.model.occ.addLine(1, 3, tag=7)
gmsh.model.occ.addLine(3, 5, tag=8)

gmsh.model.occ.addCurveLoop([6, 1], tag=1)
gmsh.model.occ.addPlaneSurface([1], tag=1)
gmsh.model.occ.addCurveLoop([7, 2, 5, -1], tag=2)
gmsh.model.occ.addPlaneSurface([2], tag=2)
gmsh.model.occ.addCurveLoop([4, -3, 8, 2], tag=3)
gmsh.model.occ.addPlaneSurface([3], tag=3)

gmsh.model.occ.synchronize()

gmsh.model.addPhysicalGroup(2, [1, 2], tag=au_tag)
gmsh.model.addPhysicalGroup(2, [3], tag=bkg_tag)
gmsh.model.addPhysicalGroup(1, [3], tag=scatt_tag)

gmsh.model.mesh.setSize([(0, 1)], size=in_sph_size)
gmsh.model.mesh.setSize([(0, 2)], size=in_sph_size)
gmsh.model.mesh.setSize([(0, 3)], size=on_sph_size)
gmsh.model.mesh.setSize([(0, 4)], size=on_sph_size)
gmsh.model.mesh.setSize([(0, 5)], size=scatt_size)
gmsh.model.mesh.setSize([(0, 6)], size=scatt_size)

gmsh.model.mesh.generate(2)

model = gmsh.model

domain, cell_tags, facet_tags = model_to_mesh(
    model, MPI.COMM_WORLD, 0, gdim=2)

gmsh.finalize()


curl_el = FiniteElement("N1curl", domain.ufl_cell(), degree)
lagr_el = FiniteElement("Lagrange", domain.ufl_cell(), degree)
V = fem.FunctionSpace(domain, MixedElement([curl_el, lagr_el]))


def f(x):
    return (x[0], x[1])


def g(x):
    return np.sin(x[1]) + 2*x[0]


Eb_m = fem.Function(V)
Eb_m.sub(0).interpolate(f)
Eb_m.sub(1).interpolate(g)
V0 = fem.FunctionSpace(domain, FiniteElement("DG", domain.ufl_cell(), 3))
curl_expr = fem.Expression(Eb_m[2].dx(1), V0.element.interpolation_points())
u0 = fem.Function(V0)
u0.interpolate(curl_expr)
u0.name = "u0"
with VTXWriter(domain.comm, f"sols/dolfinx/test_deg{degree}.bp", [u0]) as f:
    f.write(0.0)


def dgdy(domain):
    x = SpatialCoordinate(domain)
    return cos(x[1])


print(degree, fem.assemble_scalar(
    fem.form(dot(u0-dgdy(domain), u0-dgdy(domain))*dx)))

degree 2:

2 (6.431107788172592e-10+0j)

degree 3:

3 (2.5202375074077557+0j)

@mscroggs any ideas?

[jorgensd]

The last bit of the example could be further simplified to (ignore expression interpolation for outputting):

Eb_m = fem.Function(V)
Eb_m.sub(0).interpolate(f)
Eb_m.sub(1).interpolate(g)


def dgdy(domain):
    x = SpatialCoordinate(domain)
    return cos(x[1])


diff = Eb_m[2].dx(1)-dgdy(domain)
print(degree, fem.assemble_scalar(
    fem.form(dot(diff, diff)*dx)))

[mscroggs]

Failing example on a two cell mesh:

import numpy as np
from dolfinx import fem
from dolfinx import mesh
from mpi4py import MPI
from ufl import (FiniteElement, MixedElement, SpatialCoordinate, cos, dot,
                 dx, VectorElement, Mesh)

points = np.array([[0., 0.], [2., 0.], [1., 1.], [1., -1.]])
cells = [[0, 1, 2], [1, 0, 3]]

# Randomly number the points and create the mesh
order = list(range(len(points)))
domain = mesh.create_mesh(
    MPI.COMM_WORLD, cells, points,
    Mesh(VectorElement("Lagrange", "triangle", 1)))


def g(x):
    return np.sin(x[1]) + 2*x[0]


def dgdy(domain):
    x = SpatialCoordinate(domain)
    return cos(x[1])


for degree in [1, 2, 3]:
    curl_el = FiniteElement("N1curl", domain.ufl_cell(), degree)
    vlag_el = VectorElement("Lagrange", domain.ufl_cell(), degree)
    lagr_el = FiniteElement("Lagrange", domain.ufl_cell(), degree)

    V = fem.FunctionSpace(domain, MixedElement([curl_el, lagr_el]))
    Eb_m = fem.Function(V)
    Eb_m.sub(1).interpolate(g)
    diff = Eb_m[2].dx(1)-dgdy(domain)
    error = fem.assemble_scalar(fem.form(dot(diff, diff)*dx))

    V = fem.FunctionSpace(domain, MixedElement([vlag_el, lagr_el]))
    Eb_m = fem.Function(V)
    Eb_m.sub(1).interpolate(g)
    diff = Eb_m[2].dx(1)-dgdy(domain)
    error2 = fem.assemble_scalar(fem.form(dot(diff, diff)*dx))

    print(degree, error, error2)
    assert np.isclose(error, error2)

It works fine for MixedElement(vector Lagrange, Lagrange) but not MixedElement(Nedelec, Lagrange)

[mscroggs]

Update: assertion also fails with

    curl_el = FiniteElement("N1curl", domain.ufl_cell(), 1)
    vlag_el = VectorElement("Lagrange", domain.ufl_cell(), 1)
    lagr_el = FiniteElement("Lagrange", domain.ufl_cell(), degree)

So the bug is related to Lagrange being high order not Nedelec

[mscroggs]

This issue is fixed in #2347