Flow Around A Cylinder

This tutorial demonstrate how one can use cbcflow to solve a simple problem, namely a flow around a cylinder, inducing a vortex street behind the cylinder.

The source code for this can be found in FlowAroundCylinder.py.

We start by importing cbcflow and dolfin:

from cbcflow import *
from dolfin import *

Specifying the domain

The meshes for this problem is pregenerated, and is specified at the following locations:

from os import path

files = [path.join(path.dirname(path.realpath(__file__)),"../../../cbcflow-data/cylinder_0.6k.xml.gz"),

This requires that you have installed the demo data, as specified in Demos.

The domain is based on a rectangle with corners in (0,0), (0,1), (10,0) and (10,1). The cylinder is centered in (2,0.5) with radius of 0.12. The different boundaries of the domain is specified as:

class LeftBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0.0)

class RightBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 10.0)

class Cylinder(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and (sqrt((x[0]-2.0)**2+(x[1]-0.5)**2) < 0.12+DOLFIN_EPS)

class Wall(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and (near(x[1], 0.0) or near(x[1], 1.0))

Defining a NSProblem

To define a problem class recognized by cbcflow, the class must inherit from NSProblem:

class FlowAroundCylinder(NSProblem):


This class inherit from the Parameterized class, allowing for parameters in the class interface. We supply default parameters to the problem:

def default_params(cls):
    params = NSProblem.default_params()
        # Time parameters
        # Physical parameters
        # Spatial parameters
    return params

This takes the default parameters from NSProblem and replaces some parameters common for all NSProblems. We set the end time to 5.0 with a timestep of 0.1, the density \(\rho=1.0\) and dynamic viscosity \(\mu=0.001\). In addition, we add a new parameter, refinement_level, to determine which of the previously specified mesh files to use.


To initiate a FlowAroundCylinder-instance, we load the mesh and initialize the geometry:

def __init__(self, params=None):
    NSProblem.__init__(self, params)

    # Load mesh
    mesh = Mesh(files[self.params.refinement_level])

    # Create boundary markers
    facet_domains = FacetFunction("size_t", mesh)
    Wall().mark(facet_domains, 0)
    Cylinder().mark(facet_domains, 0)
    LeftBoundary().mark(facet_domains, 1)
    RightBoundary().mark(facet_domains, 2)

    # Store mesh and markers
    self.initialize_geometry(mesh, facet_domains=facet_domains)

The first call to NSProblem.__init__ updates the default parameters with any parameters passed to the constructor as a dict or ParamDict. This sets params as an attribute to self. We load the mesh from a string defined in the files-list, and define its domains. Finally, we call self.initialize_geometry to attach facet_domains to the mesh, and the mesh to self.

Initial conditions

At the initial time, the fluid is set to rest, with a zero pressure gradient. These initial conditions are prescribed by

def initial_conditions(self, spaces, controls):
    c0 = Constant(0)
    u0 = [c0, c0]
    p0 = c0
    return (u0, p0)

The argument spaces is a NSSpacePool helper object used to construct and contain the common function spaces related to the Navier-Stokes solution. This is used to limit the memory consumption and simplify the interface, so that you can, for example, call spaces.DV to get the tensor valued gradient space of the velocity regardless of velocity degree.

The argument controls is used for adjoint problems, and can be disregarded for simple forward problems such as this.

Boundary conditions

As boundary conditions, we set no-slip conditions on the cylinder, at y=0.0 and y=1.0. At the inlet we set a uniform velocity of (1.0,0.0), and zero-pressure boundary condition at the outlet.

To determine domain to apply boundary condition, we utilize the definition of facet_domains from the constructor.

def boundary_conditions(self, spaces, u, p, t, controls):
    c0 = Constant(0)
    c1 = Constant(1)

    # Create inflow and no-slip boundary conditions for velocity
    inflow = ([c1, c0], 1)
    noslip = ([c0, c0], 0)

    # Create boundary conditions for pressure
    bcp0 = (c0, 2)

    # Collect and return
    bcu = [inflow, noslip]
    bcp = [bcp0]
    return (bcu, bcp)

The way these boundary conditions are applied to the equations are determined by the scheme used to solve the equation. Note that the ordering of the boundary conditions in the lists bcu and bcp matters, in this case setting noslip last ensures the velocity is zero in the inflow corner.

Setting up the solver

Now that our FlowAroundCylinder-class is sufficiently defined, we can start thinking about solving our equations. We start by creating an instance of FlowAroundCylinder class:

problem = FlowAroundCylinder({"refinement_level": 2})

Note that we can pass a dict to the constructor to set, in this example, the desired refinement level of our mesh.

Selecting a scheme

Several schemes are implemented in cbcflow, but only a couple are properly tested and validated, and hence classified as official. Use


to list all schemes available, both official and unofficial.

In our application we select a very efficient operator-splitting scheme, IPCS,

scheme = IPCS()

Setting up postprocessing

The postprocessing is set up to determine what we want to do with our obtained solution. We start by creating a PostProcessor to handle all the logic:

casedir = "results_demo_%s_%s" % (problem.shortname(), scheme.shortname())
postprocessor = PostProcessor({"casedir": casedir})

The casedir parameter points the postprocessor to the directory where it should save the data it is being asked to save. By default, it stores the mesh, all parameters and a play log in that directory.

Then, we have to choose what we want to compute from the solution. The command


lists all available PPField to compute from the solution.

In this case, we are interested in the velocity, pressure and stream function, and we wish to both plot and save these at every timestep:

plot_and_save = dict(plot=True, save=True)
fields = [

With no saveformat prescribed, the postprocessor will choose default saveformats based on the type of data. You can use

print PPField.default_parameters()

to see common parameters of these fields.

Finally, we need to add these fields to the postprocessor:


Solving the problem

We now have instances of the classes NSProblem, NSScheme, and PostProcessor.

These can be combined in a general class to handle the logic between the classes, namely a NSSolver instance:

solver = NSSolver(problem, scheme, postprocessor)

This class has functionality to pass the solution from scheme on to the postprocessor, report progress to screen and so on. To solve the problem, simply execute