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