Add tests for non-square cells

tests/core/test_5by5.py

@@ -9,12 +9,75 @@
 from itzi.simulation_builder import SimulationBuilder
 from itzi.data_containers import SimulationConfig, SurfaceFlowParameters
 from itzi.providers.memory_output import MemoryRasterOutputProvider, MemoryVectorOutputProvider
+from itzi.providers.domain_data import DomainData
 from itzi.const import InfiltrationModelType, TemporalType
 
 if TYPE_CHECKING:
     from itzi.simulation import Simulation
 
 
+def _run_center_pulse_simulation(
+    domain_5by5,
+    helpers,
+    *,
+    dx: float,
+    dy: float,
+    duration_s: float = 1.0,
+) -> np.ndarray:
+    rows, cols = domain_5by5.domain_data.shape
+    domain_data = DomainData(
+        north=rows * dy,
+        south=0.0,
+        east=cols * dx,
+        west=0.0,
+        rows=rows,
+        cols=cols,
+        crs_wkt="",
+    )
+
+    initial_volume = float(np.sum(domain_5by5.arr_start_h) * domain_5by5.domain_data.cell_area)
+    arr_start_h = np.zeros(domain_data.shape, dtype=np.float32)
+    arr_start_h[2, 2] = initial_volume / domain_data.cell_area
+
+    sim_config = SimulationConfig(
+        start_time=datetime(2000, 1, 1, 0, 0, 0),
+        end_time=datetime(2000, 1, 1, 0, 0, 0) + timedelta(seconds=duration_s),
+        record_step=timedelta(seconds=duration_s),
+        temporal_type=TemporalType.RELATIVE,
+        input_map_names=helpers.make_input_map_names(
+            dem="z",
+            friction="n",
+            water_depth="start_h",
+        ),
+        output_map_names=helpers.make_output_map_names(
+            f"out_5by5_rect_{int(dx)}x{int(dy)}",
+            ["water_depth"],
+        ),
+        surface_flow_parameters=SurfaceFlowParameters(hmin=0.0001, dtmax=0.3, cfl=0.2),
+        infiltration_model=InfiltrationModelType.NULL,
+    )
+
+    raster_output = MemoryRasterOutputProvider({"out_map_names": sim_config.output_map_names})
+    simulation = (
+        SimulationBuilder(sim_config, domain_5by5.arr_mask, np.float32)
+        .with_domain_data(domain_data)
+        .with_raster_output_provider(raster_output)
+        .with_vector_output_provider(MemoryVectorOutputProvider({}))
+        .build()
+    )
+
+    simulation.set_array("dem", domain_5by5.arr_dem_flat.copy())
+    simulation.set_array("friction", domain_5by5.arr_n.copy())
+    simulation.set_array("water_depth", arr_start_h.copy())
+
+    simulation.initialize()
+    while simulation.sim_time < simulation.end_time:
+        simulation.update()
+    final_depth = simulation.get_array("water_depth").copy()
+    simulation.finalize()
+    return final_depth
+
+
 @pytest.fixture(scope="module")
 def sim_5by5(domain_5by5, helpers) -> Simulation:
     """Run a 5x5 simulation for 60s with 30s record step.
@@ -157,3 +220,38 @@ def test_flow_symmetry(self, sim_5by5):
             f"above={h_above:.6f}, below={h_below:.6f}, "
             f"left={h_left:.6f}, right={h_right:.6f}"
         )
+
+
+@pytest.mark.parametrize(
+    ("dx", "dy", "wetter_axis"),
+    [
+        (20.0, 10.0, "y"),
+        (10.0, 20.0, "x"),
+    ],
+)
+def test_rectangular_cells_preserve_axis_symmetry_and_bias_early_spreading(
+    domain_5by5, helpers, dx: float, dy: float, wetter_axis: str
+):
+    """Early-time spreading should stay symmetric within each axis and favor the shorter cell size."""
+    h_array = _run_center_pulse_simulation(domain_5by5, helpers, dx=dx, dy=dy)
+
+    h_above = h_array[1, 2]
+    h_below = h_array[3, 2]
+    h_left = h_array[2, 1]
+    h_right = h_array[2, 3]
+
+    assert np.isclose(h_above, h_below)
+    assert np.isclose(h_left, h_right)
+
+    if wetter_axis == "y":
+        assert min(h_above, h_below) > max(h_left, h_right)
+    else:
+        assert min(h_left, h_right) > max(h_above, h_below)
+
+
+def test_swapping_dx_and_dy_transposes_the_early_time_solution(domain_5by5, helpers):
+    """Swapping the cell sizes should swap the x/y spreading pattern."""
+    h_dx20_dy10 = _run_center_pulse_simulation(domain_5by5, helpers, dx=20.0, dy=10.0)
+    h_dx10_dy20 = _run_center_pulse_simulation(domain_5by5, helpers, dx=10.0, dy=20.0)
+
+    assert np.allclose(h_dx20_dy10, h_dx10_dy20.T)

tests/core/test_flow.py

@@ -54,6 +54,75 @@ def test_vectorizable_velocity_calculation():
         assert v_optimized == pytest.approx(expected_v)
 
 
+def test_solve_h_uses_dx_and_dy_separately_in_flow_divergence():
+    """The water-depth update must use the x and y cell sizes independently."""
+    shape = (5, 5)
+    dtype = np.float64
+
+    arr_ext = np.zeros(shape, dtype=dtype)
+    arr_qe = np.zeros(shape, dtype=dtype)
+    arr_qs = np.zeros(shape, dtype=dtype)
+    arr_bct = np.zeros(shape, dtype=dtype)
+    arr_bcv = np.zeros(shape, dtype=dtype)
+    arr_h = np.zeros(shape, dtype=dtype)
+    arr_hmax = np.zeros(shape, dtype=dtype)
+    arr_hfix = np.zeros(shape, dtype=dtype)
+    arr_herr = np.zeros(shape, dtype=dtype)
+    arr_hfe = np.zeros(shape, dtype=dtype)
+    arr_hfs = np.zeros(shape, dtype=dtype)
+    arr_v = np.zeros(shape, dtype=dtype)
+    arr_vdir = np.zeros(shape, dtype=dtype)
+    arr_vmax = np.zeros(shape, dtype=dtype)
+    arr_fr = np.zeros(shape, dtype=dtype)
+
+    center = (2, 2)
+    arr_h[center] = 1.0
+    arr_hmax[center] = 1.0
+
+    qw = 7.0
+    qe = 1.0
+    qn = 6.0
+    qs = 2.0
+    dx = 4.0
+    dy = 2.0
+    dt = 0.5
+    g = 9.81
+
+    arr_qe[2, 1] = qw
+    arr_qe[2, 2] = qe
+    arr_qs[1, 2] = qn
+    arr_qs[2, 2] = qs
+
+    flow.solve_h(
+        arr_ext=arr_ext,
+        arr_qe=arr_qe,
+        arr_qs=arr_qs,
+        arr_bct=arr_bct,
+        arr_bcv=arr_bcv,
+        arr_h=arr_h,
+        arr_hmax=arr_hmax,
+        arr_hfix=arr_hfix,
+        arr_herr=arr_herr,
+        arr_hfe=arr_hfe,
+        arr_hfs=arr_hfs,
+        arr_v=arr_v,
+        arr_vdir=arr_vdir,
+        arr_vmax=arr_vmax,
+        arr_fr=arr_fr,
+        dx=dx,
+        dy=dy,
+        dt=dt,
+        g=g,
+    )
+
+    expected_h = 1.0 + dt * (((qw - qe) / dx) + ((qn - qs) / dy))
+    swapped_dims_h = 1.0 + dt * (((qw - qe) / dy) + ((qn - qs) / dx))
+
+    assert arr_h[center] == pytest.approx(expected_h)
+    assert arr_hmax[center] == pytest.approx(expected_h)
+    assert not np.isclose(arr_h[center], swapped_dims_h)
+
+
 class TestWaterDepthFunction:
     """Integration tests for flow functions with optimized calculations"""