Using deal for Maxwell equations

[Simon Schernthanner]

Hello,

i am from TU Graz (Austria) and want to implement maxwell equations using deal. I want to evaluate deal and compare it with our in house software.

I have some basic questions:

1. Our in house software uses "2nd order edge elements" and i am not quite sure if i'm they are equal to the FE_RaviartThomas<3> elements.
   The shape functions of our edge elements are defined in that way, that the line-integral on one edge is 1 and 0 on all other edges. We have 36 edges in total.
   For example, the shape function of our 22nd enge is: (1/4)*(1-z^2)*(1-x^2)*grad(y)
   Beside the possibility that the numbering of the edges is different, is the RaviartThomas element the right one?
2. I want to implement the following functional
      Matrix:      A_ij = integral ( curl(N_i)*curl (v N_j) )
      RHS:        b_i = integral(  N_i*J  )
   where v is 1/(4*pi*10E-7) and J is the current density.
   I assemble the matrix the following way (the implementation of the RHS (current density is a little bit "hacky")):

       for (unsigned int q_point=0; q_point<n_q_points; ++q_point)

      {

        Point<dim> p = fe_values.quadrature_point (q_point);

        for (unsigned int i=0; i<dofs_per_cell; ++i)

          {

        printf("%i %f %f %f\n",i,fe_views.curl(i,q_point)[0],fe_views.curl(i,q_point)[1],fe_views.curl(i,q_point)[2]);

            for (unsigned int j=0; j<dofs_per_cell; ++j)

            {

            cell_matrix(i,j) += 795774.881*(fe_views.curl(i,q_point)*fe_views.curl(j,q_point))*fe_values.JxW(q_point);

            }

            if(p[0]<0.55&&p[0]>0.45&&p[1]<0.05&&p[1]>(-0.05)&&p[2]>-0.25&&p[2]<0.25)

            {

            cell_rhs(i) += 10000*(fe_values.shape_value_component(i,q_point,2));

            }

          }

      }

This is executed on an cube [-1,1]^3. The problem is, that this equation system does not converge. I tried many things like chaning the values of the parameters and so on

I'm planning to use deal for more complex calculations in the future (A,V-A Formulation for rotating electrical machines, muliphysics etc..).

My cc file is appendet. Thank you very much

Simon

 

 

 

[Jean-Paul Pelteret]

Quick observation: From the snippet of code in the body of your message, you are at the very least missing the fe_values.JxW(q_point) contribution for the RHS (i.e. no integration for this contribution).

[Simon Schernthanner]

Thank you very much.

Some times its difficult to see the obvious!

[Jean-Paul Pelteret]

Agreed :-)

[Simon Schernthanner]

i'd like to use ICCG to solve the system. Is this possible without using PETSc?

[Timo Heister]

Do you mean incomplete Cholesky with CG?

We don't have a Cholesky factorization inside deal.II (only through
PETSc and Trilinos). You can try SolverCG+SparseILU though.

> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to the Google Groups
> "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to dealii+un...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.



--
Timo Heister
http://www.math.clemson.edu/~heister/

[Markus Bürg]

Dear Simon,

1. Our in house software uses "2nd order edge elements" and i am not quite sure if i'm they are equal to the FE_RaviartThomas<3> elements.
   The shape functions of our edge elements are defined in that way, that the line-integral on one edge is 1 and 0 on all other edges. We have 36 edges in total.
   For example, the shape function of our 22nd enge is: (1/4)*(1-z^2)*(1-x^2)*grad(y)
   Beside the possibility that the numbering of the edges is different, is the RaviartThomas element the right one?

I doubt it. The Raviart-Thomas elements are H(div)-conforming, but, for Maxwell’s equations, you probably want H(curl)-conforming elements. This would be FE_Nedelec. Further, for higher-order elements, we do not place all degrees of freedom on the edges, but use a geometrical, hierarchical ansatz which also assigns degrees of freedom to the interior of faces and the interior of the cell.

