Add cleanup_layered_wel function
#1472
Conversation
|
…ed by adding this function, update of the unittests.
…//github.com/Deltares/imod-python into issue_#1471_add_cleanup_layered_wel_function
|
|
Besides adding a basic cleanup function (imod/mf6/wel.py, imod/prepare/cleanup.py), as was the main concern of this topic, two bugs were solved:
Solving these issues resulted in updating the well test fixtures, by adding 4 additional well locations (12 wells in total). |
There was a problem hiding this comment.
Looks good overall! Thanks for investigating the bugs and the fixes, especially updating all the tests to have 12 wells.
My main comment pertains to the code you added to the clipping functionality. I think some more info is needed here to assess if the added complexity is realIy necessary in its current form. Let's discuss that on Monday.
Furthermore we can reduce code duplication somewhat by refactoring the cleanup methods a bit. I have a suggestion how to do that which hopefully works in one go.
I also saw some TeamCity jobs are failing:
- The PipPython one will be resolved after merging the
maininto this branch. - The Examples one is flaky in the CI, I gave that another spin. Hopefully that fixes things. You could try running them locally to see if it works on your local machine by running:
pixi run examples
Furthermore there are some merge conflicts to resolve. Let me know if you run into issues with that, happy to help there.
| if (z_max is not None) or (z_min is not None): | ||
| in_bounds_z_max = None | ||
| in_bounds_z_min = None | ||
| if (z_max is not None) and (is_spatial_grid(top)): | ||
| in_bounds_z_max = np.full(ds.sizes["index"], False) | ||
| in_bounds_z_max[z_max["index"]] = True | ||
| z_max = z_max.drop_vars(lambda x: x.coords) | ||
| if (z_min is not None) and (is_spatial_grid(bottom)): | ||
| in_bounds_z_min = np.full(ds.sizes["index"], False) | ||
| in_bounds_z_min[z_min["index"]] = True | ||
| z_min = z_min.drop_vars(lambda x: x.coords) | ||
| if (in_bounds_z_max is not None) or (in_bounds_z_min is not None): | ||
| in_bounds = np.full(ds.sizes["index"], False) | ||
| if in_bounds_z_max is not None: | ||
| in_bounds = in_bounds | in_bounds_z_max | ||
| if in_bounds_z_min is not None: | ||
| in_bounds = in_bounds | in_bounds_z_min | ||
| ds = ds.loc[{"index": in_bounds}] |
There was a problem hiding this comment.
So if I'm not mistaken, this piece of code filters out values that are out of z-bounds? I expected that to be filtered out in the call to values_within_range. This adds quite some extra code complexity, where I sense the issue could be resolved in either _find_well_value_at_layer, or values_within_range. Does this code affect the tests? As in: Do they run into errors without it?
If not, I think it is best to pick this up in a separate PR, for which for my understanding I think it is good to:
- Create a separate issue describing what is going wrong, with a small reproducible example
- Sit together for 15 minutes to come up with the best solution, we have this prototype code which is very useful to see what is required.
- Add a separate unit test for the test case where this is a problem, reuse reproducible example from the issue here.
Either way, let's discuss on Monday!
| @standard_log_decorator() | ||
| def cleanup(self, dis: StructuredDiscretization | VerticesDiscretization): | ||
| """ | ||
| Clean up package inplace. This method calls | ||
| :func:`imod.prepare.cleanup_wel_layered`, see documentation of that | ||
| function for details on cleanup. | ||
|
|
||
| dis: imod.mf6.StructuredDiscretization | imod.mf6.VerticesDiscretization | ||
| Model discretization package. | ||
| """ | ||
| # Top and bottom should be forced to grids with a x, y coordinates | ||
| top, bottom = broadcast_to_full_domain(**dict(dis.dataset.data_vars)) | ||
| # Collect point variable datanames | ||
| point_varnames = list(self._write_schemata.keys()) | ||
| if "concentration" not in self.dataset.keys(): | ||
| point_varnames.remove("concentration") | ||
| point_varnames.append("id") | ||
| # Create dataset with purely point locations | ||
| point_ds = self.dataset[point_varnames] | ||
| # Take first item of irrelevant dimensions | ||
| point_ds = point_ds.isel(time=0, species=0, drop=True, missing_dims="ignore") | ||
| # Cleanup well dataframe | ||
| wells = point_ds.to_dataframe() | ||
| cleaned_wells = cleanup_wel_layered(wells, top.isel(layer=0), bottom) | ||
| # Select with ids in cleaned dataframe to drop points outside grid. | ||
| well_ids = cleaned_wells.index | ||
| dataset_cleaned = self.dataset.swap_dims({"index": "id"}).sel(id=well_ids) | ||
| # Ensure dtype of id is preserved | ||
| id_type = self.dataset["id"].dtype | ||
| dataset_cleaned = dataset_cleaned.swap_dims({"id": "index"}).reset_coords("id") | ||
| dataset_cleaned["id"] = dataset_cleaned["id"].astype(id_type) | ||
| # Override dataset | ||
| self.dataset = dataset_cleaned |
There was a problem hiding this comment.
Nice! I see you duplicated the logic in the Well package. I think we can reduce code duplication here by adding a semi-private method to the GridAgnosticWell class. Coul you try if this works?
class GridAgnosticWell(BoundaryCondition, IPointDataPackage, abc.ABC):
...
def _cleanup(self,
dis: StructuredDiscretization | VerticesDiscretization,
cleanup_func: Callable,
**cleanup_kwargs
) -> None:
# Top and bottom should be forced to grids with a x, y coordinates
top, bottom = broadcast_to_full_domain(**dict(dis.dataset.data_vars))
# Collect point variable datanames
point_varnames = list(self._write_schemata.keys())
if "concentration" not in self.dataset.keys():
point_varnames.remove("concentration")
point_varnames.append("id")
# Create dataset with purely point locations
point_ds = self.dataset[point_varnames]
# Take first item of irrelevant dimensions
point_ds = point_ds.isel(time=0, species=0, drop=True, missing_dims="ignore")
# Cleanup well dataframe
wells = point_ds.to_dataframe()
cleaned_wells = cleanup_func(wells, top.isel(layer=0), bottom, **cleanup_kwargs)
# Select with ids in cleaned dataframe to drop points outside grid.
well_ids = cleaned_wells.index
dataset_cleaned = self.dataset.swap_dims({"index": "id"}).sel(id=well_ids)
# Ensure dtype of id is preserved
id_type = self.dataset["id"].dtype
dataset_cleaned = dataset_cleaned.swap_dims({"id": "index"}).reset_coords("id")
dataset_cleaned["id"] = dataset_cleaned["id"].astype(id_type)
# Override dataset
self.dataset = dataset_cleaned
class Well(GridAgnosticWell):
...
def cleanup(self, dis: StructuredDiscretization | VerticesDiscretization):
minimum_thickness = float(self.dataset["minimum_thickness"])
self._cleanup(dis, cleanup_wel, minimum_thickness=minimum_thickness)
class LayeredWell(GridAgnosticWell):
...
def cleanup(self, dis: StructuredDiscretization | VerticesDiscretization):
self._cleanup(dis, cleanup_wel_layered)| @@ -40,6 +40,8 @@ Added | |||
| - :meth:`imod.mf6.HorizontalFlowBarrierResistance.cleanup`, | |||
| :meth:`imod.mf6.SingleLayerHorizontalFlowBarrierResistance.cleanup`, | |||
| to clean up HFB geometries crossing inactive model cells. | |||
| - :func:``imod.prepare.cleanup.cleanup_layered_wel` to clean up wells assigned | |||
| to layers. | |||
There was a problem hiding this comment.
I think you can also add under the "Fixed" header, as this is very important for users!
Fixed
~~~~~
...
- Sorting issue in :func:`imod.prepare.assign_wells`. This could cause :class:`imod.mf6.Well` to assign wells to the wrong cells.
If we also decide to include the fixes made to the well clipping (see my other comment), these should also be added to the changelog)
| @@ -61,13 +69,26 @@ def mf6wel_test_data_transient(): | |||
|
|
|||
| @pytest.fixture(scope="session") | |||
| def well_high_lvl_test_data_stationary(): | |||
There was a problem hiding this comment.
Please add a comment why we need 12 wells: To also test if there are no sorting issues because of the default assigned "id" as strings. Otherwise somebody might decide to clean this up later and reduce it to 5 wells or something ;-)
| @@ -83,13 +104,26 @@ def well_high_lvl_test_data_stationary(): | |||
|
|
|||
| @pytest.fixture(scope="session") | |||
| def well_high_lvl_test_data_transient(): | |||
| screen_top = [0.0, 0.0, 0.0, 0.0, -6.0, -6.0, -6.0, -6.0] | |||
There was a problem hiding this comment.
See my previous comment: idem dito here
Fixes #1471
Description
PR for @verkaik to push his adding of functionality.
Checklist
Issue #nr, e.g.Issue #737