diff --git a/aeolis/constants.py b/aeolis/constants.py index eae691bc..c36a727f 100644 --- a/aeolis/constants.py +++ b/aeolis/constants.py @@ -244,8 +244,15 @@ 'k' : 0.001, # [m] Bed roughness 'L' : 100., # [m] Typical length scale of dune feature (perturbation) 'l' : 10., # [m] Inner layer height (perturbation) - 'c_b' : 0.2, # [-] Slope at the leeside of the separation bubble # c = 0.2 according to Durán 2010 (Sauermann 2001: c = 0.25 for 14 degrees) - 'mu_b' : 30, # [deg] Minimum required slope for the start of flow separation + + # Separation bubble parameters + 'sep_auto_tune' : True, # [-] Boolean for automatic tuning of separation bubble parameters based on characteristic length scales + 'sep_look_dist' : 50., # [m] Flow separation: Look-ahead distance for upward curvature anticipation + 'sep_k_press_up' : 0.05, # [-] Flow separation: Press-up curvature + 'sep_k_crit_down' : 0.18, # [1/m] Flow separation: Maximum downward curvature + 'sep_s_crit' : 0.18, # [-] Flow separation: Critical bed slope below which reattachment is forced + 'sep_s_leeside' : 0.25, # [-] Maximum downward leeside slope of the streamline + 'buffer_width' : 10, # [m] Width of the bufferzone around the rotational grid for wind perturbation 'sep_filter_iterations' : 0, # [-] Number of filtering iterations on the sep-bubble (0 = no filtering) 'zsep_y_filter' : False, # [-] Boolean for turning on/off the filtering of the separation bubble in y-direction diff --git a/aeolis/separation.py b/aeolis/separation.py new file mode 100644 index 00000000..06625b03 --- /dev/null +++ b/aeolis/separation.py @@ -0,0 +1,258 @@ +""" +Streamline-based separation bubble model for AeoLiS + +This module provides a unified 1D/2D separation bubble computation based on +a streamline curvature model. The physics is contained in a single numba- +compiled 1D core. For 2D grids (FFT shear mode), the 1D streamline model +is applied to each wind-aligned row. + +Public Entry Point: + compute_separation(z_bed, dx, wind_sign, p) + +Internal: + _compute_separation_1d(...) + _compute_separation_2d(...) + _streamline_core_1d(...) # numba auto-compiled + +""" + +import numpy as np +from numba import njit + +# Wrapper for 1D vs 2D +def compute_separation(p, z_bed, dx, udir=0): + + # Get all parameters + look_dist = p['sep_look_dist'] + k_press_up = p['sep_k_press_up'] + k_crit_down = p['sep_k_crit_down'] + s_crit = p['sep_s_crit'] + s_leeside = p['sep_s_leeside'] + + # Determine size + ny_, nx = z_bed.shape + + # Flip the bed for negative wind directions (only in 1D) + if udir < 0 and ny_ == 0: + z = z[::-1] + + # Run streamline computation once in 1D + if ny_ == 0: + zsep = _streamline_core_1d( + z_bed, dx, look_dist, # Base input + k_press_up, k_crit_down, # Curviture + s_crit, s_leeside) # Leeside slopes + + + else: # Run streamline computation over every transect in 2D (Numba-compiled over 1D) + zsep = _compute_separation_2d( + z_bed, dx, look_dist, + k_press_up, k_crit_down, + s_crit, s_leeside) + + # Flip back for the 1D case + if udir < 0 and ny_ == 0: + zsep = zsep[::-1] + + return zsep + + +@njit +def _compute_separation_2d( + z_bed, dx, look_dist, + k_press_up, k_crit_down, + s_crit, s_leeside): + + ny_, nx = z_bed.shape + zsep = np.zeros_like(z_bed) + + for j in range(ny_): + + row = z_bed[j, :] + + # Quick skip: slope + curvature checks + skip = True + for i in range(1, nx-1): + s = (row[i] - row[i-1]) / dx + s_prev = (row[i-1] - row[i-2]) / dx + ds = s_prev - s # curvature approx + + if (s < -s_crit) or (ds < -k_crit_down): + skip = False + break + + if skip: + zsep[j, :] = row + else: + zsep[j, :] = _streamline_core_1d( + row, dx, look_dist, + k_press_up, k_crit_down, + s_crit, s_leeside) + + return zsep + + +# Streamline core +@njit +def _streamline_core_1d( + z_bed, dx, look_dist, + k_press_up, k_crit_down, + s_crit, s_leeside): + + # Initialize size and z (= streamline) + N = z_bed.size + z = z_bed.copy() + + # Loop over all points in windward direction + for i in range(1, N - 1): + + # Compute slope of streamline + s_str = (z[i] - z[i-1]) / dx + gap = z[i] - z_bed[i] + + # Compute slope of bed + s_bed = (z_bed[i] - z_bed[i-1]) / dx + s_bed_next = (z_bed[i+1] - z_bed[i]) / dx + ds_bed = s_bed_next - s_bed + + # Determine how far to look ahead for upward curviture + look_n = int(look_dist / dx) + i_end = min(i + look_n, N - 1) + + # Initialize maximum required curviture (k_req_max) and the resulting slope (v_z) + k_req_max = 0.0 + v_z = None + k_press_down_base = 0.05 + + # ----- 1. UPWARD CURVATURE ANTICIPATION ----- + + # Start looking forward for upward curvature anticipation + if i_end > i: + for j in range(i + 1, i_end + 1): + + # Compute difference in distance, height and slope downwind (forward = j) + dxj = (j - i) * dx # Total distance between current (i) and forward (j) + dzj = z[j] - z[i] # Total height difference between current (i) and forward (j) + sbj = (z[j] - z[j-1]) / dx # Slope at forward (j) + + # Required slope (s) to close height difference from (i) to (j) + s_req_height = dzj / dxj + + # Required curviture (k) to close slope difference for height (k_req_height) and slope (k_req_slope) + k_req_height = (s_req_height - s_str) / dxj + k_req_slope = (sbj - s_str) / dxj + + # Prevent that the streamline "overshoots" the bed level due a too steep curviture + z_est = z[i] + s_str * dxj + 0.5 * k_req_slope * dxj * dxj + if z_est > z[j]: + k_req_slope = 0.0 + + # Required distance to reach either height or slope + d_req_height = np.sqrt(2 * max(0.0, dzj - s_str * dxj) / k_press_up) if k_req_height > 0 else 0. + d_req_slope = (sbj - s_str) / k_press_up if k_req_slope > 0 else 0. + d_req = max(d_req_height, d_req_slope) + + # Check whether d_req is within reach (pass if dxj > d_req) + # i.e., if d_req < dxj; it is not necessary to bend upward yet + if d_req > dxj: + k_req_max = max(k_req_max, k_req_slope, k_req_height) + + # Apply curvature anticipation + if k_req_max >= k_press_up: + v_z = s_str + k_req_max + + # ----- 2. DOWNWARD CURVATURE BY DE- and RE-ATTACHMENT ----- + + # Don't apply downward curvature if we are doing upward anticipation + if v_z is None: + + # Check if we are at the first point of detachment + if gap < 1e-6 and (ds_bed < -k_crit_down or s_str < -s_crit): + + # The downward curvature (k_press_down) is based on an estimated attachment length + # This attachment length depends on geometry (L_attach is roughly equal to 6 * H_hill) + r_LH = 6. + s_L = 1 / r_LH + + # An addition to the dettachment length is added based on the distance its takes to bend to leeslope + L_curv = max((s_str + s_leeside) / k_press_down_base, dx) + i_curv = int(i + L_curv / dx) + + # Loop forward going downward with s_L per m and track where it intersects the bed + bed_intersect = [] + for j in range(i_curv, N): + zj = z[i] - s_L * (j - i_curv) * dx + zj_prev = z[i] - s_L * (j - i_curv - 1) * dx + if zj < z_bed[j] and zj_prev >= z_bed[j]: + bed_intersect.append(j) + + # Double crossing could be caused due to a small hill on the leeslide + # If the bed is crossed multiple times, the last one is used + if bed_intersect: + n_intersect = bed_intersect[-1] + else: + n_intersect = N-1 + + # Estimate the height of the feature by taking the max and min bed levels + h_gap = np.max(z_bed[i:n_intersect+1])- np.min(z_bed[i:n_intersect+1]) + + # Estimate attachment length and subsequent downward curvature + L_attach = h_gap * r_LH + L_curv + k_press_down = 1.5 / L_attach + + else: + k_press_down = k_press_down_base + + # Apply downward curvature + if ds_bed < -k_crit_down or s_bed < -s_crit or gap > 0.0: + + # An exponential term drives s → -s_leeside and bend steep upward slopes faster + f_converge = np.exp(1.4*(max(s_str + s_leeside, 0)**0.85)) - 1 + v_z = s_str - k_press_down * f_converge + + # Remains when no anticipation or downward curvature is applied + else: + v_z = s_bed_next + + # ----- 3. ADVANCE STREAMLINE ----- + z[i+1] = max(z[i] + v_z * dx, z_bed[i+1]) + + return z + + +def set_sep_params(p, s): + """ + Set separation–bubble parameters based on perturbation–theory + scaling (Kroy-type). Values are stored inside the main parameter + dictionary 'p'. Returns updated p. + + Inputs + ------ + p : dict + Global model parameter dictionary. + s : dict + Shear fields (contains L, l, z0 via WindShear object). + + Notes + ------ + Upward curvature (k_press_up) follows Kroy-type perturbation theory. + Downward curvature is taken as a fraction of the upward limit. + Other parameters remain simple fixed constants. + """ + + # Characteristic length scales + f_min = 0.3 # [-] Minimum ratio of tau to tau_0 + # gamma_down = 0.4 # [-] Ratio of k_down to k_up + + # --- Upward curvature limit (perturbation theory) --- + p['sep_k_press_up'] = (1 - f_min) * (2*np.pi / (np.sqrt(p['L']) * 40)) + + # # --- Downward curvature limit (softer than upward) --- + # p['sep_k_press_down'] = gamma_down * p['sep_k_press_up'] + + # print(p['sep_k_press_down']) + + # # --- Curvature threshold for detachment trigger --- + # p['sep_curv_limit_down'] = np.maximum(p['sep_s_crit'] / 0.5, 1.5 * p['sep_k_press_down']) + + return p diff --git a/aeolis/shear.py b/aeolis/shear.py index 5b2ae4e0..772083ee 100644 --- a/aeolis/shear.py +++ b/aeolis/shear.py @@ -38,7 +38,7 @@ # package modules from aeolis.utils import * - +from aeolis.separation import compute_separation # initialize logger logger = logging.getLogger(__name__) @@ -142,10 +142,8 @@ def __init__(self, x, y, z, dx, dy, L, l, z0, self.z0 = z0 - def __call__(self, x, y, z, taux, tauy, u0, udir, - process_separation, c, mu_b, - taus0, taun0, sep_filter_iterations, zsep_y_filter, - plot=False): + def __call__(self, p, x, y, z, u0, udir, + taux, tauy, taus0, taun0, plot=False): '''Compute wind shear for given wind speed and direction @@ -222,10 +220,9 @@ def __call__(self, x, y, z, taux, tauy, u0, udir, # rotate to horizontal computational grid and add to tau0 # ===================================================================== - # Compute separation bubble - if process_separation: - zsep = self.separation(c, mu_b, sep_filter_iterations, zsep_y_filter) - z_origin = gc['z'].copy() + if p['process_separation']: + z_origin = gc['z'][:].copy() + zsep = compute_separation(p, gc['z'], gc['dx']) gc['z'] = np.maximum(gc['z'], zsep) # Compute wind shear stresses on computational grid @@ -238,7 +235,7 @@ def __call__(self, x, y, z, taux, tauy, u0, udir, gc['taux'] = np.maximum(gc['taux'], 0.) # Compute the influence of the separation on the shear stress - if process_separation: + if p['process_separation']: gc['hsep'] = gc['z'] - z_origin self.separation_shear(gc['hsep']) @@ -264,7 +261,7 @@ def __call__(self, x, y, z, taux, tauy, u0, udir, gi['taux'] = self.interpolate(gc['x'], gc['y'], gc['taux'], gi['x'], gi['y'], taus0) gi['tauy'] = self.interpolate(gc['x'], gc['y'], gc['tauy'], gi['x'], gi['y'], taun0) - if process_separation: + if p['process_separation']: gi['hsep'] = self.interpolate(gc['x'], gc['y'], gc['hsep'], gi['x'], gi['y'], 0. ) # Final plots and lay-out @@ -368,158 +365,6 @@ def set_computational_grid(self, udir): return self - - - def separation(self, c, mu_b, sep_filter_iterations, zsep_y_filter): - - # Initialize grid and bed dimensions - gc = self.cgrid - - x = gc['x'] - y = gc['y'] - z = gc['z'] - - nx = len(gc['z'][1]) - ny = len(gc['z'][0]) - dx = gc['dx'] - dy = gc['dy'] - - # Initialize arrays - - dzx = np.zeros(gc['z'].shape) - - dzdx0 = np.zeros(gc['z'].shape) - dzdx1 = np.zeros(gc['z'].shape) - - stall = np.zeros(gc['z'].shape) - bubble = np.zeros(gc['z'].shape) - - k = np.array(range(0, nx)) - - zsep = np.zeros(z.shape) # total separation bubble - zsep_new = np.zeros(z.shape) # first-oder separation bubble surface - - zfft = np.zeros((ny,nx), dtype=complex) - - # Compute bed slope angle in x-dir - dzx[:,:-2] = np.rad2deg(np.arctan((z[:,2:]-z[:,:-2])/(2.*dx))) - dzx[:,-2] = dzx[:,-3] - dzx[:,-1] = dzx[:,-2] - - # Determine location of separation bubbles - '''Separation bubble exist if bed slope angle (lee side) - is larger than max angle that wind stream lines can - follow behind an obstacle (mu_b = ..)''' - - stall += np.logical_and(abs(dzx) > mu_b, dzx < 0.) - - stall[:,1:-1] += np.logical_and(stall[:,1:-1]==0, stall[:,:-2]>0., stall[:,2:]>0.) - - # Define separation bubble - bubble[:,:-1] = (stall[:,:-1] == 0.) * (stall[:,1:] > 0.) - - # Better solution for cleaner separation bubble, but no working Barchan dune (yet) - p = 1 - bubble[:,p:] = bubble[:,:-p] - bubble[:,-p:] = 0 - - bubble = bubble.astype(int) - - # Count separation bubbles - n = np.sum(bubble) - bubble_n = np.asarray(np.where(bubble == True)).T - - - # Walk through all separation bubbles and determine polynoms - j = 9999 - for k in range(0, n): - - i = bubble_n[k,1] - j = bubble_n[k,0] - - #Bart: check for negative wind direction - if np.sum(gc['taux']) >= 0: - idir = 1 - else: - idir = -1 - - ix_neg = (dzx[j, i+idir*5:] >= 0) # i + 5?? - - if np.sum(ix_neg) == 0: - zbrink = z[j,i] # z level of brink at z(x0) - else: - zbrink = z[j,i] - z[j,i+idir*5+idir*np.where(ix_neg)[0][0]] - - # Better solution and cleaner separation bubble, but no working Barchan dune (yet) - dzdx0 = (z[j,i] - z[j,i-3]) / (3.*dx) - - a = dzdx0 / c - ls = np.minimum(np.maximum((3.*zbrink/(2.*c) * (1. + a/4. + a**2/8.)), 0.1), 200.) - - a2 = -3 * zbrink/ls**2 - 2 * dzdx0 / ls - a3 = 2 * zbrink/ls**3 + dzdx0 / ls**2 - - i_max = min(i+int(ls/dx)+1,int(nx-1)) - - if idir == 1: - xs = x[j,i:i_max] - x[j,i] - else: - xs = -(x[j,i:i_max] - x[j,i]) - - zsep_new[j,i:i_max] = (a3*xs**3 + a2*xs**2 + dzdx0*xs + z[j,i]) - - # Choose maximum of bedlevel, previous zseps and new zseps - zsep[j,:] = np.maximum.reduce([z[j,:], zsep[j,:], zsep_new[j,:]]) - - for filter_iter in range(sep_filter_iterations): - - zsep_new = np.zeros(zsep.shape) - - Cut = 1.5 - dk = 2.0 * np.pi / (np.max(x)) - zfft[j,:] = np.fft.fft(zsep[j,:]) - zfft[j,:] *= np.exp(-(dk*k*dx)**2/(2.*Cut**2)) - zsep_fft = np.real(np.fft.ifft(zfft[j,:])) - - if np.sum(ix_neg) == 0: - zbrink = zsep_fft[i] - else: - zbrink = zsep_fft[i] - zsep_fft[i+idir*5+idir*np.where(ix_neg)[0][0]] - - # First order polynom - dzdx1 = (zsep_fft[i] - zsep_fft[i-3])/(3.*dx) - - a = dzdx1 / c - - ls = np.minimum(np.maximum((3.*zbrink/(2.*c) * (1. + a/4. + a**2/8.)), 0.1), 200.) - - - a2 = -3 * zbrink/ls**2 - 2 * dzdx1 / ls - a3 = 2 * zbrink/ls**3 + dzdx1 / ls**2 - - i_max1 = min(i+idir*int(ls/dx),int(nx-1)) - - if idir == 1: - xs1 = x[j,i:i_max1] - x[j,i] - else: - xs1 = -(x[j,i:i_max1] - x[j,i]) - - zsep_new[j, i:i_max1] = (a3*xs1**3 + a2*xs1**2 + dzdx1*xs1 + zbrink) - - # Pick the maximum seperation bubble hieght at all locations - zsep[j,:] = np.maximum.reduce([z[j,:], zsep[j,:], zsep_new[j,:]]) - - - # Smooth surface of separation bubbles over y direction - if zsep_y_filter: - zsep = ndimage.gaussian_filter1d(zsep, sigma=0.2, axis=0) - - #Correct for any seperation bubbles that are below the bed surface following smoothing - ilow = zsep < z - zsep[ilow] = z[ilow] - - return zsep - def compute_shear(self, u0, nfilter=(1., 2.)): '''Compute wind shear perturbation for given free-flow wind diff --git a/aeolis/wind.py b/aeolis/wind.py index db82c030..0f9b6fb2 100644 --- a/aeolis/wind.py +++ b/aeolis/wind.py @@ -36,7 +36,7 @@ # package modules import aeolis.shear from aeolis.utils import * - +from aeolis.separation import compute_separation, set_sep_params # initialize logger logger = logging.getLogger(__name__) @@ -71,6 +71,11 @@ def initialize(s, p): z0 = calculate_z0(p, s) if p['process_shear']: + + # Get separation parameters based on Perturbation theory settings + if p['process_separation'] and p['sep_auto_tune']: + p = set_sep_params(p, s) + if p['ny'] > 0: if p['method_shear'] == 'fft': s['shear'] = aeolis.shear.WindShear(s['x'], s['y'], s['zb'], @@ -219,10 +224,11 @@ def calculate_z0(p, s): return z0 +# Compute shear velocity field (including separation) def shear(s,p): - # Compute shear velocity field (including separation) + # Quasi-2D shear model (DUNA (?) approach) if 'shear' in s.keys() and p['process_shear'] and p['ny'] > 0 and p['method_shear'] == 'duna2d': shear_params = {'alfa': 3, 'beta': 1, 'c': p['c_b'], 'mu_b': p['mu_b'], 'sep_filter_iterations': p['sep_filter_iterations'], 'zsep_y_filter': p['zsep_y_filter'], 'process_separation': p['process_separation'], 'tau_sep': 0.5, 'slope': 0.2, 'rhoa': p['rhoa'], 'shear_type': p['method_shear']} @@ -235,6 +241,7 @@ def shear(s,p): s = stress_velocity(s, p) + # Quasi-2D shear model (1D analytical over 2D grid) elif 'shear' in s.keys() and p['process_shear'] and p['ny'] > 0 and p['method_shear'] == 'quasi2d': z0 = calculate_z0(p, s) @@ -248,20 +255,16 @@ def shear(s,p): s['taun'] = - s['tau'] * np.cos((-p['alfa'] + s['udir']) / 180. * np.pi) s = stress_velocity(s, p) + # 2D shear model (FFT approach) elif 'shear' in s.keys() and p['process_shear'] and p['ny'] > 0: - s['shear'](x=s['x'], y=s['y'], z=s['zb'], - taux=s['taus'], tauy=s['taun'], - u0=s['uw'][0,0], udir=s['udir'][0,0], - process_separation = p['process_separation'], - c = p['c_b'], - mu_b = p['mu_b'], - taus0 = s['taus0'][0,0], taun0 = s['taun0'][0,0], - sep_filter_iterations=p['sep_filter_iterations'], - zsep_y_filter=p['zsep_y_filter']) + # Run shear model + s['shear'](p=p, x=s['x'], y=s['y'], z=s['zb'], u0=s['uw'][0,0], udir=s['udir'][0,0], + taux=s['taus'], tauy=s['taun'], taus0 = s['taus0'][0,0], taun0 = s['taun0'][0,0]) + + # Get shear stress (tau) from module, compute magnitude and transform to shear velocity s['taus'], s['taun'] = s['shear'].get_shear() s['tau'] = np.hypot(s['taus'], s['taun']) - s = stress_velocity(s,p) # Returns separation surface @@ -269,15 +272,22 @@ def shear(s,p): s['hsep'] = s['shear'].get_separation() s['zsep'] = s['hsep'] + s['zb'] - + # 1D shear model elif p['process_shear'] and p['ny'] == 0: #NTC - Added in 1D only capabilities s = compute_shear1d(s, p) s = stress_velocity(s, p) if p['process_separation']: - zsep = separation1d(s, p) - s['zsep'] = zsep + dx = s['ds'][0, 0] + z_bed = s['zb'][0, :] # 1D bed + udir = s['udir'][0, 0] # wind-aligned direction + + # Compute separation bubble + zsep = compute_separation(p, z_bed, dx, udir) + s['zsep'] = zsep[np.newaxis, :] # shape back to (1, nx) s['hsep'] = s['zsep'] - s['zb'] + + # Compute influence of searation bubble on shear tau_sep = 0.5 slope = 0.2 # according to Durán 2010 (Sauermann 2001: c = 0.25 for 14 degrees) delta = 1. / (slope * tau_sep) @@ -286,23 +296,6 @@ def shear(s,p): s['taun'] *= zsepdelta s = stress_velocity(s, p) - # if p['process_nelayer']: - # if p['th_nelayer']: - - # ustar = s['ustar'].copy() - # ustars = s['ustars'].copy() - # ustarn = s['ustarn'].copy() - - # s['zne'][:,:] = p['ne_file'] - - # ix = s['zb'] <= s['zne'] - # s['ustar'][ix] = np.maximum(0., s['ustar'][ix] - (s['zne'][ix]-s['zb'][ix])* (1/p['layer_thickness']) * s['ustar'][ix]) - - # ix = ustar != 0. - # s['ustars'][ix] = s['ustar'][ix] * (ustars[ix] / ustar[ix]) - # s['ustarn'][ix] = s['ustar'][ix] * (ustarn[ix] / ustar[ix]) - - return s def velocity_stress(s, p):