How to correctly update solutions in a dynamic problem using MixedFunctionSpace and operator splitting ?

Question

Hello

I am trying to solve a mixed variational formulation for a dynamic FSI type of application, and I am using operator splitting techniques to set-up my problem. I have been able to set the basics (boundary conditions, initial conditions etc.) all okay. However, I think I do not really understand the way MixedFunctionSpace in FEniCS handles updates to the functions. My solution variables do not change with time at all - they seem to be remaining at the same value. To be precise only Vp seems to be getting updated - the variables for U and P and Lp are not.

I present here a minimal example (which may still be relatively long, but kindly bear with me please).

I first set up the spaces as follows:

from dolfin import *

V   = VectorFunctionSpace(mesh, "CG", 2, dim=2)
P   = FunctionSpace(mesh, "CG", 1)
Up  = VectorFunctionSpace(mesh, "R", 0, dim=2)
Lp  = VectorFunctionSpace(mesh, "R", 0, dim=2)

# for operator splitting step 1
M1          = MixedFunctionSpace([V,P])
U1, P1      = TrialFunctions(M1)  
W1, Q1      = TestFunctions(M1)

Solution_M1         = Function(M1)
(U1_sol, P1_sol)    = Solution_M1.split(True)

# for operator splitting step 2
U2          = TrialFunction(V)  
W2          = TestFunction(V) 

U2_sol = Function(V)

# for operator splitting step 3
combinedSpace     = MixedFunctionSpace([V,P,Up,Lp])
U3, P3, Vp3, Lp3  = TrialFunctions(combinedSpace) 
W3, Q3, Up3, Mp3  = TestFunctions(combinedSpace)

Solution_M3       = Function(combinedSpace)
(U3_sol, P3_sol, Vp3_sol, Lp3_sol) = Solution_M3.split()

Following this, I set up the initial values:

# additionally we need an initial space for solution
U_old   = Function(V)
P_old   = Function(P)
Vp_old  = Function(Up)

initVel   = InitialVelocity()
initPres  = InitialPressure()

U_old.interpolate(initVel)
P_old.interpolate(initPres)

Vx_old = 0.0
Vy_old = 0.0

Vp_old.vector()[0] = Vx_old
Vp_old.vector()[1] = Vy_old

Then the form definitions and the time-loop for solution:

begin("Form compilation operator splitting step 1")

# ... code to define form for B1(U1,P1; W1,Q1; U_old) ...

a1 = lhs(B1)
b1 = rhs(B1)
problem1    = LinearVariationalProblem(a1, b1, Solution_M1, bcs=M1_bcD)
solver1     = LinearVariationalSolver(problem1)

end()

begin("Form compilation operator splitting step 2")

# ... code to define form for B2(U2; W2; U1_sol) ...

a2 = lhs(B2)
b2 = rhs(B2)
problem2    = LinearVariationalProblem(a2, b2, U2_sol, bcs=U2_bcD)
solver2     = LinearVariationalSolver(problem2)

end()

while tNow <= tStop:

   solver1.solve()

   solver2.solve()

   Vx_old = Vp_old.vector()[0]
   Vy_old = Vp_old.vector()[1]

   # ... code to change Vx_old, Vy_old to Vx_new, Vy_new ...

   Vp_old.vector()[0] = Vx_new
   Vp_old.vector()[1] = Vy_new

   # ... code to define form for B3(U3,Vp3,Lp3; W3,Up3,Mp3; U2_sol,Vp_old) ...

   a3 = lhs(B3)
   b3 = rhs(B3)
   problem3    = LinearVariationalProblem(a3, b3, Solution_M3, bcs=M3_bcD)
   solver3     = LinearVariationalSolver(problem3)

   assign(U_old, Solution_M3.sub(0))

   tNow += dt

Any help, tips, or suggestions will be very helpful.

Thanks.

Comments

Sorry but this minimum working example is a) quite complex and b) not complete. eg: It doesn't state what B1, B2, etc are, so the problem could be in how you pass your '_old' variables between them.

Have you tried a simple problem (ie: only one function / form) and building up your complex model piece by piece? That would help develop your understanding...

---

Thanks a lot for your response, and apologies if the example above was too complex. I have broken down the issue into a smaller sub-issue, and further worked on this to isolate a more specific example.

I will repost the question - please ignore the previous version of the comment which had the link to another question.