From ec836c728fdddfa4176f5e22467033c144e6ab62 Mon Sep 17 00:00:00 2001 From: Laurent Courty Date: Thu, 25 Jun 2026 13:56:39 -0600 Subject: [PATCH] add tests for non-square cells --- tests/core/test_5by5.py | 98 +++++++++++++++++++++++++++++++++++++++++ tests/core/test_flow.py | 69 +++++++++++++++++++++++++++++ 2 files changed, 167 insertions(+) diff --git a/tests/core/test_5by5.py b/tests/core/test_5by5.py index 6b92950..3ea5d23 100644 --- a/tests/core/test_5by5.py +++ b/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) diff --git a/tests/core/test_flow.py b/tests/core/test_flow.py index 99ff491..e27ac37 100644 --- a/tests/core/test_flow.py +++ b/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"""