rerun CI for #383

.github/workflows/build.yml

@@ -18,7 +18,7 @@ jobs:
     runs-on: self-hosted
     # The docker container to use.
     container:
-      image: firedrakeproject/firedrake-vanilla:latest
+      image: firedrakeproject/firedrake-vanilla:2024-10
     steps:
       - uses: actions/checkout@v2
       - name: Cleanup
@@ -35,12 +35,12 @@ jobs:
         run: |
           . /home/firedrake/firedrake/bin/activate
           python $(which firedrake-clean)
-          python -m pytest -n 12 -v test
+          python -m pytest --durations=0 -n 12 -v test
       - name: Test Thetis adjoint
         run: |
           . /home/firedrake/firedrake/bin/activate
           python $(which firedrake-clean)
-          python -m pytest -n 12 -v test_adjoint
+          python -m pytest --durations=0 -n 12 -v test_adjoint
       - name: Lint
         if: ${{ always() }}
         run: |

test/examples/test_examples.py

@@ -39,7 +39,6 @@
     'tidalfarm/tidalfarm.py',
     'tidal_barrage/plotting.py',
     'channel_inversion/plot_elevation_progress.py',
-    'channel_inversion/inverse_problem.py',
     'tohoku_inversion/okada.py',
     'tohoku_inversion/plot_convergence.py',
     'tohoku_inversion/plot_elevation_initial_guess.py',

test/swe2d/test_anisotropic.py

@@ -73,7 +73,21 @@ def run(solve_adjoint=False, mesh=None, **model_options):
     options.use_grad_depth_viscosity_term = False
     options.no_exports = True
     options.update(model_options)
