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}\).

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}\).

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

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