Basic example does nothing

variational formulation
1

Hi,

I’m totally new to FeniCS and FEM (not to PDE) and try to build up a more complex application which I reduced meanwhile to the attached code. Even if so basic I don’t get the example to work, i.e. I expect to get a linear temperature decrease from bound G to bound F.

"""
Time dependent Thermal Convection-Diffusion Equation

  rho*cp*grad(T,t)+rho*cp*w*grad(T,x)=div(lambda*grad(T,x),x)
                   T = Tc    on Gamma G 
                   T = Tc/2. on Gamma F

Temperatures in degree Celsius              

Boundaries:
  x: 0 --> Thick
  y: 0 --> LengthT (y-Axis inverted)

    x=0.0    boundary  G          x = Thick
    0 0 0 0 0 0 0 0 0 0 0 0 0 ... 0
 b  4                             1
 o  4                             1  boundary M
 u  4                             1
 n  ...                           ...
 d  4                             1
 a  4                             1  y = LengthM
 r  4                             2
 y  4                             2  boundary O
    ...                           ...
 S  4                             2
    4                             2
    3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3  y = LengthT
             boundary  F
"""

from fenics import *
import matplotlib.pyplot as plt
import sys

#=============================================================
# some problem parameters

rho   = 7000           # density             kg/m**3
lmbda = 50             # conductivity        W/m K  
cp    = 500            # haéat capacity      J/kg K
thd   = lmbda/(rho*cp) # thermal diffusivity m**2/s

TG    = 1550.0         # Dirichlet Boundary G degree Celsius
TF    = TG/2.0         # Dirichlet Boundary F degree Celsius

Tini  = 1000.0         # initial temperature

Thick    = 0.2         # thickness
CenterX  = Thick/2.0
LengthM  = 0.6         # length from 0.0 to of first boundary 
npThick  = 20
LengthT  = 3.0         # total length
npLength = 300

#=============================================================
# mesh, functionspace and boundaries

p0   = Point(0.0,0.0)
p1   = Point(Thick,LengthT)
mesh = RectangleMesh(p0,p1,npThick,npLength,'left')
V    = FunctionSpace(mesh, 'P', 1)

# ALL OTHER BOUNDARIES ARE IN THIS EXAMPLE REMOVED

class BoundaryG(SubDomain): # Boundary 0 == G
    tol = 1E-8
    nr  = 0
    def inside(self, x, on_boundary):        
        if on_boundary and near(x[1],0.0,self.tol): return True
        else:                                       return False

class BoundaryF(SubDomain): # Boundary 1 == F (original Boundary 3 - all others removed)
    tol = 1E-8
    nr  = 1
    def inside(self, x, on_boundary): 
        if  on_boundary and \
            near(x[1],LengthT,self.tol): return True
        else:                            return False

boundary_markers = MeshFunction('size_t', mesh,dim=mesh.topology().dim()-1) 

boundG = BoundaryG()
boundG.mark(boundary_markers,boundG.nr)

boundF = BoundaryF()
boundF.mark(boundary_markers,boundF.nr)

# for Neumann and Robin boundary conditions (later)
ds = Measure('ds',domain=mesh,subdomain_data=boundary_markers)

boundary_file = File("t_modell/boundaries.pvd")
boundary_file << boundary_markers

# check boundary assignments
# for x in mesh.coordinates():
#     if   bG.inside(x,True): print("%s is on G" % x)
#     elif bF.inside(x,True): print("%s is on F" % x)

#=============================================================
# boundary conditions

u  = TrialFunction(V)
v  = TestFunction(V)

# default bc is Neumann symmetry

bcond = {
         boundG.nr:{'Dirichlet':Constant(TG)},
         boundF.nr:{'Dirichlet':Constant(TF)}
        }

bcs = []
for b in bcond:
    bc = DirichletBC(V,bcond[b]['Dirichlet'],boundary_markers,b)
    bcs.append(bc)

#=============================================================
# Define variational problem

fDC = Constant(thd)  # thermal diffusivity

dt  = 1.0            # time step size - here constant
fDt = Constant(dt)

# Define initial value
fIni = Constant(Tini)       # Initial Values
u_n  = interpolate(fIni,V)  # old / initial values

#===========================================================
# PDE
#
#  rho*cp*grad(T,t)+rho*cp*w*grad(T,x)=div(lambda*grad(T,x),x)
#
#===========================================================
# weak problem formulation without convection

F    = u*v*dx + fDt*fDC*dot(grad(u),grad(v))*dx - u_n*v*dx
a, L = lhs(F), rhs(F)

print("a:",a)
print("L:",L)

#=============================================================
# solve time dependent problem

tFinal = 60.0  # final simulation time

# Create VTK file for saving solution
vtkfile = File('t_modell/solution.pvd')

# Time-stepping
u = Function(V)
t = 0
while t < tFinal:

    # Update current time
    t += dt
    print(t,"s /",tFinal,'s')

    # Compute solution
    solve(a == L, u, bcs)

    # Save to file and plot solution
    vtkfile << (u, t)
    plot(u)

    # Update previous solution
    u_n.assign(u)

# Hold plot
import matplotlib.pyplot as plt
plt.show()

Could anyone with some experience give me a hint? I think it has something to do with the definition of boundary conditions.

Thank’s in advance.

2

Ok, I searched again for some hours and found, that the book example missed the following line:

boundary_markers.set_all(9999)

Therefore the inner of the domain was also assigned to boundary 0. After correction all works well.