### Assignment with mixed finite elements

I am writting a code in order to simulate a phase flow problem. To do it, I am using mexed finite elements. As the equations are time dependent, I need to assign the previous step at each iteration. I use the function assign but it is not working. Indeed, my functions of the previous step (noted with _n) are always equal to the initial fields. I hope that someone will be able to tell me how to manage this assignment issue. Thank you in advance.

My code :

```
from __future__ import print_function
from fenics import *
from dolfin import *
import numpy as np
# Define constants and known data
Tcold = Constant(-1.0)
Thot = Constant(1.0)
Pr = Constant(1.0) # Prandtl number
Ste = Constant(1.0) # Stefan number
Ra = Constant(1.0) # Rayleigh number
# Time discretization
t = 0.0
T = 1.0 # final time
num_steps = 100 # number of time steps
dt = (T-t) / num_steps # time step size
# Create mesh and define function space
N = 50 # Characteristic number of elements
mesh = RectangleMesh(Point(0,0),Point(1,1),N,N)
# Function space
V = VectorElement('P', mesh.ufl_cell(), 2) # velocity
Q = FiniteElement('P', mesh.ufl_cell(), 1) # pressure
element = MixedElement(V, Q, Q) # velocity, pressure, temperature in this order
W = FunctionSpace(mesh, element)
# Mesh related functions
norm = FacetNormal(mesh)
h = CellDiameter(mesh)
# Define boundaries
lid = 'near(x[1], 1.0)'
noslip = 'near(x[1], 0.0) || near(x[0], 0.0) || near(x[0], 1.0)'
coldwall = 'near(x[0], 1.0)'
hotwall = 'near(x[0], 0.0)'
# Boundary conditions
bc_noslip = DirichletBC(W.sub(0), Constant((0.0, 0.0)), noslip)
bc_lid = DirichletBC(W.sub(0), Constant((1.0, 0.0)), lid)
bc_hot = DirichletBC(W.sub(2), Thot, hotwall)
bc_cold = DirichletBC(W.sub(2), Tcold, coldwall)
bc = [bc_noslip, bc_lid, bc_hot, bc_cold]
# Variational formulation
w = TestFunction(W)
(v, q, theta) = split(w) # w.split() does not work here
# Define the initial conditions and previous time step values
# The previous step's solution
w_n = Function(W)
# Split the mixed function
(u_n, p_n, theta_n) = w_n.split()
# Create and assign the initial conditions
u_init = Constant((0.0,0.0))
p_init = Constant(0.0)
# T_init = Expression('Th+(Tc-Th)*x[0]', degree = 2, Tc = Tcold, Th = Thot)
T_init = Expression('x[0]<=0 ? Th : Tc', degree = 2, Tc = Tcold, Th = Thot)
u_n = interpolate(u_init, W.sub(0).collapse())
p_n = interpolate(p_init, W.sub(1).collapse())
theta_n = interpolate(T_init, W.sub(2).collapse())
# Create the current solution
w_ = Function(W)
(u_, p_, theta_) = split(w_)
# Define symmetric gradient
def epsilon(u):
return sym(nabla_grad(u))
# Define the viscosity and the phase change curve
def phi(temp):
Tr = Constant(0.01)
r = Constant(0.01)
return(1+tanh((Tr-temp)/r)/2)
def mu(temp):
muL = Constant(1.0)
muS = Constant(1.0)
return (muL + (muS - muL)*phi(temp))
# Define Boussinesq coupling term
fB = Constant((0.0, 0.0))
# fB = Ra/Pr*theta_*Constant((0.0,1.0))
# Weak formulation F == 0
# Equations terms
Eq1 = dot(u_,v)/dt*dx - dot(u_n,v)/dt*dx \
+ dot(dot(u_, nabla_grad(u_)), v)*dx \
+ 2*inner(mu(theta_)*epsilon(u_), grad(v))*dx \
- div(v)*p_*dx - dot(fB,v)*dx
Eq2 = q*div(u_)*dx
Eq3 = theta_*theta/dt*dx - theta_n*theta/dt*dx \
- 1/Ste*(phi(theta_)-phi(theta_n))*theta*dx \
+ 1/Pr*dot(nabla_grad(theta_),nabla_grad(theta))*dx \
- dot(theta_*u_,nabla_grad(theta))*dx
# Stabilization terms
jump_u = grad(u_)('+')-grad(u_)('-')
jump_v = grad(v)('+')-grad(v)('-')
jump_p = nabla_grad(p_)('+')-nabla_grad(p_)('-')
jump_q = nabla_grad(q)('+')-nabla_grad(q)('-')
Stab1 = avg(h)*avg(h)*(dot(u_('+'),norm('+'))**2)*inner(jump_u,jump_v)*dS
Stab2 = Constant(0.1)*avg(h)*avg(h)*inner(jump_u,jump_v)*dS
Stab3 = avg(h)*avg(h)*dot(jump_p,jump_q)*dS
# Sum of all contributions
F = Eq1 + Eq2 + Eq3 #+ Stab1 + Stab2 + Stab3
# Create the progress bar
progress = Progress('Time-stepping')
set_log_level(PROGRESS)
# Provide the Jacobian form
J = derivative(F,w_)
# Create the Newton solver
problem = NonlinearVariationalProblem(F, w_, bc, J)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters.newton_solver
prm.error_on_nonconvergence = False
# Save solution to file in VTK format
vtkfile1 = File('Phase change/Velocity/velocity.pvd')
vtkfile2 = File('Phase change/Pressure/pressure.pvd')
vtkfile3 = File('Phase change/Temperature/temperature.pvd')
# Form and solve the linear system
for n in range(num_steps):
# Update current time
t += dt
# Solve the variational problem
solver.solve() # for use with the NonlinearVariationalProblem structure
(u_, p_, theta_) = w_.split() # split(w_) doesn't work
# Save data
vtkfile1 << (u_, t)
vtkfile2 << (p_, t)
vtkfile3 << (theta_, t)
# Update the previous solution
w_n.assign(w_)
# Update the progress bar
progress.update(t/T)
```

