improve readibility of the Velocity class and add exceptions

features/advection_term.feature

@@ -12,3 +12,4 @@ Scenario: Generate the convection term of the convection-diffusion PDE as a line
     And the result of the matrice vector product of the flux part of the adjoint operator has the expected shape
     And the result of the matrice vector product of the flux part of the adjoint operator has the expected value
     And the linear operator is updated as expected at a given time
+    And the initialization fails if the velocity fields is not an object of the class Velocity

features/steady_velocity_field.feature

@@ -16,6 +16,4 @@ Scenario: Generate the steady velocity field as a velocity object
     And the matrices of boolean of the upwind cells has the expected values
     And the maximum of the horizontal velocity has the expected value
     And the maximum of the vertical velocity has the expected value
-    And the initialization fails if the number of dimension of the velocity field is not correct
-    Then the initialization fails if the velocity fields shape do not match
     Then the update at a given time does not change the values

features/velocity_field.feature

@@ -17,8 +17,3 @@ Scenario: Generate the velocity field as a Velocity object
     And the velocity field is correctly updated at a given time
     And the maximum of the horizontal velocity has the expected value
     And the maximum of the vertical velocity has the expected value
-    And the initialization fails if the number of dimension of the velocity field is not correct
-    And the initialization fails if the velocity fields shape do not match the shape of the mesh
-    And the initialization fails if the velocity fields shape do not match at vertical and horizontal interfaces
-    And the initialization fails if the velocity fields shape do not match with the shape of the time vector
-    And the update fails if the time is not between the lowest and largest times contained in the time vector

features/velocity_field_inputs.feature

+Feature: Velocity object
+
+Scenario: Test the validity check of the inputs of the Velocity object
+    Given a rectangular 2D mesh
+    Given a time vector
+    Given a velocity field at the vertical interfaces
+    Given a velocity field at the horizontal interfaces
+    Given an invalid velocity field at the vertical interfaces
+    Given an invalid velocity field at the horizontal interfaces
+    When initialize the velocity object
+    Then the initialization fails if the number of dimension of the velocity field is not correct
+    And the initialization fails if the velocity fields shape do not match the shape of the mesh
+    And the initialization fails if the velocity fields shape do not match at vertical and horizontal interfaces
+    And the initialization fails if the velocity fields shape do not match with the shape of the time vector
+    And the initialization fails if the velocity fields do not cover the entire time window
+    And the update fails if the time is not between the lowest and largest times contained in the time vector

src/pheromone_dispersion/advection_operator.py

@@ -2,6 +2,7 @@ import numpy as np
 from scipy.sparse.linalg import LinearOperator as LinOp
 
 from pheromone_dispersion.velocity import Velocity
+from pheromone_dispersion.velocity import _validate_velocity_field
 
 """
 Module that contains the implementation of the linear operator of the advection term of the pheromone propagation model.
@@ -36,6 +37,8 @@ def sum_advection_flux_given_U(U, x, msh):
     # for now: no income
     # change for: a given income
 
+    _validate_velocity_field(U, msh)
+
     shift_vertical = np.zeros((msh.y_horizontal_interface.size, msh.x.size))
     flux_vertical = np.zeros((msh.y_horizontal_interface.size, msh.x.size))
 
@@ -110,6 +113,7 @@ class Advection(LinOp):
         msh : ~pheromone_dispersion.geom.MeshRect2D
             The geometry of the domain.
         """
+        _validate_velocity_field(U, msh)
         self.U = U
         if U.t is None:
             self.minus_U = Velocity(msh, -U.at_vertical_interface, -U.at_horizontal_interface)

src/pheromone_dispersion/geom.py

@@ -33,7 +33,9 @@ class MeshRect2D:
     mass_cell : float
         The volume of mesh cells :math:`\Delta x\times \Delta y`.
     t : float
-        The current time :math:`t` of the modelling, initialized to :math:`t=0s`.
+        The current time :math:`t` of the modelling, initialized to :math:`t=t_0s`.
+    t_0 : float
+        The initial time :math:`t_0` of the modeling time window.
     T_final : float
         The final time :math:`T` of the modeling time window.
     dt : float
