From c088bf64234eb41b75965129def8646f9d0a98f1 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Mon, 22 Apr 2024 11:32:35 +0200 Subject: [PATCH 1/6] Compute crosscut thickness, bc_top and bc_bottom optional --- imod/prepare/topsystem/conductance.py | 34 ++++++++++++++------------- 1 file changed, 18 insertions(+), 16 deletions(-) diff --git a/imod/prepare/topsystem/conductance.py b/imod/prepare/topsystem/conductance.py index 1b3a408bc..ed45cc324 100644 --- a/imod/prepare/topsystem/conductance.py +++ b/imod/prepare/topsystem/conductance.py @@ -307,28 +307,30 @@ def _compute_layer_thickness(allocated, top, bottom): @preserve_gridtype -def _compute_crosscut_thickness(allocated, top, bottom, stage, bottom_elevation): +def _compute_crosscut_thickness(allocated, top, bottom, bc_top=None, bc_bottom=None): """ - Compute 3D grid of thicknesses crosscut by river in allocated cells. So the - upper allocated layer thickness is stage - bottom, the lower allocated layer - is top - river_bottom_elevation. + Compute 3D grid of thicknesses crosscut by boundary condition (river/drain) + in allocated cells. So the thickness in the upper allocated layer thickness + is stage - bottom, the lower allocated layer is top - + river_bottom_elevation. """ - top_layered = _enforce_layered_top(top, bottom) + if (bc_top is None) & (bc_bottom is None): + raise ValueError("`bc_top` and `bc_bottom` cannot both be None.") + top_layered = _enforce_layered_top(top, bottom) thickness = _compute_layer_thickness(allocated, top, bottom) + outside = ones_like(allocated).astype(bool) - upper_layer_allocated = (stage < top_layered) & (stage > bottom) - lower_layer_allocated = (bottom_elevation < top_layered) & ( - bottom_elevation > bottom - ) - outside = (stage < bottom) | (bottom_elevation > top_layered) + if bc_top is not None: + upper_layer_bc = (bc_top < top_layered) & (bc_top > bottom) + outside = outside | (bc_top < bottom) + thickness = thickness.where(~upper_layer_bc, thickness - (top_layered - bc_top)) + + if bc_bottom is not None: + lower_layer_bc = (bc_bottom < top_layered) & (bc_bottom > bottom) + outside = outside | (bc_bottom > top_layered) + thickness = thickness.where(~lower_layer_bc, thickness - (bc_bottom - bottom)) - thickness = thickness.where( - ~upper_layer_allocated, thickness - (top_layered - stage) - ) - thickness = thickness.where( - ~lower_layer_allocated, thickness - (bottom_elevation - bottom) - ) thickness = thickness.where(~outside, 0.0) return thickness From 14554d75996cf6cb0143793fc5dea1f58003b706 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Mon, 22 Apr 2024 13:10:13 +0200 Subject: [PATCH 2/6] Reorder arguments and add crosscut options for drain package. --- imod/prepare/topsystem/conductance.py | 161 +++++++++++++++----------- 1 file changed, 96 insertions(+), 65 deletions(-) diff --git a/imod/prepare/topsystem/conductance.py b/imod/prepare/topsystem/conductance.py index ed45cc324..1014bf1a8 100644 --- a/imod/prepare/topsystem/conductance.py +++ b/imod/prepare/topsystem/conductance.py @@ -1,4 +1,5 @@ from enum import Enum +from typing import Optional import numpy as np @@ -67,9 +68,9 @@ def distribute_riv_conductance( conductance: GridDataArray, top: GridDataArray, bottom: GridDataArray, + k: GridDataArray, stage: GridDataArray, bottom_elevation: GridDataArray, - k: GridDataArray, ) -> GridDataArray: """ Function to distribute 2D conductance over vertical layers for the RIV @@ -90,14 +91,14 @@ def distribute_riv_conductance( Model top bottom: DataArray | UgridDataArray Model layer bottoms + k: DataArray | UgridDataArray + Hydraulic conductivities stage: DataArray | UgridDataArray Planar grid with river stages, cannot contain a layer dimension. Can contain a time dimension. bottom_elevation: DataArray | UgridDataArray Planar grid with river bottom elevations, cannot contain a layer dimension. Can contain a time dimension. - k: DataArray | UgridDataArray - Hydraulic conductivities Returns ------- @@ -117,39 +118,39 @@ def distribute_riv_conductance( PLANAR_GRID.validate(conductance) match distributing_option: - case DISTRIBUTING_OPTION.by_corrected_transmissivity: - weights = _distribute_weights__by_corrected_transmissivity( - allocated, top, bottom, stage, bottom_elevation, k - ) case DISTRIBUTING_OPTION.equally: weights = _distribute_weights__equally(allocated) - case DISTRIBUTING_OPTION.by_crosscut_thickness: - weights = _distribute_weights__by_crosscut_thickness( - allocated, top, bottom, stage, bottom_elevation - ) case DISTRIBUTING_OPTION.by_layer_thickness: weights = _distribute_weights__by_layer_thickness(allocated, top, bottom) - case DISTRIBUTING_OPTION.by_crosscut_transmissivity: - weights = _distribute_weights__by_crosscut_transmissivity( - allocated, top, bottom, stage, bottom_elevation, k - ) case DISTRIBUTING_OPTION.by_layer_transmissivity: weights = _distribute_weights__by_layer_transmissivity( allocated, top, bottom, k ) case DISTRIBUTING_OPTION.by_conductivity: weights = _distribute_weights__by_conductivity(allocated, k) + case DISTRIBUTING_OPTION.by_crosscut_thickness: + weights = _distribute_weights__by_crosscut_thickness( + allocated, top, bottom, stage, bottom_elevation + ) + case DISTRIBUTING_OPTION.by_crosscut_transmissivity: + weights = _distribute_weights__by_crosscut_transmissivity( + allocated, top, bottom, k, stage, bottom_elevation + ) + case DISTRIBUTING_OPTION.by_corrected_transmissivity: + weights = _distribute_weights__by_corrected_transmissivity( + allocated, top, bottom, k, stage, bottom_elevation + ) case _: raise ValueError( "Received incompatible setting for rivers, only" - f"'{DISTRIBUTING_OPTION.by_corrected_transmissivity.name}' and" - f"'{DISTRIBUTING_OPTION.equally.name}' and" - f"'{DISTRIBUTING_OPTION.by_crosscut_thickness.name}' and" - f"'{DISTRIBUTING_OPTION.by_layer_thickness.name}' and" - f"'{DISTRIBUTING_OPTION.by_crosscut_transmissivity.name}' and" - f"'{DISTRIBUTING_OPTION.by_layer_transmissivity.name}' and" - f"'{DISTRIBUTING_OPTION.by_conductivity.name}' supported." - f"got: '{distributing_option.name}'" + f"'{DISTRIBUTING_OPTION.equally.name}', " + f"'{DISTRIBUTING_OPTION.by_layer_thickness.name}', " + f"'{DISTRIBUTING_OPTION.by_layer_transmissivity.name}', " + f"'{DISTRIBUTING_OPTION.by_conductivity.name}', " + f"'{DISTRIBUTING_OPTION.by_crosscut_thickness.name}', " + f"'{DISTRIBUTING_OPTION.by_crosscut_transmissivity.name}', and " + f"'{DISTRIBUTING_OPTION.by_corrected_transmissivity.name}' supported. " + f"Got: '{distributing_option.name}'" ) return (weights * conductance).where(allocated) @@ -162,6 +163,7 @@ def distribute_drn_conductance( top: GridDataArray, bottom: GridDataArray, k: GridDataArray, + elevation: GridDataArray, ) -> GridDataArray: """ Function to distribute 2D conductance over vertical layers for the DRN @@ -184,7 +186,9 @@ def distribute_drn_conductance( Model layer bottoms k: DataArray | UgridDataArray Hydraulic conductivities - + elevation: DataArray | UgridDataArray + Drain elevation + Returns ------- Conductances distributed over depth. @@ -197,7 +201,7 @@ def distribute_drn_conductance( ) >>> conductances_distributed = distribute_drn_conductance( DISTRIBUTING_OPTION.by_layer_transmissivity, allocated, - conductance, top, bottom, k + conductance, top, bottom, k, drain_elevation ) """ PLANAR_GRID.validate(conductance) @@ -213,14 +217,29 @@ def distribute_drn_conductance( ) case DISTRIBUTING_OPTION.by_conductivity: weights = _distribute_weights__by_conductivity(allocated, k) + case DISTRIBUTING_OPTION.by_crosscut_thickness: + weights = _distribute_weights__by_crosscut_thickness( + allocated, top, bottom, bc_top=elevation + ) + case DISTRIBUTING_OPTION.by_crosscut_transmissivity: + weights = _distribute_weights__by_crosscut_transmissivity( + allocated, top, bottom, k, bc_top=elevation + ) + case DISTRIBUTING_OPTION.by_corrected_transmissivity: + weights = _distribute_weights__by_corrected_transmissivity( + allocated, top, bottom, k, bc_top=elevation + ) case _: raise ValueError( "Received incompatible setting for drains, only" - f"'{DISTRIBUTING_OPTION.equally.name}' and" - f"'{DISTRIBUTING_OPTION.by_layer_thickness.name}' and" - f"'{DISTRIBUTING_OPTION.by_layer_transmissivity.name}' and" - f"'{DISTRIBUTING_OPTION.by_conductivity.name}' supported." - f"got: '{distributing_option.name}'" + f"'{DISTRIBUTING_OPTION.equally.name}', " + f"'{DISTRIBUTING_OPTION.by_layer_thickness.name}', " + f"'{DISTRIBUTING_OPTION.by_layer_transmissivity.name}', " + f"'{DISTRIBUTING_OPTION.by_conductivity.name}', " + f"'{DISTRIBUTING_OPTION.by_crosscut_thickness.name}', " + f"'{DISTRIBUTING_OPTION.by_crosscut_transmissivity.name}', and " + f"'{DISTRIBUTING_OPTION.by_corrected_transmissivity.name}' supported. " + f"Got: '{distributing_option.name}'" ) return (weights * conductance).where(allocated) @@ -286,17 +305,19 @@ def distribute_ghb_conductance( case _: raise ValueError( "Received incompatible setting for general head boundary, only" - f"'{DISTRIBUTING_OPTION.equally.name}' and" - f"'{DISTRIBUTING_OPTION.by_layer_thickness.name}' and" - f"'{DISTRIBUTING_OPTION.by_layer_transmissivity.name}' and" - f"'{DISTRIBUTING_OPTION.by_conductivity.name}' supported." - f"got: '{distributing_option.name}'" + f"'{DISTRIBUTING_OPTION.equally.name}', " + f"'{DISTRIBUTING_OPTION.by_layer_thickness.name}', " + f"'{DISTRIBUTING_OPTION.by_layer_transmissivity.name}', and " + f"'{DISTRIBUTING_OPTION.by_conductivity.name}' supported. " + f"Got: '{distributing_option.name}'" ) return (weights * conductance).where(allocated) @preserve_gridtype -def _compute_layer_thickness(allocated, top, bottom): +def _compute_layer_thickness( + allocated: GridDataArray, top: GridDataArray, bottom: GridDataArray +): """ Compute 3D grid of thicknesses in allocated cells """ @@ -307,12 +328,17 @@ def _compute_layer_thickness(allocated, top, bottom): @preserve_gridtype -def _compute_crosscut_thickness(allocated, top, bottom, bc_top=None, bc_bottom=None): +def _compute_crosscut_thickness( + allocated: GridDataArray, + top: GridDataArray, + bottom: GridDataArray, + bc_top: Optional[GridDataArray] = None, + bc_bottom: Optional[GridDataArray] = None, +): """ Compute 3D grid of thicknesses crosscut by boundary condition (river/drain) - in allocated cells. So the thickness in the upper allocated layer thickness - is stage - bottom, the lower allocated layer is top - - river_bottom_elevation. + in allocated cells. So the thickness in the upper allocated layer is bc_top + - bottom and the lower allocated layer is top - bc_bottom. """ if (bc_top is None) & (bc_bottom is None): raise ValueError("`bc_top` and `bc_bottom` cannot both be None.") @@ -340,9 +366,9 @@ def _distribute_weights__by_corrected_transmissivity( allocated: GridDataArray, top: GridDataArray, bottom: GridDataArray, - stage: GridDataArray, - bottom_elevation: GridDataArray, k: GridDataArray, + bc_top: Optional[GridDataArray], + bc_bottom: Optional[GridDataArray], ): """ Distribute conductances according to default method in iMOD 5.6, as @@ -351,27 +377,28 @@ def _distribute_weights__by_corrected_transmissivity( Furthermore it corrects distribution weights for the mismatch between the midpoints of crosscut areas and layer midpoints. """ - PLANAR_GRID.validate(stage) - PLANAR_GRID.validate(bottom_elevation) - crosscut_thickness = _compute_crosscut_thickness( - allocated, top, bottom, stage, bottom_elevation + allocated, top, bottom, bc_top=bc_top, bc_bottom=bc_bottom ) transmissivity = crosscut_thickness * k top_layered = _enforce_layered_top(top, bottom) - - upper_layer_allocated = (stage < top_layered) & (stage > bottom) - lower_layer_allocated = (bottom_elevation < top_layered) & ( - bottom_elevation > bottom - ) - layer_thickness = _compute_layer_thickness(allocated, top, bottom) midpoints = (top_layered + bottom) / 2 + Fc = midpoints.copy() + + if bc_top is not None: + PLANAR_GRID.validate(bc_top) + upper_layer_bc = (bc_top < top_layered) & (bc_top > bottom) + # Computing vertical midpoint of river crosscutting layers. + Fc = Fc.where(~upper_layer_bc, (bottom + bc_top) / 2) + + if bc_bottom is not None: + PLANAR_GRID.validate(bc_bottom) + lower_layer_bc = (bc_bottom < top_layered) & (bc_bottom > bottom) + # Computing vertical midpoint of river crosscutting layers. + Fc = Fc.where(~lower_layer_bc, (top_layered + bc_bottom) / 2) - # Computing vertical midpoint of river crosscutting layers. - Fc = midpoints.where(~upper_layer_allocated, (bottom + stage) / 2) - Fc = Fc.where(~lower_layer_allocated, (top_layered + bottom_elevation) / 2) # Correction factor for mismatch between midpoints of crosscut layers and # layer midpoints. F = 1.0 - np.abs(midpoints - Fc) / (layer_thickness * 0.5) @@ -400,15 +427,17 @@ def _distribute_weights__by_crosscut_thickness( allocated: GridDataArray, top: GridDataArray, bottom: GridDataArray, - stage: GridDataArray, - bottom_elevation: GridDataArray, + bc_top: Optional[GridDataArray] = None, + bc_bottom: Optional[GridDataArray] = None, ): """Compute weights based on crosscut thickness""" - PLANAR_GRID.validate(stage) - PLANAR_GRID.validate(bottom_elevation) + if bc_top is not None: + PLANAR_GRID.validate(bc_top) + if bc_bottom is not None: + PLANAR_GRID.validate(bc_bottom) crosscut_thickness = _compute_crosscut_thickness( - allocated, top, bottom, stage, bottom_elevation + allocated, top, bottom, bc_top, bc_bottom ).where(allocated) return crosscut_thickness / crosscut_thickness.sum(dim="layer") @@ -431,16 +460,18 @@ def _distribute_weights__by_crosscut_transmissivity( allocated: GridDataArray, top: GridDataArray, bottom: GridDataArray, - stage: GridDataArray, - bottom_elevation: GridDataArray, k: GridDataArray, + bc_top: Optional[GridDataArray] = None, + bc_bottom: Optional[GridDataArray] = None, ): """Compute weights based on crosscut transmissivity""" - PLANAR_GRID.validate(stage) - PLANAR_GRID.validate(bottom_elevation) + if bc_top is not None: + PLANAR_GRID.validate(bc_top) + if bc_bottom is not None: + PLANAR_GRID.validate(bc_bottom) crosscut_thickness = _compute_crosscut_thickness( - allocated, top, bottom, stage, bottom_elevation + allocated, top, bottom, bc_top=bc_top, bc_bottom=bc_bottom ) transmissivity = crosscut_thickness * k From ea9a440c9d686a8eec46b5a890efd8e40fc47c81 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Mon, 22 Apr 2024 13:32:16 +0200 Subject: [PATCH 3/6] Add tests --- imod/prepare/topsystem/conductance.py | 14 +++---- imod/tests/test_prepare/test_topsystem.py | 2 +- .../test_prepare/test_topsystem_cases.py | 42 +++++++++++++++++++ 3 files changed, 50 insertions(+), 8 deletions(-) diff --git a/imod/prepare/topsystem/conductance.py b/imod/prepare/topsystem/conductance.py index 1014bf1a8..a1c7bb174 100644 --- a/imod/prepare/topsystem/conductance.py +++ b/imod/prepare/topsystem/conductance.py @@ -6,7 +6,7 @@ from imod.prepare.topsystem.allocation import _enforce_layered_top from imod.schemata import DimsSchema from imod.typing import GridDataArray -from imod.typing.grid import ones_like, preserve_gridtype +from imod.typing.grid import ones_like, preserve_gridtype, zeros_like from imod.util.dims import enforced_dim_order @@ -219,15 +219,15 @@ def distribute_drn_conductance( weights = _distribute_weights__by_conductivity(allocated, k) case DISTRIBUTING_OPTION.by_crosscut_thickness: weights = _distribute_weights__by_crosscut_thickness( - allocated, top, bottom, bc_top=elevation + allocated, top, bottom, bc_bottom=elevation ) case DISTRIBUTING_OPTION.by_crosscut_transmissivity: weights = _distribute_weights__by_crosscut_transmissivity( - allocated, top, bottom, k, bc_top=elevation + allocated, top, bottom, k, bc_bottom=elevation ) case DISTRIBUTING_OPTION.by_corrected_transmissivity: weights = _distribute_weights__by_corrected_transmissivity( - allocated, top, bottom, k, bc_top=elevation + allocated, top, bottom, k, bc_bottom=elevation ) case _: raise ValueError( @@ -345,7 +345,7 @@ def _compute_crosscut_thickness( top_layered = _enforce_layered_top(top, bottom) thickness = _compute_layer_thickness(allocated, top, bottom) - outside = ones_like(allocated).astype(bool) + outside = zeros_like(allocated).astype(bool) if bc_top is not None: upper_layer_bc = (bc_top < top_layered) & (bc_top > bottom) @@ -367,8 +367,8 @@ def _distribute_weights__by_corrected_transmissivity( top: GridDataArray, bottom: GridDataArray, k: GridDataArray, - bc_top: Optional[GridDataArray], - bc_bottom: Optional[GridDataArray], + bc_top: Optional[GridDataArray]=None, + bc_bottom: Optional[GridDataArray]=None, ): """ Distribute conductances according to default method in iMOD 5.6, as diff --git a/imod/tests/test_prepare/test_topsystem.py b/imod/tests/test_prepare/test_topsystem.py index ed31e5ae6..f8b43cdc6 100644 --- a/imod/tests/test_prepare/test_topsystem.py +++ b/imod/tests/test_prepare/test_topsystem.py @@ -136,7 +136,7 @@ def test_distribute_drn_conductance( conductance = zeros_like(elevation) + 1.0 actual_da = distribute_drn_conductance( - option, allocated, conductance, top, bottom, k + option, allocated, conductance, top, bottom, k, elevation ) actual = take_first_cell_in_xy_plane(actual_da) diff --git a/imod/tests/test_prepare/test_topsystem_cases.py b/imod/tests/test_prepare/test_topsystem_cases.py index ed8db2a2a..12d8d6a7d 100644 --- a/imod/tests/test_prepare/test_topsystem_cases.py +++ b/imod/tests/test_prepare/test_topsystem_cases.py @@ -165,6 +165,48 @@ def distribution_by_corrected_transmissivity__TFTF(): return option, allocated_layer, expected +@case(tags=["drn"]) +def distribution_by_corrected_transmissivity__drn(): + option = DISTRIBUTING_OPTION.by_corrected_transmissivity + allocated_layer = xr.DataArray( + [False, True, True, False], coords={"layer": [1, 2, 3, 4]}, dims=("layer",) + ) + expected = [np.nan, 1.0, 0.0, np.nan] + return option, allocated_layer, expected + + +@case(tags=["drn"]) +def distribution_by_corrected_transmissivity__first_active__drn(): + """First active cell allocated, while stage in second layer""" + option = DISTRIBUTING_OPTION.by_corrected_transmissivity + allocated_layer = xr.DataArray( + [True, True, True, False], coords={"layer": [1, 2, 3, 4]}, dims=("layer",) + ) + expected = [(2 / 3), (1 / 3), 0.0, np.nan] + return option, allocated_layer, expected + + +@case(tags=["drn"]) +def distribution_by_corrected_transmissivity__third_only__drn(): + """Third layer active only, while stage in second layer""" + option = DISTRIBUTING_OPTION.by_corrected_transmissivity + allocated_layer = xr.DataArray( + [False, False, True, False], coords={"layer": [1, 2, 3, 4]}, dims=("layer",) + ) + expected = [np.nan, np.nan, 0.0, np.nan] + return option, allocated_layer, expected + + +@case(tags=["drn"]) +def distribution_by_corrected_transmissivity__TFTF__drn(): + option = DISTRIBUTING_OPTION.by_corrected_transmissivity + allocated_layer = xr.DataArray( + [True, False, True, False], coords={"layer": [1, 2, 3, 4]}, dims=("layer",) + ) + expected = [1.0, np.nan, 0.0, np.nan] + return option, allocated_layer, expected + + @case(tags=["riv"]) def distribution_by_crosscut_transmissivity(): option = DISTRIBUTING_OPTION.by_crosscut_transmissivity From 8b7fb53d566c16ee46d6c5594ee03e217c6ca5a4 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Mon, 22 Apr 2024 14:20:42 +0200 Subject: [PATCH 4/6] Add test cases for new drain options --- .../test_prepare/test_topsystem_cases.py | 91 ++++++++++++++++++- 1 file changed, 90 insertions(+), 1 deletion(-) diff --git a/imod/tests/test_prepare/test_topsystem_cases.py b/imod/tests/test_prepare/test_topsystem_cases.py index 12d8d6a7d..7ddab24b8 100644 --- a/imod/tests/test_prepare/test_topsystem_cases.py +++ b/imod/tests/test_prepare/test_topsystem_cases.py @@ -193,7 +193,8 @@ def distribution_by_corrected_transmissivity__third_only__drn(): allocated_layer = xr.DataArray( [False, False, True, False], coords={"layer": [1, 2, 3, 4]}, dims=("layer",) ) - expected = [np.nan, np.nan, 0.0, np.nan] + # xarray turns 0.0 / 0.0 into np.nan, so third element should also be np.nan + expected = [np.nan, np.nan, np.nan, np.nan] return option, allocated_layer, expected @@ -250,6 +251,50 @@ def distribution_by_crosscut_transmissivity__TFTF(): return option, allocated_layer, expected +@case(tags=["drn"]) +def distribution_by_crosscut_transmissivity__drn(): + option = DISTRIBUTING_OPTION.by_crosscut_transmissivity + allocated_layer = xr.DataArray( + [False, True, True, False], coords={"layer": [1, 2, 3, 4]}, dims=("layer",) + ) + expected = [np.nan, 1.0, 0.0, np.nan] + return option, allocated_layer, expected + + +@case(tags=["drn"]) +def distribution_by_crosscut_transmissivity__first_active__drn(): + """First active cell allocated, while drain elevation in second layer""" + option = DISTRIBUTING_OPTION.by_crosscut_transmissivity + allocated_layer = xr.DataArray( + [True, True, True, False], coords={"layer": [1, 2, 3, 4]}, dims=("layer",) + ) + expected = [0.5, 0.5, 0.0, np.nan] + return option, allocated_layer, expected + + +@case(tags=["drn"]) +def distribution_by_crosscut_transmissivity__third_only__drn(): + """Third layer active only, while drain elevation in second layer""" + option = DISTRIBUTING_OPTION.by_crosscut_transmissivity + allocated_layer = xr.DataArray( + [False, False, True, False], coords={"layer": [1, 2, 3, 4]}, dims=("layer",) + ) + # xarray turns 0.0 / 0.0 into np.nan, so third element should also be np.nan + expected = [np.nan, np.nan, np.nan, np.nan] + return option, allocated_layer, expected + + +@case(tags=["drn"]) +def distribution_by_crosscut_transmissivity__TFTF__drn(): + """First and third layer active, while drain elevation in second layer""" + option = DISTRIBUTING_OPTION.by_crosscut_transmissivity + allocated_layer = xr.DataArray( + [True, False, True, False], coords={"layer": [1, 2, 3, 4]}, dims=("layer",) + ) + expected = [1.0, np.nan, 0.0, np.nan] + return option, allocated_layer, expected + + @case(tags=["riv"]) def distribution_by_crosscut_thickness(): option = DISTRIBUTING_OPTION.by_crosscut_thickness @@ -293,6 +338,50 @@ def distribution_by_crosscut_thickness__TFTF(): return option, allocated_layer, expected +@case(tags=["drn"]) +def distribution_by_crosscut_thickness__drn(): + option = DISTRIBUTING_OPTION.by_crosscut_thickness + allocated_layer = xr.DataArray( + [False, True, True, False], coords={"layer": [1, 2, 3, 4]}, dims=("layer",) + ) + expected = [np.nan, 1.0, 0.0, np.nan] + return option, allocated_layer, expected + + +@case(tags=["drn"]) +def distribution_by_crosscut_thickness__first_active__drn(): + """First active cell allocated, while drain elevation in second layer""" + option = DISTRIBUTING_OPTION.by_crosscut_thickness + allocated_layer = xr.DataArray( + [True, True, True, False], coords={"layer": [1, 2, 3, 4]}, dims=("layer",) + ) + expected = [0.5, 0.5, 0.0, np.nan] + return option, allocated_layer, expected + + +@case(tags=["drn"]) +def distribution_by_crosscut_thickness__third_only__drn(): + """Third layer active only, while drain elevation in second layer""" + option = DISTRIBUTING_OPTION.by_crosscut_thickness + allocated_layer = xr.DataArray( + [False, False, True, False], coords={"layer": [1, 2, 3, 4]}, dims=("layer",) + ) + # xarray turns 0.0 / 0.0 into np.nan, so third element should also be np.nan + expected = [np.nan, np.nan, np.nan, np.nan] + return option, allocated_layer, expected + + +@case(tags=["drn"]) +def distribution_by_crosscut_thickness__TFTF__drn(): + """First and third layer active, while drain elevation in second layer""" + option = DISTRIBUTING_OPTION.by_crosscut_thickness + allocated_layer = xr.DataArray( + [True, False, True, False], coords={"layer": [1, 2, 3, 4]}, dims=("layer",) + ) + expected = [1.0, np.nan, 0.0, np.nan] + return option, allocated_layer, expected + + @case(tags=["riv", "drn", "ghb"]) def distribution_by_conductivity(): option = DISTRIBUTING_OPTION.by_conductivity From fe0ce3c6f3adec18b8028371bcfbbba7a9c7b0e4 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Mon, 22 Apr 2024 14:20:54 +0200 Subject: [PATCH 5/6] Format --- imod/prepare/topsystem/conductance.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/imod/prepare/topsystem/conductance.py b/imod/prepare/topsystem/conductance.py index a1c7bb174..f8d5ad92a 100644 --- a/imod/prepare/topsystem/conductance.py +++ b/imod/prepare/topsystem/conductance.py @@ -188,7 +188,7 @@ def distribute_drn_conductance( Hydraulic conductivities elevation: DataArray | UgridDataArray Drain elevation - + Returns ------- Conductances distributed over depth. @@ -367,8 +367,8 @@ def _distribute_weights__by_corrected_transmissivity( top: GridDataArray, bottom: GridDataArray, k: GridDataArray, - bc_top: Optional[GridDataArray]=None, - bc_bottom: Optional[GridDataArray]=None, + bc_top: Optional[GridDataArray] = None, + bc_bottom: Optional[GridDataArray] = None, ): """ Distribute conductances according to default method in iMOD 5.6, as From d839a360e540bbb788f265acbcc357d0a89991c0 Mon Sep 17 00:00:00 2001 From: Joeri van Engelen Date: Mon, 22 Apr 2024 16:31:25 +0200 Subject: [PATCH 6/6] Fix wrong arg order --- imod/tests/test_prepare/test_topsystem.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/imod/tests/test_prepare/test_topsystem.py b/imod/tests/test_prepare/test_topsystem.py index f8b43cdc6..55bd10a2f 100644 --- a/imod/tests/test_prepare/test_topsystem.py +++ b/imod/tests/test_prepare/test_topsystem.py @@ -111,7 +111,7 @@ def test_distribute_riv_conductance( conductance = zeros_like(bottom_elevation) + 1.0 actual_da = distribute_riv_conductance( - option, allocated, conductance, top, bottom, stage, bottom_elevation, k + option, allocated, conductance, top, bottom, k, stage, bottom_elevation ) actual = take_first_cell_in_xy_plane(actual_da)