Generalize inlet coordinates and flow direction

pyDeltaRCM/default.yml

@@ -186,4 +186,13 @@ clobber_netcdf:
   default: False
 legacy_netcdf:
   type: 'bool'
-  default: False
+  default: False
+inlet_x:
+  type: ['list', 'None']
+  default: null
+inlet_y:
+  type: ['list', 'None']
+  default: null
+inlet_flow_dir:
+  type: ['list', 'None']
+  default: null

pyDeltaRCM/init_tools.py

@@ -631,7 +631,44 @@ def create_domain(self) -> None:
         self.cell_type[: self.L0, :] = cell_land
         self.cell_type[: self.L0, channel_inds:y_channel_max] = cell_channel
 
-        self.inlet = np.array(np.unique(np.where(self.cell_type == 1)[1]))
+        has_x = hasattr(self, 'inlet_x') and self.inlet_x is not None
+        has_y = hasattr(self, 'inlet_y') and self.inlet_y is not None
+        if has_x or has_y:
+            if not (has_x and has_y):
+                missing = 'inlet_y' if has_x else 'inlet_x'
+                provided = 'inlet_x' if has_x else 'inlet_y'
+                val = getattr(self, provided)
+                raise ValueError(
+                    f"Both `inlet_x` and `inlet_y` must be provided if custom inlet coordinates are used. "
+                    f"Specified `{provided}`={val}, but `{missing}` is missing or None."
+                )
+            if len(self.inlet_x) != len(self.inlet_y):
+                raise ValueError(
+                    f"`inlet_x` and `inlet_y` must have the same length, "
+                    f"but got len(inlet_x)={len(self.inlet_x)} ({self.inlet_x}) and len(inlet_y)={len(self.inlet_y)} ({self.inlet_y})."
+                )
+            inlet_x_arr = np.array(self.inlet_x)
+            inlet_y_arr = np.array(self.inlet_y)
+            invalid_x = inlet_x_arr[(inlet_x_arr < 0) | (inlet_x_arr >= self.L)]
+            if len(invalid_x) > 0:
+                raise ValueError(
+                    f"inlet_x values must be within domain length (0 to L-1 = {self.L - 1}), "
+                    f"but got invalid values: {invalid_x.tolist()}."
+                )
+            invalid_y = inlet_y_arr[(inlet_y_arr < 0) | (inlet_y_arr >= self.W)]
+            if len(invalid_y) > 0:
+                raise ValueError(
+                    f"inlet_y values must be within domain width (0 to W-1 = {self.W - 1}), "
+                    f"but got invalid values: {invalid_y.tolist()}."
+                )
+            self.inlet = np.ravel_multi_index((inlet_x_arr, inlet_y_arr), self.cell_type.shape)
+            self.cell_type[inlet_x_arr, inlet_y_arr] = cell_channel
+        else:
+            inlet_y = np.array(np.unique(np.where(self.cell_type[0, :] == 1)[0]))
+            self.inlet = np.ravel_multi_index((np.zeros_like(inlet_y), inlet_y), self.cell_type.shape)
+
+        if not hasattr(self, 'inlet_flow_dir') or self.inlet_flow_dir is None:
+            self.inlet_flow_dir = [1, 0]
         self.eta[:] = self.stage - self.depth
 
         # update eta trackers with initial bed elevation

pyDeltaRCM/iteration_tools.py

@@ -95,13 +95,13 @@ def finalize_timestep(self) -> None:
 
         # apply bed elevation boundary condition at inlet
         #   first, calc the change in eta at inlet
-        _eta_change = (self.stage[0, self.inlet] - self._h0) - self.eta[0, self.inlet]
+        _eta_change = (self.stage.flat[self.inlet] - self._h0) - self.eta.flat[self.inlet]
         self._Vp_inletbc = (
             np.sum(_eta_change) * self._dx * self._dx
         )  # for mass cons checks
         #   now apply boundary condition
-        self.eta[0, self.inlet] += _eta_change
-        self.depth[0, self.inlet] = self._h0
+        self.eta.flat[self.inlet] += _eta_change
+        self.depth.flat[self.inlet] = self._h0
 
         self.hook_compute_sand_frac()
         self.compute_sand_frac()

pyDeltaRCM/model.py

@@ -361,6 +361,81 @@ def L0_meters(self, L0_meters: float) -> None:
             raise ValueError("L0_meters must be a greater than or equal to 0.")
         self._L0_meters = L0_meters
 
