# Burgers equation
# ================
#
# The Burgers equation is a non-linear equation for the advection and
# diffusion of momentum. Here we choose to write the Burgers equation in
# two dimensions to demonstrate the use of vector function spaces:
#
# .. math::
#
# \frac{\partial u}{\partial t} + (u\cdot\nabla) u - \nu\nabla^2 u = 0
#
# (n\cdot \nabla) u = 0 \ \textrm{on}\ \Gamma
#
# where :math:`\Gamma` is the domain boundary and :math:`\nu` is a
# constant scalar viscosity. The solution :math:`u` is sought in some
# suitable vector-valued function space :math:`V`. We take the inner
# product with an arbitrary test function :math:`v\in V` and integrate
# the viscosity term by parts:
#
# .. math::
#
# \int_\Omega\frac{\partial u}{\partial t}\cdot v +
# ((u\cdot\nabla) u)\cdot v + \nu\nabla u\cdot\nabla v \ \mathrm d x = 0.
#
# The boundary condition has been used to discard the surface
# integral. Next, we need to discretise in time. For simplicity and
# stability we elect to use a backward Euler discretisation:
#
# .. math::
#
# \int_\Omega\frac{u^{n+1}-u^n}{dt}\cdot v +
# ((u^{n+1}\cdot\nabla) u^{n+1})\cdot v + \nu\nabla u^{n+1}\cdot\nabla v \ \mathrm d x = 0.
#
# We can now proceed to set up the problem. We choose a resolution and set up a square mesh::
from firedrake import *
n = 30
mesh = UnitSquareMesh(n, n)
# We choose degree 2 continuous Lagrange polynomials. We also need a
# piecewise linear space for output purposes::
V = VectorFunctionSpace(mesh, "CG", 2)
V_out = VectorFunctionSpace(mesh, "CG", 1)
# We also need solution functions for the current and the next
# timestep. Note that, since this is a nonlinear problem, we don't
# define trial functions::
u_ = Function(V, name="Velocity")
u = Function(V, name="VelocityNext")
v = TestFunction(V)
# For this problem we need an initial condition::
x = SpatialCoordinate(mesh)
ic = project(as_vector([sin(pi*x[0]), 0]), V)
# We start with current value of u set to the initial condition, but we
# also use the initial condition as our starting guess for the next
# value of u::
u_.assign(ic)
u.assign(ic)
# :math:`\nu` is set to a (fairly arbitrary) small constant value::
nu = 0.0001
# The timestep is set to produce an advective Courant number of
# around 1. Since we are employing backward Euler, this is stricter than
# is required for stability, but ensures good temporal resolution of the
# system's evolution::
timestep = 1.0/n
# Here we finally get to define the residual of the equation. In the advection
# term we need to contract the test function :math:`v` with
# :math:`(u\cdot\nabla)u`, which is the derivative of the velocity in the
# direction :math:`u`. This directional derivative can be written as
# ``dot(u,nabla_grad(u))`` since ``nabla_grad(u)[i,j]``:math:`=\partial_i u_j`.
# Note once again that for a nonlinear problem, there are no trial functions in
# the formulation. These will be created automatically when the residual
# is differentiated by the nonlinear solver::
F = (inner((u - u_)/timestep, v)
+ inner(dot(u,nabla_grad(u)), v) + nu*inner(grad(u), grad(v)))*dx
# We now create an object for output visualisation::
outfile = File("burgers.pvd")
# Output only supports visualisation of linear fields (either P1, or
# P1DG). In this example we project to a linear space by hand. Another
# option is to let the :class:`~.File` object manage the decimation. It
# supports both interpolation to linears (the default) or projection (by
# passing ``project_output=True`` when creating the :class:`~.File`).
# Outputting data is carried out using the :meth:`~.File.write` method
# of :class:`~.File` objects::
outfile.write(project(u, V_out, name="Velocity"))
# Finally, we loop over the timesteps solving the equation each time and
# outputting each result::
t = 0.0
end = 0.5
while (t <= end):
solve(F == 0, u)
u_.assign(u)
t += timestep
outfile.write(project(u, V_out, name="Velocity"))
# A python script version of this demo can be found `here `__.