2. I want to implement the following functional
      Matrix:      A_ij = integral ( curl(N_i)*curl (v N_j) )
      RHS:        b_i = integral(  N_i*J  )

This is executed on an cube [-1,1]^3. The problem is, that this equation system does not converge. I tried many things like chaning the values of the parameters and so on

This has probably two reasons. First, try with FE_Nedelec, then things should improve. Second, your problem is not uniquely solvable, since the null space of your matrix is quite large.

Best,

Markus

[Che Gl]

Dear alls and Dear Simon ,

hope i can post my question here, because it is more or less general question based on your problem Simon.
I'm  planning to do some calculation for time-harmonic electromagnetic problem in my thesis (i.e rotating electrical machine, maybe also multiphysic. let see how far i can go ..). After searching many c++ library i have found two most advanced libraries (for me) : libmesh and deal II (my first impression). But Deal II seem to be more suitable with high order nedelec element, which is not included in the current version of libmesh.
So before start, i just have some question:
- Simon : are you successful in implementing Deal II for any application you planned?
- Did anyone implement successfully the Deal II code for any time-harmonic electromagnetic problem?, If not what is about the effort to realize it from the sight of you as experts (i know stupid question, but time for the thesis is limited, and my background is: just successfully implemented my own FEM-c++ code included matrix creation, direct() and iterativ(GMRES) solver for structure mechanic with only quad4 element in former project work ).

I would appreciate any suggestion.

Best regards

To

[Guido Kanschat]

Dear all,

As an elaboration of Markus' answer, the curl-curl system indeed has no unique solution,  since all gradients of functions in H1 have zero curl.

If any of this is new to you, I  strongly recommend browsing through the book by Peter Monk. It has all the relevant information on finite elements for Maxwell equations.

Best,
Guido

--

[Che Gl]

Dear all and dear Guido,

thanks for your answer and your recommendation. Could you please explain more the about the problem with nedelec-element implemented in DEALII for electromagnetic? or more detail what need to be implemented additionally on the existing nedelec-element in DEALII, so that i.e. time harmonic problem can be solved?.
I know, i have to read the book as well as learning the dealII examples first. But as someone with engineer background the book from peter monk seems to be too much expensive for the first thought: which tools can be used to reach the aims (calculation of some time harmonic problem) ?

Thanks in advance for your answer.
To

[Guido Kanschat]

Dear To,

there are a lot of lecture notes on finite elements for Maxwell equations online; just ask your favorite search engine. It is a topic that requires some mathematical investigation if you do not want to run into possibly completely wrong approximations. Using Nedelec elements is the first step into the right direction, but there are more issues which are too complex for this mailing list.

Best,

Guido

[Che Gl]

Dear Guido,
many thanks for your effort in explaination.
After searching around and reading some other explanation, i would have some further questions:
+ in respect to to the dissertation
http://www.hpfem.jku.at/publications/szthesis.pdf

there are some explanation:

"Due to the non-trivial, large kernel of the curl-operator - the gradient fields of H1- not only the convergence analysis but also the iterative solution of discretized Maxwell problems becomes very challenging. The main difficulty stems from the different scaling of solenoidal and irrotational fields in the curl-curl problem. This leads to very ill-conditioned system matrices and standard Schwarz-type preconditioners like multigrid/multilevel techniques yield only bad convergence behavior. In this context, we want to mention the pioneering works on robust preconditioning in H(curl) and H(div) by Arnold and Hiptmair"

"deals with parameter-robust Schwarz-type preconditioners for H(curl). We prove that parameter-robust solvers are obtained also even by simple Schwarz-type smoothers, if the presented conforming high-order FE-basis is used. In the second part the concept of reduced basis gauging is introduced. Finally, we present numerical tests for magneto-static and magneto-quasi-static problems, illustrating the benefits of our methods."

