### AMR for a phase change problem

190

views

-1

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 :

The error I get :

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_)
# 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.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 \
+ 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 + 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 \
+ 1/Pr*dot(nabla_grad(theta_),nabla_grad(theta))*dx \
- dot(theta_*u_,nabla_grad(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)
solver = AdaptiveNonlinearVariationalSolver(problem, M)
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 Answer

1

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 (

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.

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 :

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_)
# 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.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 \
+ 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
# 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)
solver = AdaptiveNonlinearVariationalSolver(problem, M)
# 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
***
*** fenics-support@googlegroups.com
***
*** 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

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

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?

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

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.