@@ -51,7 +53,7 @@ class MeshRect2D:
         Coordinates of the origin of the mesh :math:`(x_0, y_0)`. By default set to :math:`(x_0, y_0)=(0, 0)`.
     """
 
-    def __init__(self, L_x, L_y, dx, dy, T_final, X_0=None):
+    def __init__(self, L_x, L_y, dx, dy, T_final, X_0=None, t_0=None):
         r"""
         Constructor method.
 
@@ -69,9 +71,12 @@ class MeshRect2D:
             The final time :math:`T` of the modeling time window.
         X_0 : tuple of float, default: None
             Coordinates of the origin of the mesh :math:`(x_0, y_0)`. If `None`, set to :math:`(x_0, y_0)=(0, 0)`.
+        t_0 : float, default: None
+            The initial time :math:`t_0` of the modeling time window. If `None`, set to :math:`t_0=0`
         """
         # To do:
         # Add exceptions to ensure all inputs are float.
+        # Add exceptions to make sure that t_0<T_f
 
         self.L_x = L_x
         self.L_y = L_y
@@ -85,13 +90,15 @@ class MeshRect2D:
             ny += 1
         if X_0 is None:
             X_0 = (0.0, 0.0)
+        if t_0 is None:
+            self.t_0 = 0.0
         self.x_vertical_interface = np.arange(0, (nx + 0.5) * dx, dx) + X_0[0]
         self.x = 0.5 * (self.x_vertical_interface[:-1] + self.x_vertical_interface[1:])
         self.y_horizontal_interface = np.arange(0, (ny + 0.5) * dy, dy) + X_0[1]
         self.y = 0.5 * (self.y_horizontal_interface[:-1] + self.y_horizontal_interface[1:])
         self.mass_cell = dx * dy
 
-        self.t = 0
+        self.t = self.t_0
         self.T_final = T_final
         self.dt = None
         self.t_array = None
@@ -118,7 +125,7 @@ class MeshRect2D:
         # To do:
         # add exceptions to make sure that all the inputs are float
         self.dt = np.min([1.0 / (1.2 * U.max_horizontal_U / self.dx + 1.2 * U.max_vertical_U / self.dy), dt_max])
-        self.t_array = np.arange(0, self.T_final + self.dt, self.dt)
+        self.t_array = np.arange(self.t_0, self.T_final + self.dt, self.dt)
         self.T_final = self.t_array[-1]
 
     def calc_dt_implicit_solver(self, dt):

src/pheromone_dispersion/velocity.py

@@ -17,6 +17,12 @@ class Velocity:
 
     Attributes
     ----------
+    t : ~numpy.ndarray
+        The array of times :math:`t` at which the velocity fields maps :math:`\vec{u}(x,y)` are given.
+    msh_dx : float
+        Space step of the mesh along the :math:`x`-axis :math:`\Delta x`.
+    msh_dy : float
+        Space step of the mesh along the :math:`y`-axis :math:`\Delta y`.
     at_vertical_interface: ~numpy.ndarray
         The wind velocity field :math:`\vec{u}(x,y)` at the vertical interfaces of each cells of the mesh at the current time :math:`t`.
     at_horizontal_interface: ~numpy.ndarray
@@ -25,8 +31,6 @@ class Velocity:
         The divergence of the wind velocity field
         :math:`div~\vec{u}(x,y)=\partial_xu(x,y)+\partial_yv(x,y)`
         at the center of each cells of the mesh at the current time :math:`t`.
-    t : ~numpy.ndarray
-        The array of times :math:`t` at which the velocity fields maps :math:`\vec{u}(x,y)` are given.
         If :math:`\vec{u}` is stationary, set to `None`.
     time_interpolation_at_vertical_interface : scipy.interpolate.interp1d
         Callable function used to compute the wind velocity field :math:`\vec{u}(x,y)`
@@ -36,9 +40,6 @@ class Velocity:
         Callable function used to compute the wind velocity field :math:`\vec{u}(x,y)`
         at the horizontal interfaces of each cells of the mesh given a time :math:`t`.
         Computation using the linear interpolation of the wind fields provided on the given time array.
-    time_interpolation_div : scipy.interpolate.interp1d
-        Callable function used to compute the divergence of the wind velocity field :math:`div~\vec{u}(x,y)` given a time :math:`t`.
-        Computation using the linear interpolation of the divergence of the wind fields originally computed on the given time array.
     cell_above_upwind : ~numpy.ndarray of boolean
         For each cell, is `True` if the cell above is the upwind cell
         for the :math:`y`-component of the wind velocity :math:`v`
@@ -86,75 +87,67 @@ class Velocity:
 
         """
 
