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

545
views
2
12 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
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 12 months ago by Jan Blechta

2
12 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_)

written 12 months ago by Duy Duc NGUYEN
I do not have this error with direct substitution of Hernán's edits.
written 12 months ago by pf4d
1
After restarting my terminal, and applying your idea again, it works correctly now. Thank all.
written 12 months ago by Duy Duc NGUYEN
2
12 months ago by
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
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 12 months ago by Adam Janecka
1
12 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 12 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 12 months ago by Martin Řehoř
This is what happens when you call assign, I believe.
written 12 months ago by pf4d
This is what happens when you call assign, I believe.
written 12 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 12 months ago by Jan Blechta
That's what I thought. But why does it give different result?
written 12 months ago by Adam Janecka
It certainly should not....  You want a [:] at the end of that though.
written 12 months ago by pf4d