# 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"),
path.join(path.dirname(path.realpath(__file__)),"../../../cbcflow-data/cylinder_2k.xml.gz"),
path.join(path.dirname(path.realpath(__file__)),"../../../cbcflow-data/cylinder_8k.xml.gz"),
path.join(path.dirname(path.realpath(__file__)),"../../../cbcflow-data/cylinder_32k.xml.gz"),
path.join(path.dirname(path.realpath(__file__)),"../../../cbcflow-data/cylinder_129k.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):
```

### Parameters¶

This class inherit from the `Parameterized`

class,
allowing for parameters in the class interface. We supply default parameters to
the problem:

```
@classmethod
def default_params(cls):
params = NSProblem.default_params()
params.replace(
# Time parameters
T=5.0,
dt=0.1,
# Physical parameters
rho=1.0,
mu=1.0/1000.0,
)
params.update(
# Spatial parameters
refinement_level=0,
)
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.

### Constructor¶

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)
facet_domains.set_all(4)
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

```
show_schemes()
```

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

```
show_fields()
```

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 = [
Pressure(plot_and_save),
Velocity(plot_and_save),
StreamFunction(plot_and_save),
]
```

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:

```
postprocessor.add_fields(fields)
```

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

```
solver.solve()
```