Helmholtz eigenvalue with different domains

1

Hello,

I’m new to FEM, as a first step I’d like to compute the magnetic field around a superconducting (SC) sphere. It is described by a Helmholtz equation (London equations - Wikipedia see 2 equation in “London penetration depth”). As “warmup” I wrote scripts for 1D SC with fields at both edges, and 2D SC with constant fields at both edges, in both cases I got the exponential solution I was looking for.

from dolfin import *

mesh = UnitIntervalMesh( 16 )
V = FunctionSpace( mesh, "CG", 2 )

u = TrialFunction( V )
v = TestFunction( V )

bc1 = DirichletBC( V, Constant(1), lambda x, _: near( x[0], 0, 1e-5) )
bc2 = DirichletBC( V, Constant(1), lambda x, _: near( x[0], 1, 1e-5) )
bc = [bc1,bc2]

a = ( u.dx(0)*v.dx(0) + 10*u*v ) *dx
L = v*Constant(0)*dx

sol =Function( V )
solve( a==L, sol, bc )

plt.figure()
plot(sol)

and

from dolfin import *

mesh = UnitSquareMesh( 16, 16 )
V = VectorFunctionSpace( mesh, "CG", 2 )

u = TrialFunction( V )
v = TestFunction( V )

bc1 = DirichletBC( V, Constant((0,1)), lambda x, _: near( x[0], 0, 1e-5) )
bc2 = DirichletBC( V, Constant((0,1)), lambda x, _: near( x[0], 1, 1e-5) )
bc = [bc1,bc2]

a = ( inner( nabla_grad(u), nabla_grad(v) ) + 100*inner( u, v ) ) *dx
L = inner( v, Constant((0,0)) )*dx

sol = Function( V )
solve( a==L, sol, bc )

plt.figure()
a = plot(sol)
plt.colorbar(a)

For the 2D case, where I introduced a region with the SC (n_s = 10) and a normal region (n_s=0) I do not get the results show for example here. The code I use right now is

from dolfin import *

mesh = UnitSquareMesh( 16, 16 )
V = VectorFunctionSpace( mesh, "CG", 2 )

u = TrialFunction( V )
v = TestFunction( V )

bc1 = DirichletBC( V, Constant((0,1)), lambda x, _: near( x[1], 0, 1e-5) )
bc2 = DirichletBC( V, Constant((0,1)), lambda x, _: near( x[1], 1, 1e-5) )
bc = [bc1,bc2]

k = Expression( "pow( pow( x[0]-.5, 2 )+pow( x[1]-.5, 2 ), .5 )<.2 ? 50 : 0", degree=0 )

a = ( inner( nabla_grad(u), nabla_grad(v) ) + k*inner( u, v ) ) *dx
L = inner( v, Constant((0,0)) )*dx

sol = Function( V )
solve( a==L, sol, bc )

plt.figure()
patch = patches.Circle((.5,.5),.2,edgecolor="r",facecolor='none')
plt.gca().add_patch(patch)
a = plot(sol)
plt.colorbar(a)

I’ve seen that similar problems canbe handled with eigenvalue solvers, but I don’t see how I can incorporate the different domains into that.

2

Is \vec{u} supposed to be a complex phasor \vec{u}: \Omega \rightarrow \mathbb{C}^d? Or \vec{u} : \Omega \rightarrow \mathbb{R}^d? If \vec{u} is real valued then I can’t see how you can achieve a coupling between the x and y components to get the effect you’re looking for.

Also, is the problem well posed even if k isn’t an eigenvalue of the system? There’s a demo for the generalised eigenvalue problem which may help.