+        _validate_inputs(msh, U_at_vertical_interface, U_at_horizontal_interface, t)
         self.t = t
-        if self.t is not None:
-            if U_at_horizontal_interface.ndim != 4 or U_at_vertical_interface.ndim != 4:
-                raise ValueError(
-                    "The number dimension of the velocity field is incorrect,"
-                    + " the velocity field is unsteady,"
-                    + " the number of dimension shoud be 3"
-                )
-            if self.t.size != U_at_vertical_interface.shape[0] or self.t.size != U_at_horizontal_interface.shape[0]:
-                raise ValueError("The shape of the velocity field does not coincide with the shape of the time vector")
-            div_U_all_t = ((U_at_horizontal_interface[:, 1:, :, 1] - U_at_horizontal_interface[:, :-1, :, 1]) / msh.dy) + (
-                (U_at_vertical_interface[:, :, 1:, 0] - U_at_vertical_interface[:, :, :-1, 0]) / msh.dx
-            )
-            self.at_vertical_interface = U_at_vertical_interface[0, :, :, :]
-            self.at_horizontal_interface = U_at_horizontal_interface[0, :, :, :]
-            self.div = div_U_all_t[0, :, :]
+        self.msh_dx = msh.dx
+        self.msh_dy = msh.dy
+
+        is_unsteady = self.t is not None
+
+        if is_unsteady:
             self.time_interpolation_at_vertical_interface = interp1d(t, U_at_vertical_interface, axis=0)
             self.time_interpolation_at_horizontal_interface = interp1d(t, U_at_horizontal_interface, axis=0)
-            self.time_interpolation_div = interp1d(t, div_U_all_t, axis=0)
-            self.max_horizontal_U = np.max(
-                (np.max(np.abs(U_at_vertical_interface[:, :, :, 0])), np.max(np.abs(U_at_horizontal_interface[:, :, :, 0])))
-            )
-            self.max_vertical_U = np.max(
-                (np.max(np.abs(U_at_vertical_interface[:, :, :, 1])), np.max(np.abs(U_at_horizontal_interface[:, :, :, 1])))
-            )
-
-        else:
-            if U_at_horizontal_interface.ndim != 3 or U_at_vertical_interface.ndim != 3:
-                raise ValueError(
-                    "The number dimension of the velocity field is incorrect,"
-                    + " the velocity field is unsteady,"
-                    + " the number of dimension shoud be 2"
-                )
-            self.at_vertical_interface = U_at_vertical_interface
-            self.at_horizontal_interface = U_at_horizontal_interface
-            self.div = ((self.at_horizontal_interface[1:, :, 1] - self.at_horizontal_interface[:-1, :, 1]) / msh.dy) + (
-                (self.at_vertical_interface[:, 1:, 0] - self.at_vertical_interface[:, :-1, 0]) / msh.dx
-            )
-            self.max_horizontal_U = np.max((np.max(U_at_vertical_interface[:, :, 0]), np.max(U_at_horizontal_interface[:, :, 0])))
-            self.max_vertical_U = np.max((np.max(U_at_vertical_interface[:, :, 1]), np.max(U_at_horizontal_interface[:, :, 1])))
 
