Calculating Jacobian with a mixed FunctionSpace of CG and R

Question

My goal is to solve a nonlinear problem with one additional degree of freedom. In order to demonstrate my problem I will use the Helmholtz equation in 1D

$$ (\nabla^2 + n^2 k^2) u = 0 $$

defined on the unit interval with Neumann bcs on the left and Dirichlet bcs on the right. Given an initial guess for both u and k, I want to use a Newton solver to find a solution of the Helmholtz equation, where k is my additional degree of freedom, i.e. in the weak form,

$$ F = \int_\Omega \nabla u \nabla v dx - n^2 k^2 \int_\Omega u v= 0 $$

In order to make this problem well-posed an additional constraint is necessary to cope with the additional degree of freedom, for which I choose u(0) = 1.

How can I include this constraint correctly in fenics such that the Jacobian of F can be calculated automatically, the code below fails with

ufl.log.UFLException: All terms in form must have same rank.

which is to be expected as F at the moment does not include the constraint as an additional equation.

#!/usr/bin/env python

from dolfin import *

ka = 10.0
gammaT = 4.0

mesh = IntervalMesh(50, 0.0, 1.0)

V = FunctionSpace(mesh, "CG", 1)
R = FunctionSpace(mesh, "R", 0)
W = V * R

k = 3.0*pi/2.0
n = Constant(1.0)
amplitude = Constant(1.0)

u_mixed = interpolate(Expression((("cos(k*x[0])"),("k")), k=k), W)
uw, kw = u_mixed.split()

v = TestFunction(V)

class LeftBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0],0.0)
class RightBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0],1.0)

boundaries = FacetFunction('size_t', mesh)
LeftBoundary().mark(boundaries, 1)
RightBoundary().mark(boundaries, 2)
# plot(boundaries, interactive=True)

bc = DirichletBC(V, Constant(0.0), boundaries, 2)
F = (inner(nabla_grad(uw), nabla_grad(v)) - uw * v * n**2 * kw**2) * dx + (uw - amplitude) * v * ds(1)

# Check if F is correctly defined and almost zero for the analytic solution
Fv = assemble(F, bcs = bc, exterior_facet_domains=boundaries)
print(Fv)
print(Fv.max())

J = derivative(F, u_mixed)
assemble(J, bcs = bc, exterior_facet_domains=boundaries)

Thanks for your help

Answer 1

Hi,

I think you need to use

uw, kw = split(u_mixed)

and TestFunctions for the mixed space, like

v, c = TestFunctions(W)

Your additional constraint should then use c, and not v

(uw - amplitude) * c * ds(1)

I have not seen Dirichlet and Neuman mixed this way though, so I can't really tell whether it will work or not. Best of luck:-)

Comments

Thanks, it finally works perfectly. The nonlinear problem wouldn't converge at first, but that was due to the fact, that I used assemble(F, bcs = bcs,...) instead of applying the boundary conditions myself, which is necessary here.

I have one more question: Is there a way to generalize the additional constraint to two dimensions? Since ds defines a boundary integral, I can no longer use it in 2d for this specific case. Is there an equivalent for "integral over one vertex" in 2d, such that again the corresponding row for the constraint in the Jacobian has only a single entry.