+    @property
+    def inlet_x(self) -> list:
+        """
+        inlet_x specifies the x-coordinates of the inlet cells.
+
+        This should be a list of integers representing the x-coordinates (rows)
+        where the inlet is located. If set to `None`, the default boundary 
+        behavior is used.
+        """
+        return self._inlet_x
+
+    @inlet_x.setter
+    def inlet_x(self, inlet_x: list) -> None:
+        if inlet_x is not None:
+            invalid_vals = [x for x in inlet_x if x < 0]
+            if invalid_vals:
+                raise ValueError(
+                    f"inlet_x values must be greater than or equal to 0, but got negative values: {invalid_vals}."
+                )
+            if hasattr(self, "_inlet_y") and self._inlet_y is not None:
+                if len(inlet_x) != len(self._inlet_y):
+                    raise ValueError(
+                        f"inlet_x and inlet_y must have the same length, but got len(inlet_x)={len(inlet_x)} and len(inlet_y)={len(self._inlet_y)}."
+                    )
+        self._inlet_x = inlet_x
+
+    @property
+    def inlet_y(self) -> list:
+        """
+        inlet_y specifies the y-coordinates of the inlet cells.
+
+        This should be a list of integers representing the y-coordinates (columns)
+        where the inlet is located. If set to `None`, the default boundary 
+        behavior is used. Must be the same length as `inlet_x`.
+        """
+        return self._inlet_y
+
+    @inlet_y.setter
+    def inlet_y(self, inlet_y: list) -> None:
+        if inlet_y is not None:
+            invalid_vals = [y for y in inlet_y if y < 0]
+            if invalid_vals:
+                raise ValueError(
+                    f"inlet_y values must be greater than or equal to 0, but got negative values: {invalid_vals}."
+                )
+            if hasattr(self, "_inlet_x") and self._inlet_x is not None:
+                if len(inlet_y) != len(self._inlet_x):
+                    raise ValueError(
+                        f"inlet_x and inlet_y must have the same length, but got len(inlet_x)={len(self._inlet_x)} and len(inlet_y)={len(inlet_y)}."
+                    )
+        self._inlet_y = inlet_y
+
+    @property
+    def inlet_flow_dir(self) -> list:
+        """
+        inlet_flow_dir sets the primary flow direction of water parcels at the inlet.
+
+        This should be a list containing two elements representing the flow vector
+        [dx, dy] at the inlet. If `None`, defaults to [1, 0].
+        """
+        return self._inlet_flow_dir
+
+    @inlet_flow_dir.setter
+    def inlet_flow_dir(self, inlet_flow_dir: list) -> None:
+        if inlet_flow_dir is not None:
+            if len(inlet_flow_dir) != 2:
+                raise ValueError(
+                    f"inlet_flow_dir must be a list of length 2, but got {inlet_flow_dir} of length {len(inlet_flow_dir)}."
+                )
+            if not all(isinstance(v, (int, float)) for v in inlet_flow_dir):
+                raise ValueError(
+                    f"inlet_flow_dir elements must be numeric, but got {inlet_flow_dir}."
+                )
+        self._inlet_flow_dir = inlet_flow_dir
+
     @property
     def S0(self) -> float:
         """

pyDeltaRCM/water_tools.py

@@ -114,9 +114,9 @@ def run_water_iteration(self) -> None:
 
         # flux from ghost node
         start_inlets, start_counts = np.unique(start_indices, return_counts=True)
-        self.qxn.flat[start_inlets] += start_counts
-        self.qyn.flat[start_indices] += 0  # this could be omitted...
-        self.qwn.flat[start_indices] += self.Qp_water / self._dx / 2
+        self.qxn.flat[start_inlets] += start_counts * self.inlet_flow_dir[0]
+        self.qyn.flat[start_inlets] += start_counts * self.inlet_flow_dir[1]
+        self.qwn.flat[start_inlets] += start_counts * (self.Qp_water / self._dx / 2)
 
         # load the initial indices into the walk indices
         self.free_surf_walk_inds[:, _step] = start_indices
@@ -514,9 +514,9 @@ def update_flow_field(self, iteration: int) -> None:
 
         self.qw = (self.qx**2 + self.qy**2) ** (0.5)
 
-        self.qx[0, self.inlet] = self.qw0
-        self.qy[0, self.inlet] = 0
-        self.qw[0, self.inlet] = self.qw0
+        self.qx.flat[self.inlet] = self.qw0 * self.inlet_flow_dir[0]
+        self.qy.flat[self.inlet] = self.qw0 * self.inlet_flow_dir[1]
+        self.qw.flat[self.inlet] = self.qw0
 
     def update_velocity_field(self) -> None:
         """Update flow velocity fields.

tests/test_model.py

@@ -201,6 +201,93 @@ def test_negative_L0_meters(self, tmp_path: Path) -> None:
         with pytest.raises(ValueError):
             _ = DeltaModel(input_file=p)
 