+        self.at_vertical_interface = U_at_vertical_interface[(0,) if is_unsteady else slice(None)]
+        self.at_horizontal_interface = U_at_horizontal_interface[(0,) if is_unsteady else slice(None)]
+        self.compute_divergence()
+        self.compute_upwind_cells()
+        self.compute_max_velocity(U_at_vertical_interface, U_at_horizontal_interface)
+
+    def compute_divergence(self):
+        r"""
+        Compute the divergence of the velocity field at the current time
+        and store it in the attribute :attr:`div`.
+        """
+        dudx = (self.at_vertical_interface[:, 1:, 0] - self.at_vertical_interface[:, :-1, 0]) / self.msh_dx
+        dvdy = (self.at_horizontal_interface[1:, :, 1] - self.at_horizontal_interface[:-1, :, 1]) / self.msh_dy
+        self.div = dudx + dvdy
+
+    def compute_upwind_cells(self):
+        r"""
+        Compute the boolean arrays to determine which cells are upwind cells
+        and store it in the attributes :attr:`cell_above_upwind`, :attr:`cell_under_upwind`,
+        :attr:`cell_right_upwind` and :attr:`cell_left_upwind`.
+        """
         self.cell_above_upwind = self.at_horizontal_interface[:, :, 1] < 0
         self.cell_under_upwind = self.at_horizontal_interface[:, :, 1] > 0
         self.cell_right_upwind = self.at_vertical_interface[:, :, 0] < 0
         self.cell_left_upwind = self.at_vertical_interface[:, :, 0] >= 0
 
-        if self.at_vertical_interface.shape != (msh.y.size, msh.x_vertical_interface.size, 2):
-            raise ValueError(
-                "The shape of the velocity field at the vertical and horizontal interfaces do not match with the shape of the msh."
-            )
-        if (
-            self.at_vertical_interface.shape[0] + 1 != self.at_horizontal_interface.shape[0]
-            or self.at_horizontal_interface.shape[1] + 1 != self.at_vertical_interface.shape[1]
-        ):
-            raise ValueError(
-                "The shape of the velocity field at the vertical interfaces"
-                + "(that is (nb of cells along the y axis, nb of cells along the x axis + 1)) "
-                + "does not coincide with the shape of the velocity field at the horizontal interfaces"
-                + "(that is (nb of cells along the y axis + 1, nb of cells along the x axis))."
-            )
+    def compute_max_velocity(self, U_v, U_h):
+        r"""
+        Compute the maximum of each components of the velocity field
+        and store it in the attributes :attr:`max_horizontal_U` and :attr:`max_vertical_U`.
+
+        Parameters
+        ----------
+        U_v : ~numpy.ndarray
+            The wind velocity fields :math:`\vec{u}(x,y)` at one or multiple times at the vertical interfaces of each cells of the mesh.
+        U_h : ~numpy.ndarray
+            The wind velocity fields :math:`\vec{u}(x,y)` at one or multiple times at the horizontal interfaces of each cells of the mesh.
+        """
+        self.max_horizontal_U = np.max((np.max(np.abs(U_v[..., 0])), np.max(np.abs(U_h[..., 0]))))
+        self.max_vertical_U = np.max((np.max(np.abs(U_v[..., 1])), np.max(np.abs(U_h[..., 1]))))
 
     def at_current_time(self, tc):
         r"""
-        Update the attributes :attr:`at_vertical_interface`,
-        :attr:`at_horizontal_interface` and :attr:`div` at a given time
+        Update the attributes :attr:`at_vertical_interface` and
+        :attr:`at_horizontal_interface` at a given time
         using the linear interpolation callable attributes resp
-        :attr:`time_interpolation_at_vertical_interface`,
-        :attr:`time_interpolation_at_horizontal_interface` and
-        :attr:`time_interpolation_div`.
+        :attr:`time_interpolation_at_vertical_interface` and
+        :attr:`time_interpolation_at_horizontal_interface`.
+        Update the attributes :attr:`div`, :attr:`cell_above_upwind`, :attr:`cell_under_upwind`,
+        :attr:`cell_right_upwind` and :attr:`cell_left_upwind` using the methods to compute them.
 
         Parameters
         ----------
@@ -169,15 +162,65 @@ class Velocity:
         """
 
         if self.t is not None:
