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.

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
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  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »