Issue #1052 log filtered wells (#1067)

Fixes #1052 

# Description
when importing wells from imod5: If not all wells could be assigned to
the modflow6 model (this can be due to inactive cells,
or due to constraints of minimum thickness/permeability) then a log
message will be generated to document in the log
file that not all wells were successfully imported. 

# Checklist
- [X] Links to correct issue
- [X] Update changelog, if changes affect users
- [X] PR title starts with ``Issue #nr``, e.g. ``Issue #737``
- [X] Unit tests were added
- [ ] **If feature added**: Added/extended example

---------

Co-authored-by: Joeri van Engelen <[email protected]>

docs/api/changelog.rst

@@ -14,6 +14,7 @@ Changed
 ~~~~~~~
 - :class:`imod.mf6.Well` now also validates that well filter top is above well filter bottom
 - :function:`open_projectfile_data` now also imports well filter top and bottom.
+- :class:`imod.mf6.Well` now logs a warning if any wells are removed during writing.
 
 [Unreleased]
 ------------

imod/mf6/wel.py

@@ -1,5 +1,6 @@
 from __future__ import annotations
 
+import textwrap
 import warnings
 from typing import Any, Optional, Tuple, Union
 
@@ -11,7 +12,8 @@
 import xugrid as xu
 
 import imod
