From f9bbefa1404d7afa6755a6d66b6c3f0fb4e4b8ac Mon Sep 17 00:00:00 2001 From: Sierd Date: Thu, 23 Oct 2025 19:38:19 +0200 Subject: [PATCH 01/10] Avalanching numba patch Tested for 2D cases, seems to work fine now. Need to test for 1D cases. --- aeolis/avalanching.py | 285 +++++++++++++++++++----------------------- 1 file changed, 131 insertions(+), 154 deletions(-) diff --git a/aeolis/avalanching.py b/aeolis/avalanching.py index e77fa56f..7957a351 100644 --- a/aeolis/avalanching.py +++ b/aeolis/avalanching.py @@ -33,6 +33,8 @@ # package modules from aeolis.utils import * +from numba import njit + # initialize logger logger = logging.getLogger(__name__) @@ -90,178 +92,153 @@ def avalanche(s, p): Spatial grids ''' - if p['process_avalanche']: + nx = p['nx'] + 1 + ny = p['ny'] + 1 - nx = p['nx']+1 - ny = p['ny']+1 - - #parameters - + # parameters tan_stat = np.tan(np.deg2rad(s['theta_stat'])) tan_dyn = np.tan(np.deg2rad(s['theta_dyn'])) - - - - E = 0.2 - grad_h_down = np.zeros((ny,nx,4)) - flux_down = np.zeros((ny,nx,4)) - slope_diff = np.zeros((ny,nx)) - grad_h = np.zeros((ny,nx)) + E = 0.1 max_iter_ava = p['max_iter_ava'] - - max_grad_h, grad_h, grad_h_down = calc_gradients(s['zb'], nx, ny, s['ds'], s['dn'], s['zne']) - - s['gradh'] = grad_h.copy() - - initiate_avalanche = (max_grad_h > tan_stat) - - if initiate_avalanche: - - for i in range(0,max_iter_ava): - - grad_h_down *= 0. - flux_down *= 0. - slope_diff *= 0. - grad_h *= 0. - - max_grad_h, grad_h, grad_h_down = calc_gradients(s['zb'], nx, ny, s['ds'], s['dn'], s[ 'zne']) - - if max_grad_h < tan_dyn: - break - - # Calculation of flux - - grad_h_nonerod = (s['zb'] - s['zne']) / s['ds'] # HAS TO BE ADJUSTED! - - ix = np.logical_and(grad_h > tan_dyn, grad_h_nonerod > 0) - slope_diff[ix] = np.tanh(grad_h[ix]) - np.tanh(0.9*tan_dyn) - - ix = grad_h_nonerod < grad_h - tan_dyn - slope_diff[ix] = np.tanh(grad_h_nonerod[ix]) - - ix = grad_h != 0 - - if ny == 1: - #1D interpretation - flux_down[:,:,0][ix] = slope_diff[ix] * grad_h_down[:,:,0][ix] / grad_h[ix] - flux_down[:,:,2][ix] = slope_diff[ix] * grad_h_down[:,:,2][ix] / grad_h[ix] - - # Calculation of change in bed level - - q_in = np.zeros((ny,nx)) - - q_out = 0.5*np.abs(flux_down[:,:,0]) + 0.5*np.abs(flux_down[:,:,2]) - - q_in[0,1:-1] = 0.5*(np.maximum(flux_down[0,:-2,0],0.) \ - - np.minimum(flux_down[0,2:,0],0.) \ - + np.maximum(flux_down[0,2:,2],0.) \ - - np.minimum(flux_down[0,:-2,2],0.)) - else: - # 2D interpretation - flux_down[:,:,0][ix] = slope_diff[ix] * grad_h_down[:,:,0][ix] / grad_h[ix] - flux_down[:,:,1][ix] = slope_diff[ix] * grad_h_down[:,:,1][ix] / grad_h[ix] - flux_down[:,:,2][ix] = slope_diff[ix] * grad_h_down[:,:,2][ix] / grad_h[ix] - flux_down[:,:,3][ix] = slope_diff[ix] * grad_h_down[:,:,3][ix] / grad_h[ix] - - # Calculation of change in bed level - - q_in = np.zeros((ny,nx)) - - q_out = 0.5*np.abs(flux_down[:,:,0]) + 0.5* np.abs(flux_down[:,:,1]) + 0.5*np.abs(flux_down[:,:,2]) + 0.5* np.abs(flux_down[:,:,3]) - - q_in[1:-1,1:-1] = 0.5*(np.maximum(flux_down[1:-1,:-2,0],0.) \ - - np.minimum(flux_down[1:-1,2:,0],0.) \ - + np.maximum(flux_down[:-2,1:-1,1],0.) \ - - np.minimum(flux_down[2:,1:-1,1],0.) \ - - + np.maximum(flux_down[1:-1,2:,2],0.) \ - - np.minimum(flux_down[1:-1,:-2,2],0.) \ - + np.maximum(flux_down[2:,1:-1,3],0.) \ - - np.minimum(flux_down[:-2,1:-1,3],0.)) - - s['zb'] += E * (q_in - q_out) + + # call the njit-compiled loop for performance + zb, grad_h = avalanche_loop( + s['zb'].copy(), s['zne'], s['ds'], s['dn'], nx, ny, E, max_iter_ava, tan_dyn + ) # Ensure water level is up-to-date with bed level - s['zs'] = s['SWL'].copy() + s['zb'] = zb + s['gradh'] = grad_h + s['zs'] = s['SWL'] ix = (s['zb'] > s['zs']) s['zs'][ix] = s['zb'][ix] - return s + return s -def calc_gradients(zb, nx, ny, ds, dn, zne): - '''Calculates the downslope gradients in the bed that are needed for - avalanching module - - Parameters - ---------- - +@njit(cache=True) +def avalanche_loop(zb, zne, ds, dn, nx, ny, E, max_iter_ava, tan_dyn): + # Rewritten to use explicit loops and avoid numpy boolean indexing + for it in range(max_iter_ava): + # temporaries + grad_h_down = np.zeros((ny, nx, 2)) + flux_down = np.zeros((ny, nx, 2)) + slope_diff = np.zeros((ny, nx)) + grad_h = np.zeros((ny, nx)) + + # first calculate the downslope gradients to see if there is avalanching + # Compute downslope gradients grad_h_down (ny,nx,2), grad_h (ny,nx), and max_grad_h + # Initialize + max_grad_h = 0.0 + + # Directions: 0 => +X, 1 => +Y + for i in range(ny): + for j in range(nx): + center = zb[i, j] + + # +X direction + g0 = 0.0 + if 1:#(j > 0) and (j < nx - 1): + right = zb[i, (j + 1) % nx] + left = zb[i, (j - 1) % nx] + if not ((right > center) and (left > center)): + if right > left: + g0 = left - center + else: + g0 = center - right + + grad_h_down[i, j, 0] = g0 + + # +Y direction + g1 = 0.0 + if 1:#(i > 0) and (i < ny - 1): + down = zb[(i + 1) % ny, j] + up = zb[(i - 1) % ny, j] + if not ((down > center) and (up > center)): + if down > up: + g1 = up - center + else: + g1 = center - down + grad_h_down[i, j, 1] = g1 + + # normalize by grid spacing (assume ds, dn are 2D fields) + 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): + gh2 = grad_h_down[i, j, 0] * grad_h_down[i, j, 0] + grad_h_down[i, j, 1] * grad_h_down[i, j, 1] + # optional suppression near zne disabled + gh = np.sqrt(gh2) + grad_h[i, j] = gh + # derive maximum slope + if gh > max_grad_h: + max_grad_h = gh + + # ok now the gradients are calculated + # these are max_grad_h, grad_h, grad_h_down + # check for stopping criterion + if max_grad_h < tan_dyn: + break - Returns - ------- - np.ndarray - Downslope gradients in 4 different directions (nx*ny, 4) - ''' - - grad_h_down = np.zeros((ny,nx,4)) - - # Calculation of slope (positive x-direction) - grad_h_down[:,1:-1,0] = zb[:,1:-1] - zb[:,2:] - ix = zb[:,2:] > zb[:,:-2] - grad_h_down[:,1:-1,0][ix] = - (zb[:,1:-1][ix] - zb[:,:-2][ix]) - ix = np.logical_and(zb[:,2:]>zb[:,1:-1], zb[:,:-2]>zb[:,1:-1]) - grad_h_down[:,1:-1,0][ix] = 0. - - # Calculation of slope (positive y-direction) - grad_h_down[1:-1,:,1] = zb[1:-1,:] - zb[2:,:] - ix = zb[2:,:] > zb[:-2,:] - grad_h_down[1:-1,:,1][ix] = - (zb[1:-1,:][ix] - zb[:-2,:][ix]) - ix = np.logical_and(zb[2:,:]>zb[1:-1,:], zb[:-2,:]>zb[1:-1,:]) - grad_h_down[1:-1,:,1][ix] = 0. - - # Calculation of slope (negative x-direction) - grad_h_down[:,1:-1,2] = zb[:,1:-1] - zb[:,:-2] - ix = zb[:,:-2] > zb[:,2:] - grad_h_down[:,1:-1,2][ix] = - (zb[:,1:-1][ix] - zb[:,2:][ix]) - ix = np.logical_and(zb[:,:-2]>zb[:,1:-1], zb[:,2:]>zb[:,1:-1]) - grad_h_down[:,1:-1,2][ix] = 0. - - # Calculation of slope (negative y-direction) - grad_h_down[1:-1,:,3] = zb[1:-1,:] - zb[:-2,:] - ix = zb[:-2,:] > zb[2:,:] - grad_h_down[1:-1,:,3][ix] = - (zb[1:-1,:][ix] - zb[2:,:][ix]) - ix = np.logical_and(zb[:-2,:]>zb[1:-1,:], zb[2:,:]>zb[1:-1,:]) - grad_h_down[1:-1,:,3][ix] = 0. - - if ny == 1: - #1D interpretation - grad_h_down[:,0,:] = 0 - grad_h_down[:,-1,:] = 0 - - else: - # 2D interpretation - grad_h_down[:,0,:] = 0 - grad_h_down[:,-1,:] = 0 - grad_h_down[0,:,:] = 0 - grad_h_down[-1,:,:] = 0 + # if we continue we compute fluxes and update zb - grad_h_down[:,:,0] /= ds - grad_h_down[:,:,1] /= dn - grad_h_down[:,:,2] /= ds - grad_h_down[:,:,3] /= dn + # compute grad_h_nonerod and slope_diff per cell using explicit loops + for i in range(ny): + for j in range(nx): + grad_h_nonerod = (zb[i, j] - zne[i, j]) / ds[i, j] + if grad_h[i, j] > tan_dyn and grad_h_nonerod > 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) + + for i in range(ny): + for j in range(nx): + if grad_h[i, j] != 0.0: + 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] + + # Build q_in and q_out from 2-component flux representation + f_x = flux_down[:, :, 0] + f_y = flux_down[:, :, 1] - grad_h2 = 0.5*grad_h_down[:,:,0]**2 + 0.5*grad_h_down[:,:,1]**2 + 0.5*grad_h_down[:,:,2]**2 + 0.5*grad_h_down[:,:,3]**2 + # preserve sign: compute positive outgoing components per direction + out_east = np.maximum(f_x, 0.0) + out_south = np.maximum(f_y, 0.0) - if 0: #Sierd_com; to be changed in future release - ix = zb < zne + 0.005 - grad_h2[ix] = 0. + # average with neighbor contributions at faces (keeps sign info) + out_north = np.maximum(-f_y, 0.0) + out_west = np.maximum(-f_x, 0.0) - grad_h = np.sqrt(grad_h2) + q_out = out_east + out_west + out_south + out_north - max_grad_h = np.max(grad_h) + inc_west = np.zeros_like(f_x) + # from west neighbor (positive f_x of west cell) with periodic wrap + inc_west[:, 1:] = np.maximum(f_x[:, :-1], 0.0) + inc_west[:, 0] = np.maximum(f_x[:, -1], 0.0) - return max_grad_h, grad_h, grad_h_down - + inc_east = np.zeros_like(f_x) + # from east neighbor (negative f_x of east cell) with periodic wrap + inc_east[:, :-1] = np.maximum(-f_x[:, 1:], 0.0) + inc_east[:, -1] = np.maximum(-f_x[:, 0], 0.0) + + inc_north = np.zeros_like(f_y) + # from north neighbor (positive f_y of north cell) with periodic wrap + inc_north[1:, :] = np.maximum(f_y[:-1, :], 0.0) + inc_north[0, :] = np.maximum(f_y[-1, :], 0.0) + + inc_south = np.zeros_like(f_y) + # from south neighbor (negative f_y of south cell) with periodic wrap + inc_south[:-1, :] = np.maximum(-f_y[1:, :], 0.0) + inc_south[-1, :] = np.maximum(-f_y[0, :], 0.0) + + q_in = (inc_west + inc_east + inc_north + inc_south) + + zb += E * (q_in - q_out) + return zb, grad_h \ No newline at end of file From eea5f1fdd4c132e404f7495cbf4594687e855cef Mon Sep 17 00:00:00 2001 From: Sierd Date: Thu, 23 Oct 2025 23:03:27 +0200 Subject: [PATCH 02/10] updated avalanching.py to include non erodible layers. Not tested yet. --- aeolis/avalanching.py | 64 ++++++++++++++++++++++++------------------- 1 file changed, 36 insertions(+), 28 deletions(-) diff --git a/aeolis/avalanching.py b/aeolis/avalanching.py index 7957a351..a1d4381c 100644 --- a/aeolis/avalanching.py +++ b/aeolis/avalanching.py @@ -108,7 +108,7 @@ def avalanche(s, p): zb, grad_h = avalanche_loop( s['zb'].copy(), s['zne'], s['ds'], s['dn'], nx, ny, E, max_iter_ava, tan_dyn ) - + # Ensure water level is up-to-date with bed level s['zb'] = zb s['gradh'] = grad_h @@ -136,11 +136,13 @@ def avalanche_loop(zb, zne, ds, dn, nx, ny, E, max_iter_ava, tan_dyn): # Directions: 0 => +X, 1 => +Y for i in range(ny): for j in range(nx): - center = zb[i, j] - - # +X direction - g0 = 0.0 - if 1:#(j > 0) and (j < nx - 1): + # disable avalanching where zne >= zb + if zne[i, j] >= zb[i, j]: + continue + else: + center = zb[i, j] + # +X direction + g0 = 0.0 right = zb[i, (j + 1) % nx] left = zb[i, (j - 1) % nx] if not ((right > center) and (left > center)): @@ -148,12 +150,10 @@ def avalanche_loop(zb, zne, ds, dn, nx, ny, E, max_iter_ava, tan_dyn): g0 = left - center else: g0 = center - right + grad_h_down[i, j, 0] = g0 / ds[i, j] - grad_h_down[i, j, 0] = g0 - - # +Y direction - g1 = 0.0 - if 1:#(i > 0) and (i < ny - 1): + # +Y direction + g1 = 0.0 down = zb[(i + 1) % ny, j] up = zb[(i - 1) % ny, j] if not ((down > center) and (up > center)): @@ -161,17 +161,17 @@ def avalanche_loop(zb, zne, ds, dn, nx, ny, E, max_iter_ava, tan_dyn): g1 = up - center else: g1 = center - down - grad_h_down[i, j, 1] = g1 + grad_h_down[i, j, 1] = g1 / dn[i, j] # normalize by grid spacing (assume ds, dn are 2D fields) - 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] + # 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): + # for i in range(ny): + # for j in range(nx): gh2 = grad_h_down[i, j, 0] * grad_h_down[i, j, 0] + grad_h_down[i, j, 1] * grad_h_down[i, j, 1] # optional suppression near zne disabled gh = np.sqrt(gh2) @@ -186,22 +186,22 @@ def avalanche_loop(zb, zne, ds, dn, nx, ny, E, max_iter_ava, tan_dyn): if max_grad_h < tan_dyn: break - # if we continue we compute fluxes and update zb + # we continue to compute fluxes and update zb # compute grad_h_nonerod and slope_diff per cell using explicit loops for i in range(ny): for j in range(nx): - grad_h_nonerod = (zb[i, j] - zne[i, j]) / ds[i, j] - if grad_h[i, j] > tan_dyn and grad_h_nonerod > 0.0: + # 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) + # elif grad_h_nonerod < (grad_h[i, j] - tan_dyn): + # slope_diff[i, j] = np.tanh(grad_h_nonerod) - for i in range(ny): - for j in range(nx): + # for i in range(ny): + # for j in range(nx): if grad_h[i, j] != 0.0: - 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] # Build q_in and q_out from 2-component flux representation f_x = flux_down[:, :, 0] @@ -239,6 +239,14 @@ def avalanche_loop(zb, zne, ds, dn, nx, ny, E, max_iter_ava, tan_dyn): q_in = (inc_west + inc_east + inc_north + inc_south) - zb += E * (q_in - q_out) + # # check mass balance in the presence of non-erodible layer + 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) return zb, grad_h \ No newline at end of file From 24c527263da81b39683e63d52050680496058a04 Mon Sep 17 00:00:00 2001 From: Sierd de Vries Date: Sat, 25 Oct 2025 10:40:25 +0200 Subject: [PATCH 03/10] Remove theta_stat and update theta_dyn description --- aeolis/constants.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/aeolis/constants.py b/aeolis/constants.py index 4f8a7b62..6bb0fb75 100644 --- a/aeolis/constants.py +++ b/aeolis/constants.py @@ -99,8 +99,7 @@ 'ustar0', # [m/s] Initial shear velocity (without perturbation) 'zsep', # [m] Z level of polynomial that defines the separation bubble 'hsep', # [m] Height of separation bubbel = difference between z-level of zsep and of the bed level zb - 'theta_stat', # [degrees] Updated, spatially varying static angle of repose - 'theta_dyn', # [degrees] Updated, spatially varying dynamic angle of repose + 'theta_dyn', # [degrees] spatially varying dynamic angle of repose for avalanching 'rhoveg', # [-] Vegetation cover 'drhoveg', # Change in vegetation cover 'hveg', # [m] height of vegetation From 29a6d909b4ad634565d85b5044cdb4c238745164 Mon Sep 17 00:00:00 2001 From: Sierd Date: Sat, 25 Oct 2025 10:42:45 +0200 Subject: [PATCH 04/10] removed static angle of repose. Could be added back later if needed. --- aeolis/avalanching.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/aeolis/avalanching.py b/aeolis/avalanching.py index a1d4381c..b23fd3e5 100644 --- a/aeolis/avalanching.py +++ b/aeolis/avalanching.py @@ -62,10 +62,10 @@ def angele_of_repose(s,p): # comment Lisa: dependence on moisture content is not yet implemented # Can we do something with theta dependent on vegetation cover (larger rhoveg = larger theta?) - theta_stat = p['theta_stat'] + # theta_stat = p['theta_stat'] theta_dyn = p['theta_dyn'] - s['theta_stat'] = theta_stat + # s['theta_stat'] = theta_stat s['theta_dyn'] = theta_dyn return s @@ -92,12 +92,15 @@ def avalanche(s, p): Spatial grids ''' + + print('testing avalanching process2') if p['process_avalanche']: nx = p['nx'] + 1 ny = p['ny'] + 1 - # parameters - tan_stat = np.tan(np.deg2rad(s['theta_stat'])) + # parameters - only dynamic angle used in loop for now. + # Static angle can be used for more complex criterions in later + # tan_stat = np.tan(np.deg2rad(s['theta_stat'])) tan_dyn = np.tan(np.deg2rad(s['theta_dyn'])) E = 0.1 From e83ea56ed1731656195369431fa77b49abde311c Mon Sep 17 00:00:00 2001 From: Sierd de Vries Date: Sat, 25 Oct 2025 11:12:52 +0200 Subject: [PATCH 05/10] Update aeolis/avalanching.py copilot suggested changes on memory allocation Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- aeolis/avalanching.py | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/aeolis/avalanching.py b/aeolis/avalanching.py index b23fd3e5..d9e08e18 100644 --- a/aeolis/avalanching.py +++ b/aeolis/avalanching.py @@ -124,12 +124,17 @@ def avalanche(s, p): @njit(cache=True) def avalanche_loop(zb, zne, ds, dn, nx, ny, E, max_iter_ava, tan_dyn): # Rewritten to use explicit loops and avoid numpy boolean indexing + # Allocate temporaries once + grad_h_down = np.zeros((ny, nx, 2)) + flux_down = np.zeros((ny, nx, 2)) + slope_diff = np.zeros((ny, nx)) + grad_h = np.zeros((ny, nx)) for it in range(max_iter_ava): - # temporaries - grad_h_down = np.zeros((ny, nx, 2)) - flux_down = np.zeros((ny, nx, 2)) - slope_diff = np.zeros((ny, nx)) - grad_h = np.zeros((ny, nx)) + # Reset temporaries to zero + grad_h_down.fill(0) + flux_down.fill(0) + slope_diff.fill(0) + grad_h.fill(0) # first calculate the downslope gradients to see if there is avalanching # Compute downslope gradients grad_h_down (ny,nx,2), grad_h (ny,nx), and max_grad_h From b7c17046d2f42c4374f302f22df9ecfb3384e37f Mon Sep 17 00:00:00 2001 From: Sierd de Vries Date: Sat, 25 Oct 2025 11:14:19 +0200 Subject: [PATCH 06/10] Update aeolis/avalanching.py added boundary definition. Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- aeolis/avalanching.py | 39 +++++++++++++++++++++++---------------- 1 file changed, 23 insertions(+), 16 deletions(-) diff --git a/aeolis/avalanching.py b/aeolis/avalanching.py index d9e08e18..cfc32d63 100644 --- a/aeolis/avalanching.py +++ b/aeolis/avalanching.py @@ -151,25 +151,32 @@ def avalanche_loop(zb, zne, ds, dn, nx, ny, E, max_iter_ava, tan_dyn): center = zb[i, j] # +X direction g0 = 0.0 - right = zb[i, (j + 1) % nx] - left = zb[i, (j - 1) % nx] - if not ((right > center) and (left > center)): - if right > left: - g0 = left - center - else: - g0 = center - right - grad_h_down[i, j, 0] = g0 / ds[i, j] + # Handle boundaries: set gradient to zero at edges + if j == 0 or j == nx - 1: + grad_h_down[i, j, 0] = 0.0 + else: + right = zb[i, j + 1] + left = zb[i, j - 1] + if not ((right > center) and (left > center)): + if right > left: + g0 = left - center + else: + g0 = center - right + grad_h_down[i, j, 0] = g0 / ds[i, j] # +Y direction g1 = 0.0 - down = zb[(i + 1) % ny, j] - up = zb[(i - 1) % ny, j] - if not ((down > center) and (up > center)): - if down > up: - g1 = up - center - else: - g1 = center - down - grad_h_down[i, j, 1] = g1 / dn[i, j] + if i == 0 or i == ny - 1: + grad_h_down[i, j, 1] = 0.0 + else: + down = zb[i + 1, j] + up = zb[i - 1, j] + if not ((down > center) and (up > center)): + if down > up: + g1 = up - center + else: + g1 = center - down + grad_h_down[i, j, 1] = g1 / dn[i, j] # normalize by grid spacing (assume ds, dn are 2D fields) # for i in range(ny): From d37d6d090f22051c5e91673b53ce5345cf668e53 Mon Sep 17 00:00:00 2001 From: Sierd de Vries Date: Sat, 25 Oct 2025 11:43:16 +0200 Subject: [PATCH 07/10] Update aeolis/avalanching.py copilot suggestion to filter for ne Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- aeolis/avalanching.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/aeolis/avalanching.py b/aeolis/avalanching.py index cfc32d63..58603e74 100644 --- a/aeolis/avalanching.py +++ b/aeolis/avalanching.py @@ -255,7 +255,7 @@ def avalanche_loop(zb, zne, ds, dn, nx, ny, E, max_iter_ava, tan_dyn): q_in = (inc_west + inc_east + inc_north + inc_south) # # check mass balance in the presence of non-erodible layer - if np.any((E * (q_in - q_out)) < (zne-zb)): + 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 From 9baab2418b4ece480d90744f7b9542c141eeb21e Mon Sep 17 00:00:00 2001 From: Sierd Date: Sat, 25 Oct 2025 12:25:50 +0200 Subject: [PATCH 08/10] Update to work with NE layer in avalanche module Removed requirement for ne_file when avalanching is enabled --- aeolis/avalanching.py | 31 ++++--------------------------- aeolis/inout.py | 7 ------- 2 files changed, 4 insertions(+), 34 deletions(-) diff --git a/aeolis/avalanching.py b/aeolis/avalanching.py index 58603e74..2ce3772e 100644 --- a/aeolis/avalanching.py +++ b/aeolis/avalanching.py @@ -178,17 +178,8 @@ def avalanche_loop(zb, zne, ds, dn, nx, ny, E, max_iter_ava, tan_dyn): g1 = center - down grad_h_down[i, j, 1] = g1 / dn[i, j] - # normalize by grid spacing (assume ds, dn are 2D fields) - # 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 gh2 = grad_h_down[i, j, 0] * grad_h_down[i, j, 0] + grad_h_down[i, j, 1] * grad_h_down[i, j, 1] - # optional suppression near zne disabled gh = np.sqrt(gh2) grad_h[i, j] = gh # derive maximum slope @@ -206,15 +197,8 @@ def avalanche_loop(zb, zne, ds, dn, nx, ny, E, max_iter_ava, tan_dyn): # compute grad_h_nonerod and slope_diff per cell using explicit loops for i in range(ny): for j in range(nx): - # 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: + if grad_h[i, j] > tan_dyn: 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) - - # for i in range(ny): - # for j in range(nx): - if grad_h[i, j] != 0.0: 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] @@ -254,14 +238,7 @@ def avalanche_loop(zb, zne, ds, dn, nx, ny, E, max_iter_ava, tan_dyn): q_in = (inc_west + inc_east + inc_north + inc_south) - # # 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 without non-erodible layer + zb += E * (q_in - q_out) return zb, grad_h \ No newline at end of file diff --git a/aeolis/inout.py b/aeolis/inout.py index a79cd3ad..e29dadd3 100644 --- a/aeolis/inout.py +++ b/aeolis/inout.py @@ -118,13 +118,6 @@ def read_configfile(configfile, parse_files=True, load_defaults=True): if 'nsavetimes' in p and not p['nsavetimes']: p['nsavetimes'] = int(p['dzb_interval']/p['dt']) - # catch some incompatible parameter combinations. - if (p['ne_file'] is None) and (p['process_avalanche'] == True): - print('Please provide a valid ne_file path in the configuration file when using Avalanching.') - print('Code does not proceed until this is provided.') - print('Hint: If you do not have a ne_file, you can create one with all zeros, with the same dimensions as your grid.') - exit("Exiting due to error") - return p From 29231b647280668af0ab08e018a672a185e1f6a7 Mon Sep 17 00:00:00 2001 From: Sierd de Vries Date: Sat, 25 Oct 2025 12:50:28 +0200 Subject: [PATCH 09/10] Update aeolis/avalanching.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- aeolis/avalanching.py | 1 - 1 file changed, 1 deletion(-) diff --git a/aeolis/avalanching.py b/aeolis/avalanching.py index 58603e74..4a2cdc99 100644 --- a/aeolis/avalanching.py +++ b/aeolis/avalanching.py @@ -93,7 +93,6 @@ def avalanche(s, p): ''' - print('testing avalanching process2') if p['process_avalanche']: nx = p['nx'] + 1 ny = p['ny'] + 1 From dc239ab174c2a15ddafb8dda5c86e643e824b193 Mon Sep 17 00:00:00 2001 From: Sierd de Vries Date: Sat, 25 Oct 2025 12:52:39 +0200 Subject: [PATCH 10/10] Update aeolis/constants.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- aeolis/constants.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/aeolis/constants.py b/aeolis/constants.py index 6bb0fb75..eae691bc 100644 --- a/aeolis/constants.py +++ b/aeolis/constants.py @@ -98,7 +98,7 @@ 'ustarn', # [m/s] Component of shear velocity in y-direction by wind 'ustar0', # [m/s] Initial shear velocity (without perturbation) 'zsep', # [m] Z level of polynomial that defines the separation bubble - 'hsep', # [m] Height of separation bubbel = difference between z-level of zsep and of the bed level zb + 'hsep', # [m] Height of separation bubble = difference between z-level of zsep and of the bed level zb 'theta_dyn', # [degrees] spatially varying dynamic angle of repose for avalanching 'rhoveg', # [-] Vegetation cover 'drhoveg', # Change in vegetation cover