### AMR + time dependent problem issue with assignment

198

views

1

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.

Thank you in advance.

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.

Thank you in advance.

```
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_)
# Define symmetric gradient
def epsilon(u):
return sym(nabla_grad(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 \
+ dot(dot(grad(u_), u_), v)*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 \
+ 1/Pr*dot(nabla_grad(theta_), nabla_grad(theta))*dx \
- dot(theta_*u_, nabla_grad(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)
solver = AdaptiveNonlinearVariationalSolver(problem, M)
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)
# Adapted mesh refinement saved
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

Please login to add an answer/comment or follow this question.

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

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'.

2 main changes :

`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.