Update channel inversion example (#366)

- Correct file paths in example for plotting
- Correct VTKFile paths for import
- Load timesteps for the optimisation fields (one timestep less than existing files)
- Import from collections.abc (deprecation update)
- Correct Makefile

examples/channel_inversion/Makefile

@@ -1,6 +1,7 @@
 all: invert plot
 
 CONTROLS	= Bathymetry Manning InitialElev
+STATIONS    = stationA stationB stationC stationD stationE
 
 invert:
 	for controls in $(CONTROLS); do \
@@ -9,7 +10,9 @@ invert:
 
 plot:
 	for controls in $(CONTROLS); do \
-		python3 plot_elevation_progress.py -c $$controls
+		for station in $(STATIONS); do \
+			python3 plot_elevation_progress.py -c $$controls -s $$station; \
+		done; \
 	done
 
 clean:

examples/channel_inversion/inverse_problem.py

@@ -1,4 +1,5 @@
 from thetis import *
+import firedrake as fd
 from firedrake.adjoint import *
 import numpy
 import thetis.inversion_tools as inversion_tools
@@ -133,4 +134,4 @@
     oc.rename(name)
     print_function_value_range(oc, prefix='Optimal')
     if not no_exports:
-        VTKFile(f'{options.output_directory}/{name}_optimised.pvd').write(oc)
+        fd.VTKFile(f'{options.output_directory}/{name}_optimised.pvd').write(oc)

examples/channel_inversion/plot_elevation_progress.py

@@ -2,7 +2,7 @@
 import matplotlib
 import matplotlib.pyplot as plt
 import argparse
-from collections import Iterable
+from collections.abc import Iterable
 
 matplotlib.rc('font', size=10)
 
@@ -26,7 +26,7 @@
 
 ctrl_str = '-'.join(controls)
 station_str = '-'.join(station_names)
-inv_dir = f'outputs_{ctrl_str}-opt'
+inv_dir = f'outputs_new_{ctrl_str}-opt'
 
 nplots = len(station_names)
 
@@ -45,14 +45,15 @@
 
     g = f'{inv_dir}/diagnostic_timeseries_progress_{sta}_elev.hdf5'
     with h5py.File(g, 'r') as h5file:
+        iter_times = h5file['time'][:].flatten()
         iter_vals = h5file['elev'][:]
 
     niter = iter_vals.shape[0] - 1
 
     ax = next(ax_iter)
     ax.plot(time, vals, 'k:', zorder=3, label='observation', lw=1.3)
     for i, v in enumerate(iter_vals):
-        ax.plot(time, v, label=f'iteration {i}', lw=0.5)
+        ax.plot(iter_times, v, label=f'iteration {i}', lw=0.5)
     ax.set_title(f'{sta}', size='small')
     ax.set_ylabel('Elevation (m)')
     ax.grid(True)