AMR for a phase change problem

190
views
-1
3 months ago by
Hello, I am trying to solve a phase change problem. Because of the sharpness of my functions I want to implement an AMR in the process. However, it seems that it doesn't match very well with the non linear solver based on the Newton Method. How to correct this bug ? I precise it was well working without AMR.
Thank you in advance for helping me !

My code :

from __future__ import print_function
from fenics import *
from dolfin import *
import numpy as np

# Define constants and known data
Tcold_ = -0.01
Thot_ = 2.
Tcold = Constant(Tcold_)
Thot = Constant(Thot_)
Pr = Constant(1.0)  # Prandtl number
Ste = Constant(1.0) # Stefan number
Ra = Constant(10**6)  # Rayleigh number

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

# Create mesh and define function space
N = 20 # 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) # w.split() does not work here

# 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.1 ? 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)
(u_, p_, theta_) = split(w_)

def epsilon(u):

# Define the viscosity and the phase change curve
def phi(temp):
Tr = Constant(0.01)
r = Constant(0.025)
return(1+tanh((Tr-temp)/r)/2)

def mu(temp):
muL = Constant(1.0)
muS = Constant(10**8)
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 \
- div(v)*p_*dx - dot(fB,v)*dx

Eq2 = q*div(u_)*dx + Constant(10**(-7))*p_*q*dx

Eq3 = theta_*theta/dt*dx - theta_n*theta/dt*dx \
- 1/Ste*(phi(theta_)-phi(theta_n))*theta*dx \

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

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

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

# Define goal functional for AMR
M = phi(theta_)*dx

# Define error tolerance for AMR
tol = Constant(10**(-5))

# Create the Newton solver
problem = NonlinearVariationalProblem(F, w_, bc, J)
prm = solver.parameters
prm['newton_solver']['absolute_tolerance'] = 1E-8
prm['newton_solver']['relative_tolerance'] = 1E-7
prm['newton_solver']['relaxation_parameter'] = 0.8

# Save solution to file in VTK format
vtkfile1 = File('Sharp transition AMR/Velocity/velocity.pvd')
vtkfile2 = File('Sharp transition AMR/Pressure/pressure.pvd')
vtkfile3 = File('Sharp transition AMR/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(tol)
(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)​

The error I get :

Generating forms required for error control, this may take some time...
Traceback (most recent call last):
File "PhaseChange.py", line 134, in <module>
prm['newton_solver']['absolute_tolerance'] = 1E-8
File "/usr/lib/python3/dist-packages/dolfin/cpp/common.py", line 2453, in __getitem__
raise KeyError("'%s'"%key)
KeyError: "'newton_solver'"
Aborted (core dumped)
​
Community: FEniCS Project
1

What do you mean by "it seems that it doesn't match very well with the non linear solver based on the Newton Method"? What "bug" are you referring to? Have you tried using AMR for any other problem? Perhaps start with something more simple, and then turn this into a minimum working example (MWE) with which this community can help.

Also we should definitely compare notes if you're going to be spending much time on this. It looks like you're repeating the same process I went through over the past year, with the same model that I implemented using FEniCS. Your notation seems to almost exactly follow a combination of my proceedings paper and the 2014 paper by Danaila et al. How do I contact you privately?

written 3 months ago by Alexander G. Zimmerman
Indeed, I am doing an internship and my point is how to stabilize the FEM using continuous interior penalty stabilization method. S.Claus recommended me to look at your paper and the paper by Danaila in order to compute a phase change problem. Your paper helped me by giving me the system of equations that I implemented. I had a look at your codes too. The main difference with the one I computed is that it is written in object-oriented programming.

When I was saying bug I was refering to the error I just printed below in the comment.

I am not sure that your exact model will be the one my tutor wants to use in the end, but it is certain it is a good step to introduce me to FEniCS and to introduce me to the type of problems she wants me to solve.

You can easily contact me by mail : matthieu.diaz@ens-paris-saclay.fr
written 3 months ago by Matthieu Diaz
Thanks for clarifying the background.
written 12 weeks ago by Alexander G. Zimmerman
written 3 months ago by Alexander G. Zimmerman
I am sorry I didn't see that we could. That is done now.
written 3 months ago by Matthieu Diaz

I'll just repeat a small part of my other comment which might have been lost in all the rest of it:

Perhaps start with something more simple, and then turn this into a minimum working example (MWE) with which this community can help.

written 12 weeks ago by Alexander G. Zimmerman

1
12 weeks ago by
I realized that my paper and the Phaseflow python package might not be the optimal resources for learning how to do this for the first time.

So I made a series of Jupyter notebooks of increasing complexity, each using mixed elements ( fenics.MixedElement ) and goal-oriented adaptive mesh refinement ( fenics.AdaptiveNonlinearVariationalSolver).

This notebook solves the convection-coupled melting benchmark. There are some minor differences between this and what's shown in the script in your question, some of which are likely errors in your script. For example, your definition of $phi$ looks incorrect. Try plotting it in the temperature range and see if it's what you expect.

For each notebook, I also put a Python script ( .py ) with otherwise the same name in the same folder, where the comments and plotting have been stripped.
I took note from your tutorials and tried to used them on my code. The simulation perfectly works without AMR but I want to add it to avoid memory limitations.
Here is my code :

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

parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["optimize"] = True

# 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(10000.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 = 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((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)
(u_, p_, theta_) = split(w_)

def epsilon(u):

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

def mu(temp):
muL = Constant(1)
muS = Constant(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 \

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

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

# Provide the Jacobian form
J = derivative(F, w_, TrialFunction(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)

# 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('Phase Change/Concentration/concentration.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_.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_.leaf_node().vector()

# Update the progress bar
progress.update(t/T)​
And here is the error I get :

Summary of adaptive data:

Level  |  functional_value  error_estimate  tolerance  num_cells  num_dofs
--------------------------------------------------------------------------
0      |         -0.201431        0.000072   0.000010       5000     25604
1      |         -0.201362        0.000004   0.000010       5724     29236

Time spent for adaptive solve (in seconds):

Level  |  solve_primal  estimate_error  compute_indicators  mark_mesh  adapt_mesh    update
-------------------------------------------------------------------------------------------
0      |      5.929581        2.379780            2.659883   0.000755    0.005996  0.823443
1      |      3.652143        2.846782                   0          0           0         0

Traceback (most recent call last):
File "PhaseChangeAMR.py", line 156, in <module>
w_n.vector()[:] = w_.leaf_node().vector()
File "/usr/lib/python3/dist-packages/dolfin/cpp/la.py", line 1488, in __setitem__
self.__setslice__(0, len(self), values)
File "/usr/lib/python3/dist-packages/dolfin/cpp/la.py", line 1405, in __setslice__
self._assign(values)
RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to assign one vector to another.
*** Reason:  Vectors must be of the same length when assigning. Consider using the copy constructor instead.
*** Where:   This error was encountered inside PETScVector.cpp.
​
Could you help me with that ?
written 12 weeks ago by Matthieu Diaz
Try AMR first on a more simple problem, and if you're still having trouble, post a minimal working example (MWE) which shows the version working without AMR, and the version not working with AMR.
written 11 weeks ago by Alexander G. Zimmerman