Avalanching numba patch

[Sierd]

Code still needs to be tested. Especially the NE_layer.

Also, 1D compatibility needs to be tested.

[Copilot]

Pull Request Overview

This PR refactors the avalanching module to use Numba's JIT compilation for performance optimization. The main change extracts the iterative avalanche computation into a separate @njit-decorated function that uses explicit loops instead of NumPy boolean indexing, which is more compatible with Numba.

Key changes:

• Extracted avalanche iteration logic into a Numba-compiled avalanche_loop function
• Simplified flux computation from 4-directional to 2-directional (X and Y only)
• Modified parameter E from 0.2 to 0.1

---

_{Tip: Customize your code reviews with copilot-instructions.md. Create the file or learn how to get started.}

[Copilot]

Commented-out conditions suggest incomplete implementation. The original code had logic to handle grad_h_nonerod when it was less than grad_h - tan_dyn (lines 197-198 also commented). These conditions should either be properly implemented or removed with explanation in commit message if intentionally excluded.

Suggested change

	# grad_h_nonerod = (zb[i, j] - zne[i, j]) / (ds[i, j]*dn[i, j])
	if grad_h[i, j] > tan_dyn: # and (zb[i, j] - zne[i, j]) > 0.0:
	slope_diff[i, j] = np.tanh(grad_h[i, j]) - np.tanh(0.9 * tan_dyn)
	# elif grad_h_nonerod < (grad_h[i, j] - tan_dyn):
	# slope_diff[i, j] = np.tanh(grad_h_nonerod)
	# Only the main slope condition is currently implemented.
	# Handling for non-erodible gradients (grad_h_nonerod) is not included.
	if grad_h[i, j] > tan_dyn:
	slope_diff[i, j] = np.tanh(grad_h[i, j]) - np.tanh(0.9 * tan_dyn)

[Copilot]

The division by grad_h[i, j] is commented out, changing the flux calculation significantly from the original implementation. The original code normalized flux components by the total gradient magnitude. This change should be documented as it affects the physical behavior of the avalanching model.

Suggested change

	flux_down[i, j, 0] = slope_diff[i, j] * grad_h_down[i, j, 0]# / grad_h[i, j]
	flux_down[i, j, 1] = slope_diff[i, j] * grad_h_down[i, j, 1]# / grad_h[i, j]
	flux_down[i, j, 0] = slope_diff[i, j] * grad_h_down[i, j, 0] / grad_h[i, j]
	flux_down[i, j, 1] = slope_diff[i, j] * grad_h_down[i, j, 1] / grad_h[i, j]

[Copilot]

This unconditionally sets zb = zne for all cells when any cell violates the non-erodible layer condition. This is incorrect - only cells where zb + E * (q_in - q_out) < zne[i, j] should be limited to zne[i, j]. The update should be element-wise with proper masking or a per-cell conditional.

Suggested change

	if np.any((E * (q_in - q_out)) < (zne-zb)):
	# update bed level with non-erodible layer limit
	# this will effectively shut down further avalanching from the cells concerned
	# because zb will equal zne there in the next iteration
	zb += (zne - zb)
	else:
	# update bed level without non-erodible layer
	zb += E * (q_in - q_out)
	# Compute proposed new bed level
	zb_new = zb + E * (q_in - q_out)
	# Mask for cells where new bed level would fall below non-erodible layer
	mask = zb_new < zne
	# For masked cells, set zb to zne; otherwise, apply normal update
	zb[mask] = zne[mask]
	zb[~mask] = zb_new[~mask]

[Copilot]

Pull Request Overview

Copilot reviewed 2 out of 2 changed files in this pull request and generated 7 comments.

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[Copilot]

Large blocks of commented-out code should be removed. The normalization by grid spacing is already performed inline at lines 165 and 179, making this commented block redundant. If this alternative approach is needed for reference, it should be documented in version control history instead.

Suggested change

	# for i in range(ny):
	# for j in range(nx):
	# grad_h_down[i, j, 0] = grad_h_down[i, j, 0] / ds[i, j]
	# grad_h_down[i, j, 1] = grad_h_down[i, j, 1] / dn[i, j]
	# gradient magnitude and maximum
	# for i in range(ny):
	# for j in range(nx):
	# gradient magnitude and maximum

[Copilot]

Commented-out conditions and alternative logic should either be implemented or removed. The partial condition # and (zb[i, j] - zne[i, j]) > 0.0 and the entire elif block (lines 212-213) suggest incomplete implementation. Either complete this logic or remove the comments if the simplified version is intended.

Suggested change

	if grad_h[i, j] > tan_dyn: # and (zb[i, j] - zne[i, j]) > 0.0:
	slope_diff[i, j] = np.tanh(grad_h[i, j]) - np.tanh(0.9 * tan_dyn)
	# elif grad_h_nonerod < (grad_h[i, j] - tan_dyn):
	# slope_diff[i, j] = np.tanh(grad_h_nonerod)
	if grad_h[i, j] > tan_dyn:
	slope_diff[i, j] = np.tanh(grad_h[i, j]) - np.tanh(0.9 * tan_dyn)

[Copilot]

The commented-out division by grad_h[i, j] suggests uncertainty about the correct formula. Either implement the division with proper explanation, or remove the comment if the current implementation is correct. The presence of this comment indicates the flux calculation may need verification.

Suggested change

	flux_down[i, j, 0] = slope_diff[i, j] * grad_h_down[i, j, 0]# / grad_h[i, j]
	flux_down[i, j, 1] = slope_diff[i, j] * grad_h_down[i, j, 1]# / grad_h[i, j]
	# Normalize the direction vector to ensure flux is in the direction of the gradient with correct magnitude
	flux_down[i, j, 0] = slope_diff[i, j] * grad_h_down[i, j, 0] / grad_h[i, j]
	flux_down[i, j, 1] = slope_diff[i, j] * grad_h_down[i, j, 1] / grad_h[i, j]

[Copilot]

The implementation applies periodic boundary conditions in the X direction (lines 237-238, 242-243) and Y direction (lines 247-248, 252-253), but the gradient calculations use fixed zero-gradient boundaries (lines 155-156, 169-170). This inconsistency between gradient and flux boundary treatments will cause mass conservation issues at domain boundaries.

[Copilot]

When the non-erodible layer constraint is violated, all cells are updated to zb = zne (line 262) rather than only the cells that violate the constraint. This incorrectly affects cells where zb > zne. The update should be applied selectively: zb = np.maximum(zb + E * (q_in - q_out), zne).

Suggested change

	# # check mass balance in the presence of non-erodible layer
	if np.any(zb + E * (q_in - q_out) < zne):
	# update bed level with non-erodible layer limit
	# this will effectively shut down further avalanching from the cells concerned
	# because zb will equal zne there in the next iteration
	zb += (zne - zb)
	else:
	# update bed level without non-erodible layer
	zb += E * (q_in - q_out)
	# update bed level, ensuring it does not fall below the non-erodible layer
	zb = np.maximum(zb + E * (q_in - q_out), zne)