Heat-conduction with neumann-boundaries only

Question

tl;dr:
fenics keeps telling me "An Integral without a Domain is now illegal" and I dont know why.
Long version:
Hi guys, i want to solve the heat conduction for a cube with neumann-condition on the boundaries (convection and radiation on top and bottom and source of heat on bottom which might be constant in this example). Therefore I installed FEniCS on my MacBook and wrote the following code:

mesh = UnitCubeMesh(10,10,10)
xmax, ymax, zmax = 1, 1, 1
V = FunctionSpace(mesh, 'Lagrange', 1)
boundary_parts = FacetFunction("size_t", mesh)

class Bottom(SubDomain):
   def inside(self, x, on_boundary):
      tol = 1E-14
      return on_boundary and abs(x[2]) < tol
bottom = Bottom()
bottom.mark(boundary_parts, 0)

class Top(SubDomain):
   def inside(self, x, on_boundary):
      tol = 1E-14
      return on_boundary and abs(x[2] - zmax) < tol
top = Top()
top.mark(boundary_parts, 0)

Classes Left, Right, Front and Back analog.
I continue as described in http://fenicsproject.org/documentation/tutorial/timedep.html and and http://fenicsproject.org/documentation/tutorial/materials.html#multiple-neumann-robin-and-dirichlet-condition (both examples work with my version of FEniCS, so I guess i installed it correctly):

T = TrialFunction(V)
v = TestFunction(V)
T0 = Constant("293")
T_1= interpolate(T0, V)
del_t=1
t_ges=10
source=Constant("5")
a = (T*v + inner(nabla_grad(T), nabla_grad(v)))*dx - del_t*(T-T0)*ds(0) - del_t*(T-T0)*ds(1)
L= T_1*v*dx - del_t*source*ds(0)

Now I need to combine both tutorial chapters:

A = assemble(a, exterior_facet_domains = boundary_parts)
b = None
T = Function(V)
t = del_t
while t <= t_ges:
   b = assemble(L, tensor=b, exterior_face_domains=boundary_parts)
   solve(A, T.vector(), b, 'lu')
   t += dt
   plot(T)
   T_1.assign(T)
interactive()

But if I try to run the program now, it tells me Integrals without a Domain would be illegal now - the domains were declared above.
The error seems to be in the line where i declared L = T_1*v*dx -del_t*source*dx; when i delete the part with the source the error is replaced with another error message:
Now it seems to find an error in the line where i wrote A=assemble(...) and says (after a few lines of directories): "ufl.algorithms.check_arities.ArityMismatch: Multiplying expressions with overlapping form argument number 1, argument is v_1"

Please help me, I am using FEniCS for the first time and cant find my mistake in the code above.

Comments

The domains are the classes Bottom, Top etc., arent they?

Answer 1

Hi,

When I try to run the provided code I get the error "All terms in form must have same rank" which means that your weak formula is not well-posed. The generic weak formulation is a(u,v)=L(v) and it seems you forgot to multiply by the test function v, here's the fix:

a = (T*v + inner(nabla_grad(T), nabla_grad(v)))*dx() - del_t*(T-T0)*v*ds(0) - del_t*(T-T0)*v*ds(1)
L = T_1*v*dx() - del_t*source*v*ds(0)

If you have trouble defining the formula you can just add all the terms in a single functionnal then call rhs and lhs (right hand side/left hand side) to sort it:

Fa = (T*v + inner(nabla_grad(T), nabla_grad(v)))*dx() - del_t*(T-T0)*ds(0) - del_t*(T-T0)*ds(1)
FL= T_1*v*dx() - del_t*source*ds(0)
F=Fa+FL
a=lhs(F)
L=rhs(F)

Which automatically builds the form. Then the code runs fine but you didn't define dt in the while loop so it stops.

If you define dt it will run. Be careful you marked both subdomains as 0:

bottom.mark(boundary_parts, 0)
top.mark(boundary_parts, 0)

in order to distinguish the top of the box, the bottom of the box and the interior you should do:

