Nonlinear solve for one component in MixedFunctionSpace

Question

I would like to solve a problem in the MixedFunctionSpace([W, P]), and instead of allocating separate functions in W and P, I start off with up in WP (see also this question about stitching together functions).

At one point, I need to solve a nonlinear equation for u only, though (the tentative velocity update in Chorin's method for Navier-Stokes). The corresponding call

solver.solve(step_problem, ui.vector())

won't work since .vector() isn't valid for subfunctions. Also, I cannot allocate a new function ui from W and fill in the values back into up (see bug #210).

What is a better way to deal with this than

from dolfin import *

mesh = UnitSquareMesh(20, 20)

W = VectorFunctionSpace(mesh, 'CG', 2)
P = FunctionSpace(mesh, 'CG', 1)

WP = MixedFunctionSpace([W, P])

up = Function(WP)
u, p = up.split()

f = Expression(('sin(x[0])', 'cos(x[1])'))

v = TestFunction(W)
F = dot(u, v) * dx - dot(f, v) * dx

solve(F == 0, u)

which yields

*** Error:   Unable to access vector of degrees of freedom.
*** Reason:  Cannot access a non-const vector from a subfunction.
*** Where:   This error was encountered inside Function.cpp.

Answer 1

Hi,

Here is a suggested solution:

from dolfin import *

mesh = UnitSquareMesh(20, 20)

W = VectorFunctionSpace(mesh, 'CG', 2)
P = FunctionSpace(mesh, 'CG', 1)

WP = MixedFunctionSpace([W, P])

up = Function(WP)
# Make a deep copy to create two new Functions u and p (not subfunctions of WP)
u, p = up.split(deepcopy=True) 

f = Expression(('sin(x[0])', 'cos(x[1])'))

v = TestFunction(W)
F = dot(u, v) * dx - dot(f, v) * dx

solve(F == 0, u)

# Now use black magic to put newly computed vector back in mixed space.
# WP.sub(0).dofmap().collapse(mesh) returns the collapsed dofmap and a 
# dictionary (->[1]) that maps the collapsed dofmap to the dofs in the mixed 
# space, i.e. the final values() 
d = WP.sub(0).dofmap().collapse(mesh)[1].values()
up.vector()[d] = u.vector()

Comments

@mikael-mortensen Could you edit your answer to give more detail on why a deepcopy fixes the subfunction space issues, and more detail on the black magic that is d = WP.sub(0).dofmap().collapse(mesh)[1].values()?