# 2D channel example¶

This example demonstrates a depth-averaged 2D simulation in a closed rectangular domain, where the flow is forced by an initial pertubation in the water elevation field.

We begin by importing Thetis and creating a rectangular mesh with `RectangleMesh()`

.
The domain is 40 km long and 2 km wide.
We generate 25 elements in the along-channel direction and 2 in the
cross-channel direction:

```
from thetis import *
lx = 40e3
ly = 2e3
nx = 25
ny = 2
mesh2d = RectangleMesh(nx, ny, lx, ly)
```

Next we define a bathymetry function in the 2D mesh, using continuous linear elements. In this example we set the bathymetry to constant 20 m depth:

```
P1_2d = FunctionSpace(mesh2d, 'CG', 1)
bathymetry_2d = Function(P1_2d, name='Bathymetry')
depth = 20.0
bathymetry_2d.assign(depth)
```

Note

See Firedrake manual for more information on mesh generation, functions and function spaces.

We are now ready to create a 2D solver object, and set some options:

```
# total duration in seconds
t_end = 2 * 3600
# export interval in seconds
t_export = 100.0
solver_obj = solver2d.FlowSolver2d(mesh2d, bathymetry_2d)
options = solver_obj.options
options.simulation_export_time = t_export
options.simulation_end_time = t_end
```

Here we simply define the total duration of the run, and the
export interval. See `ModelOptions`

for more information about the
available options.

Next we define the used time integrator for the shallow water equations and set the time step:

```
options.swe_timestepper_type = 'CrankNicolson'
options.timestep = 50.0
```

Because Crank-Nicolson is an uncondionally stable method, we can set the time step freely.

We then define the initial condition for elevation. We begin by creating a function (in the same linear continous function space):

```
elev_init = Function(P1_2d, name='initial elevation')
```

We then need to define an analytical expression the the x,y coordinates of the
mesh. To this end, we use
`SpatialCoordinate`

and define a UFL expression (see
Firedrake’s interpolation manual
for more information):

```
xy = SpatialCoordinate(mesh2d)
gauss_width = 4000.
gauss_ampl = 2.0
gauss_expr = gauss_ampl * exp(-((xy[0]-lx/2)/gauss_width)**2)
```

This defines a 2 m tall Gaussian hill in the x-direction in the middle on the domain. We can then interpolate this expression on the function:

```
elev_init.interpolate(gauss_expr)
```

and set this function as an initial condition to the elevation field:

```
solver_obj.assign_initial_conditions(elev=elev_init)
```

Model setup is now complelete. We run the model by issuing:

```
solver_obj.iterate()
```

While the model is running, Thetis prints some statistics on the command line:

```
0 0 T= 0.00 eta norm: 6251.2574 u norm: 0.0000 0.00
1 2 T= 100.00 eta norm: 5905.0262 u norm: 1398.1128 0.76
2 4 T= 200.00 eta norm: 5193.5227 u norm: 2377.8512 0.03
3 6 T= 300.00 eta norm: 4656.5334 u norm: 2856.5165 0.03
...
```

The first column is the export index, the second one the number of executed
time steps, followed by the simulation time. `eta norm`

and `u norm`

are
the L2 norms of the elevation and depth averaged velocity fields, respectively.
The last column stands for the (approximate) wall-clock time between exports.

The simulation terminates once the end time is reached. See outputs and visualization page on how to visualize the results.

This tutorial can be dowloaded as a Python script here.