-    solver_obj.create_function_spaces()
+    solver_obj.create_equations()
+
+    # recover Hessian
+    if not solve_adjoint:
+        if 'hessian_2d' in field_metadata:
+            field_metadata.pop('hessian_2d')
+        P1t_2d = get_functionspace(mesh2d, 'CG', 1, tensor=True)
+        hessian_2d = Function(P1t_2d)
+        u_2d = solver_obj.fields.uv_2d[0]
+        hessian_calculator = thetis.diagnostics.HessianRecoverer2D(u_2d, hessian_2d)
+        solver_obj.add_new_field(hessian_2d,
+                                 'hessian_2d',
+                                 'Hessian of x-velocity',
+                                 'Hessian2d',
+                                 preproc_func=hessian_calculator)
 
     # apply boundary conditions
     solver_obj.bnd_functions['shallow_water'] = {
@@ -121,36 +135,12 @@ def bump(mesh, locs, scale=1.0):
     farm_options.turbine_options.diameter = D
     C_T = thrust_coefficient * correction
     farm_options.turbine_options.thrust_coefficient = C_T
-    solver_obj.options.tidal_turbine_farms['everywhere'] = [farm_options]
-    solver_obj.create_equations()
+    solver_obj.options.tidal_turbine_farms['everywhere'] = farm_options
 
-    # recover Hessian
-    if not solve_adjoint:
-        if 'hessian_2d' in field_metadata:
-            field_metadata.pop('hessian_2d')
-        P1t_2d = get_functionspace(mesh2d, 'CG', 1, tensor=True)
-        hessian_2d = Function(P1t_2d)
-        u_2d = solver_obj.fields.uv_2d[0]
-        hessian_calculator = thetis.diagnostics.HessianRecoverer2D(u_2d, hessian_2d)
-        solver_obj.add_new_field(hessian_2d,
-                                 'hessian_2d',
-                                 'Hessian of x-velocity',
-                                 'Hessian2d',
-                                 preproc_func=hessian_calculator)
-
-    # Apply initial guess of inflow velocity and solve
+    # apply initial guess of inflow velocity, solve and return number of nonlinear solver iterations
     solver_obj.assign_initial_conditions(uv=inflow_velocity)
     solver_obj.iterate()
-
     if not solve_adjoint:
-        # Check that turbines have been picked up
-        assert not np.allclose(solver_obj.fields.uv_2d.dat.data[0], 0.5, atol=1e-3)
-
-        # Check the turbine density has been set up correctly
-        total_turbine_density = assemble(solver_obj.tidal_farms[0].turbine_density * dx)
-        assert np.isclose(total_turbine_density, 2, atol=0.01), f"Expected 2, but got {total_turbine_density}"
-
-        # Return number of nonlinear solver iterations
         return solver_obj.timestepper.solver.snes.getIterationNumber()
 
     # quantity of interest: power output

test_adjoint/examples/test_examples.py

@@ -15,7 +15,7 @@
 # list of all adjoint examples to run
 adjoint_files = [
     'tidalfarm/tidalfarm.py',
-    'channel_inversion/inverse_problem.py',
+    # 'channel_inversion/inverse_problem.py',  # FIXME requires obs time series
 ]
 
 cwd = os.path.abspath(os.path.dirname(__file__))
@@ -34,8 +34,6 @@ def example_file(request):
 
 def test_examples(example_file, tmpdir, monkeypatch):
     assert os.path.isfile(example_file), 'File not found {:}'.format(example_file)
-    if 'examples/channel_inversion/inverse_problem.py' in example_file:
-        pytest.xfail("Known issue with Firedrake and mixed function spaces. See Firedrake issue #3368.")
     # copy mesh files
     source = os.path.dirname(example_file)
     for f in glob.glob(os.path.join(source, '*.msh')):

thetis/callback.py

@@ -529,9 +529,6 @@ def __init__(self, solver_obj,
         self.field_names = field_names
         self._name = name
 
-        # initialise interpolation functions using vom_interpolator_functions
-        self.interp_functions = vom_interpolator_functions(solver_obj, field_names, detector_locations)
-
     @property
     def name(self):
         return self._name
@@ -542,8 +539,7 @@ def variable_names(self):
 
     def _values_per_field(self, values):
         """
-        Given all values evaulated in a detector location, return the values per field
-        """
+        Given all values evaulated in a detector location, return the values per field"""
         i = 0
         result = []
         for dim in self.field_dims:
@@ -558,11 +554,7 @@ def message_str(self, *args):
             for name, values in zip(self.detector_names, args))
 
     def _evaluate_field(self, field_name):
-        field = self.solver_obj.fields[field_name]
-        f_at_points, f_at_input_points = self.interp_functions[field_name]
-        f_at_points.interpolate(field)
-        f_at_input_points.interpolate(f_at_points)
-        return f_at_input_points.dat.data_ro[:]
+        return self.solver_obj.fields[field_name](self.detector_locations)
 
     def __call__(self):
         """
@@ -692,16 +684,10 @@ def _initialize(self):
         xyz = (self.x, self.y, self.z) if self.on_sphere else (self.x, self.y)
         self.xyz = numpy.array([xyz])
 
-        # initialise interpolation functions using vom_interpolator_functions
-        self.interp_functions = vom_interpolator_functions(self.solver_obj, self.fieldnames, self.xyz)
-
         # test evaluation
         try:
             if self.eval_func is None:
-                field = self.solver_obj.fields[self.fieldnames[0]]
-                f_at_points, f_at_input_points = self.interp_functions[self.fieldnames[0]]
-                f_at_points.interpolate(field)
-                f_at_input_points.interpolate(f_at_points)
+                self.solver_obj.fields.bathymetry_2d.at(self.xyz, tolerance=self.tolerance)
             else:
                 self.eval_func(self.solver_obj.fields.bathymetry_2d, self.xyz, tolerance=self.tolerance)
         except PointNotInDomainError as e:
@@ -721,10 +707,7 @@ def __call__(self):
             try:
                 field = self.solver_obj.fields[fieldname]
                 if self.eval_func is None:
-                    f_at_points, f_at_input_points = self.interp_functions[fieldname]
-                    f_at_points.interpolate(field)
-                    f_at_input_points.interpolate(f_at_points)
-                    val = f_at_input_points.dat.data_ro[:]
+                    val = field.at(self.xyz, tolerance=self.tolerance)
                 else:
                     val = self.eval_func(field, self.xyz, tolerance=self.tolerance)
                 arr = numpy.array(val)

thetis/exporter.py

@@ -22,11 +22,11 @@ def get_visu_space(fs):
     """
     mesh = fs.mesh()
     family = 'Lagrange' if is_cg(fs) else 'Discontinuous Lagrange'
-    if len(fs.value_shape) == 1:
-        dim = fs.value_shape[0]
+    if len(fs.ufl_element().value_shape) == 1:
+        dim = fs.ufl_element().value_shape[0]
         visu_fs = get_functionspace(mesh, family, 1, family, 1,
                                     vector=True, dim=dim)
-    elif len(fs.value_shape) == 2:
+    elif len(fs.ufl_element().value_shape) == 2:
         visu_fs = get_functionspace(mesh, family, 1, family, 1,
                                     tensor=True)
     else:

thetis/utility.py

@@ -328,7 +328,7 @@ def get_facet_mask(function_space, facet='bottom'):
         Here we assume that the mesh has been extruded upwards (along positive
         z axis).
     """
-    from finat.element_factory import create_element as create_finat_element
+    from tsfc.finatinterface import create_element as create_finat_element
 
     # get base element
     elem = get_extruded_base_element(function_space.ufl_element())
@@ -605,8 +605,8 @@ def compute_elem_height(zcoord, output):
                 }
             }
         }""" % {'nodes': zcoord.cell_node_map().arity,
-                'func_dim': zcoord.function_space().block_size,
-                'output_dim': output.function_space().block_size},
+                'func_dim': zcoord.function_space().value_size,
+                'output_dim': output.function_space().value_size},
         'my_kernel')
     op2.par_loop(
         kernel, fs_out.mesh().cell_set,
@@ -1154,38 +1154,3 @@ def form2indicator(F):
         },
     )
     return indicator
-
-
-@PETSc.Log.EventDecorator("thetis.vom_interpolator_functions")
-def vom_interpolator_functions(solver_obj, field_names, locations):
-    r"""
-    Creates function spaces and associated Functions for interpolation
-    on a VertexOnlyMesh (VOM) and returns them for reuse.
-
-    :arg solver_obj: Thetis solver object
-    :arg field_names: List of field names to create functions for.
-    :arg locations: List of locations for interpolation.
-    :return: A dictionary mapping field names to a tuple of (f_at_points, f_at_input_points)
-             which are Functions for interpolation.
-    """
-    vom = VertexOnlyMesh(solver_obj.mesh2d, locations, redundant=True)
-
-    functions_dict = {}
-
-    for field_name in field_names:
-        field = solver_obj.fields[field_name]
-
-        if isinstance(field.function_space().ufl_element(), VectorElement):
-            P0DG = VectorFunctionSpace(vom, "DG", 0)
-            P0DG_input_ordering = VectorFunctionSpace(vom.input_ordering, "DG", 0)
-        else:
-            P0DG = FunctionSpace(vom, "DG", 0)
-            P0DG_input_ordering = FunctionSpace(vom.input_ordering, "DG", 0)
-
-        f_at_points = Function(P0DG)
-        f_at_input_points = Function(P0DG_input_ordering)
-
-        # Store the Functions in the dictionary keyed by field name
-        functions_dict[field_name] = (f_at_points, f_at_input_points)
-
-    return functions_dict

thetis/utility3d.py

@@ -527,8 +527,8 @@ def __init__(self, input_2d, output_3d, elem_height=None):
                     }
                 }
             }""" % {'nodes': self.fs_2d.finat_element.space_dimension(),
-                    'func2d_dim': self.input_2d.function_space().block_size,
-                    'func3d_dim': self.fs_3d.block_size,
+                    'func2d_dim': self.input_2d.function_space().value_size,
+                    'func3d_dim': self.fs_3d.value_size,
                     'v_nodes': n_vert_nodes},
             'my_kernel')
 