### 1 Answer

The problem is coming from a Python mistake that I was also making frequently when first using FEniCS.

When you do call `interpolate`

for `u_n`

for example, this returns a new object `u_n`

which replaces the previous one (which was previously returned by `w_n.split()`

) in the namespace. Now later when you update `w_n`

, the updates values aren't actually associated with your variational form.

I ran your example and saw the behavior you explained. I made some changes which perhaps aren't the minimal required (so you can change it differently, and I would indeed recommend making one change at a time and seeing what happens); but in this version the initial values are set in a way where you keep the correct references. When I run this, the initial values do appear to be updating appropriately:

```
from __future__ import print_function
from fenics import *
from dolfin import *
import numpy as np
# Define constants and known data
Tcold_ = -1.
Thot_ = 1.
Tcold = Constant(Tcold_)
Thot = Constant(Thot_)
Pr = Constant(1.0) # Prandtl number
Ste = Constant(1.0) # Stefan number
Ra = Constant(1.0) # Rayleigh number
# Time discretization
t = 0.0
T = 0.03 # final time
num_steps = 3 # number of time steps
dt = (T-t) / num_steps # time step size
# Create mesh and define function space
N = 50 # Characteristic number of elements
mesh = RectangleMesh(Point(0,0),Point(1,1),N,N)
# Function space
V = VectorElement('P', mesh.ufl_cell(), 2) # velocity
Q = FiniteElement('P', mesh.ufl_cell(), 1) # pressure
element = MixedElement(V, Q, Q) # velocity, pressure, temperature in this order
W = FunctionSpace(mesh, element)
# Mesh related functions
norm = FacetNormal(mesh)
h = CellDiameter(mesh)
# Define boundaries
lid = 'near(x[1], 1.0)'
noslip = 'near(x[1], 0.0) || near(x[0], 0.0) || near(x[0], 1.0)'
coldwall = 'near(x[0], 1.0)'
hotwall = 'near(x[0], 0.0)'
# Boundary conditions
bc_noslip = DirichletBC(W.sub(0), Constant((0.0, 0.0)), noslip)
bc_lid = DirichletBC(W.sub(0), Constant((1.0, 0.0)), lid)
bc_hot = DirichletBC(W.sub(2), Thot, hotwall)
bc_cold = DirichletBC(W.sub(2), Tcold, coldwall)
bc = [bc_noslip, bc_lid, bc_hot, bc_cold]
# Variational formulation
w = TestFunction(W)
(v, q, theta) = split(w) # w.split() does not work here
# Create and assign the initial conditions
T_init = 'x[0]<=0 ? Th : Tc'.replace('Tc', str(Tcold_)).replace('Th', str(Thot_))
w_n = interpolate(Expression(("0.", "0.", "0.", T_init), element=element), W)
u_n, p_n, theta_n = w_n.split()
# Create the current solution
w_ = Function(W)
(u_, p_, theta_) = split(w_)
# Define symmetric gradient
def epsilon(u):
return sym(nabla_grad(u))
# Define the viscosity and the phase change curve
def phi(temp):
Tr = Constant(0.01)
r = Constant(0.01)
return(1+tanh((Tr-temp)/r)/2)
def mu(temp):
muL = Constant(1.0)
muS = Constant(1.0)
return (muL + (muS - muL)*phi(temp))
# Define Boussinesq coupling term
fB = Constant((0.0, 0.0))
# fB = Ra/Pr*theta_*Constant((0.0,1.0))
# Weak formulation F == 0
# Equations terms
Eq1 = dot(u_,v)/dt*dx - dot(u_n,v)/dt*dx \
+ dot(dot(u_, nabla_grad(u_)), v)*dx \
+ 2*inner(mu(theta_)*epsilon(u_), grad(v))*dx \
- div(v)*p_*dx - dot(fB,v)*dx
Eq2 = q*div(u_)*dx
Eq3 = theta_*theta/dt*dx - theta_n*theta/dt*dx \
- 1/Ste*(phi(theta_)-phi(theta_n))*theta*dx \
+ 1/Pr*dot(nabla_grad(theta_),nabla_grad(theta))*dx \
- dot(theta_*u_,nabla_grad(theta))*dx
# Stabilization terms
jump_u = grad(u_)('+')-grad(u_)('-')
jump_v = grad(v)('+')-grad(v)('-')
jump_p = nabla_grad(p_)('+')-nabla_grad(p_)('-')
jump_q = nabla_grad(q)('+')-nabla_grad(q)('-')
Stab1 = avg(h)*avg(h)*(dot(u_('+'),norm('+'))**2)*inner(jump_u,jump_v)*dS
Stab2 = Constant(0.1)*avg(h)*avg(h)*inner(jump_u,jump_v)*dS
Stab3 = avg(h)*avg(h)*dot(jump_p,jump_q)*dS
# Sum of all contributions
F = Eq1 + Eq2 + Eq3 #+ Stab1 + Stab2 + Stab3
# Create the progress bar
progress = Progress('Time-stepping')
set_log_level(PROGRESS)
# Provide the Jacobian form
J = derivative(F,w_)
# Create the Newton solver
problem = NonlinearVariationalProblem(F, w_, bc, J)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters.newton_solver
prm.error_on_nonconvergence = False
# Save solution to file in VTK format
vtkfile1 = File('Phase change/Velocity/velocity.pvd')
vtkfile2 = File('Phase change/Pressure/pressure.pvd')
vtkfile3 = File('Phase change/Temperature/temperature.pvd')
# Form and solve the linear system
for n in range(num_steps):
# Update current time
t += dt
# Solve the variational problem
solver.solve() # for use with the NonlinearVariationalProblem structure
(u_, p_, theta_) = w_.split() # split(w_) doesn't work
# Save data
vtkfile1 << (u_, t)
vtkfile2 << (p_, t)
vtkfile3 << (theta_, t)
# Update the previous solution
w_n.vector()[:] = w_.vector()
# Update the progress bar
progress.update(t/T)
```

Note that I reduced `T`

and `num_steps`

because, after several, the Newton solver beings to hang and the loop just keeps going.

I'll also take this opportunity to shamelessly plug [Phaseflow](https://github.com/geo-fluid-dynamics/phaseflow-fenics), since it looks like you're trying to solve a **very** similar problem. We even seem to use mostly the same notation, so it is hopefully easy to read. At a glance it seems that the difference is that you've added some stabilization terms (those these are commented out, i.e. it looks like you're not actually adding them to the variational form yet).

This approach has some robustness issues that I discuss in my proceedings paper. I spent the better part of the last year figuring out how to implement this approach in FEniCS. My lessons learned are mostly captured in that paper and in the the Phaseflow repository; so I suggest reviewing these materials.

Also you should ask a new question about why the Newton solver hangs, since perhaps other people could contribute to that discussion and it's really a separate question. Please do accept my answer if it solved your problem.