+    def test_negative_inlet_x(self, tmp_path: Path) -> None:
+        file_name = "user_parameters.yaml"
+        p, f = utilities.create_temporary_file(tmp_path, file_name)
+        utilities.write_parameter_to_file(f, "out_dir", tmp_path / "out_dir")
+        utilities.write_parameter_to_file(f, "inlet_x", [-1])
+        f.close()
+        with pytest.raises(ValueError):
+            _ = DeltaModel(input_file=p)
+
+    def test_negative_inlet_y(self, tmp_path: Path) -> None:
+        file_name = "user_parameters.yaml"
+        p, f = utilities.create_temporary_file(tmp_path, file_name)
+        utilities.write_parameter_to_file(f, "out_dir", tmp_path / "out_dir")
+        utilities.write_parameter_to_file(f, "inlet_y", [-1])
+        f.close()
+        with pytest.raises(ValueError):
+            _ = DeltaModel(input_file=p)
+
+    def test_bad_length_inlet_flow_dir(self, tmp_path: Path) -> None:
+        file_name = "user_parameters.yaml"
+        p, f = utilities.create_temporary_file(tmp_path, file_name)
+        utilities.write_parameter_to_file(f, "out_dir", tmp_path / "out_dir")
+        utilities.write_parameter_to_file(f, "inlet_flow_dir", [1, 0, 0])
+        f.close()
+        with pytest.raises(ValueError):
+            _ = DeltaModel(input_file=p)
+
+    def test_mismatched_length_inlet_x_y(self, tmp_path: Path) -> None:
+        file_name = "user_parameters.yaml"
+        p, f = utilities.create_temporary_file(tmp_path, file_name)
+        utilities.write_parameter_to_file(f, "out_dir", tmp_path / "out_dir")
+        utilities.write_parameter_to_file(f, "inlet_x", [0, 0])
+        utilities.write_parameter_to_file(f, "inlet_y", [10])
+        f.close()
+        with pytest.raises(ValueError):
+            _ = DeltaModel(input_file=p)
+
+    def test_inlet_x_only_raises(self, tmp_path: Path) -> None:
+        file_name = "user_parameters.yaml"
+        p, f = utilities.create_temporary_file(tmp_path, file_name)
+        utilities.write_parameter_to_file(f, "out_dir", tmp_path / "out_dir")
+        utilities.write_parameter_to_file(f, "inlet_x", [0, 0])
+        f.close()
+        with pytest.raises(ValueError):
+            _ = DeltaModel(input_file=p)
+
+    def test_inlet_y_only_raises(self, tmp_path: Path) -> None:
+        file_name = "user_parameters.yaml"
+        p, f = utilities.create_temporary_file(tmp_path, file_name)
+        utilities.write_parameter_to_file(f, "out_dir", tmp_path / "out_dir")
+        utilities.write_parameter_to_file(f, "inlet_y", [10, 11])
+        f.close()
+        with pytest.raises(ValueError):
+            _ = DeltaModel(input_file=p)
+
+    def test_out_of_bounds_inlet_x(self, tmp_path: Path) -> None:
+        file_name = "user_parameters.yaml"
+        p, f = utilities.create_temporary_file(tmp_path, file_name)
+        utilities.write_parameter_to_file(f, "out_dir", tmp_path / "out_dir")
+        utilities.write_parameter_to_file(f, "inlet_x", [9999])
+        utilities.write_parameter_to_file(f, "inlet_y", [10])
+        f.close()
+        with pytest.raises(ValueError):
+            _ = DeltaModel(input_file=p)
+
+    def test_out_of_bounds_inlet_y(self, tmp_path: Path) -> None:
+        file_name = "user_parameters.yaml"
+        p, f = utilities.create_temporary_file(tmp_path, file_name)
+        utilities.write_parameter_to_file(f, "out_dir", tmp_path / "out_dir")
+        utilities.write_parameter_to_file(f, "inlet_x", [0])
+        utilities.write_parameter_to_file(f, "inlet_y", [9999])
+        f.close()
+        with pytest.raises(ValueError):
+            _ = DeltaModel(input_file=p)
+
+    def test_custom_inlet_valid_initialization(self, tmp_path: Path) -> None:
+        file_name = "user_parameters.yaml"
+        p, f = utilities.create_temporary_file(tmp_path, file_name)
+        utilities.write_parameter_to_file(f, "out_dir", tmp_path / "out_dir")
+        utilities.write_parameter_to_file(f, "inlet_x", [2, 2, 2])
+        utilities.write_parameter_to_file(f, "inlet_y", [10, 11, 12])
+        utilities.write_parameter_to_file(f, "inlet_flow_dir", [0, 1])
+        f.close()
+        delta = DeltaModel(input_file=p)
+        assert np.all(delta.cell_type[2, [10, 11, 12]] == 1)
+        assert delta.inlet_flow_dir == [0, 1]
+
     def test_negative_itermax(self, tmp_path: Path) -> None:
         file_name = "user_parameters.yaml"
         p, f = utilities.create_temporary_file(tmp_path, file_name)