-            if tc < min(self.t) or tc > max(self.t):
+            if not (min(self.t) <= tc <= max(self.t)):
                 raise ValueError("The given time must be between the lowest and largest times contained in the time vector.")
             self.at_vertical_interface = self.time_interpolation_at_vertical_interface(tc)
             self.at_horizontal_interface = self.time_interpolation_at_horizontal_interface(tc)
-            self.div = self.time_interpolation_div(tc)
-            self.cell_above_upwind = self.at_horizontal_interface[:, :, 1] < 0
-            self.cell_under_upwind = self.at_horizontal_interface[:, :, 1] >= 0
-            self.cell_right_upwind = self.at_vertical_interface[:, :, 0] < 0
-            self.cell_left_upwind = self.at_vertical_interface[:, :, 0] >= 0
+            self.compute_divergence()
+            self.compute_upwind_cells()
+
+
+def _validate_inputs(msh, U_vi, U_hi, t_U):
+    # Validate the inputs of the Velocity class and ensure correct shapes.
+
+    expected_dims = {True: 4, False: 3}  # True if unsteady, False if stationary
+    is_unsteady = t_U is not None
+
+    # Check that the velocity field has the expected dimension
+    if U_vi.ndim != expected_dims[is_unsteady] or U_hi.ndim != expected_dims[is_unsteady]:
+        raise ValueError(f"Velocity field should be {expected_dims[is_unsteady]}D when {'unsteady' if is_unsteady else 'stationary'}.")
+
+    # Check that the shape of the velocity field matches the size of the time vector
+    if is_unsteady and any(t_U.size != arr.shape[0] for arr in (U_vi, U_hi)):
+        raise ValueError("The shape of the velocity field does not match the size of the time vector.")
+
+    dim_y, dim_x = (1, 2) if is_unsteady else (0, 1)
+    # Check that the shape of the velocity field at the horizontal and vertical interfaces match
+    if U_vi.shape[dim_y] + 1 != U_hi.shape[dim_y] or U_hi.shape[dim_x] + 1 != U_vi.shape[dim_x]:
+        raise ValueError(
+            f"Mismatch between velocity field shapes:\n"
+            f"  - Vertical interfaces: expected ({U_hi.shape[dim_y] - 1}, {U_vi.shape[dim_x]}) but got {U_vi.shape[dim_y:]}\n"
+            f"  - Horizontal interfaces: expected ({U_vi.shape[dim_y]}, {U_hi.shape[dim_x] - 1}) but got {U_hi.shape[dim_y:]}"
+        )
+
+    # Check that the shape of the shape velocity field match with the mesh
+    U_vi_slice, U_hi_slice = U_vi[(0,) if is_unsteady else slice(None)], U_hi[(0,) if is_unsteady else slice(None)]
+    _validate_velocity_field_shape_wrt_mesh(msh, t_U, U_vi_slice, U_hi_slice)
+
+
+def _validate_velocity_field_shape_wrt_mesh(msh, t_U, U_vi, U_hi):
+    # Check that the shape of the velocity field at a given time match with the mesh.
+
+    if U_vi.shape != (msh.y.size, msh.x_vertical_interface.size, 2):
+        raise ValueError(
+            "The shape of the velocity field at the vertical and horizontal interfaces do not match with the shape of the msh."
+            f" Expected respectively {msh.y.size, msh.x_vertical_interface.size, 2} and {msh.y_horizontal_interface.size, msh.x.size, 2}"
+            f" but got respectively {U_vi.shape} and {U_hi.shape}"
+        )
+    if t_U is not None:
+        if np.min(t_U) > msh.t_0:
+            raise ValueError("The given velocity field lacks of coverage at the beginning of the time window of msh.")
+        if np.max(t_U) < msh.T_final:
+            raise ValueError("The given velocity field lacks of coverage at the end of the time window of msh.")
+
+
+def _validate_velocity_field(U, msh):
+    # Check that the velocity field is of class Velocity
+    # and that the shape of the velocity field match with the mesh.
+
+    if not isinstance(U, Velocity):
+        raise TypeError('The provided velocity field is not of type Velocity.')
+    _validate_velocity_field_shape_wrt_mesh(msh, U.t, U.at_vertical_interface, U.at_vertical_interface)
 
 
 def velocity_field_from_meteo_data(path_data, file_name_data, msh):