As far as i understand, to solve electromagnetic  problem two main things need to be considered
+ to improve the convegence behavior (obviously the problem Simon had as described in the first thread), the appropriated preconditioner (Schwarz-type) for H(curl) need to be applied.
+ Also using element of high-order FE-basis could lead to the improvement of convegence.
Since DEAL II supports hp-refinement and Nedelec elements of any order, so from my point of view i would see two challenges at the moment:
              - find the right appropriated preconditoner to solve the problem with DEAL II?
              - and apply the high order elemnent with conforming shape funtion as stated in the dissertation. Are those conformed shape function implemented for high order Nedelec element in DEAL II? or do you mean , it is a challenge to implement those elements in DEAL II?

I would appreciate any suggestion, means i 'm really at the start ... and thanks in advance for your answer.

Best regards.
To

[Wolfgang Bangerth]

To,

> thanks for your answer and your recommendation. Could you please explain more
> the about the problem with nedelec-element implemented in DEALII for
> electromagnetic? or more detail what need to be implemented additionally on
> the existing nedelec-element in DEALII, so that i.e. time harmonic problem can
> be solved?.

Nothing, rally. The building blocks are all essentially there.

What Guido was pointing out is that if you want to solve the equation

curl curl B = 0

then that by itself is difficult because there are many fields B that satisfy
this. You need additional conditions or equations, such as for example
div B = 0
Getting this right is not a question of what deal.II provides, but what
exactly the problem it is you want to solve and in particular how it is
formulated. You will need to read papers solving similar problems to the ones
you are interested in.

Now, you say that you want to solve time harmonic problems, which I assume is
-omega^2 B + curl curl B = 0
I do not know whether you need additional conditions for this (I suspect that
you do not need them from a mathematical perspective, but that you want them
for physical correctness). The literature on similar problems may help you
find what the proper weak formulation for this problem is.

Best
Wolfgang

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

[Wolfgang Bangerth]

> As far as i understand, to solve electromagnetic problem two main things need
> to be considered
> + to improve the convegence behavior (obviously the problem Simon had as
> described in the first thread), the appropriated preconditioner (Schwarz-type)
> for H(curl) need to be applied.
> + Also using element of high-order FE-basis could lead to the improvement of
> convegence.
> Since DEAL II supports hp-refinement and Nedelec elements of any order, so
> from my point of view i would see two challenges at the moment:
> - find the right appropriated preconditoner to solve the
> problem with DEAL II?
> - and apply the high order elemnent with conforming shape
> funtion as stated in the dissertation. Are those conformed shape function
> implemented for high order Nedelec element in DEAL II? or do you mean , it is
> a challenge to implement those elements in DEAL II?

You confuse two separate things. The first "convergence" is about how fast an
iterative solver of the linear system you will get on every mesh converges. In
other words, we are talking about the number of iterations of, say, the CG or
GMRES method you will need to perform to solve the linear system.

The second "convergence" is about how fast the discretization error decreases
with the mesh size (e.g., whether the error goes like O(h^2) or O(h^3) where h
is the size of each cell).

These things are independent. You address the first by using a better
preconditioner, and the second by using higher order elements.

[Simon Schernthanner]

Hello,

fortunately someone other has taken up my idea to use deal.ii for Maxwell equations.

Im quite aware that the problem is not unique solvable, but after a short review the ungauged A formulation seems to be state of the art.

In our in house software uses edge-elements with 36 DoFs like described in [1]. In [2] there is a review about the type of problem i want to solve

and there are also some remarks about non unique solvable equation systems. I would be very interested in implementing capabilities for maxwell equations

into deal.ii, but as i am a electrical engineer and not an mathematician it's very hard for me to get in.

Thank You

[1] A. Kameari, "Calculation of transient 3D eddy current using edge 

elements," IEEE Trans. Magn., vol. 26, pp. 466-469, March 1990

[1] O. Biro, K. Preis, K. R. Richter  “On the Use of the Magnetic Vector 

Potential in the Nodal and Edge Finite Element Analysis of 3D 

Magnetostatic Problems”. IEEE Transactions on Magnetics, Vol 

32, No. 5, 1996.