broken curl-curl formulation of magnetostatic maxwell equations

4 months ago by
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:
  1. 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 ```
  2. 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)
Regarding boundary conditions I use zero Dirichlet Boundary condations at the domain boundary. For the nedelec elements I also tried to set only the tangential component to zero, but without any significant changes.

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
Community: FEniCS Project
Just as an addition (the source files are also included within the curl_curl.tar.gz file from above) here a combined sourcefile:

import numpy as np
from dolfin import *

parameters["allow_extrapolation"] = True

# Permeability
mu0 = 4.*np.pi*1e-7
mu_values = [1., # coil
1000., # magnet
1. # air

# Domain structure
dom_c = 1 # Coil
dom_s = 2 # Soft magnet
dom_a = 3 # Air

mesh = Mesh("mesh.xml")
subdomains = MeshFunction("size_t", mesh, "mesh_physical_region.xml")

dx = Measure('dx', mesh)(subdomain_data=subdomains)
V = VectorFunctionSpace(mesh, "CG", 1)
#V = FunctionSpace(mesh, "N1curl", 1)
V0 = FunctionSpace(mesh, "DG", 0)

# Define current
class CurrentExpression(Expression):
def __init__(self, j=-1., **kwargs):
self._j = j

def eval(self, value, x):
value[2] = 0.0

rx, ry = x[0], x[1]
r = np.sqrt(rx**2 + ry**2)
value[0] = -ry/r
value[1] = rx/r

value *= self._j

def value_shape(self):
return (3,)
J = CurrentExpression(j=1./(4.*20.)*mu0*1.e6, degree=1)

# Define material
mu = Function(V0)
mu.vector()[:] = np.array([mu_values[int(ID-1)] for ID in subdomains.array()])

# Define variational problem
A = TrialFunction(V)
v = TestFunction(V)

a = inner(1./(mu*mu0)*curl(v), curl(A))*dx + 1/(mu*mu0)*div(v)*div(A)*dx
#a = inner(1./(mu*mu0)*curl(v), curl(A))*dx
L = inner(v, J)*dx(dom_c)
bc = DirichletBC(V, Constant((0.0, 0.0, 0.0)), 'on_boundary')

# Solve variational problem
A, b = assemble_system(a, L, bc)

#solver = KrylovSolver('gmres')
#prm = solver.parameters
#prm.absolute_tolerance = 1E-7
#prm.relative_tolerance = 1E-4
#prm.maximum_iterations = 1000
solver = LUSolver('mumps')

u = Function(V)
solver.solve(A, u.vector(), b)
B = project(curl(u), V)
#B = project(curl(u), VectorFunctionSpace(mesh, "CG", 1))
File("data/A.pvd") << u
File('data/B.pvd') << B

print "B(42.5, 0, 0):", B(42.5, 0, 0)​
written 4 months ago by Florian Bruckner  
Why not keep things simple and test with a manufactured solution?

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.
written 4 months ago by Nate  
We also simulated a simple cylindical coil with a cylindical iron core. If setting mu_iron = 1 the numerical solution fits to the analytical solution of an air-coil very well. If using mu > 1 we see similar behaviour like in the example above (some solvers do not converge, some show wrong scaling compared with COMSOL and analytic estimations). I will summarize those results and post it here.

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
written 4 months ago by Florian Bruckner  
Nedelec elements are $H\left(\text{curl};\Omega\right)$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.
written 4 months ago by Nate  
Dear Nate,

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

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
written 4 months ago by Florian Bruckner  
Is this the same thing as Gauss's law? You need to enforce it as shown in the fenics-qa link you posted via an involution.
written 4 months ago by Nate  
Regarding the Nedelec elements and the Coulomb gauge conditon:
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.
written 4 months ago by wagnandr  

1 Answer

3 months ago by
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:

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

Similar posts:
Search »