-from imod.logging import init_log_decorator
+from imod.logging import init_log_decorator, logger
+from imod.logging.loglevel import LogLevel
 from imod.logging.logging_decorators import standard_log_decorator
 from imod.mf6.boundary_condition import (
     BoundaryCondition,
@@ -414,7 +416,7 @@ def __create_dataset_vars(
         """
         Create dataset with all variables (rate, concentration), with a similar shape as the cellids.
         """
-        data_vars = ["rate"]
+        data_vars = ["id", "rate"]
         if "concentration" in wells_assigned.columns:
             data_vars.append("concentration")
 
@@ -548,9 +550,6 @@ def to_mf6_pkg(
 
         Parameters
         ----------
-        is_partitioned: bool
-        validate: bool
-            Run validation before converting
         active: {xarry.DataArray, xugrid.UgridDataArray}
             Grid with active cells.
         top: {xarry.DataArray, xugrid.UgridDataArray}
@@ -559,6 +558,11 @@ def to_mf6_pkg(
             Grid with bottom of model layers.
         k: {xarry.DataArray, xugrid.UgridDataArray}
             Grid with hydraulic conductivities.
+        validate: bool
+            Run validation before converting
+        is_partitioned: bool
+            Set to true if model has been partitioned
+
         Returns
         -------
         Mf6Wel
@@ -574,11 +578,11 @@ def to_mf6_pkg(
         minimum_thickness = self.dataset["minimum_thickness"].item()
 
         wells_df = self.__create_wells_df()
+        nwells_df = len(wells_df["id"].unique())
         wells_assigned = self.__create_assigned_wells(
             wells_df, active, top, bottom, k, minimum_k, minimum_thickness
         )
 
-        nwells_df = len(wells_df["id"].unique())
         nwells_assigned = (
             0 if wells_assigned.empty else len(wells_assigned["id"].unique())
         )
@@ -602,8 +606,31 @@ def to_mf6_pkg(
         ds["print_flows"] = self["print_flows"].values[()]
         ds["print_input"] = self["print_input"].values[()]
 
+        filtered_wells = [
+            id for id in wells_df["id"].unique() if id not in ds["id"].values
+        ]
+        if len(filtered_wells) > 0:
+            message = self.to_mf6_package_information(filtered_wells)
+            logger.log(loglevel=LogLevel.WARNING, message=message, additional_depth=2)
+
+        ds = ds.drop_vars("id")
+
         return Mf6Wel(**ds.data_vars)
 
+    def to_mf6_package_information(self, filtered_wells):
+        message = textwrap.dedent(
+            """Some wells were not placed in the MF6 well package. This 
+            can be due to inactive cells or permeability/thickness constraints.\n"""
+        )
+        if len(filtered_wells) < 10:
+            message += "The filtered wells are: \n"
+        else:
+            message += " The first 10 unplaced wells are: \n"
+
+        for i in range(min(10, len(filtered_wells))):
+            message += f" id = {filtered_wells[i]} x = {self.dataset['x'][int(filtered_wells[i])].values[()]}  y = {self.dataset['y'][int(filtered_wells[i])].values[()]} \n"
+        return message
+
     def mask(self, domain: GridDataArray) -> Well:
         """
         Mask wells based on two-dimensional domain. For three-dimensional

imod/tests/test_mf6/test_mf6_wel.py

@@ -1,4 +1,5 @@
 import pathlib
+import sys
 import tempfile
 import textwrap
 from contextlib import nullcontext as does_not_raise
@@ -10,6 +11,8 @@
 from pytest_cases import parametrize_with_cases
 
 import imod
+from imod.logging.config import LoggerType
+from imod.logging.loglevel import LogLevel
 from imod.mf6.utilities.grid import broadcast_to_full_domain
 from imod.mf6.write_context import WriteContext
 from imod.schemata import ValidationError
@@ -179,6 +182,110 @@ def test_to_mf6_pkg__high_lvl_transient(basic_dis, well_high_lvl_test_data_trans
     np.testing.assert_equal(mf6_ds["rate"].values, rate_expected)
 
 
+def test_to_mf6_pkg__logging_with_message(
+    tmp_path, basic_dis, well_high_lvl_test_data_transient
+):
+    # Arrange
+    logfile_path = tmp_path / "logfile.txt"
+    idomain, top, bottom = basic_dis
+    modified_well_fixture = list(well_high_lvl_test_data_transient)
+
+    # create an idomain where layer 1 is active and layer 2 and 3 are inactive.
+    # layer 1 has a bottom at -6, layer 2 at -35 and layer 3 at -120
+    # so only wells that have a filter top above -6 will end up in the simulation
+    active = idomain == 1
+
+    active.loc[1, :, :] = True
+    active.loc[2:, :, :] = False
+
+    # modify the well filter top and filter bottoms so that
+    # well 0 is not placed
+    # well 1 is partially placed
+    # well 2 is fully placed
+    # well 3 is partially placed
+    # wells 4 to 7 are nog placed
+    modified_well_fixture[2] = [
+        -6.0,
+        -3.0,
+        0.0,
+        0.0,
+        -6.0,
+        -6.0,
+        -6.0,
+        -6.0,
+    ]
+    modified_well_fixture[3] = [
+        -102.0,
+        -102.0,
+        -6,
+        -102.0,
+        -1020.0,
+        -1020.0,
+        -1020.0,
+        -1020.0,
+    ]
+
+    with open(logfile_path, "w") as sys.stdout:
+        imod.logging.configure(
+            LoggerType.PYTHON,
+            log_level=LogLevel.DEBUG,
+            add_default_file_handler=False,
+            add_default_stream_handler=True,
+        )
+
+        wel = imod.mf6.Well(*modified_well_fixture)
+
+        k = xr.ones_like(idomain)
+        _ = wel.to_mf6_pkg(active, top, bottom, k)
+
+    # the wells that were fully or partially placed should not appear in the log message
+    # but all the wells that are completely left out should be listed
+    with open(logfile_path, "r") as log_file:
+        log = log_file.read()
+        assert "Some wells were not placed" in log
+        assert "id = 1" not in log
+        assert "id = 2" not in log
+        assert "id = 3" not in log
+        assert "id = 0" in log
+        assert "id = 4" in log
+        assert "id = 5" in log
+        assert "id = 6" in log
+        assert "id = 7" in log
+
+
+def test_to_mf6_pkg__logging_without_message(
+    tmp_path, basic_dis, well_high_lvl_test_data_transient
+):
+    # This test activates logging, and then converts a high level well package to
+    # an MF6 package, in such a way that all the wells can be placed.
+    # Logging is active, and the log file should not include the "Some wells were not placed"
+    # message
+    logfile_path = tmp_path / "logfile.txt"
+    idomain, top, bottom = basic_dis
+
+    with open(logfile_path, "w") as sys.stdout:
+        imod.logging.configure(
+            LoggerType.PYTHON,
+            log_level=LogLevel.DEBUG,
+            add_default_file_handler=False,
+            add_default_stream_handler=True,
+        )
+
+        wel = imod.mf6.Well(*well_high_lvl_test_data_transient)
+        active = idomain == 1
+        k = xr.ones_like(idomain)
+
+        active.loc[1, :, :] = True
+        active.loc[2:, :, :] = True
+
+        # Act
+        _ = wel.to_mf6_pkg(active, top, bottom, k)
+
+    with open(logfile_path, "r") as log_file:
+        log = log_file.read()
+        assert "Some wells were not placed" not in log
+
+
 @pytest.mark.parametrize("save_flows", [True, False])
 @pytest.mark.parametrize("print_input", [True, False])
 @pytest.mark.parametrize("print_flows", [True, False])