boundary_parts = FacetFunction("size_t", mesh)
boundary_parts.set_all(0) #marks whole cube as domain 0

class Bottom(SubDomain):
   def inside(self, x, on_boundary):
      tol = 1E-14
      return on_boundary and abs(x[2]) < tol
bottom = Bottom()
bottom.mark(boundary_parts, 1) #marks bottom of the cube as 1

class Top(SubDomain):
   def inside(self, x, on_boundary):
      tol = 1E-14
      return on_boundary and abs(x[2] - zmax) < tol
top = Top()
top.mark(boundary_parts, 2) #marks the top of the cube as 2

then you can define a measure on each of the subdomains:

dx=Measure('dx')[boundary_parts]
ds=Measure('ds')[boundary_parts]

With all this, you can get a running code, still, I'm not much into heat conduction and I don't know if the output is what you were looking for, anyway here's the full code, hoping this could help you:

from dolfin import *
import numpy as np

mesh = UnitCubeMesh(10,10,10)
xmax, ymax, zmax = 1, 1, 1
V = FunctionSpace(mesh, 'Lagrange', 1)
boundary_parts = FacetFunction("size_t", mesh)
boundary_parts.set_all(0) #marks whole cube as domain 0

class Bottom(SubDomain):
   def inside(self, x, on_boundary):
      tol = 1E-14
      return on_boundary and abs(x[2]) < tol
bottom = Bottom()
bottom.mark(boundary_parts, 1) #marks bottom of the cube as 1

class Top(SubDomain):
   def inside(self, x, on_boundary):
      tol = 1E-14
      return on_boundary and abs(x[2] - zmax) < tol
top = Top()
top.mark(boundary_parts, 2) #marks the top of the cube as 2

dx=Measure('dx')[boundary_parts]
ds=Measure('ds')[boundary_parts]

plot(boundary_parts)

T = TrialFunction(V)
v = TestFunction(V)
T0 = Constant("293")
T_1= interpolate(T0, V)
del_t=1
t_ges=10
source=Constant("5")
Fa = (T*v + inner(nabla_grad(T), nabla_grad(v)))*dx() - del_t*(T-T0)*v*ds(1) - del_t*(T-T0)*v*ds(2)
FL= T_1*v*dx() - del_t*source*v*ds(1)

F=Fa+FL

a=lhs(F)
L=rhs(F)

A = assemble(a, exterior_facet_domains = boundary_parts)
b = None
T = Function(V)
t = del_t
dt = 0.5
while t <= t_ges:
   b = assemble(L, tensor=b, exterior_facet_domains=boundary_parts)
   solve(A, T.vector(), b, 'lu')
   t += dt
   T_1.assign(T)

plot(T_1)

Comments

Thank you so much! I was searching everywhere and couldnt find this mistake.

---

Your solution is quite good, but you made a little mistake when you wrote:

F = Fa + FL

It was correct to write:

F = Fa - FL

Otherwise the results were very weird, because a or L would have the wrong sign.

---

Well, indeed it looked very odd because of this mistake... Thank you for the update!

---

Hello Maxmeier,

I am solving a heat equation for the first time and was looking at the code above.

Here, you say you are implying convection at top and radiation at bottom in your bilinear form through this:

- del_t*(T-T0)*v*ds(1) - del_t*(T-T0)*v*ds(2)

Also, your initial temperature of cube is 293 while the source of heat at bottom is 5.

Could you pls. spare a few minutes to explain the extended variational equation (as above) and physical meaning behind what is happening here. It would be very helpful for me.

What does T-T0 mean? Your heat source is less than your initial temperature and what are you implying by radiation at bottom and convection at top. Also, I thought that convection required 'movement of a medium', so what does convection at top mean?

---

I'm kind of in a hurry right now, so I will answer shortly:
The source is not a temperature with value 5, its power with for example 5 Watt. You will need to add some constants like density, thermal conductivity and capacity. T-T0 means convection to the surrounding air. The air is moving, thats why its called convection. In the equation above there is no radiation, I added it later because the solution algorithm becomes a little bit more tricky with nonlinear boundary conditions.