Turn a Function into ufl form

Question

The operations between Function and ufl Form seem not possible. I guess its determined by the assemble and other structures that are needed in ufl Form.

But in my formulation I need to solve a Nonlinear Problem. I have found two alternatives to formulate the problem one is to supply the linear and bilinear form in NonlinearVariationalProblem(), another approach is to subclass the NonlinearProblem, where I need to define the user-F and user-J. Both of them fail.

1. In my formulation I need to add Function to a Form in the linear form, then it gives me an error
2. If I subclass NonlinearProblem like the following, it throws the nasty bug Segmentation Error.

class ExtendedNonlinear(NonlinearProblem): def __init__(self, a, L, vec): NonlinearProblem.__init__(self) self.L = L self.a = a self.vec = vec

def F(self, b, x):
    # print('F assemble')
    # assemble(self.L) + self.vec.vector()
    b = PETScVector()
    # assemble(self.L, tensor=b) + self.vec.vector()

    assemble(self.L, tensor=b)
    # b.axpy(1., self.vec.vector())

def J(self, A, x):
    # print('J assemble')
    assemble(self.a, tensor=A)

Is there a trick to turn a Function to Form or resolve the Segmentation Fault error?

Thanks!

Answer 1

You can't add Function and Form, it's not an operation that makes sense mathematically. A (linear) form is a functional, not a function.

What you can do is something similar to

f = Function(V)
v = TestFunction(V)
form += inner(f,v) * dx()

In a sense this adds f to form, but does so in mathematically meaningful way by describing how f acts on the testfunction v.

For examples of subclassing NonlinearProblem, see the Cahn-Hilliard demo or the example below.

class ExtendedNonlinear(NonlinearProblem): 
    def __init__(self, F, u, d, bcs = []):
        NonlinearProblem.__init__(self)
        self.F_form = F
        self.u = u
        self.d = d
        self.J_form = derivative(F, u, TrialFunction(u.function_space()))
        self.bcs = bcs

    def F(self, b, x):
        self.u.assign(Function(u.function_space(), x))
        # assemble then add the vector d
        assemble(self.F_form, tensor = b)
        b += self.d
        for bc in self.bcs: bc.apply(b)

    def J(self, A, x):
        self.u.assign(Function(u.function_space(), x))
        assemble(self.J_form, tensor=A)
        for bc in self.bcs: bc.apply(A)

mesh = UnitSquareMesh(64, 64)
V = FunctionSpace(mesh, "CG", 1)

u = Function(V)
v = TestFunction(V)

F = inner(grad(u), grad(v)) * dx
d = - assemble(Constant(1.) * v * dx)
bcs = [DirichletBC(V, 0., "on_boundary")]

nlp = ExtendedNonlinear(F, u, d, bcs)
solver = NewtonSolver()
solver.solve(nlp, u.vector())

Comments

Thx for the answer!

why do you do the assign in F and J,

self.u.assign(Function(u.function_space(), x))

I have viewed the Cahn Hilliard demo, it did not set the dependent function u as in your example.

class CahnHilliardEquation(NonlinearProblem):
    def __init__(self, a, L):
        NonlinearProblem.__init__(self)
        self.L = L
        self.a = a
    def F(self, b, x):
        assemble(self.L, tensor=b)
    def J(self, A, x):
        assemble(self.a, tensor=A)

---

The assignment makes sure the forms are evaluated at x. Consider for example

nlp = ExtendedNonlinear(F, u, d, bcs)
x = u.vector()
y = Function(V).vector()
f = Function(V).vector()
g = Function(V).vector()

# randomize y
from numpy.random import randn
y[:] = randn(y.local_size())

nlp.F(f, x)
nlp.F(g, y)
print "||F(x) - F(y)|| = {}".format((f-g).norm("l2"))

F(x) and F(y) should be different, but if you comment out the assignment lines you will see that this is not the case.

---

Very informative!

It seems that the update of a Function or a Form is sometimes not very clear.