### broken curl-curl formulation of magnetostatic maxwell equations

360

views

0

Dear FEniCS Users,

I am experiencing problems with the solution of the magnetostatic Maxwell equations using the curl-curl formulation ```\text{curl} (\nu \text{curl} A) = j```

I currently have two implementations which both show some strange behaviour:

For both systems I tried a direct solver (mumps) as well as an iterative method (CG+AMG) and finally get the following strange behaviour:

Here the solutions using Lagrange elements with CG+AMG (looks correct by shows wrong scaling):

File attached: curl_curl.tar.gz (4.83 MB)

Does anyone have experience using the curl-curl formulation?

It seems that there is some problem with a gauging when using Lagrange elements.

I have no explaination for the behaviour of the nedelec formulation (if CG works, why should the direct method fail?)

I would appreciate any advice

best wishes

Florian

I am experiencing problems with the solution of the magnetostatic Maxwell equations using the curl-curl formulation ```\text{curl} (\nu \text{curl} A) = j```

I currently have two implementations which both show some strange behaviour:

- using Lagrange elements and adding an additional gauging term to the weak form:

``` int \text{curl}(u) \text{curl}(v) + \text{div}(u)*\text{div}(v) = int j * v ``` - using Nedelec elements with implicit gauging:

``` int \text{curl}(u) \text{curl}(v) = int j * v ```

For both systems I tried a direct solver (mumps) as well as an iterative method (CG+AMG) and finally get the following strange behaviour:

- Lagrange elements with CG solver looks like a correct solution, but the field amplitude is off by a factor of 10 (lagrange_cg.png)
- Lagrange elements with MUMPS show the same behaviour
- Nedelec elements with CG solver seem to work, correctly (we also verfied these results via some COMSOL simulations)
- BUT Nedelec elements with MUMPS give completely wrong results (nedelec_mumps.png)

Here the solutions using Lagrange elements with CG+AMG (looks correct by shows wrong scaling):

Here how the broken MUMPS solution using Nedelec elements looks like:

Finally here are the used codes for the 4 different implementation as well as the used mesh:

File attached: curl_curl.tar.gz (4.83 MB)

And here a short summary of the calculated fieldstrength within the airgap of the magnetic yoke:

```
lagrange_cg.out B(42.5, 0, 0): [ 3.13337283e-10 -4.51505482e-10 -8.50192787e-08]
lagrange_mumps.out B(42.5, 0, 0): [ -7.76810904e-11 -5.23699844e-10 -1.32724796e-07]
nedelec_cg.out B(42.5, 0, 0): [ 1.07924577e-08 -6.51315344e-09 -1.00436369e-06] >> correct
nedelec_mumps.out B(42.5, 0, 0): [ 3.09678597e-07 5.34212960e-07 -4.00908910e-06]
```

Does anyone have experience using the curl-curl formulation?

It seems that there is some problem with a gauging when using Lagrange elements.

I have no explaination for the behaviour of the nedelec formulation (if CG works, why should the direct method fail?)

I would appreciate any advice

best wishes

Florian

Community: FEniCS Project

### 1 Answer

0

Thanks for all contributions. I think our curl-curl problem is now solved using edge elements in combination with a regularization term.

Alternatively a special preconditioner could be used as mentioned in the following post:

solverpreconditioner-for-the-gauged-vector-potential-formulation-curl-curl-equation

best wishes

Florian

Alternatively a special preconditioner could be used as mentioned in the following post:

solverpreconditioner-for-the-gauged-vector-potential-formulation-curl-curl-equation

best wishes

Florian

Please login to add an answer/comment or follow this question.

After glancing at your code, what are your boundary conditions? In the Nedelec examples you seem to have $n\times u=0$

n×u=0 and in the CG examples you have $u=0$u=0 . Perhaps I'm missing the details of your formulation.Regarding the boundary conditions the current implementations always use u=0 (also for nedelec elements), because i am not sure how to force n x u = 0 on arbitrary boundaries!?

For testing reasons I tried n x u = 0 for our cubic airbox, but it does not show significant differences from the u=0 results.

best wishes

Florian

H(curl;Ω) conforming. Setting the DirichletBC to zero on the boundary implies strongly enforcing $n\times u=0$n×u=0 . Perhaps examine the great demo here to be sure you understand what you're doing numerically.thanks for your comment about the boundary conditions. What I would like to have if ```n \times u = 0```. I thought setting DirichletBC((0,0,0)) would mean ```u=0```. If I understand you correctly, this should be equivalent to ```n \times u = 0``` since I am using Nedelec Elements. (Just to be save, i also tried to explicitly set the transversal components of u to zero for a simple box geometry, which gave similar results).

Regarding our problem I additionally want to mention that the type of boundary conditions (Neumann, Dirichlet, ...) should not have a large influence on the solution, if the airbox is large enough. The analytical solution should be 0 at infinity. Therefor it should lead to similar results if ```u=0``` or ```{du}/{dn}=0``` is used.

We already found a workaround for the problem with the curl-curl formulation:

Using the following regularization term seems to solve the problem with the non-converging direct method (MUMPS):

```

a = 1/{mu*mu0}*curl(v)*curl(A))*dx + 1e-6 * 1/{mu*mu0}*v*A*dx

```

This seems to be related with the following post

https://fenicsproject.org/qa/11451/nedelec-conforming-elements-converge-simple-test-problem/

As already asked in this previous post: "I always thought, the Coulomb gauge condition is already intrinsically embedded in the Nedelec elements, where the ansatz functions have zero divergence per construction?"

Could someone explain to me, why this additional regularization is needed?

Additionally the original problem using Lagrange elements is still an open problem!

thanks again

Florian

The main problem is, that the normal component of a Nedelec function does not have to be continuous at the element boundaries (e.g. the edges in 2D).

Hence these boundaries can act as (very small) "sources" and "drains".

Thus while you still have zero divergence locally on every element, the total flux can be different from zero.

Therefore just using Nedelec elements, does not guarantee to find a good approximation for the function with zero divergence and your solution is not unique.

And this is most likely the reason why the direct solver fails in your original questions and you need a stabilizer.

Just to clear up:

Using DirichletBC((0,0,0)) for Nedelec elements means - as far as I know - that you set the dofs at the boundary to zero. The dofs for Nedelec elements are on the edges. That means, that you just set the tangential component of the function to zero, i.e. $u\times n=0$

u×n=0.That means that enforcing $u=0$

u=0 is very hard with Nedelec elements.So, that's the reason why you should get same results if you set the tangential component manually.

That's all I can contribute to your issue (for the other things I do not know enough theory (yet) :)).

Still hope it helps.