The curl curl equation for vector potential A with the Coulomb gauge

variational formulation
1

Hello

I am trying to solve a simple magnetostatic problem through a curl curl equation with Coulomb gauge.
\nabla \times \nabla \times \mathbf{A} = \mu_0 \mathbf{J} and \nabla \cdot \mathbf{A} = 0 on \Omega.
with homogeneous essential conditions on the exterior boundary \hat{\mathbf{n}} \times \mathbf{A} = 0 on \partial \Omega.

I have a simple domain, which has a circular coil (rectangular cross-section) in a cube. (I could not upload pictures here. Sorry.)

I adopted Lagrange multipliers for the Coulomb gauge and a mixed element like following. (I guess this combination of element orders is stable.)

el_u = basix.ufl.element('N1curl', domain.basix_cell(), 2)
el_p = basix.ufl.element('CG', domain.basix_cell(), 1)
el = basix.ufl.mixed_element([el_u, el_p])
W = dolfinx.fem.functionspace(domain, el)

u, p = ufl.TrialFunctions(W)
v, q = ufl.TestFunctions(W)

I wrote a bilinear form and a rhs linear form as following. Earlier, I create a ufl expression j_expr for the current density for the coil.

a = (
    # The stiffness matrix for the curl curl equation.
    ufl.dot(ufl.curl(u), ufl.curl(v)) * ufl.dx + 
    # Lagrangian multipliers p for the Coulomb gauge.
    ufl.dot(ufl.grad(p), v) * ufl.dx +
    # Enforce the Coulomb gauge weakly.
    ufl.dot(u, ufl.grad(q)) * ufl.dx
)

L = ufl.dot(j_expr, v) * ufl.dx

One thing I learned is that I might not want to use ufl.div(u) for N1curl since the normal components may not be continuous and its (weak) divergence may not be defined on facets. So, I used ufl.dot(ufl.grad(p), v) instead of p * ufl.div(v).

I applied homogeneous boundary conditions for both of vector potential u and Lagrange multipliers p like following and tried the preonly solver option (to make sure my variational formulation has some wrong with it.).

domain.topology.create_connectivity(2, 3)
boundary_facets = dolfinx.mesh.exterior_facet_indices(domain.topology)
boundary_dofs = dolfinx.fem.locate_dofs_topological(
    W,
    domain.topology.dim - 1,
    boundary_facets
)

w_D = dolfinx.fem.Function(W)
w_D.x.array[:] = 0.

bc = dolfinx.fem.dirichletbc(w_D, boundary_dofs)

opts = {
    "ksp_type": "preonly",
    "pc_type": "lu",
    "pc_factor_mat_solver_type": "mumps"
}

prob = dolfinx.fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options=opts)

Unfortunately, this formulation gives an unphysical answer. (Actually, all infs). What could be a problem in this approach?

Best,

2

I am not the one who can explain but you probably need both function spaces of same polynomial order.

1 Like
3

Actually, I realized that \nabla CG(order=3) \subset N1curl(order=2). Since grad(q) in CG(3) is tangentially continuous and curl(q) is also integrable. I compared these two since I have a term dot(u, grad(q)) * dx in my bilinear form.

4

Thanks sbhasan. :folded_hands:
You are absolutely right. It works with the same order! Actually, I tried the second before I read you reply. Now, I need to figure out why it works.

1 Like
5

If you are using the H1curl basis for \mathbf{A}, there is no need to enforce a constraint on \nabla\cdot\mathbf{A}. The H1curl basis is divergence free by definition, so the Coulomb Gauge is enforced, by default. Explicit inclusion of the Coulomb Gauge is redundant.

1 Like
6

Hello BillS,
I don’t think this is not the case.
If you mean H(\mathrm{curl}), its definition does not say about the divergence at all (its value or existence).
If you mean N1curl, yes, its divergence may vanish for special cases (the first order N1curl element, affine cell (no curved surface), interior of each cell). But in general, its divergence does not vanish. Its surface normal is not even continuous (its (weak) divergence is not even L^2 integrable.)

1 Like
7

Sorry. Yes I was referring to the classic Nedelec elements. Within each element, the divergence should always vanish for the bog-standard affine simplex elements. On the interelement boundaries, jumps may occur in the normal field component (this is required on material boundaries). The surface normal should exhibit continuity across element boundaries (within numerical limits) for homogeneous material parameters.

If we are talking about global divergence, this is another matter. Spurious solutions, which arise as the result of \nabla\phi\ne 0 and are solutions of the Laplace equation (fictitious charge appearing on interelement boundaries) are also divergence free within elements, but not globally. (The curl operator introduces a null space that corresponds to the Laplace solution.) The weak divergence \int \mathbf{v}\cdot \nabla\phi \, dV in this case will not vanish, but is of course integrable within elements (finite polynomials will always be integrable). Globally integrable, is as you say, may not be the case, as singularities may occur (jumps in direction along sharp boundaries).

One way to handle charge distributions is to introduce the scalar electric potential alongside the vector potential (e.g. see Silvester & Ferrari. Finite Elements for Electrical Engineers or similar for some discussion of eddy current problems).

8

I must correct you here. This is a misconception I’ve heard more often. Compatible discretization does not mean that the basis functions are constructed to satisfy a condition (e.g. \nabla\cdot \boldsymbol{v} = 0 ) . Rather, in mixed form, their compatibility ensures that the condition is satisfied post solving.

Take a look at the classic example of the mixed poisson: Mixed formulation of the Poisson equation with a block-preconditioner/solver # noqa — DOLFINx 0.10.0.0 documentation
where you’ll see that RT elements are also combined with DG elements to enforce the constraint. Intrinsically, RT elements are not divergence free.

1 Like
9

@Stein - Thanks for that! I’ll have a look at the example. I’ll admit, I fell for the misconception.

10

I will checkout Silverster & Ferrari. Thanks :folded_hands:
I also used the magnetic scalar potential to solve a magnetostatic problem for a given \mathbf{H}_{\mathrm{ext}}. This involved solving the Laplace equation, where the sources are the volume and surface magnetic ‘charges’ created by the divergence and the jumps in normal components of the magnetization \mathbf{M}.
\langle \nabla u, \nabla v \rangle = \langle \mathbf{M}(\mathbf{H}_{\mathrm{ext}} - \nabla u ),\nabla v \rangle
I think the solution \nabla u just looks like a projection of \mathbf{M} onto the space \nabla H^1(\Omega) or (\nabla u - \mathbf{M}) \perp \nabla H^1(\Omega).