Apply non-erodible layer constraint only at domain boundaries

Per feedback, the problem is specifically at downwind boundaries, not throughout
the domain. Modified implementation to apply constraint only at domain edges:

- Added boundary detection logic in all quadrant solvers
- Constrain pickup only for cells at n=0, n=max, s=0, or s=max
- Added post-processing for boundary cells with copied pickup values
- Updated tests to verify boundary-only constraint behavior
- Updated documentation to reflect targeted boundary approach

This prevents the issue at downwind boundaries without affecting interior cells.

Co-authored-by: Sierd <14054272+Sierd@users.noreply.github.com>

Author: Copilot

SOLUTION_SUMMARY.md

@@ -30,9 +30,9 @@ At the downwind boundary:
 
 ## Solution Implementation
 
-A dual-layer defensive approach was implemented:
+A targeted boundary-focused approach was implemented:
 
-### Layer 1: Prevention in Sweep Solver
+### Layer 1: Prevention in Sweep Solver (Boundary-Only)
 
 Modified `aeolis/utils.py`:
 
@@ -43,13 +43,16 @@ Modified `aeolis/utils.py`:
    - `porosity`: Bed porosity
 
 2. Updated all quadrant solvers and generic stencil to:
-   - Calculate total pickup for each cell
+   - **Apply constraint ONLY at domain boundaries** (edges of grid)
+   - Calculate total pickup for each boundary cell
    - Convert pickup to bed level change: `dz = total_pickup / (rhog * (1 - porosity))`
    - Check if `zb - dz < zne`
    - If true, limit pickup to: `max_pickup = (zb - zne) * rhog * (1 - porosity)`
    - Scale down pickup proportionally for all fractions
 
-This prevents the solver from calculating excessive pickup in the first place.
+3. Added post-processing constraint for boundary cells that receive copied pickup values
+
+This prevents excessive pickup at domain boundaries where the problem occurs.
 
 ### Layer 2: Safety Net in Bed Update
 

aeolis/tests/test_nonerodible.py

@@ -16,10 +16,10 @@ class TestNonErodibleLayerConstraint:
     non-erodible layer constraint to prevent erosion below zne.
     """
 
-    def test_sweep_respects_nonerodible_layer(self):
+    def test_sweep_respects_nonerodible_layer_at_boundaries(self):
         """
-        Test that the sweep function limits pickup to prevent bed from
-        dropping below the non-erodible layer.
+        Test that the sweep function limits pickup at domain boundaries to prevent bed from
+        dropping below the non-erodible layer. The constraint is only applied at boundaries.
         """
         # Setup a simple 5x5 grid with 1 fraction
         ny, nx, nf = 4, 4, 1
@@ -66,15 +66,20 @@ def test_sweep_respects_nonerodible_layer(self):
         dz = total_pickup / mass_per_meter
         new_zb = zb - dz
 
-        # Check that bed level doesn't drop below non-erodible layer
-        # Allow small numerical tolerance
-        assert np.all(new_zb >= zne - 1e-6), \
-            f"Bed level dropped below non-erodible layer. Min new_zb: {new_zb.min()}, zne: {zne.min()}"
+        # Check that bed level at BOUNDARIES doesn't drop below non-erodible layer
+        # Define boundary cells (edges of domain)
+        boundary_mask = np.zeros_like(zb, dtype=bool)
+        boundary_mask[0, :] = True  # Top edge
+        boundary_mask[-1, :] = True  # Bottom edge
+        boundary_mask[:, 0] = True  # Left edge
+        boundary_mask[:, -1] = True  # Right edge
 
-        # Check that pickup was limited (should be less than available mass)
-        max_allowed_pickup = mass_per_meter * (zb - zne).max()
-        assert np.all(total_pickup <= max_allowed_pickup + 1e-6), \
-            f"Pickup exceeded maximum allowed. Max pickup: {total_pickup.max()}, max allowed: {max_allowed_pickup}"
+        # Check only boundary cells
+        assert np.all(new_zb[boundary_mask] >= zne[boundary_mask] - 1e-6), \
+            f"Bed level at boundaries dropped below non-erodible layer. Min new_zb: {new_zb[boundary_mask].min()}, zne: {zne[boundary_mask].min()}"
+        
+        # Interior cells may drop below zne (constraint not applied there)
+        # This is expected behavior as the problem is specifically at boundaries
 
     def test_sweep_without_nonerodible_layer_backward_compatible(self):
         """

aeolis/utils.py

