### AMR + time dependent problem issue with assignment

198
views
1
11 weeks ago by
Hello,
I am solving under FEniCS a time dependent non linear variational problem of phase change. I am using mixed finite elements to do it.
I solve at each time step the problem. It provides me the velocity field, the temperature and the pressure in the domain I defined.
The code is running but it does not actualize the problem at each iteration. I mean I always observe the first iteration because the solution at time step n is not well assigned.
I copied below the essential parts of my code. I hope someone will be able to help me.

from __future__ import print_function
from fenics import *
from dolfin import *

# 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(100000.0)  # Rayleigh number

# Time discretization
t = 0.0
T = 0.5 # final time
dt = 0.1
num_steps = int((T-t)/dt)   # number of time steps

# Create mesh and define function space
N = 10 # 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)

# 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((0.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)

# Create and assign the initial conditions
# theta_init = 'Th+(Tc-Th)*x[0]'.replace('Tc', str(Tcold_)).replace('Th', str(Thot_))
theta_init = 'x[0]<=0.0 ? Th : Tc'.replace('Tc', str(Tcold_)).replace('Th', str(Thot_))
# u_init and p_init are null
w_n = interpolate(Expression(("0.", "0.", "0.", theta_init), element=element), W)
u_n, p_n, theta_n = w_n.split()

# Create the current solution
w_ = Function(W, name="Solution")
(u_, p_, theta_) = split(w_)

def epsilon(u):

# Define the viscosity and the phase change curve
def phi(temp):
Tr = 0.01
r = 0.05
return(0.5+0.5*tanh((Tr-temp)/r))

def mu(temp):
muL = 1
muS = 100000
return(muL + (muS - muL)*phi(temp))

# Define Boussinesq coupling term
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 \
+ 2*mu(theta_)*inner(epsilon(u_), epsilon(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 \

gamma = Constant(1E-7)
Stab_pressure = gamma*p_*q*dx

# Sum of all contributions
F = Eq1 + Eq2 + Eq3 + Stab_pressure

# Create the progress bar
progress = Progress('Time-stepping')
set_log_level(PROGRESS)

# Provide the Jacobian form
J = derivative(F, w_)

# AMR parameters
M = ( phi(theta_)-phi(theta_n) )*dx() # dphi/dt goal function
epsilon_M = 1.e-5

# Create the Newton solver
problem = NonlinearVariationalProblem(F, w_, bc, J)

prm = solver.parameters
prm["nonlinear_variational_solver"]["newton_solver"]['absolute_tolerance'] = 1E-7
prm["nonlinear_variational_solver"]["newton_solver"]['relative_tolerance'] = 1E-7
prm["nonlinear_variational_solver"]["newton_solver"]['relaxation_parameter'] = 1.0

# Save solution to file in VTK format
vtkfile1 = File('AMR/Velocity/velocity.pvd')
vtkfile2 = File('AMR/Pressure/pressure.pvd')
vtkfile3 = File('AMR/Temperature/temperature.pvd')
vtkfile4 = File('AMR/Mesh/mesh.pvd')

# Form and solve the linear system
for n in range(num_steps):

# Update current time
t += dt

# Solve the variational problem
solver.solve(epsilon_M)
solver.summary()

(u_, p_, theta_) = w_.leaf_node().split()

# Save data
vtkfile1 << (u_, t)
vtkfile2 << (p_, t)
vtkfile3 << (theta_, t)

leaf_mesh = w_.leaf_node().function_space().mesh()
vtkfile4 << leaf_mesh

# Update the previous solution
w_n.assign(w_.leaf_node()) ###### Issue here I think

# Update the progress bar
progress.update(t/T)​
Community: FEniCS Project
1
Can you reproduce this issue with something more simple, e.g. the heat equation in 1D?

Otherwise maybe try finding the issue by comparing to this tutorial: https://github.com/geo-fluid-dynamics/phaseflow-fenics/blob/master/tutorials/FEniCS/03-ConvectionCoupledMelting-MixedElement-AMR.ipynb
written 11 weeks ago by Alexander G. Zimmerman
Amazing tutorial, thank you so much for the contribution!
written 11 weeks ago by Miguel
Thanks and you're welcome. Jupyter notebooks seem to be quite nice for this stuff.
written 11 weeks ago by Alexander G. Zimmerman
What is surprising to me is that I wrote
print('Infos')
info(w_n.leaf_node().vector())
info(w_.leaf_node().vector())

w_n.leaf_node.vector()[:] = w_.leaf_node.vector()​

The info function returns that each one of the vectors is a <PETScVector of size ...> but when the update returns the error 'Function' object has no attribute 'vector'.

written 10 weeks ago by Matthieu Diaz
Problem solved :

2 main changes :

u_n, p_n, theta_n = split(w_n)

...

w_n.leaf_node().assign(w_.leaf_node())​
written 10 weeks ago by Matthieu Diaz
Note that you changed a couple of things here. With the other approach which does not use assign, you should have written

w_n.leaf_node().vector()[:] = w_.leaf_node().vector()​

, i.e. you must call Function.leaf_node. Anyway, I think it's best to use assign if that actually works.
written 10 weeks ago by Alexander G. Zimmerman