@@ -664,8 +664,8 @@ def __init__(self, input_3d, output_2d,
                         }
                     }
                 }""" % {'nodes': self.output_2d.cell_node_map().arity,
-                        'func2d_dim': self.output_2d.function_space().block_size,
-                        'func3d_dim': self.fs_3d.block_size},
+                        'func2d_dim': self.output_2d.function_space().value_size,
+                        'func3d_dim': self.fs_3d.value_size},
                 'my_kernel')
         else:
             self.kernel = op2.Kernel("""
@@ -676,8 +676,8 @@ def __init__(self, input_3d, output_2d,
                         }
                     }
                 }""" % {'nodes': self.output_2d.cell_node_map().arity,
-                        'func2d_dim': self.output_2d.function_space().block_size,
-                        'func3d_dim': self.fs_3d.block_size},
+                        'func2d_dim': self.output_2d.function_space().value_size,
+                        'func3d_dim': self.fs_3d.value_size},
                 'my_kernel')
 
         if self.do_hdiv_scaling:
@@ -771,8 +771,8 @@ def __init__(self, solver):
                     }
                 }
             }""" % {'nodes': self.fs_2d.finat_element.space_dimension(),
-                    'func2d_dim': self.fs_2d.block_size,
-                    'func3d_dim': self.fs_3d.block_size,
+                    'func2d_dim': self.fs_2d.value_size,
+                    'func3d_dim': self.fs_3d.value_size,
                     'v_nodes': n_vert_nodes},
             'my_kernel')
 
@@ -790,8 +790,8 @@ def __init__(self, solver):
                     }
                 }
             }""" % {'nodes': self.fs_2d.finat_element.space_dimension(),
-                    'func2d_dim': self.fs_2d.block_size,
-                    'func3d_dim': self.fs_3d.block_size,
+                    'func2d_dim': self.fs_2d.value_size,
+                    'func3d_dim': self.fs_3d.value_size,
                     'v_nodes': n_vert_nodes},
             'my_kernel')