From a7d7f990a4be7d367ba614bb1e068c20ec88e683 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 12 Mar 2026 15:50:13 +0000 Subject: [PATCH 1/2] Initial plan From aeda8cc7713781f78d32c1d7d9f7d7d852a0ea4c Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 12 Mar 2026 16:25:36 +0000 Subject: [PATCH 2/2] Add spatial sweep alternative to FFT-based wind shear perturbation - shear.py: add compute_shear_spatial + _cauchy_hilbert_sweep (numba JIT) - shear.py: add method parameter to WindShear.__init__ ('fft'|'spatial') - wind.py: pass method_shear='spatial' to WindShear when selected - constants.py: document 'spatial' as valid method_shear option - tests/test_shear.py: 12 new tests for helper and WindShear Co-authored-by: Sierd <14054272+Sierd@users.noreply.github.com> --- aeolis/constants.py | 2 +- aeolis/shear.py | 184 +++++++++++++++++++++++++++++++++++- aeolis/tests/test_shear.py | 186 +++++++++++++++++++++++++++++++++++++ aeolis/wind.py | 5 +- 4 files changed, 370 insertions(+), 7 deletions(-) create mode 100644 aeolis/tests/test_shear.py diff --git a/aeolis/constants.py b/aeolis/constants.py index eae691bc..386b94ee 100644 --- a/aeolis/constants.py +++ b/aeolis/constants.py @@ -320,7 +320,7 @@ 'method_transport' : 'bagnold', # Name of method to compute equilibrium sediment transport rate 'method_roughness' : 'constant', # Name of method to compute the roughness height z0, note that here the z0 = k, which does not follow the definition of Nikuradse where z0 = k/30. 'method_grainspeed' : 'windspeed', # Name of method to assume/compute grainspeed (windspeed, duran, constant) - 'method_shear' : 'fft', # Name of method to compute topographic effects on wind shear stress (fft, quasi2d, duna2d (experimental)) + 'method_shear' : 'fft', # Name of method to compute topographic effects on wind shear stress (fft, spatial, quasi2d, duna2d (experimental)) 'max_error' : 1e-6, # [-] Maximum error at which to quit iterative solution in implicit numerical schemes 'max_iter' : 1000, # [-] Maximum number of iterations at which to quit iterative solution in implicit numerical schemes 'max_iter_ava' : 1000, # [-] Maximum number of iterations at which to quit iterative solution in avalanching calculation diff --git a/aeolis/shear.py b/aeolis/shear.py index 5b2ae4e0..32fa0f4c 100644 --- a/aeolis/shear.py +++ b/aeolis/shear.py @@ -35,6 +35,7 @@ #import scipy.interpolate as spint #import scipy.spatial.qhull as qhull import time +from numba import njit # package modules from aeolis.utils import * @@ -90,7 +91,7 @@ class WindShear: def __init__(self, x, y, z, dx, dy, L, l, z0, - buffer_width, buffer_relaxation=None): + buffer_width, buffer_relaxation=None, method='fft'): '''Class initialization Parameters @@ -120,6 +121,10 @@ def __init__(self, x, y, z, dx, dy, L, l, z0, Height of inner layer (default: 10) z0 : float, optional Aerodynamic roughness (default: .001) + method : str, optional + Method to compute wind shear perturbation. Options are + ``'fft'`` (default, 2D FFT-based Weng et al. approach) and + ``'spatial'`` (loop-based Cauchy sweep that avoids FFT). ''' @@ -140,6 +145,11 @@ def __init__(self, x, y, z, dx, dy, L, l, z0, self.L = L self.l = l self.z0 = z0 + + # Method for shear perturbation computation + if method not in ('fft', 'spatial'): + raise ValueError("method must be 'fft' or 'spatial', got %r" % method) + self.method = method def __call__(self, x, y, z, taux, tauy, u0, udir, @@ -229,7 +239,10 @@ def __call__(self, x, y, z, taux, tauy, u0, udir, gc['z'] = np.maximum(gc['z'], zsep) # Compute wind shear stresses on computational grid - self.compute_shear(u0, nfilter=(1., 2.)) + if self.method == 'spatial': + self.compute_shear_spatial(u0, nfilter=(1., 2.)) + else: + self.compute_shear(u0, nfilter=(1., 2.)) # Add shear self.add_shear() @@ -627,8 +640,126 @@ def compute_shear(self, u0, nfilter=(1., 2.)): gc['dtauy'] = np.real(np.fft.ifft2(dtauy_t)) - - + + def compute_shear_spatial(self, u0, nfilter=(1., 2.)): + '''Compute wind shear perturbation using a spatial sweep algorithm. + + Implements the same Weng et al. (1991) perturbation theory as + :func:`compute_shear` but avoids the 2-D FFT entirely. Instead the + computation is broken into a row-by-row (y-constant) sweep in the + wind (x) direction using a discrete Cauchy / Hilbert-transform kernel, + analogous to the sweeping algorithm used for the sediment-concentration + field ``Ct``. + + For each y-row the streamwise shear perturbation ``dtaux`` is + computed as the Cauchy principal-value integral of the streamwise bed + slope, and the cross-wind perturbation ``dtauy`` is computed as the + same integral applied to the cross-wind bed slope. The amplitude + coefficients are derived directly from the inner-layer parameters of + the Weng et al. theory (``L``, ``l``, ``z0``). + + High-frequency topographic features are suppressed with a spatial + Gaussian filter (replacing the FFT-based filter used in + :func:`compute_shear`). + + .. note:: + The computational cost scales as O(ny × nx²) because for each + of the ny × nx cells the kernel sums over all nx other cells in + the same row. For large grids this is more expensive than the + FFT method (O(ny × nx × log nx)). The inner loop is compiled + with Numba to keep the constant factor small, but for grids + significantly larger than ~200 × 200 cells the FFT method may + be preferable on performance grounds. + + Parameters + ---------- + u0 : float + Free-flow wind speed + nfilter : 2-tuple + Width of the Gaussian pre-filter applied to the topography + (in grid cells). Corresponds to the wavenumber range used in + :func:`filter_highfrequenies`. + ''' + gc = self.cgrid + gi = self.igrid + + if u0 == 0.: + gc['dtaux'] = np.zeros(gc['z'].shape) + gc['dtauy'] = np.zeros(gc['z'].shape) + return + + ny, nx = gc['z'].shape + dx = gc['dx'] + dy = gc['dy'] + + z0 = self.z0 + L = self.L / 4. + l = self.l + + # Interpolate roughness length z0 to computational grid (same as FFT method) + if np.size(z0) > 1: + z0new = self.interpolate(gi['x'], gi['y'], z0, gc['x'], gc['y'], 0) + else: + z0new = z0 + + # Inner layer height (iterative, same as FFT method) + for _ in range(5): + l = 2 * 0.41**2 * L / np.log(l / z0new) + + # Middle layer height (iterative, same as FFT method) + hm = 1.0 + for _ in range(5): + hm = L / np.sqrt(np.log(hm / z0new)) + + # Non-dimensional velocity (same as FFT method) + ul = np.log(l / z0new) / np.log(hm / z0new) + + # Scalar effective values when z0 is spatially varying + if np.size(z0new) > 1: + log_l_z0 = float(np.mean(np.log(l / z0new))) + ul_eff = float(np.mean(ul)) + else: + log_l_z0 = float(np.log(l / z0new)) + ul_eff = float(ul) + + # Apply Gaussian spatial filter to topography. + # This replaces the FFT-based high-frequency filter: the sigma of the + # Gaussian is taken as the upper bound of nfilter (in grid cells). + if nfilter is not None: + sigma_filter = float(np.max(nfilter)) + z_filtered = ndimage.gaussian_filter(gc['z'], sigma=sigma_filter) + else: + z_filtered = gc['z'].copy() + + # Amplitude coefficients derived from the Weng et al. leading-order + # perturbation theory. The dominant contribution to the streamwise + # perturbation scales as 2*log(l/z0)/ul^2, and the cross-wind + # contribution scales as 2/ul^2. + alpha_s = 2.0 * log_l_z0 / ul_eff**2 # streamwise amplitude + alpha_n = 2.0 / ul_eff**2 # cross-wind amplitude + + # Compute bed gradients by central differences + dzds = np.zeros_like(z_filtered) # streamwise (x) slope + dzdn = np.zeros_like(z_filtered) # cross-wind (y) slope + + dzds[:, 1:-1] = (z_filtered[:, 2:] - z_filtered[:, :-2]) / (2. * dx) + dzds[:, 0] = (z_filtered[:, 1] - z_filtered[:, 0]) / dx + dzds[:, -1] = (z_filtered[:, -1] - z_filtered[:, -2]) / dx + + if ny > 1: + dzdn[1:-1, :] = (z_filtered[2:, :] - z_filtered[:-2, :]) / (2. * dy) + dzdn[0, :] = (z_filtered[1, :] - z_filtered[0, :]) / dy + dzdn[-1, :] = (z_filtered[-1, :] - z_filtered[-2, :]) / dy + + # Compute shear perturbations via row-by-row Cauchy sweep. + # For each y-row the Hilbert-transform kernel 1/(k*pi) is applied + # to the streamwise (or cross-wind) bed slope, giving an antisymmetric + # response: increased shear on the windward side and reduced shear in + # the lee – consistent with the Weng et al. perturbation theory. + gc['dtaux'] = _cauchy_hilbert_sweep(dzds, alpha_s, ny, nx) + gc['dtauy'] = _cauchy_hilbert_sweep(dzdn, alpha_n, ny, nx) + + def separation_shear(self, hsep): '''Reduces the computed wind shear perturbation below the separation surface to mimic the turbulence effects in the @@ -918,4 +1049,49 @@ def interpolate(self, x, y, z, xi, yi, z0): + +@njit(cache=True) +def _cauchy_hilbert_sweep(field, alpha, ny, nx): + '''Compute wind shear perturbation via row-by-row Cauchy sweep. + + For each y-row the discrete Cauchy principal-value integral of ``field`` + is evaluated using the antisymmetric kernel ``1/(k * pi)``, which + implements the Hilbert transform of the input field in the x (wind) + direction. This gives an antisymmetric perturbation around each + topographic feature: positive (increased shear) on the windward side and + negative (reduced shear) on the leeward side – consistent with the + Weng et al. perturbation theory. + + The sweep processes all y-rows independently and can therefore be + thought of as a 1-D sweep in the wind direction applied row-by-row, + analogous to the directional sweep used for the sediment-concentration + field ``Ct``. + + Parameters + ---------- + field : numpy.ndarray + 2-D input field of shape (ny, nx). Typically the streamwise or + cross-wind bed-slope gradient. + alpha : float + Amplitude coefficient derived from the Weng et al. inner-layer + parameters. + ny, nx : int + Grid dimensions. + + Returns + ------- + result : numpy.ndarray + Wind shear perturbation field of shape (ny, nx). + ''' + result = np.zeros((ny, nx)) + for j in range(ny): + for i in range(nx): + integ = 0.0 + for k in range(1, nx): + im = (i - k) % nx # upstream index (periodic) + ip = (i + k) % nx # downstream index (periodic) + # Antisymmetric Cauchy kernel: upstream minus downstream + integ += (field[j, im] - field[j, ip]) / (k * np.pi) + result[j, i] = alpha * integ + return result diff --git a/aeolis/tests/test_shear.py b/aeolis/tests/test_shear.py new file mode 100644 index 00000000..3c226997 --- /dev/null +++ b/aeolis/tests/test_shear.py @@ -0,0 +1,186 @@ +""" +Tests for aeolis.shear – specifically the spatial sweep alternative to the +FFT-based compute_shear method. + +Run with:: + + pytest aeolis/tests/test_shear.py +""" + +import numpy as np +import pytest + +from aeolis.shear import WindShear, _cauchy_hilbert_sweep + + +# --------------------------------------------------------------------------- +# Helpers +# --------------------------------------------------------------------------- + +def _make_hill_grid(ny=10, nx=40, dx=2.0, dy=2.0, hill_height=3.0, hill_sigma=10.0): + """Return a small 2-D grid with a centred Gaussian hill.""" + x1d = np.arange(nx) * dx - nx * dx / 2.0 + y1d = np.arange(ny) * dy - ny * dy / 2.0 + xx, yy = np.meshgrid(x1d, y1d) + z = hill_height * np.exp(-(xx**2 + yy**2) / (2.0 * hill_sigma**2)) + return xx, yy, z + + +def _run_windshear(method, xx, yy, z, dx=2.0, dy=2.0, + L=100.0, l=10.0, z0=0.001, buffer_width=20.0, + u0=10.0, udir=270.0): + """Instantiate WindShear with the given method and run it.""" + ws = WindShear(xx, yy, z, dx=dx, dy=dy, L=L, l=l, z0=z0, + buffer_width=buffer_width, method=method) + taux0 = np.ones_like(z) * 0.1 + tauy0 = np.zeros_like(z) + ws(x=xx, y=yy, z=z, + taux=taux0, tauy=tauy0, + u0=u0, udir=udir, + process_separation=False, c=0.5, mu_b=30.0, + taus0=0.1, taun0=0.0, + sep_filter_iterations=0, zsep_y_filter=False) + return ws.get_shear() + + +# --------------------------------------------------------------------------- +# Tests for _cauchy_hilbert_sweep +# --------------------------------------------------------------------------- + +class TestCauchyHilbertSweep: + """Unit tests for the module-level numba-compiled helper.""" + + def test_flat_bed_gives_zero(self): + """A flat bed has no slope, so the shear perturbation must be zero.""" + ny, nx = 5, 30 + field = np.zeros((ny, nx)) + result = _cauchy_hilbert_sweep(field, 1.0, ny, nx) + np.testing.assert_allclose(result, 0.0, atol=1e-12) + + def test_output_shape(self): + """Output shape must match input shape.""" + ny, nx = 7, 25 + field = np.random.rand(ny, nx) + result = _cauchy_hilbert_sweep(field, 1.0, ny, nx) + assert result.shape == (ny, nx) + + def test_antisymmetric_response_around_symmetric_hill(self): + """For a symmetric Gaussian hill, the perturbation must be + antisymmetric: positive on the windward side and negative on the + leeward side of the crest.""" + nx = 60 + dx = 1.0 + x = np.arange(nx) * dx - nx * dx / 2.0 + z = 5.0 * np.exp(-x**2 / (2.0 * 8.0**2)) + dzdx = np.gradient(z, dx) + field = dzdx.reshape(1, nx) + + result = _cauchy_hilbert_sweep(field, 1.0, 1, nx) + row = result[0] + + # The hill is centred at x=0, which falls at index nx//2. + # Indices 0..(nx//2 - 1) correspond to x < 0 (windward / upwind side). + # Indices nx//2..nx-1 correspond to x >= 0 (leeward / downwind side). + # The windward half should contain the peak perturbation (increased shear) + # and the leeward half should contain the trough (reduced shear). + windward = row[:nx // 2] # x < 0, upwind of crest + leeward = row[nx // 2:] # x >= 0, at/downwind of crest + + assert windward.max() > 0.0, "Expected positive perturbation on windward side" + assert leeward.min() < 0.0, "Expected negative perturbation on leeward side" + + def test_amplitude_scales_with_alpha(self): + """Output must scale linearly with alpha.""" + ny, nx = 4, 20 + field = np.random.rand(ny, nx) * 0.1 + r1 = _cauchy_hilbert_sweep(field, 1.0, ny, nx) + r2 = _cauchy_hilbert_sweep(field, 3.0, ny, nx) + np.testing.assert_allclose(r2, 3.0 * r1, rtol=1e-12) + + def test_row_independence(self): + """Each y-row must be processed independently.""" + ny, nx = 4, 20 + # Make a field where every row is identical + row = np.random.rand(nx) * 0.1 + field = np.tile(row, (ny, 1)) + result = _cauchy_hilbert_sweep(field, 1.0, ny, nx) + for j in range(1, ny): + np.testing.assert_allclose(result[j], result[0], rtol=1e-12) + + +# --------------------------------------------------------------------------- +# Tests for WindShear with method='spatial' +# --------------------------------------------------------------------------- + +class TestWindShearSpatial: + """Integration tests for WindShear.compute_shear_spatial.""" + + @pytest.fixture + def hill_grid(self): + return _make_hill_grid() + + def test_init_accepts_spatial_method(self, hill_grid): + """WindShear must accept method='spatial' without raising.""" + xx, yy, z = hill_grid + ws = WindShear(xx, yy, z, dx=2.0, dy=2.0, L=100.0, l=10.0, + z0=0.001, buffer_width=20.0, method='spatial') + assert ws.method == 'spatial' + + def test_init_rejects_invalid_method(self, hill_grid): + """WindShear must raise ValueError for unknown method strings.""" + xx, yy, z = hill_grid + with pytest.raises(ValueError, match="method must be"): + WindShear(xx, yy, z, dx=2.0, dy=2.0, L=100.0, l=10.0, + z0=0.001, buffer_width=20.0, method='unknown') + + def test_zero_wind_gives_zero_perturbation(self, hill_grid): + """With u0=0 the shear perturbation must be exactly zero.""" + xx, yy, z = hill_grid + ws = WindShear(xx, yy, z, dx=2.0, dy=2.0, L=100.0, l=10.0, + z0=0.001, buffer_width=20.0, method='spatial') + taux0 = np.ones_like(z) * 0.1 + tauy0 = np.zeros_like(z) + ws(x=xx, y=yy, z=z, + taux=taux0, tauy=tauy0, + u0=0.0, udir=270.0, + process_separation=False, c=0.5, mu_b=30.0, + taus0=0.1, taun0=0.0, + sep_filter_iterations=0, zsep_y_filter=False) + taux, tauy = ws.get_shear() + # With no wind the stresses are returned unchanged (no perturbation) + np.testing.assert_allclose(taux, taux0, atol=1e-10) + + def test_spatial_produces_nonzero_perturbation(self, hill_grid): + """The spatial method must produce a non-trivial shear perturbation + for a Gaussian hill and non-zero wind.""" + xx, yy, z = hill_grid + taux, tauy = _run_windshear('spatial', xx, yy, z) + # There should be some cells above and below the flat-bed reference + assert taux.max() > 0.1, "Expected positive taux somewhere" + assert taux.min() >= 0.0, "taux must be clipped to non-negative" + + def test_spatial_and_fft_produce_correlated_results(self, hill_grid): + """The spatial and FFT methods should produce qualitatively similar + taux fields (correlation > 0.7).""" + xx, yy, z = hill_grid + taux_fft, _ = _run_windshear('fft', xx, yy, z) + taux_sp, _ = _run_windshear('spatial', xx, yy, z) + corr = np.corrcoef(taux_fft.flatten(), taux_sp.flatten())[0, 1] + assert corr > 0.7, ( + f"Expected correlation > 0.7 between FFT and spatial methods, " + f"got {corr:.4f}" + ) + + def test_spatial_output_shape_matches_input(self, hill_grid): + """The output shear fields must have the same shape as the input grid.""" + xx, yy, z = hill_grid + taux, tauy = _run_windshear('spatial', xx, yy, z) + assert taux.shape == z.shape + assert tauy.shape == z.shape + + def test_fft_method_still_works(self, hill_grid): + """The original FFT method must continue to work after the refactor.""" + xx, yy, z = hill_grid + taux, tauy = _run_windshear('fft', xx, yy, z) + assert taux.max() > 0.1 + assert taux.shape == z.shape diff --git a/aeolis/wind.py b/aeolis/wind.py index db82c030..dacacf9b 100644 --- a/aeolis/wind.py +++ b/aeolis/wind.py @@ -72,11 +72,12 @@ def initialize(s, p): if p['process_shear']: if p['ny'] > 0: - if p['method_shear'] == 'fft': + if p['method_shear'] in ('fft', 'spatial'): s['shear'] = aeolis.shear.WindShear(s['x'], s['y'], s['zb'], dx=p['dx'], dy=p['dy'], L=p['L'], l=p['l'], z0=z0, - buffer_width=p['buffer_width']) + buffer_width=p['buffer_width'], + method=p['method_shear']) else: s['shear'] = aeolis.rotation.rotationClass(s['x'], s['y'], s['zb'], dx=p['dx'], dy=p['dy'],