# Heat equation in moving media¶

Find approximate solution to following linear PDE

\begin{align}\begin{aligned}u_t + \mathbf{b}\cdot\nabla{u} - \operatorname{div}(K \nabla u) &= f &&\quad\text{ in }\Omega\times(0, T),\\u &= u_\mathrm{D} &&\quad\text{ in }\Omega_\mathrm{D}\times(0, T),\\\tfrac{\partial u}{\partial\mathbf{n}} &= g &&\quad\text{ on }\Gamma_\mathrm{N}\times(0, T),\\u &= u_0 &&\quad\text{ on }\Omega\times\{0\}\end{aligned}\end{align}

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}$$.

Discretize the equation in time and write variational formulation of the problem.

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))
right  = AutoSubDomain(lambda x: near(x, 1.0))
bottom = AutoSubDomain(lambda x: near(x, 0.0))
left  .mark(boundary_parts, 1)
right .mark(boundary_parts, 2)
bottom.mark(boundary_parts, 2)


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.5
vy = x - 0.5
value = -vy
value =  vx
def value_shape(self):
return (2,)

b = B()

# oneline C++
b = Expression(("-(x - 0.5)", "x - 0.5"))


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)


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.

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

Download 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))
right  = AutoSubDomain(lambda x: near(x, 1.0))
bottom = AutoSubDomain(lambda x: near(x, 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.25, 2) + pow(x - 0.25, 2) < 0.2*0.2
? -25.0 * ((pow(x - 0.25, 2) + pow(x - 0.25, 2)) - 0.2*0.2)
: 0.0""", degree=1)
f = Expression("""pow(x - 0.75, 2) + pow(x - 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 - 0.5)", "x - 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):

# 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)