Hi all,
I am sure I am making a silly mistake with my time-stepped heat equation problem, but I cannot for the life of me find the error. I am a bit of newbie with using time stepping in Fenics.
I am solving the usual 3D heat equation in time within a cylinder of conductivity K with no internal sources
\frac{\partial T}{\partial t} - \nabla\cdot K\nabla T = 0,
where I have a constant flux source on the top (Neumann boundary) and a constant temperature on the bottom (Dirichlet boundary):
\begin{array}{ll}
K\frac{\partial T}{\partial n} = 10 & \textrm{on top} \\
T = 0 & \textrm{on bottom} \end{array}
The variational form for implicit Euler is
a = \int\limits _\Omega(\Delta t K \nabla v\cdot\nabla T^{n+1} + v\cdot T^{n+1} - v\cdot T^n)dV - \Delta t \int\limits_{\partial\Omega_{top}} v \cdot \Phi ds
where \Phi is a fixed flux value, v is the test function, \Delta t is the time increment and n is the time step number.
The code runs without error. The problem is that the solution is always zero. The solution does not advance with the time step.
from mshr import *
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
import os, sys, traceback
tol = 1.0e-8
K = Constant((1.0))
# time stepping params
dt = 0.1
NumSteps = 2
C1 = Cylinder(dolfin.Point(0,0,0), dolfin.Point(0, 0, 2), 1.0, 1.0)
mesh = generate_mesh(C1, 45)
class Free(SubDomain):
def inside(self, x, on_boundary):
return on_boundary
class FixedTemp(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[2], 0.0, tol)
class HeatFlux(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[2], 2.0, tol)
info(mesh)
# Output mesh for inspection
File("CylMesh.pvd").write(mesh)
# Mark boundaries
sub_domains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
sub_domains.set_all(3)
FREE = Free()
FREE.mark(sub_domains, 0)
FIXEDTEMP = FixedTemp()
FIXEDTEMP.mark(sub_domains, 1)
HEATFLUX = HeatFlux()
HEATFLUX.mark(sub_domains,2)
File("BoxSubDomains.pvd").write(sub_domains)
#Input heat flux
Qs = Expression(('10.0'), degree=0)
# Set up function spaces
cell = tetrahedron
ele = FiniteElement('Lagrange', cell, 1)
V = FunctionSpace(mesh, ele)
u = TrialFunction(V)
v = TestFunction(V)
un = Function(V)
u1 = Function(V)
#Initial condition
g = Expression('0.0', degree = 0)
pp0 = interpolate(g, V)
un.assign(pp0)
u1.assign(pp0)
#Boundary condition dictionary
boundary_conditions = {0: {'FREE' : 0.0},
1: {'FIXEDTEMP': 0.0},
2: {'HEATFLUX': Qs}}
#Build boundary conditions
bcs = []
integrals_source = []
# Fixed temp to zero on edges
for i in boundary_conditions:
if 'FIXEDTEMP' in boundary_conditions[i]:
bcs.append(DirichletBC(V, Constant((0.0)), sub_domains, i))
# Build constant power heat source on surface
for i in boundary_conditions:
if 'HEATFLUX' in boundary_conditions[i]:
Q = boundary_conditions[i]['HEATFLUX']
bb = v * Q * ds(i)
integrals_source.append(bb)
a = (dt * K * inner(grad(v), grad(u)) + v * u) * dx
L = dt * sum(integrals_source) + v * un * dx
vdim = u1.vector().size()
print ('Solution vector size = ',vdim)
#A, b = assemble_system(a, L, bcs)
#solver = PETScLUSolver("mumps")
#solver.parameters['symmetric'] = True
fpp = File('Temp.pvd')
# Do time stepping
t = 0
for nn in range(NumSteps):
t += dt
#Solve
#solver.solve(A, u1.vector(), b)
solve(a == L, u1, bcs, solver_parameters = {'linear_solver' : 'mumps'})
#Update solution
un.assign(u1)
fpp << (u1, t)
sys.exit(0)