Fix time lag issue for inversion (#370)

* Fix time lag issue for inversion

- t passed to cost_fn(t) was out of sync because update_forcings was being used so the optimal result was not correct for channel inversion
- Now need to pass the cost function through the export_func rather than through update_forcings
- Typo fixed in Tohoku Makefile as well

* Update inversion tools for parallelisation

- Previously would hang because you could have a processor with no Station coordinates within its partition of the mesh, causing issues

Author: cpjordan

examples/channel_inversion/inverse_problem.py

@@ -113,7 +113,7 @@
 cost_function = inv_manager.get_cost_function(solver_obj)
 
 # Solve and setup reduced functional
-solver_obj.iterate(update_forcings=cost_function)
+solver_obj.iterate(export_func=cost_function)  # note that export time should equal dt if not using a custom callback
 inv_manager.stop_annotating()
 
 # Run inversion

examples/tohoku_inversion/Makefile

@@ -9,7 +9,7 @@ invert:
 
 plot:
 	python3 plot_elevation_progress.py
-	python3 plot_elevation_optimised.py
+	python3 plot_elevation_optimized.py
 	python3 plot_convergence.py
 
 clean:

examples/tohoku_inversion/inverse_problem.py

@@ -101,7 +101,7 @@
 cost_function = inv_manager.get_cost_function(solver_obj, weight_by_variance=True)
 
 # Solve and setup the reduced functional
-solver_obj.iterate(update_forcings=cost_function)
+solver_obj.iterate(export_func=cost_function)
 inv_manager.stop_annotating()
 
 # Run inversion

thetis/inversion_tools.py

@@ -252,7 +252,8 @@ def get_cost_function(self, solver_obj, weight_by_variance=False):
                 var.dat.data[i] = numpy.var(self.sta_manager.observation_values[j])
             self.sta_manager.station_weight_0d.interpolate(1/var)
 
-        def cost_fn(t):
+        def cost_fn():
+            t = solver_obj.simulation_time
             misfit = self.sta_manager.eval_cost_function(t)
             self.J_misfit += misfit
             self.J += misfit
@@ -487,26 +488,27 @@ def construct_evaluator(self):
         # Construct timeseries interpolator
         self.station_interpolators = []
         self.local_station_index = []
-        for i in range(self.fs_points_0d.dof_dset.size):
-            # loop over local DOFs and match coordinates to observations
-            # NOTE this must be done manually as VertexOnlyMesh reorders points
-            x_mesh, y_mesh = mesh0d.coordinates.dat.data[i, :]
-            xy_diff = xy - numpy.array([x_mesh, y_mesh])
-            xy_dist = numpy.sqrt(xy_diff[:, 0]**2 + xy_diff[:, 1]**2)
-            j = numpy.argmin(xy_dist)
-            self.local_station_index.append(j)
-
-            x, y = xy[j, :]
-            t = self.observation_time[j]
-            v = self.observation_values[j]
-            x_mesh, y_mesh = mesh0d.coordinates.dat.data[i, :]
-
-            msg = 'bad station location ' \
-                f'{j} {i} {x} {x_mesh} {y} {y_mesh} {x-x_mesh} {y-y_mesh}'
-            assert numpy.allclose([x, y], [x_mesh, y_mesh]), msg
-            # create temporal interpolator
-            ip = interp1d(t, v, **interp_kw)
-            self.station_interpolators.append(ip)
+        if len(mesh0d.coordinates.dat.data[:]) > 0:
+            for i in range(self.fs_points_0d.dof_dset.size):
+                # loop over local DOFs and match coordinates to observations
+                # NOTE this must be done manually as VertexOnlyMesh reorders points
+                x_mesh, y_mesh = mesh0d.coordinates.dat.data[i, :]
+                xy_diff = xy - numpy.array([x_mesh, y_mesh])
+                xy_dist = numpy.sqrt(xy_diff[:, 0]**2 + xy_diff[:, 1]**2)
+                j = numpy.argmin(xy_dist)
+                self.local_station_index.append(j)
+
+                x, y = xy[j, :]
+                t = self.observation_time[j]
+                v = self.observation_values[j]
+                x_mesh, y_mesh = mesh0d.coordinates.dat.data[i, :]
+
+                msg = 'bad station location ' \
+                    f'{j} {i} {x} {x_mesh} {y} {y_mesh} {x-x_mesh} {y-y_mesh}'
+                assert numpy.allclose([x, y], [x_mesh, y_mesh]), msg
+                # create temporal interpolator
+                ip = interp1d(t, v, **interp_kw)
+                self.station_interpolators.append(ip)
 
         # Process start and end times for observations
         self.obs_start_times = numpy.array([