Preconditioner for 3D electromagnetic problem

Dear all,

I am using NGSolve to validate a few preconditioners for the 3D electromagnetic scattering using Hcurl-conforming elements. Right now I am experimenting with Gauss Seidel iterations on block based preconditioners.

The problem I am trying to solve is the propagation of the dominant mode (obtained by modal analysis) in a step-index optical fiber.

I am evaluating if GMRES converges in a decent rate and what is the initial preconditioned residual (preconditioner applied to rhs vec).

For lowest order elements, the scheme seems to be reasonable when using element blocks. Edge-based blocks are also fine. Zaglmayr’s thesis mentions that vertex-based patches (AFW precond) are supposed to be robust for this equation. However I get a super high initial residual (10**14)

By any chance, have you ever seen such a behaviour?

In case it is relevant, the code is available at this repository and the AFW blocks are created as

for v in mesh.vertices:
     vdofs = set()
     for edge in mesh[v].edges:
         vdofs |= set(d for d in fes.GetDofNrs(edge)
                      if freedofs[d])

Thank you very much in advance.

Hi Francisco,

we see similar behaviour already for the Helmholtz equation:

For the Poisson equation the preconditioner becomes better if you invest more into larger blocks. However, for the Helmholtz equation, if the blocks are geometrically so large that you get a standing wave (with Dirichlet b.c.) in a block, the sub-space correction may go into the wrong direction (and can get arbitrarily bad).

Let me advertise a hybrid method with double facet variables from Martin Huber’s dissertation:
Here iterative solvers work, the wave is propagating further through the domain with every iteration. NGSolve-Python examples for this method are going to come very soon.


Dear Joachim,

Thank you for the insight!

I will definitely take a look at this dissertation and see what could be applied to my problem.

For this particular example, I don’t really understand what is the difference between building the blocks as

blocks = []
freedofs = fes.FreeDofs()
for v in mesh.vertices:
     vdofs = set()
     for edge in mesh[v].edges:
         vdofs |= set(d for d in fes.GetDofNrs(edge)
                      if freedofs[d])


blocks = fes.CreateSmoothingBlocks(subassembled=True)

the latter gives me 125800 blocks and 22 colors. The initial residual is 1.6915 and after 20 iterations the residual is already at 0.0843.

The former, however, results in the same number of blocks, but 81 colors and the initial residual is 373765086541229.

Checking the HCurlHighOrderFESpace::CreateSmoothingBlocks source code, I had the impression that they would result in the same preconditioner. Regarding dirichlet BCs, from what I’ve understood, you still have the vertex patch for a vertex at the BC, but only including edges that are not at the boundary surface.

Hi Francsico,
I found the problem:

returned all edges of elements containing v.
It’s fixed now, the coming nightly will only return those edges containing the vertex v.

Have you found this tutorial: You can also write a C++ snippet to build the blocks:


Ah, now it makes sense! Thank you very much for looking into this.

Wow, it is quite nice to be able to build the blocks directly in c++. I am definitely checking this out.

If I may ask one more question:

The AFW blocks (high order or not) won’t have any entries associated with face/internal dofs, correct? Therefore, I think these remaining dofs won’t be taken into account when using this precond in the left-preconditioned GMRES, as

q = rhs.CreateVector() = A * Q[k] = pre * tmp

will result in having null entries for these positions.
If that is correct, which approach you would recommend for creating the remaining blocks? Simply creating one additional block for each face/volume?

Sorry if the question is dumb, but preconditioners are something that I am still getting acquainted with.

yes, as you said, one block per face (and condensing out the element dofs via static condensation) is a good choice

Thank you very much, Joachim!

Preliminary tests indicate good results for this setting.

If I may ask, would you see any reasons on why the same strategy wouldn’t work on meshes with hexahedra? On my other code, not using NGSolve, I have observed this weird behaviour with high order AFW + face blocks + static condensation of element dofs.