@@ -582,6 +582,30 @@ def sweep(Ct, Cu, mass, dt, Ts, ds, dn, us, un, w, zb=None, zne=None, rhog=2650.
         # print(np.shape(visited[0,:]==False))
         pickup[0,:,0] = pickup[1,:,0].copy() 
         pickup[-1,:,0] = pickup[-2,:,0].copy() 
+        
+        # Apply non-erodible layer constraint to boundary cells that were set by copying
+        if zb is not None and zne is not None:
+            for boundary_row in [0, Ct.shape[0] - 1]:
+                for s in range(Ct.shape[1]):
+                    # Calculate total pickup for this cell
+                    total_pickup = 0.0
+                    for f in range(nf):
+                        total_pickup += pickup[boundary_row, s, f]
+                    
+                    if total_pickup > 0:
+                        # Calculate bed level change
+                        dz = total_pickup / (rhog * (1.0 - porosity))
+                        
+                        # Check if bed would drop below non-erodible layer
+                        if zb[boundary_row, s] - dz < zne[boundary_row, s]:
+                            # Limit pickup
+                            max_dz = max(0.0, zb[boundary_row, s] - zne[boundary_row, s])
+                            max_pickup_total = max_dz * rhog * (1.0 - porosity)
+                            
+                            if total_pickup > 0.0:
+                                scale_factor = max_pickup_total / total_pickup
+                                for f in range(nf):
+                                    pickup[boundary_row, s, f] *= scale_factor
 
         k+=1
 
@@ -657,8 +681,12 @@ def _solve_quadrant1(Ct, Cu, mass, pickup, dt, Ts, ds, dn, ufs, ufn, w, visited,
 
                         Ct[n, s, f] = num_limited / den_limited
 
-                # Apply non-erodible layer constraint if zb and zne are provided
-                if zb is not None and zne is not None:
+                # Apply non-erodible layer constraint at domain boundaries only
+                # Check if cell is at a boundary (any edge of domain)
+                is_boundary = (n == 1 or n == Ct.shape[0] - 1 or 
+                              s == 1 or s == Ct.shape[1] - 1)
+                
+                if zb is not None and zne is not None and is_boundary:
                     # Calculate total pickup mass for this cell
                     total_pickup = 0.0
                     for f in range(nf):
@@ -736,8 +764,12 @@ def _solve_quadrant2(Ct, Cu, mass, pickup, dt, Ts, ds, dn, ufs, ufn, w, visited,
 
                         Ct[n, s, f] = num_limited / den_limited
 
-                # Apply non-erodible layer constraint if zb and zne are provided
-                if zb is not None and zne is not None:
+                # Apply non-erodible layer constraint at domain boundaries only
+                # Check if cell is at a boundary (edges of domain)
+                is_boundary = (n == 1 or n == Ct.shape[0] - 1 or 
+                              s == 1 or s == Ct.shape[1] - 2)
+                
+                if zb is not None and zne is not None and is_boundary:
                     # Calculate total pickup mass for this cell
                     total_pickup = 0.0
                     for f in range(nf):
@@ -815,8 +847,12 @@ def _solve_quadrant3(Ct, Cu, mass, pickup, dt, Ts, ds, dn, ufs, ufn, w, visited,
 
                         Ct[n, s, f] = num_limited / den_limited
 
-                # Apply non-erodible layer constraint if zb and zne are provided
-                if zb is not None and zne is not None:
+                # Apply non-erodible layer constraint at domain boundaries only
+                # Check if cell is at a boundary (edges of domain)
+                is_boundary = (n == 0 or n == Ct.shape[0] - 2 or 
+                              s == 0 or s == Ct.shape[1] - 2)
+                
+                if zb is not None and zne is not None and is_boundary:
                     # Calculate total pickup mass for this cell
                     total_pickup = 0.0
                     for f in range(nf):
@@ -894,8 +930,12 @@ def _solve_quadrant4(Ct, Cu, mass, pickup, dt, Ts, ds, dn, ufs, ufn, w, visited,
 
                         Ct[n, s, f] = num_limited / den_limited
 
-                # Apply non-erodible layer constraint if zb and zne are provided
-                if zb is not None and zne is not None:
+                # Apply non-erodible layer constraint at domain boundaries only
+                # Check if cell is at a boundary (edges of domain)
+                is_boundary = (n == 0 or n == Ct.shape[0] - 2 or 
+                              s == 1 or s == Ct.shape[1] - 1)
+                
+                if zb is not None and zne is not None and is_boundary:
                     # Calculate total pickup mass for this cell
                     total_pickup = 0.0
                     for f in range(nf):
@@ -1004,8 +1044,12 @@ def _solve_generic_stencil(Ct, Cu, mass, pickup, dt, Ts, ds, dn, ufs, ufn, w, vi
 
                         Ct[n, s, f] = num_lim / den_lim
 
-                # Apply non-erodible layer constraint if zb and zne are provided
-                if zb is not None and zne is not None:
+                # Apply non-erodible layer constraint at domain boundaries only
+                # Check if cell is at a boundary (edges of domain)  
+                is_boundary = (n == 0 or n == Ct.shape[0] - 2 or 
+                              s == 1 or s == Ct.shape[1] - 1)
+                
+                if zb is not None and zne is not None and is_boundary:
                     # Calculate total pickup mass for this cell
                     total_pickup = 0.0
                     for f in range(nf):