Preconditioner for the time harmonic Maxwell equation

[Abbas Ballout]

This isn't a dealii problem and I am sorry if this isn't the right place to post to, but if there is someone who could give me some hints on how I would approach this then it would be you guys. 

I have been trying to solve the time harmonic maxwell equation discritized with the use of Nedelec elements. Details are in This stack overflow post. 

I have an issue with the linear system, since using any black box Trilinos preconditioner yields no convergence or worse convergence than a simple Jaccobi or an identity preconditioner and even in that case the restarted GMRES solver needs more than 400 

vectors to converge and it gets worse with a refined mesh. My system is broken into 2x2 block matrices much like step-29  by the way. 

Creating a multigrid solver like step-16 might solve this issue but I would rather use that as plan B. 

Any hint is appreciated.  

Thanx :3 
Abbas 

[Sebastian Kinnewig]

Hi Abbas,

as you already noticed, standard preconditioners do not yield very good convergence for Maxwell's equations, for more details see Why it is Difficult to Solve Helmholtz Problems with Classical Iterative Methods. Since the scalar Helmholtz problem is closely related to the Maxwell's equations, the same arguments from the Helmholtz problem also hold true for the Maxwell's equations. So not even your plan B would work (see Geometric Multigrid Methods for Maxwell’s Equations, Chapter 6.2).

Therefore you need to apply specially suited preconditioners for the Maxwell's equations. For the problem you are considering, I would recommend the block-preconditioner, which is presented in this paper: Large-scale 3D geoelectromagnetic modeling using parallel adaptive high-order finite element method. As a starting point for your implementation of that block-preconditioner you can use step-22. 

An other approach to solve the Maxwell's equations is via domain decomposition method (for example see: A quasi-optimal domain decomposition algorithm for the time-harmonic Maxwell's equations).

I hope this helps you, best regards
Sebastian

[Abbas]

Thank you Sebestian, 

In case someone comes across this thread in the future 
It seems that the problem that is being solved in the "Large scale 3D geomag..." paper is different than mine. 
The paper's: curl x curl x E + sigma E = J   , sigma > 0 
Mine is        : curl x curl x E  - k^2 E  = 0 
and because of that it looks Hypre's AMS solver can only work with problems as given in the paper

Am I correct here? 

best,
Abbas 

[Alexander]

Hi Abbas

AMS requires the matrix to be (semi-)positive definite. The problem you stated is not going to satisfy this if I understand it right.

Alexander

[Sebastian Kinnewig]

Hello Abbas,

okay, I see, previously there was just a sign error in your stack overflow post, but the sign is very important in this case. So let's sum it up:

The problem you previously stated in your stack overflow post, i.e.
    curl ( curl ( E ) ) + k^2 E = 0
is sometimes called the well-posed Maxwell's equations or the definite Maxwell's equations. As the name suggests this leads to a positive definite system and is therefore easier to solve. 
For example, the definite Maxwell's equations can be solved efficiently via an iterative solver (e.g. GMRES) together with the block-preconditioner given in the in the "Large scale 3D geomag..." paper.
Also AMS is only a solver for the definite Maxwell's equations.

The problem you are considering, i.e.
    curl ( curl ( E ) ) - k^2 E = 0
is respectively sometimes called the ill-posed Maxwell's equations or the indefinite Maxwell's equations.
When you are considering only the 2D case a direct solver should be enough.
But for the 3D case you will need more advanced methods. The most common approach to solve these equations is via domain decomposition method (ddm). The basic idea of the ddm is to divide the the domain into smaller subdomains, where each subdomain becomes small enough so it can be handled by a direct solver. The ddm for Maxwell is explained in greater detail in this paper: A quasi-optimal domain decomposition algorithm for the time-harmonic Maxwell's equations. But it is not easy to build a Maxwell solver based on ddm in deal.II (I know this, since I am currently working on such a solver by myself). 
I am sorry, there is no easy way to do this.

Best,
Sebastian

[Wolfgang Bangerth]

