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

417

views

2

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

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

### 3 Answers

2

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
9 months ago by
Duy Duc NGUYEN

I do not have this error with direct substitution of Hernán's edits.

written
9 months ago by
pf4d

1

After restarting my terminal, and applying your idea again, it works correctly now. Thank all.

written
9 months ago by
Duy Duc NGUYEN

2

While the answer by Hernán Mella is correct, there is no need to make deep copies.

Simply,

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
9 months ago by
Adam Janecka

1

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

Any idea why it happens?

`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
9 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
9 months ago by
Martin Řehoř

This is what happens when you call assign, I believe.

written
9 months ago by
pf4d

This is what happens when you call assign, I believe.

written
9 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
9 months ago by
Jan Blechta

That's what I thought. But why does it give different result?

written
9 months ago by
Adam Janecka

It certainly should not.... You want a [:] at the end of that though.

written
9 months ago by
pf4d

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

`Function.assign()`

goes against Pythonicexplicit 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()`

.