tests/test_advection_term.py

 from pathlib import Path
 
 import numpy as np
+import pytest
 from pytest_bdd import given
 from pytest_bdd import scenario
 from pytest_bdd import then
@@ -102,3 +103,11 @@ def the_linear_operator_is_updated_as_expected_at_a_given_time(advection):
     assert (advection.minus_U.at_vertical_interface[:, :, 0] == -1.5).all() and (
         advection.minus_U.at_vertical_interface[:, :, 1] == 0.5
     ).all()
+
+
+@then("the initialization fails if the velocity fields is not an object of the class Velocity")
+def the_initialization_fails_if_the_velocity_fields_is_not_an_object_of_the_class_Velocity(msh):
+    """the initialization fails if the velocity fields is not an object of the class Velocity"""
+    with pytest.raises(TypeError) as e:
+        Advection(3.0, msh)
+    assert e.type == TypeError

tests/test_steady_velocity_field.py

 from pathlib import Path
 
 import numpy as np
-import pytest
 from pytest_bdd import given
 from pytest_bdd import scenario
 from pytest_bdd import then
@@ -119,25 +118,6 @@ def the_maximum_of_the_vertical_velocity_has_the_expected_value(U):
     assert U.max_vertical_U == 1
 
 
-@then("the initialization fails if the number of dimension of the velocity field is not correct")
-def the_initialization_fails_if_the_number_of_dimension_of_the_velocity_field_is_not_correct(U_vi, U_hi, msh):
-    """the initialization fails if the number of dimension of the velocity is not correct"""
-    with pytest.raises(ValueError) as e:
-        Velocity(msh, U_vi[:, :, 0], U_hi)
-    assert e.type == ValueError
-    with pytest.raises(ValueError) as e:
-        Velocity(msh, U_vi, U_hi[:, :, 0])
-    assert e.type == ValueError
-
-
-@then("the initialization fails if the velocity fields shape do not match")
-def the_initialization_fails_if_the_velocity_fields_shape_do_not_match(msh, U_vi, U_hi):
-    """the initialization fails if the velocity fields shape do not match"""
-    with pytest.raises(ValueError) as e:
-        Velocity(msh, U_hi, U_vi)
-    assert e.type == ValueError
-
-
 @then("the update at a given time does not change the values")
 def the_update_at_a_given_time_does_not_change_the_values(U):
     """the update at a given time does not change the values"""

tests/test_velocity_field.py

 from pathlib import Path
 
 import numpy as np
-import pytest
 from pytest_bdd import given
 from pytest_bdd import scenario
 from pytest_bdd import then
@@ -132,47 +131,3 @@ def the_maximum_of_the_horizontal_velocity_has_the_expected_value(U):
 def the_maximum_of_the_vertical_velocity_has_the_expected_value(U):
     """the maximum of the vertical velocity has the expected value"""
     assert U.max_vertical_U == 2