On 5/21/21 2:23 AM, Sebastian Kinnewig wrote:
>
> The problem you are considering, i.e.
>     curl ( curl ( E ) ) - k^2 E = 0
> is respectively sometimes called the ill-posed Maxwell's equations or the
> indefinite Maxwell's equations.
> When you are considering only the 2D case a direct solver should be enough.
> But for the 3D case you will need more advanced methods. The most common
> approach to solve these equations is via domain decomposition method (ddm).
> The basic idea of the ddm is to divide the the domain into smaller subdomains,
> where each subdomain becomes small enough so it can be handled by a direct
> solver. The ddm for Maxwell is explained in greater detail in this paper: A
> quasi-optimal domain decomposition algorithm for the time-harmonic Maxwell's
> equations

> <https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.sciencedirect.com%2Fscience%2Farticle%2Fpii%2FS0021999115001965&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7Cae145f9e73094301870808d91c31afa9%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637571823692274784%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=ZE74xHiT8FLjtTNhQCA68VeVYp3BjnWxDl3o%2Brrn8nM%3D&reserved=0>.

> But it is not easy to build a Maxwell solver based on ddm in deal.II (I know
> this, since I am currently working on such a solver by myself).
> I am sorry, there is no easy way to do this.

This is somewhat outside of my realm of experience, but I'll chime in anyway:

Domain decomposition methods conceptually originated from the desire to use
serial codes that were developed in the 1990s and re-use them in a parallel
context. The question was how to connect them into a parallel code. deal.II
does not follow this paradigm: it just distributes *all* data structures
(mesh, matrices, vectors, etc). As a consequence, it is difficult to press
deal.II into service for DD methods because they philosophically just don't
quite fit into the deal.II world view.

But, I believe you also don't have to because essentially everything that has
been developed in the DD world is also available in the fully-distributed
world deal.II lives in. In the case you describe above, the key piece is the
direct solver that can be run on each subdomain sequentially. But there are
good parallel direct solvers that one could use -- specifically, I would
suggest to explore MUMPS, which you can access via its PETSc interfaces. (I
think I recall that there is also a Trilinos parallel direct solver, but I've
forgotten its name and we don't have an interface to it.) If you were to
consider something like step-40, where the domain is split into subdomains,
the matrix is split into corresponding sub-matrices, then a parallel direct
solver such as MUMPS can conceptually be thought of as inverting the diagonal
blocks of the matrix corresponding to each subdomain -- i.e., not so different
from what a direct solver applied to one subdomain in the DD context does. The
difference is that in the DD context, you then have to think of how to stitch
the individual subdomains back together (the boundary or "transmission"
conditions between subdomains) whereas in the fully-distributed,
one-global-linear-system world this is unnecessary.

Best
Wolfgang

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Thomas WICK]

Dear Wolfgang,


Am 21.05.21 um 14:48 schrieb Wolfgang Bangerth:

Thanks for your email concerning DD and deal.II. In general, we agree
that DD in the deal.II context
is not really necessary. Just using a parallel sparse direct solver
(e.g. MUMPS) is not an option because
of memory.

Maxwell is indeed very special and one cannot simply apply well-known
concepts. Of course, Hiptmair proposed in 1998 a multigrid method for
Maxwell. The key is a special
construction of the smoother (Gauß-Seidel-like / multiplicative Schwarz)
and a decomposition of the function spaces.

Of course this paper is highly cited and used, but it does not seem that
idea resolves all difficulties one may have
in Maxwell because of the indefinite structure.

Recently, over the last years, DD seems to be more competitive (e.g.,
Bouajaj et al. JCP 2015) and for exactly this reason, we started
developing a DD solution within deal.II. Our current results are promising.

We need to couple this later with Navier-Stokes and/or solids and for
this reason,
we decided to go with deal.II, despite that other packages (ngsolve for
instance) could
have been better for Maxwell alone.

Any further comments are specifically welcome of course.

Best regards,

Thomas

---
Thomas Wick
Leibniz University Hannover
https://thomaswick.org







> Best
> Wolfgang
>

--
---
++++--------------------------------------------++++
Prof. Dr. Thomas Wick
Leibniz Universität Hannover (LUH)
Institut für Angewandte Mathematik (IfAM)
Arbeitsgruppe Wissenschaftliches Rechnen (GWR)

Welfengarten 1
30167 Hannover, Germany

Tel.: +49 511 762 3360
Email: thoma...@ifam.uni-hannover.de
www: https://ifam.uni-hannover.de/wick
www: https://thomaswick.org
++++--------------------------------------------++++
---
--