Hello,
Not sure if this is the best place to ask this question, but I think I am
not understanding the "MixedFunctionSpace" concept in FEniCS well, and I
did not find any existing discussions on something like this in the Q & A
forum. I have posted this as a question also in Q&A (link:
http://fenicsproject.org/qa/6380/correctly-solutions-mixedfunctionspace-operator-splitting
).
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 Vpseems 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 helpful. If more details are
required, then I can provide additional code segments, and other details.
Thanks again, and apologies for a rather long email to the thread.
Warm regards.
Dr. Debanjan Mukherjee
_______________________________________________
fenics mailing list
[email protected]
http://fenicsproject.org/mailman/listinfo/fenics