-
-
-@then("the initialization fails if the number of dimension of the velocity field is not correct")
-def the_initialization_fails_if_the_number_of_dimension_of_the_velocity_field_is_not_correct(U_vi, U_hi, t, msh):
-    """the initialization fails if the number of dimension of the velocity is not correct"""
-    with pytest.raises(ValueError) as e:
-        Velocity(msh, U_vi[:, :, :, 0], U_hi, t=t)
-    assert e.type == ValueError
-    with pytest.raises(ValueError) as e:
-        Velocity(msh, U_vi, U_hi[:, :, :, 0], t=t)
-    assert e.type == ValueError
-
-
-@then("the initialization fails if the velocity fields shape do not match the shape of the mesh")
-def the_initialization_fails_if_the_velocity_fields_shape_do_not_match_the_shape_of_the_mesh(msh, U_hi_invalid, U_vi_invalid):
-    """the initialization fails if the velocity fields shape do not match the shape of the mesh"""
-    with pytest.raises(ValueError) as e:
-        Velocity(msh, U_vi_invalid, U_hi_invalid)
-    assert e.type == ValueError
-
-
-@then("the initialization fails if the velocity fields shape do not match at vertical and horizontal interfaces")
-def the_initialization_fails_if_the_velocity_fields_shape_do_not_match_at_vertical_and_horizontal_interfaces(msh, t, U_vi, U_hi):
-    """the initialization fails if the velocity fields shape do not match at vertical and horizontal interfaces"""
-    with pytest.raises(ValueError) as e:
-        Velocity(msh, U_hi, U_vi, t=t)
-    assert e.type == ValueError
-
-
-@then("the initialization fails if the velocity fields shape do not match with the shape of the time vector")
-def the_initialization_fails_if_the_velocity_fields_shape_do_not_match_with_the_shape_of_the_time_vector(msh, U_vi, U_hi):
-    """the initialization fails if the velocity fields shape do not match with the shape of the time vector"""
-    with pytest.raises(ValueError) as e:
-        t_fail = np.array([0])
-        Velocity(msh, U_vi, U_hi, t=t_fail)
-    assert e.type == ValueError
-
-
-@then("the update fails if the time is not between the lowest and largest times contained in the time vector")
-def the_update_fails_if_the_time_is_not_between_the_lowest_and_largest_times_contained_in_the_time_vector(U):
-    """the update fails if the time is not between the lowest and largest times contained in the time vector"""
-    with pytest.raises(ValueError) as e:
-        U.at_current_time(max(U.t) + 1)
-    assert e.type == ValueError

tests/test_velocity_field_inputs.py

