Heat equation in moving media¶
Find approximate solution to following linear PDE
using \(\theta\)-scheme discretization in time and arbitrary FE discretization in space with data
- \(\Omega = (0, 1)^2\)
- \(T = 10\)
- \(\Gamma_\mathrm{N} = \left\{ x = 0 \right\}\)
- \(\Gamma_\mathrm{D} = \left\{ x = 1 \right\} \cup \left\{ y = 0 \right\}\)
- \(g = 0.1\)
- \(K = 0.01\)
- \(\mathbf{b} = \left( -(y-\tfrac{1}{2}), x-\tfrac{1}{2} \right)\)
- \(f = \chi_{ B_{1/5}\left(\left[\frac{3}{4}, \frac{3}{4}\right]\right) }\)
- \(u_0(\mathbf{x}) = \left( 1 - 25 \operatorname{dist}\left(\mathbf{x}, \left[\frac{1}{4}, \frac{1}{4}\right]\right) \right) \chi_{ B_{1/5}\left(\left[\frac{1}{4}, \frac{1}{4}\right]\right) }\)
- \(u_\mathrm{D} = 0\)
where \(\chi_X\) is a characteristic function of set \(X\), \(B_R(\mathbf{z})\) is a ball of radius \(R\) and center \(\mathbf{z}\) and \(\operatorname{dist}(\mathbf{p}, \mathbf{q})\) is Euclidian distance between points \(\mathbf{p}\), \(\mathbf{q}\).
Task 1
Discretize the equation in time and write variational formulation of the problem.
Task 2
Build mesh, prepare facet function marking \(\Gamma_\mathrm{N}\) and \(\Gamma_\mathrm{D}\) and plot it to check its correctness.
Hint
You can follow the procedure from subdomains-poisson demo.
(Follow a construction of boundaries
object therein.)
mesh = UnitSquareMesh(10, 10, 'crossed')
# Create boundary markers
tdim = mesh.topology().dim()
boundary_parts = MeshFunction('size_t', mesh, tdim-1)
left = AutoSubDomain(lambda x: near(x[0], 0.0))
right = AutoSubDomain(lambda x: near(x[0], 1.0))
bottom = AutoSubDomain(lambda x: near(x[1], 0.0))
left .mark(boundary_parts, 1)
right .mark(boundary_parts, 2)
bottom.mark(boundary_parts, 2)
Task 3
Define expressions \(\mathbf{b}\), \(f\), \(u_0\) and plot them.
Hint
According to your personal taste, get hint at Expression class documentation
or any documented demo.
Use any kind of expression you wish (subclassing
Python Expression
, oneline C++, subclassing C++ Expression
).
# python Expression subclass
class B(Expression):
def eval(self, value, x):
vx = x[0] - 0.5
vy = x[1] - 0.5
value[0] = -vy
value[1] = vx
def value_shape(self):
return (2,)
b = B()
# oneline C++
b = Expression(("-(x[1] - 0.5)", "x[0] - 0.5"))
Task 4
Use facet markers from Task 2 to define DirichletBC
object
and Measure
for integration along \(\Gamma_\mathrm{N}\).
Hint
dsN = Measure("ds", subdomain_id=1, subdomain_data=boundary_parts)
Task 5
Now proceed to variational formulation and time-stepping loop. Write bilinear and linear form representing PDE. How is solution at previous time-step represented therein?
Hint
Use LinearVariationalProblem
and LinearVariationalSolver
classes so that solve
method of an instance of the latter is called
every time-step while nothing else is touched excepted updating value
of solution from previous time-step figuring in variational form. You
can use for instance Function.assign
method to do that.
Task 5
Add solution output for external visualisation, like Paraview.
Hint
# Create file for storing results
f = XDMFFile("results/u.xdmf")
u.rename("u", "temperature")
f.write(u, t)
Reference solution¶
Show/Hide Code
from dolfin import *
import matplotlib.pyplot as plt
# Create mesh and build function space
mesh = UnitSquareMesh(40, 40, 'crossed')
V = FunctionSpace(mesh, "Lagrange", 1)
# Create boundary markers
tdim = mesh.topology().dim()
boundary_parts = MeshFunction('size_t', mesh, tdim-1)
left = AutoSubDomain(lambda x: near(x[0], 0.0))
right = AutoSubDomain(lambda x: near(x[0], 1.0))
bottom = AutoSubDomain(lambda x: near(x[1], 0.0))
left .mark(boundary_parts, 1)
right .mark(boundary_parts, 2)
bottom.mark(boundary_parts, 2)
# Initial condition and right-hand side
ic = Expression("""pow(x[0] - 0.25, 2) + pow(x[1] - 0.25, 2) < 0.2*0.2
? -25.0 * ((pow(x[0] - 0.25, 2) + pow(x[1] - 0.25, 2)) - 0.2*0.2)
: 0.0""", degree=1)
f = Expression("""pow(x[0] - 0.75, 2) + pow(x[1] - 0.75, 2) < 0.2*0.2
? 1.0
: 0.0""", degree=1)
# Equation coefficients
K = Constant(1e-2) # thermal conductivity
g = Constant(0.01) # Neumann heat flux
b = Expression(("-(x[1] - 0.5)", "x[0] - 0.5"), degree=1) # convecting velocity
# Define boundary measure on Neumann part of boundary
dsN = Measure("ds", subdomain_id=1, subdomain_data=boundary_parts)
# Define steady part of the equation
def operator(u, v):
return ( K*inner(grad(u), grad(v)) - f*v + dot(b, grad(u))*v )*dx - K*g*v*dsN
# Define trial and test function and solution at previous time-step
u = TrialFunction(V)
v = TestFunction(V)
u0 = Function(V)
# Time-stepping parameters
t_end = 10
dt = 0.1
theta = Constant(0.5) # Crank-Nicolson scheme
# Define time discretized equation
F = (1.0/dt)*inner(u-u0, v)*dx + theta*operator(u, v) + (1.0-theta)*operator(u0, v)
# Define boundary condition
bc = DirichletBC(V, Constant(0.0), boundary_parts, 2)
# Prepare solution function and solver
u = Function(V)
problem = LinearVariationalProblem(lhs(F), rhs(F), u, bc)
solver = LinearVariationalSolver(problem)
# Prepare initial condition
u0.interpolate(ic)
# Create file for storing results
f = XDMFFile("results/u.xdmf")
# Time-stepping
t = 0.0
u.rename("u", "temperature")
u.interpolate(ic)
# Save initial solution
f.write(u, t)
# Open figure for plots
fig = plt.figure()
plt.show(block=False)
while t < t_end:
# Solve the problem
solver.solve()
# Store solution to file and plot
f.write(u, t)
p = plot(u, title='Solution at t = %g' % t)
if p is not None:
plt.colorbar(p)
fig.canvas.draw()
plt.clf()
# Move to next time step
u0.assign(u)
t += dt
# Report flux
n = FacetNormal(mesh)
flux = assemble(K*dot(grad(u), n)*dsN)
info('t = %g, flux = %g' % (t, flux))