Fix some warnings

Add pytest.ini that errors for any unknown warning.

pytest.ini

@@ -0,0 +1,7 @@
+[pytest]
+filterwarnings =
+    error
+    ignore:.*arg_id_to_.*
+    ignore:.*row_stack.*:DeprecationWarning
+    ignore::pytools.persistent_dict.RecommendedHashNotFoundWarning
+    ignore:.*pkg_resources.*:DeprecationWarning

thetis/coordsys.py

@@ -6,6 +6,7 @@
 import pyproj
 import numpy
 from abc import ABC, abstractmethod
+import warnings
 
 LL_WGS84 = pyproj.Proj(proj='latlong', datum='WGS84', errcheck=True)
 
@@ -147,8 +148,12 @@ def get_vector_rotation_matrix(trans, x, y, delta=None):
     """
     if delta is None:
         delta = 1e-6  # ~1 m in LL_WGS84
-    x1, y1 = trans.transform(x, y)
-    x2, y2 = trans.transform(x, y + delta)
+    with warnings.catch_warnings():
+        # supress spurious DeprecationWarning in numpy 1.25
+        # (see https://github.com/pyproj4/pyproj/issues/1307)
+        warnings.simplefilter("ignore", category=DeprecationWarning)
+        x1, y1 = trans.transform(x, y)
+        x2, y2 = trans.transform(x, y + delta)
     dxdl = (x2 - x1) / delta
     dydl = (y2 - y1) / delta
     theta = numpy.arctan2(-dxdl, dydl)

thetis/diagnostics.py

@@ -4,6 +4,7 @@
 from .utility import *
 from .configuration import *
 from abc import ABCMeta, abstractmethod
+from firedrake.ufl_expr import extract_unique_domain
 
 __all__ = ["VorticityCalculator2D", "HessianRecoverer2D", "KineticEnergyCalculator",
            "ShallowWaterDualWeightedResidual2D", "TracerDualWeightedResidual2D"]
@@ -259,8 +260,9 @@ def __init__(self, solver_obj, dual):
         if mesh2d.topological_dimension() != 2:
             dim = mesh2d.topological_dimension()
             raise ValueError(f"Expected a mesh of dimension 2, not {dim}")
-        if mesh2d != dual.ufl_domain():
-            raise ValueError(f"Mismatching meshes ({mesh2d} vs {func.ufl_domain()})")
+        dual_mesh = extract_unique_domain(dual)
+        if mesh2d != dual_mesh:
+            raise ValueError(f"Mismatching meshes ({mesh2d} vs {dual_mesh})")
         self.F = replace(self.form, {TestFunction(self.space): dual})
 
     @abstractmethod

thetis/interpolation.py

@@ -42,7 +42,7 @@ def set_fields(self, time):
 import os
 from .timezone import *
 from .log import *
-import scipy.spatial.qhull as qhull
+from scipy.spatial import Delaunay, QhullError
 import netCDF4
 from abc import ABC, abstractmethod
 from firedrake import *
@@ -190,7 +190,7 @@ def get_norm_params(x, scale=None):
         self.cannot_interpolate = False
         try:
             d = ngrid_xyz.shape[1]
-            tri = qhull.Delaunay(ngrid_xyz)
+            tri = Delaunay(ngrid_xyz)
             # NOTE this becomes expensive in 3D for npoints > 10k
             simplex = tri.find_simplex(ntarget_xyz)
             vertices = numpy.take(tri.simplices, simplex, axis=0)
@@ -208,7 +208,7 @@ def get_norm_params(x, scale=None):
                 from scipy.spatial import cKDTree
                 dist, ix = cKDTree(ngrid_xyz).query(ntarget_xyz[self.outside])
                 self.outside_to_nearest = ix
-        except qhull.QhullError as e:
+        except QhullError as e:
             if not dont_raise:
                 raise e
             self.cannot_interpolate = True
@@ -312,7 +312,7 @@ def _get_subset_nodes(grid_x, grid_y, target_x, target_y):
     orig_shape = grid_x.shape
     grid_xy = numpy.array((grid_x.ravel(), grid_y.ravel())).T
     target_xy = numpy.array((target_x.ravel(), target_y.ravel())).T
-    tri = qhull.Delaunay(grid_xy)
+    tri = Delaunay(grid_xy)
     simplex = tri.find_simplex(target_xy)
     vertices = numpy.take(tri.simplices, simplex, axis=0)
     nodes = numpy.unique(vertices.ravel())

thetis/options.py

@@ -881,11 +881,11 @@ class ModelOptions2d(CommonModelOptions):
 
         Note this is only relevant if `use_automatic_wetting_and_drying_alpha` is set to ``True``.
         """).tag(config=True)
-    tidal_turbine_farms = Dict(trait=List(TidalTurbineFarmOptions()), default_value={},
+    tidal_turbine_farms = Dict(value_trait=List(TidalTurbineFarmOptions()), default_value={},
                                help='Dictionary mapping subdomain ids to a list of TidalTurbineFarmOptions instances '
                                     'corresponding to one or more farms.')
 
-    discrete_tidal_turbine_farms = Dict(trait=List(DiscreteTidalTurbineFarmOptions()), default_value={},
+    discrete_tidal_turbine_farms = Dict(value_trait=List(DiscreteTidalTurbineFarmOptions()), default_value={},
                                         help='Dictionary mapping subdomain ids to a list of DiscreteTidalTurbineFarmOptions '
                                              'instances corresponding to one or more farms.')
 

thetis/shallowwater_eq.py

@@ -225,9 +225,9 @@ def __init__(self, space,
         p = self.function_space.ufl_element().degree()
         self.quad_degree = 2*p + 1
         self.dx = dx(degree=self.quad_degree,
-                     domain=self.function_space.ufl_domain())
+                     domain=self.function_space.mesh())
         self.dS = dS(degree=self.quad_degree,
-                     domain=self.function_space.ufl_domain())
+                     domain=self.function_space.mesh())
 
     def get_bnd_functions(self, eta_in, uv_in, bnd_id, bnd_conditions):
         """

thetis/utility.py

@@ -9,6 +9,7 @@
 import ufl  # NOQA
 from firedrake import *
 from firedrake.petsc import PETSc
+from firedrake.ufl_expr import extract_unique_domain
 from mpi4py import MPI  # NOQA
 from pyop2.profiling import timed_stage  # NOQA
 import numpy
@@ -1125,7 +1126,7 @@ def form2indicator(F):
     """
     if len(F.arguments()) > 0:
         raise ValueError("Input form should be 0-form")
-    mesh = F.ufl_domain()
+    mesh = extract_unique_domain(F)
     P0 = FunctionSpace(mesh, "DG", 0)
     p0test = TestFunction(P0)
     h = ufl.CellVolume(mesh)