+from pathlib import Path
+
+import numpy as np
+import pytest
+from pytest_bdd import given
+from pytest_bdd import scenario
+from pytest_bdd import then
+from pytest_bdd import when
+
+from pheromone_dispersion.geom import MeshRect2D
+from pheromone_dispersion.velocity import Velocity
+
+
+@scenario("velocity_field_inputs.feature", "Test the validity check of the inputs of the Velocity object")
+def test_the_validity_check_of_the_inputs_of_the_Velocity_object():
+    """Test the validity check of the inputs of the Velocity object."""
+
+
+@given("a rectangular 2D mesh", target_fixture="msh")
+def a_rectangular_2D_mesh():
+    "a rectangular 2D mesh"
+    Lx = np.load(Path("tests/test_data") / "Lx.npy")
+    Ly = np.load(Path("tests/test_data") / "Ly.npy")
+    dx = np.load(Path("tests/test_data") / "dx.npy")
+    dy = np.load(Path("tests/test_data") / "dy.npy")
+    Tfinal = np.load(Path("tests/test_data") / "Tfinal.npy")
+    return MeshRect2D(Lx, Ly, dx, dy, Tfinal)
+
+
+@given("a time vector", target_fixture='t')
+def a_time_vector():
+    "a time vector"
+    t = np.load(Path("tests/test_data") / "t.npy")
+    return t
+
+
+@given("a velocity field at the vertical interfaces", target_fixture='U_vi')
+def a_velocity_field_at_the_vertical_interfaces():
+    "a velocity field at the vertical interfaces"
+    U_vi = np.load(Path("tests/test_data") / "U_vi.npy")
+    return U_vi
+
+
+@given("a velocity field at the horizontal interfaces", target_fixture='U_hi')
+def a_velocity_field_at_the_horizontal_interfaces():
+    "a velocity field at the horizontal interfaces"
+    U_hi = np.load(Path("tests/test_data") / "U_hi.npy")
+    return U_hi
+
+
+@given("an invalid velocity field at the vertical interfaces", target_fixture='U_vi_invalid')
+def an_invalid_velocity_field_at_the_vertical_interfaces():
+    "an invalid velocity field at the vertical interfaces"
+    U_vi_invalid = np.load(Path("tests/test_data") / "U_vi_invalid.npy")
+    return U_vi_invalid
+
+
+@given("an invalid velocity field at the horizontal interfaces", target_fixture='U_hi_invalid')
+def an_invalid_velocity_field_at_the_horizontal_interfaces():
+    "an invalid velocity field at the horizontal interfaces"
+    U_hi_invalid = np.load(Path("tests/test_data") / "U_hi_invalid.npy")
+    return U_hi_invalid
+
+
+@when("initialize the velocity object", target_fixture='U')
+def initialize_the_velocity_object(msh, t, U_vi, U_hi):
+    """initialize the velocity object"""
+    return Velocity(msh, U_vi, U_hi, t=t)
+
+
+@then("the initialization fails if the number of dimension of the velocity field is not correct")
+def the_initialization_fails_if_the_number_of_dimension_of_the_velocity_field_is_not_correct(U_vi, U_hi, t, msh):
+    """the initialization fails if the number of dimension of the velocity is not correct"""
+    with pytest.raises(ValueError) as e:
+        Velocity(msh, U_vi[:, :, :, 0], U_hi, t=t)
+    assert e.type == ValueError
+    with pytest.raises(ValueError) as e:
+        Velocity(msh, U_vi, U_hi[:, :, :, 0], t=t)
+    assert e.type == ValueError
+    with pytest.raises(ValueError) as e:
+        Velocity(msh, U_vi[:, :, :, 0], U_hi)
+    assert e.type == ValueError
+    with pytest.raises(ValueError) as e:
+        Velocity(msh, U_vi, U_hi[:, :, :, 0])
+    assert e.type == ValueError
+
+
+@then("the initialization fails if the velocity fields shape do not match the shape of the mesh")
+def the_initialization_fails_if_the_velocity_fields_shape_do_not_match_the_shape_of_the_mesh(msh, U_hi_invalid, U_vi_invalid):
+    """the initialization fails if the velocity fields shape do not match the shape of the mesh"""
+    with pytest.raises(ValueError) as e:
+        Velocity(msh, U_vi_invalid, U_hi_invalid)
+    assert e.type == ValueError
+
+
+@then("the initialization fails if the velocity fields shape do not match at vertical and horizontal interfaces")
+def the_initialization_fails_if_the_velocity_fields_shape_do_not_match_at_vertical_and_horizontal_interfaces(msh, t, U_vi, U_hi):
+    """the initialization fails if the velocity fields shape do not match at vertical and horizontal interfaces"""
+    with pytest.raises(ValueError) as e:
+        Velocity(msh, U_hi, U_vi, t=t)
+    assert e.type == ValueError
+    with pytest.raises(ValueError) as e:
+        Velocity(msh, U_hi[0, :, :, :], U_vi[0, :, :, :])
+    assert e.type == ValueError
+
+
+@then("the initialization fails if the velocity fields shape do not match with the shape of the time vector")
+def the_initialization_fails_if_the_velocity_fields_shape_do_not_match_with_the_shape_of_the_time_vector(msh, U_vi, U_hi):
+    """the initialization fails if the velocity fields shape do not match with the shape of the time vector"""
+    with pytest.raises(ValueError) as e:
+        t_fail = np.array([0])
+        Velocity(msh, U_vi, U_hi, t=t_fail)
+    assert e.type == ValueError
+
+
+@then("the initialization fails if the velocity fields do not cover the entire time window")
+def the_initialization_fails_if_the_velocity_fields_do_not_cover_the_entire_time_window(msh, t, U_vi, U_hi):
+    """the initialization fails if the velocity fields do not cover the entire time window"""
+    with pytest.raises(ValueError) as e:
+        t_fail = t.astype(np.float32)
+        t_fail[0] = msh.t_0 + 0.1
+        Velocity(msh, U_vi, U_hi, t=t_fail)
+    assert e.type == ValueError
+    with pytest.raises(ValueError) as e:
+        t_fail = np.copy(t)
+        t_fail[-1] = msh.T_final - 0.1
+        Velocity(msh, U_vi, U_hi, t=t_fail)
+    assert e.type == ValueError
+
+
+@then("the update fails if the time is not between the lowest and largest times contained in the time vector")
+def the_update_fails_if_the_time_is_not_between_the_lowest_and_largest_times_contained_in_the_time_vector(U):
+    """the update fails if the time is not between the lowest and largest times contained in the time vector"""
+    with pytest.raises(ValueError) as e:
+        U.at_current_time(max(U.t) + 1)
+    assert e.type == ValueError