Heat transfer model does not work correctly in case of multi-physics


239
views
0
3 months ago by
Hi everyone 

I am new to Fenics and I am trying to improve the example t04_heat_gaussian.py in the tutorial of Fenics a little bit by adding another physic, but the result seems not to be correct. The initial condition in this example gives a certain value of variable u in the beginning and it should decrease with time, see pictures. However, when I add another variable, namely p, which is belong to another function space, the variable u takes the initial value and keeps constantly with time.

If anyone has any idea, please let me know. 
Thank you very much

from __future__ import print_function
from fenics import *
import numpy as np
from pylab import show, triplot
from mshr import *
import time

#========Parameters==============

width = 2.0                 # width of Omega
height = 2.0                # height of Omega

T = 2.0
num_steps = 50              # Number of time steps
dt = T / num_steps          # Time step size
t = 0.0                     # Time variale

# Create domain 
Omega = Rectangle(Point(0.0, 0.0), Point(width, height))

# Create mesh
mesh = generate_mesh(Omega, 30)

# Define function space
Q = FiniteElement('Lagrange', triangle, 1) # p space
V = FiniteElement('Lagrange', triangle, 1) # u space
element = Q*V
W = FunctionSpace(mesh, element)

# Define boundary conditions
bc1 = DirichletBC(W.sub(1), Constant(0), 'on_boundary')
bc2 = DirichletBC(W.sub(0), Constant(0), 'on_boundary')
bcs = [bc1, bc2]


# Variational problem
unk = Function(W)
test = TestFunction(W)
(p, u) = split(unk)
(w, v) = split(test)

# Define initial value
initial = Expression(('0.0', 'exp(-a*pow(x[0]-1, 2) - a*pow(x[1]-1, 2))'), degree=2, a=5) 
unk0 = interpolate(initial, W)
(p0, u0) = unk0.split()

f = Constant(0)

# Dynamic variation form 
F1 = u*v*dx + dt*dot(grad(u), grad(v))*dx - (u0 + dt*f)*v*dx 
F2 = dot(grad(p), grad(w))*dx - u*w*dx
F =  F1 + F2

# Create VTK file for visualization output
vtkfile_p = File('postprocessing/solution_p.pvd')
vtkfile_u = File('postprocessing/solution_u.pvd') 

# Dynamic Solver
for n in range(num_steps):
	
	# Compute solution
	solve(F==0, unk, bcs)
	p_, u_ = unk.split()	
	
	vtkfile_p << (p_, t)
	vtkfile_u << (u_, t)
	
	# Update previous solution
	unk0.assign(unk)
	
	# Update current time
	t += dt
	
# Hold plot
interactive()
Community: FEniCS Project
2
Function.assign() goes against Pythonic explicit is better than implicit. Its problems are reported, see https://bitbucket.org/fenics-project/dolfin/issues/918.

Because you know that unk and unk0 live on the same space, it is natural to do explicit in-place assignment of vectors: unk0.vector()[:] = unk.vector().
written 3 months ago by Jan Blechta  

3 Answers


2
3 months ago by

The functions that you get by using (p0, u0) = unk0.split(True) are not related with unk0 after this step, and hence you have to update p0 and u0 instead of unk0.

Try something like this:

# Compute solution
solve(F==0, unk, bcs)
p_, u_ = unk.split(True)	

# Update previous solution
u0.assign(u_)
p0.assign(p_)
Hi Mella. Thank for your reply. I tried you suggestion and I got an error saying the function assign in /dolfin/functions/function.py had a problem with self._assign(rhs). 
written 3 months ago by Duy Duc NGUYEN  
I do not have this error with direct substitution of Hernán's edits.
written 3 months ago by pf4d  
1
After restarting my terminal, and applying your idea again, it works correctly now. Thank all. 
written 3 months ago by Duy Duc NGUYEN  
2
3 months ago by
pf4d  
While the answer by Hernán Mella is correct, there is no need to make deep copies.

Simply,

from __future__ import print_function
from fenics import *
import numpy as np
from pylab import show, triplot
from mshr import *
import time

#========Parameters==============

width     = 2.0             # width of Omega
height    = 2.0             # height of Omega
T         = 2.0             # end time
num_steps = 50              # Number of time steps
dt        = T / num_steps   # Time step size
t         = 0.0             # Time variale

# Create domain 
Omega = Rectangle(Point(0.0, 0.0), Point(width, height))

# Create mesh
mesh  = generate_mesh(Omega, 30)

# Define function space
Q       = FiniteElement('Lagrange', triangle, 1) # p space
V       = FiniteElement('Lagrange', triangle, 1) # u space
element = Q*V
W       = FunctionSpace(mesh, element)

# Define boundary conditions
bc1     = DirichletBC(W.sub(1), Constant(0), 'on_boundary')
bc2     = DirichletBC(W.sub(0), Constant(0), 'on_boundary')
bcs     = [bc1, bc2]


# Variational problem
unk     = Function(W)
test    = TestFunction(W)
(p, u)  = split(unk)
(w, v)  = split(test)

# Define initial value
initial = Expression(('0.0', 'exp(-a*pow(x[0]-1, 2) - a*pow(x[1]-1, 2))'),
                     degree=2, a=5)
unk0    = interpolate(initial, W)
p0, u0  = unk0 # NOTE: the only change

f       = Constant(0)

# Dynamic variation form 
F1 = u*v*dx + dt*dot(grad(u), grad(v))*dx - (u0 + dt*f)*v*dx 
F2 = dot(grad(p), grad(w))*dx - u*w*dx
F  =  F1 + F2

# Create VTK file for visualization output
vtkfile_p = File('postprocessing/solution_p.pvd')
vtkfile_u = File('postprocessing/solution_u.pvd')

# Dynamic Solver
for n in range(num_steps):

  # Compute solution
  solve(F==0, unk, bcs)
  p_, u_ = unk.split()     
  
  vtkfile_p << (p_, t)
  vtkfile_u << (u_, t)
  
  # Update previous solution
  unk0.assign(unk)
  
  # Update current time
  t += dt
​
2
Just for the sake of clarity, the above notation is equivalent to
p0, u0 = split(unk0)​
written 3 months ago by Adam Janecka  
1
3 months ago by
Leaving your code as it is and just updating the previous solution using slicing, i.e.,
unk0.vector()[:] = unk.vector()​

also works, even though I'm not sure why.

1
The following example demonstrates that the problem is related with the call of assign method. It seems that it somehow breaks the link between w0 and (u0, p0).
Any idea why it happens?

from dolfin import *

mesh = UnitSquareMesh(1, 1)
Q = FiniteElement('CG', triangle, 1) # p space
V = FiniteElement('CG', triangle, 1) # u space

W = FunctionSpace(mesh, MixedElement([Q, V]))
w, w0 = Function(W), Function(W)

# Create expression for easy setup of components in w, w0
expr = Expression(('p', 'u'), degree=2, p=0.0, u=1.0)

# Create functions for viewing components of w0
_p = Function(FunctionSpace(mesh, Q))
_u = Function(FunctionSpace(mesh, V))

# Initialize w0
w0.interpolate(expr)
p0, u0 = w0.split(deepcopy=False) # False is default
_p.interpolate(p0)
_u.interpolate(u0)

# Print w0 and its components
info(w0.vector(), True)
info(_p.vector(), True) # p0 = 0
info(_u.vector(), True) # u0 = 1

# Initialize w with some other values
expr.p = 2.0
expr.u = 3.0
w.interpolate(expr)

# Update w0 with values from w (p0, u0 should be updated as well)
w0.assign(w) # THIS CALL seems to "break a link" between w0 and (p0, u0)
#w0.vector()[:] = w.vector() # works fine if we comment out w0.assign(w)
#w0.interpolate(w)           # works fine if we comment out w0.assign(w)

# View w0 and its components again
_p.interpolate(p0)
_u.interpolate(u0)
info(w0.vector(), True)
info(_p.vector(), True) # should be p0 = 2
info(_u.vector(), True) # should be u0 = 3​
written 3 months ago by Martin Řehoř  
1
Look what happens when w0.assign(w) is called (the following snippet is taken from dolfin/function/Function.cpp):

// Make a copy of all the data, or if v is a sub-function, then we
// collapse the dof map and copy only the relevant entries from the
// vector of v.
if (v._vector->size() == v._function_space->dim())
{
  // Copy function space
  _function_space = v._function_space;

  // Copy vector
  _vector = v._vector->copy();

  // Clear subfunction cache
  _sub_functions.clear();
}
else
...  ​

The problem is probably caused by calling  _sub_functions.clear() which breaks the previously mentioned link between w0 and (p0, u0). Nothing like this happens when we use w0.vector()[:] = w.vector().

written 3 months ago by Martin Řehoř  
This is what happens when you call assign, I believe.
written 3 months ago by pf4d  
This is what happens when you call assign, I believe.
written 3 months ago by pf4d  
1
`Function.assign()` does too much magic, see https://bitbucket.org/fenics-project/dolfin/issues/918. `unk0.vector()[:] = unk.vector()` essentially happens in this case although not that efficiently - one temporary vector is created in fact.
written 3 months ago by Jan Blechta  
That's what I thought. But why does it give different result?
written 3 months ago by Adam Janecka  
It certainly should not....  You want a [:] at the end of that though.
written 3 months ago by pf4d  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »